A mesh-refinement scheme for finite element computations

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 7 (1976) 93-105 0 NORTH-HOLLAND PUBLISHING COMPANY

A MESH-REFINEMENT

SCHEME FOR FINITE

ELEMENT

COMPUTATIONS

Graham F. CAREY Center for Quantitative Science in Forestry, Fisheries, and Widlife University of Washington, Seattle, Washington 98195,

USA

Received 6 February 1975 The paper considers the use of unorthodox grids where rapid transition from refined zones to coarser zones is effected, thus introducing exposed nodal freedoms at the zone interfaces. A technique for automated mesh enrichment of finite element discretizations is devised. Suitable convergence criteria are designed to delineate those subregions of a domain where rermement is necessary. Enriched and unaltered regions are separated by a fringe of semirefined elements. A valid finite element theory is provided for the fringe elements. A pilot numerical study illustrates the practical implementation of the automated refinement strategy. The principal features of the program are described.

1. Introduction Certain practical considerations arise during finite element applications in connection with the choice of an appropriate idealization. A related question concerns local mesh point adequacy, solution accuracy, and mode of refinemknt. There are two basic approaches possible to improve solution accuracy for a given finite element idealization: (1) Higher order elements can be introduced and (2) the mesh can be refined. The former scheme is less flexible from the viewpoint of programmed implementation and leads to difficulties similar to those encountered in direct refinement. Therefore, we shall be concerned here with the mesh enrichment approach. Directly refining the mesh introduces new node points increasing the size of the subsequent algebraic problem and inevitably refining the mesh in regions where enrichment is unnecessary. Manual description of the refined idealization can avoid this to some extent, but as problem size increases and three-dimensional problems become increasingly important, this approach becomes impractical. Thus the established trend of automatic mesh generation and verification pre- and post-processers must continue to explore suitable refinement schemes. Ideally, one wishes to enrich the mesh only in those regions where the solution is not adequate. In one-dimensional finite element or finite difference methods this is achieved simply by including additional nodes in the appropriate intervals. However, in higher dimensions this cannot be so readily achieved, as we observed above. A simple scheme, which is not generally satisfactory for obvious reasons, is to “window-in” on a region where mesh enrichment is necessary (such as a crack tip). The size of the domain is reduced with each refinement step, using the solution on the previous mesh to generate boundary data for the new “window”. The aim of this present investigation is to devise a suitable, readily automated, local enrichment strategy for two and higher dimensions. In principle, this can be attempted in two ways. One method is to refine the solution locally and restrict its behavior at the edge of the refined zone by introducing constraints enforced through Lagrange multipliers in the variational problem (somewha

G.F. Carey, A mesh-refinement

94

scheme for finite element computations

analogously to some hybrid element formulations). This increases the number of unknowns to include the undetermined multipliers. Further, it is difficult to incorporate in existing programs without extensive reprogramming. Nevertheless, it is a viable technique [ 11. The second approach, developed herein, is to imbed the constraint directly in the element basis for those elements at the edge of the enriched zone and so introduce a semirefined element. The method is straightforward and can be readily incorporated in existing programs. The following investigation then examines criteria for solution adequacy and proposes a suitable automated refinement strategy. For clarity the study treats two-dimensional domains and a piecewise-linear basis on triangular elements applied to the Dirichlet integral. However, the extension to higher dimensions and alternative element types is immediate. The important features of this investigation concern both the underlying analysis and the programmed implementation. Refinement is advocated by subdivision of a triangular “parent” element into four congruent triangular subelements. fn the subdivision scheme, internal interfacing boundaries confine locally enriched regions, outside of which the former unaltered mesh is retained. The fully subdivided and unaltered zones are separated by a fringe of semirefined elements with exposed midside solution components. removed by linear solution constraints on that side. A pilot numerical study is devised to illustrate the automated refinement strategy and employs nested confined zones of subdivision for successive mesh-refinement iterations. A description of the principal program design aspects is included to highlight the necessary distinctive steps involved. Several immediate applications - such as stress analysis of cutouts in panels and, generally, problems with regions of marked solution variation - are evident. On this basis, the scheme can be implemented at considerable advantage into special purpose and general sophisticated finite element computer systems of the types presently existing in aerospace and other industries.

2. Discussion Most computations utilize the analyst’s experience of expected physical behavior in the problem and rely on judgment to devise an appropriate mesh. This is often sufficient, particularly if a series of problems are very similar and solutions are well-behaved. Boundary layer and other phenomena give rise to problems in which the mesh distribution is critical and cannot be determined a priori. Refinement by trial and error, if successful, will usually lead to a solution involving many more nodal equations than necessary. As the solution of linear systems by elimination is of order n3 (or w2n operations for bandwidth w), where IZis the number of equations, the computation time is adversely affected. The teal-and-error iteration itself constitutes a major concern with respect to Bow times for successive data preparations, verifications, and executions. Several accuracy criteria are possible for the nodal solution + in a specific finite element idealization: (I) The (weighted) mean-square residual norm on an element is

i/z J’w{~)~(#)’ 1 dR,

,

w(x)

> 0 >

&

where w(x) is the weight function and N(G) is the governing differential operator. The solution may be judged adequate if IIN,liZis less than some prescribed tolerance e. (This norm is particularly

G.F. Carey, A mesh-refinement scheme for jkite element computations

95

appropriate for weighted-residual finite element applications such as the Galerkin method.) (2) In another way, the variational (energy) functional itself induces a norm in

111,ll=

[.I’Wb#x,#+.)d~, I

>

L positive definite.

R,

Here L is the Lagrangian and e denotes the general element. The stationary condition on I is insufficient alone, and to test the adequacy of the mesh we require comparisons of ll1,ll for successive solutions with refined meshes. (3) A pointwise error control on successive solution refinements can be particularly useful, as it directly indicates which regions require successive refinement. (4) Finally, possibly the simplest though least precise scheme is to locate refinement zones as those regions where the solution is most rapidly changing. This is, in fact, the criterion applied in the numerical example considered later. Engineering experience has shown that idealizations with slender elements may often yield poor results. This situation has subsequently been mathematically examined, and the conclusion is that use of slender elements is generally to be avoided. This requires particular note with respect to automatic mesh refinement schemes. To enlarge on this, the reasons for a good subdivision scheme such as the one presently described are two-fold: 1. A desirable feature of any discretization is that the basis be uniform [ 21. This is effectively a geometrical constraint on the element, whereby the smallest angle of each element is bounded below (see the discussion of the subsequent inequality). 2. Slender triangles arising from various subdivisions cause deterioration of the numerical condition of the algebraic system. Both criteria are desirable design features. Some analysts [2, p. 1391 recommend triangulation with element angles greater than w/S, and Synge [ 3, p, 2111 has shown that if the finite element approximation $(x,v) coincides with the exact solution $(x,v) at the nodes,

where K is a bound on I$,,[, d is the diameter of the element, and 8 is the largest angle formed by the sides of the triangle. The need for caution in constructing successive refinements is evident from the inequality, which suggests that slender elements be avoided. Of course, uniformity and numerical conditioning are both problem-dependent properties (note that K above is undetermined). In practice one can construct meshes with slender elements that perform adequately provided the solution is well-behaved in the regions so idealized. However, the concern here is to determine an automatic refinement scheme, and thus some form of the two guidelines above now becomes essential. For this purpose, consider the simplest and best-known element, a triangle with a linear solution approximation - this appears for instance in plane stress and potential flow analyses. The natural refinement of a three-node triangle is to introduce new nodes at the side midpoints, subdividing a triangle into four sub-triangles, as in fig. 1. The additional midside nodes are shared with adjacent elements, necessitating similar partitions in these to ensure compatibility across the interface. Then, for an initial idealization of N, elements and iV,, nodes, the local refinement pattern propagates through the entire network, and the refined region has 4N, elements and precisely

G.F. Carey, A mesh-refinement scheme for finite element computations

96

2N, + N, - 1 nodes. This is proved from an important property of polygonal subdivisions of a plane surface, proposed by Descartes and first proved by Euler: viz. V - E + F = 1, where V is the number of vertices, E is the number of edges, and F is the number of faces. Let subscripts 1 and 2 denote the discretization before and after subdivision respectively. Then Vi - Ej + Fi = 1, i = 1,2. The following relations hold for this triangular subdivision: F, = 4F,, E, = 2E, + 3 F,, so that V,=1tE’L--F2=1+2E1-F1=1t2(P’1+F1-1)-F1=2V’1+F1--1. Suppose now that subdivisions are made repeatedly so that the region assumes a sequence of partitions. If #i designates the value of #(x) at xi in the j-th partition, then the size of I&- @‘I indicates the accuracy of the solution. Certain regions are thus distinguished in which considerable change occurs between refinement iteratesj and jt 1. It is desirable to enrich the mesh in these regions.

Y \ t

IX

Fig. 1. Subdivision of triangle l-2-3 to four subtriangles.

Fig. 2. Piecewise linear solution profile corresponding to mesh gradation across an interface. The broken lines are intended to show solution variation along element edges, not the unrefined solution.

Again consider the subdivided triangle in fig. 1 as part of a local enrichment scheme. Initially, interpolation is linear through nodal solutions Cpi,&, and (pa corresponding to the “parent” element vertices. On subdivision the solution approximation is of higher order over the “parent” element in the sense of piecewise approximation spaces. The central and simple idea is to avoid mesh propagation across an exposed side, say 2-3, by constraining the solution to vary linearly along that edge with #S ‘=$($, + (p3).With this constraint a linear solution is permissible in the neighboring element to the right of side 2-3, and the solution approximation across the interface is continuous. In fig. 2 the approximation profile is indicated for three elements, showing the mesh gradation across such an interface. For a general domain an internal boundary confines locally enriched regions outside of which the former unaltered mesh is retained (fig. 3).

G.F. Carey, A mesh-refinement scheme for finite element computations

97

Unalterec ext me

Fig. 3. Local refinement by subdivision within an internal confining interface.

3. Finite element theory

There have been several instances, initially in research developments of higher-order compatible elements and subsequently in industry, where some form of linear nodal-solution averaging is employed. Such a modified model is most readily obtained by accumulating subelement matrix contributions and then applying the nodal averaging constraint. However, in some cases the application of the constraint condition has been carried out in an imprecise manner and violates the variational principle. That is, the finite element model is not a consistent one in the sense that the final system does not determine a solution that minimizes the energy integral in the constrained approximation subspace. Thus some care should be exercised in incorporating the constraint condition during subelement accumulation. An alternative approach, avoiding this difficulty, is to imbed the constraint in the element basis prior to applying the variational requirement or, equivalently, to use the correct elementary matrix operations that are described later. The correct derivation is outlined below to clarify the issue. Consider a typical subelement i-j-k in fig. 1 with nodal vector I$~= {I#+#,. &}, Then in triangular coordinates 5 the linear element interpolant is

where the triangular (or area) coordinates are local coordinates for the element and are defined as area ratios or equivalently by

with subelement area a and nodal vertices (xn, u,), n = i, j, k. In the case of the Dirichlet integral for the Laplacian the stationary functional contribution

at

G.F. Carey, A mesh-refinement

98

scheme for finite element computations

node IZfor i, =

$

s N#~>~+ (#yj211dR, He

is

K? Ji :K =Re

dx dy .

(3)

Substitution for #J(S) from eq. (1) yields the usual 3 X 3 matrix for this element. Each subelement contribution is now developed for the first integrand expression in eq. (3), but employing the linearity constraint on side l-2 a priori. The general element expressions for variation with respect to 9, are

with a = iA, the subelement area, and p = i, j, k. In order, the subelement matrices are now formed. For element 1, designated l-4-6 with Q4 = k(#, + &), the matrix terms are

with p = 1,2,6. The subelement matrix is then (&)Z +(c11),(~2)X+:($& ~~)z$

5(‘!,),(&+&2);

symmet~c

(rl),(~3)x+~(~2)x(~3)x D&),(&)X

&,)E

(63);

I where $A(E,l, M&l,

=Yj-Yf

=yYi-yk =Yl-Yq

=y4-y6

= k(v2-y3),

kA(E,),

=Yk-yi

=y6-Y1

= $(Y~-Y~),

1

.3

(3

1 and

= k(Yl--Y,).

Repetition of the analysis leading to eq, (5) determines similar matrices for subelements 2, 3 and 4. Accumulating subelement contributions nodally, the fringe element matrix is Kr = K{ + KT with KT again determined from hrfy on replacement of (xi} for Cyi}, i = 1, . . . . 6. The final simplified result for Kf is

99

G.F. Carey, A m~~~-re~~ernent scheme for finite element Gorn~ut~ti~ns

(YZ-Y3)(YZ-Yd &tY3-Y1)2+

cps

@3

4%

@l

cc!-Y3)(Y3-Y1) 2(Yl-Y2)2

+;iCv3-Yd2

+ w-Y2)2

0

(YZ-Y3)

til -Y2)

-(YrY2)2

+(Y2-Y3)21

+(Y2-Y3121

(Y2-Y3)(YZ-Y1)

0

+$[fY3-Y*)2+2(YrY2)*

-(Y1-Y212

(Yl-Y2)

(Y3-Y’

(Y2-Y3)

(Yl -Y:

+tY2-Y3)*l

(Y3-Yl)(Yi-Y2)

(Yl-Y2)2

symmetric

(YZ-Y3j2

+ fy3-Yd2

%Yz-Y3)

(Ys-Y

+ (Y1-Y212

(YZ-Y3)* + (Y3-Yd2 + (Y,-Y212

The consistency of this model can be verified by considering additional constraints to linear variation along sides 2-3 and 3- 1. This yields on simplification the correct element matrix for the three-node triangle l-2-3. The consistent model of eq. (6) can be constructed directly by elementary matrix operations as follows. Letting $4 = ~(c#I,+ c#,), we define the elementary matrix

(7)

If K = (kji) is the 6 X 6 symmetric matrix for the unconstrained is simply -kontracted from E *KE :

K,- =

k,, + %I,

k,, + ik,,

k,, + k&s

k,, + f k,,

k 33

symmetric.

model, the consistent fringe matrix

&,+ %

k 35

k 36

k 55

k 56 k 66

100

G.F. Carey, A mesh-refinement

scheme for finite element computations

which on evaluation yields precisely the result of eq. (6). Some earlier formulations have included the constraint in the assembled matrix by direct substitution of 4, = $(#, + (PJ. This corresponds to the matrix operation EIE rather than EtKE. The resulting element matrix is too flexible (in the sense of structural applications). It is easy to show that constraining all three midside nodes in this way gives Kf = K E, E, E, = f Kf for the 3-node element. A numerical example, illustrating the mesh refinement strategy and use of the element, are taken up subsequently. Prior to this, the distinctive features of computations involving automated refinement are briefly described.

4. Computer analysis Assume we have an existing idealization for a given region and the associated nodal solution. Consider a refinement in some automatically detected region. The variational functional I retains the same form over all “parent” elements. For interior subdivided elements the quadrature is performed over subelements. Newly introduced midside nodes are included in the nodal solution vector for the discretized domain. The interfacing fringe elements are most conveniently treated in a manner similar to the interior subdivided elements. Accumulate subelement contributions to K and then form EtKE to determine Kf. This will necessitate the use of a vector of “dummy” midside nodes at the interface. Only those elements that are refined will introduce new nodal contributions to the difference system. Where new subelements contribute to an old node, the coefficient values will be changed. New nodes will set up new equations increasing the size of the system. In unrefined regions the finite element equations will be unchanged and need not be regenerated. Introduction of new nodes in refinement is very easily accommodated by iterative solution schemes, such as Gauss-Seidel with relaxation, as nodal numbering does not influence solution computations. This favors iterative schemes for finite element computations where automatic repeated refinements occur. On the other hand, most current finite element applications, particularly in structural mechanics, exploit the sparse, banded symmetric, positive-definite coefficient matrix by means of banded Cholesky decomposition. Economic computation hinges on judicious nodal numbering which induces a narrow dense bandwidth. Inclusion of new nodes poses a considerable difficulty, which has been overcome by programming to allow flexible bandwidth, or by automatic nodal renumbering schemes to reoptimize the bandwidth for successive iterations [ 4,Sl. These approaches become even less appealing when mesh refinement is moderate. Regeneration of finite element contributions may be simply avoided in regions where no subdivision is desired through use of pointer arrays relating former to new nodal and element numbering schemes. The reordering problem in enriched zones can be handled in like manner using pointer arrays. One such scheme, very easily implemented, works as follows: the enrichment zones are located as sequences of elements by applying an appropriate mesh point adequacy test of the type previously mentioned. Progressing sequentially in nodal order through the idealization, new nodes are inserted and numbered in succession from adjacent “old” nodes of the former unrefined idealization. As each new node is inserted, the existing nodes of larger index are correspondingly incremented, the pointer arrays recording all such transitions. For example, let several elements represent part of a small zone to be subdivided, with (say) node 7 1 the largest nodal index for elements in this zone of the existing well-numbered partition. Introduce new midside nodes 72, 73, . .. . 86 (say)

G.F. Carey, A yeah-re~ne~ent scheme for finite element computations

101

for the subdivision in this local zone. All former node indices greater than 7 1 are simultaneously incremented by 15 so that previous node 72 is now node 87, and so on. The confined nature of the refinement strategy implies that an initial well-numbered nodal set will be maintained under such a scheme. A great variety of similar approaches of varying degree of sophistication and computational expense can be constructed. The general sequence of automated operations is indicated below: 1. Generate the initial idealization. 2. Compute the solution vector. 3. Scan the mesh and use an appropriate test to delineate zones of refinement. Exit if the solution is acceptable. 4. Refine the local zones by subdivision, inserting new nodes. 5. Compute element matrices in the subdivided zones, including the fringe elements. 6. Assemble contributions of the refined elements into the gross (enlarged) coefficient matrix. 7. Return to 2. One important point is that extensive application of mesh refinement techniques may lead us to re-evaluate the role of iterative solution algorithms, at least to a limited extent. Most linear structural problems require solutions for a variety of loading conditions, that is many right-hand sides. Thus, iterative methods are quite inappropriate. However, for problems with, say, a single right-hand side and requiring extensive refinement, the solution for a given intermediate mesh provides an excellent starting guess for the iterative algorithm at the next refinement. Conceivably a combination of the iterative and direct methods would be a suitable arrangement, the former being used to prepare an idealization to which the direct method can be applied.

5. Numerical investigation Conceptually the automatic mesh refinement scheme is quite straightforward, but programming for very general applications is intricate. The principal features of such general programmed implementation are elaborated in the preceding discussion. To indicate the procedure involved, without becoming immersed in the programming details necessary for production programs in industry, an illustrative “pilot” study is devised. The problem is posed on a square domain with the initial 3 X 3 regular partition indicated in fig. 4. The initial ide~ization is program-generated, and for expedience the refinement strategy is restricted to subdivide regular squares in the partition, guided by a simplified mesh-point adequacy test. The problem is to determine an approximate solution to the Poisson equation A$=S(x-i)6(y-i)

in

(0,1)X(0,1),

(9)

with Cp(x,y) = 0 on the boundary. Physically, $J(x, y) can be interpreted as the electrostatic potential in a unit square, corresponding to a line charge at the center of the square and with boundary at zero potential. Mathematically, an interpretation is that cp(x,y) is the Green’s function g(x, y; & Q) evaluated for a source at (F, rl) = ($, i). The associated low-order variational functional is (see [ 61)

102

G.F. Carey, A mesh-refinement

I= f

scheme for finite element computations

~td~,,dR - jbtx, VI 6(x - 5)S(v R

t> dR ,

R

while, from the fundamental property of a delta function, I = 5 s(#J2

dR - &,;)

.

(10)

R

Introducing the finite element discretization, let n be the nodal index for the singular point (i , f), and let i be the nodal index for any nodal point in the discretization. The approximated functional is denoted I”,and the variational statement Sl”= 0 then implies that for the nodal vector $,

$=F%= CA,(L,L:+L,L:)Q~,-~,=~.

(11)

e

Here a/ad, is the gradient (with components a/&#+), Te is the element functional contribution, A, is the element area, and Cp,(x, y) = L(x, y)‘$, with subscript e the element index; e, is a vector with nth component unity and all other components zero. The mesh-improvement scheme is automated as follows. Solve the problem with the coarse initial mesh. Over the initial 3 X 3 partition of fig. 4 the nodal solution is scanned, comparing neighboring nodal values. Only the central nodes of each square are interrogated, and those where most solution change is observed locate the refinement zones. In the first pass the square subregion [ $, $1 X [k, $1 centered at ($, 4) is selected for enrichment. Following enrichment, the finite element solution is repeated with the new mesh, the new nodal solution examined, and the mesh again enriched where necessary. Again, a central square about the s~g~a~ty at (4, i) is subdivided and the solution is repeated. The refinement ceases at a specified difference tolerated in successive solution at the critical node(s). The successive nested zones are shown in fii. 4. An analytic expression for the Green’s function in the square domain can be obtained directly by reflection (the method of images of potential theory) and yields an infinite series corresponding to the accumulated potential of the resulting distribution of sources and sinks. Successive reflections in the four side-planes generate a full, doubly infinite, two-dimensional lattice for potential summations. In addition to the complexity of this summation is the problem of accurate numerical evaluation of the infinite series. An alternative and more satisfactory procedure, particularly in view of the central location of the source delta function, utilizes conformal mapping. For a square domain the SchwarzChristoffel theory applies, and the mapping w

Z=

sJzF

(12)

0

(an elliptic integral of the first kind [ 7 1) conformally maps the unit circle onto the interior of the square with diagonal

G.F. Carey, A mesh-refinement scheme for jInite element computations

103

Fig. 4. Nested refinements of 3 X 3 partition for numerical example with singularly at (i, i).

The inverse mapping follows as w =

1 sn(flzIm)

x

(13)

dn(flzlm)’

withm =$. The special functions I’, sn and dn are the gamma function and Jacobian elliptic functions and are extensively tabulated. A second, simpler conformal transformation c = z/w maps the square of diagonal l?2(f)/fi onto the given domain of diagonal dF, where cy= r2(&‘(2ti-,). Then the final transformation is (14) which maps the given square onto the circle. The boundary condition @= 0 holds on the circle boundary, and the source point maps onto the circle center. The Green’s function in the unit circle for a source at (0,O) is simply G(u, u; 0,O) = -+oglwl,

w=u+iv.

Substituting w = f(c), the desired Green’s function on the square with center at (0, 0) is

11.9

G.F.Carey, A mesh-rejinement

104

g(x,y;O,O)=-&log

scheme for fmite element computations

1 sntfil~3.I f>

t/z dn (fl$,li)

(16)

In particular, along one diagonal in the r-plane, c is real, and the analytic result for g(x, y; 0,O) is sketched in fig. 5 interpolated at discrete { points with sn and dn computed from tabulated theta functions [8, pp. 58 l-5831. The solution is graphed radially along y = x with i G x < 1 in this figure for each of the three iterated solution meshes of fig. 4. The linear system is solved in each case by Gauss-Seidel iteration with the stiffness matrix completely recalculated for each mesh iteration.

***+hh=.*= Analytic solution -----Second refinement --First refine~nt Enitiol mesh

0.6

0.5

0.4 , ; 0.3 k 0.2

0.

I 0 x-

Fig. 5. Solution @(x,y) alongy =x, i
The intention in the numerical example is not to imply that problems with singularities be attempted by repeated mesh enrichment with simple element models. There is a recent and rapidly growing literature on the use of singular basis functions for these problems [9,10]. Of course, this does not exclude incorporation of automated enrichment schemes in these formulations. Rather, the example presents a convenient simple situation to illustrate the enrichment theory and implementation. In practice the method might be employed to some advantage for such problems as stress analysis of “cutouts” (for example, in aircraft window design) where accurate stress predictions are warranted. 6. Conclusion

A simple and easily implemented technique suited to automated mesh refinements of general finite element discretizations has been proposed. The scheme is based on subdivision of elements of the current mesh in regions where enrichment is necessary. Without loss of generality, triangular subdivision of a two-dimensional region with analysis by a piecewise linear basis describes the

G.F. Carey, A mesh-refinement

scheme for finite element computations

105

process. The fully subdivided and unaltered zones are separated by a fringe of semirefined elements with exposed midside solution component removed by linear solution constraint on that side. A consistent model analysis provides a valid finite element theory for the fringe elements of the refinement scheme as well as an approach for the formation of correct models in other unrelated applications which require nodally constrained elements. A numerical example with simplified mesh-point adequacy and enrichment criteria illustrates the concepts of automated mesh refinement discussed. Poisson’s equation on a unit square is considered with solution singularity at the central point ($, i). The numerical solution employs nested confined zones of subdivision for successive mesh refinement iterations. Important industrial implementations are conceptually the same but demand more sophisticated program design and tactics. Computational features which differ from conventional nonrefining analyses are discussed. It is inferred that difficulties concerning inclusion of the automatic refinement scheme in existing finite element programs will be quite minimal. One cautionary note must be added. The concept presumes the usual demand of a reasonable initial idealization, avoiding very slender elements and too rapid gradation in element size. Further constraints on the number of iterations may be programmed for such problems as a cracked panel in plane stress or boundary layers in fluids. In problems of this type an extrapolatory refinement scheme can be imbedded in the analysis. Then refinement will not proceed as a single “parent” element partition, but instead as multiple subdivisions of that type within nested subelements prior to each solution iterate. In this last point problems of round-off and ill-conditioning with excessive refinement must be considered, and a possible safequard is to underrelax the mesh extrapolation.

Acknowledgements

The author expresses his appreciation to Carl Pearson, Jim Tocher and Carlos Felippa for contributing through discussion of the topic.

References [ 1] I. J. SomervailIe and A.P. Kabaila, Mesh grading techniques for compatible and equilibrium elements, in: Proc. Int. Conf. on Finite Element Methods in Engineering, V.A. Pulmano and A.P. Kabaila (eds.) (Clarendon Press, Sydney, Aug. 1974) 251-271. [2] G. Strang and G. Fix, An analysis of the finite element method (Prentice Hall, 1973). [3] J.L. Synge, The hypercircle in mathematical physics (Cambridge Univ. Press, 1957). [4] M. Hellborg, Design of a CDC 6600 computer program for the solution of very large systems of equations, Boeing Dot. D6-29662 TN, Dot-Ads, Renton (1968). [S] B.M. Irons, A frontal solution program for finite element analysis, Int. J. Num. Meth. Eng. 2 (1970) 5-32. [6] H.C. Martin and G.F. Carey, Introduction to finite element analysis (McGraw-Hill, New York, 1973). [7] P.F. Byrd and M.D. Friedman, Handbook of elliptic integrals for engineers and scientists, 2nd ed. (Springer, 1971). [8] M. Abramowitz and I. Stegun (eds.), Handbook of mathematical functions (AMS 55) (National Bureau of Standards, 1964). [9] G. Fix, S. Gulati and I. Wakoff, On the use of singular functions with the finite element method, J. Comp. Phys. 2 (1973) 209-228. [IO] G.F. Carey, Basis functions in finite element theory and applications, in: Proc. Int. Conf. on Finite Element Methods in Engineering, V.A. Pulmano and A.P. Kabaila (eds.) (Clarendon Press, Sydney, 1974) 55-74.