UNIVERS1TV0F „HSH5S- The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. To renew call Telephone Center, 333-8400 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN MAR22W3 JAN 2 mi DEC 15 1917 $£P 8 19BB L161— O-1096 /)0-/&Jl/ UIUCDCS-R-80-1021 UILU-ENG 80 1712 Uf. RAPID SOLUTION OF FINITE ELEMENT EQUATIONS ON LOCALLY REFINED GRIDS BY MULTI-LEVEL METHODS by John R. Van Rosendale [tHE LIBRARY OF THE JUW 3 U isou UNIVERSITY OF ILLINOIS URBANA-CHAMPAIGN May 1980 UIUCDCS-R-80-1021 RAPID SOLUTION OF FINITE ELEMENT EQUATIONS ON LOCALLY REFINED GRIDS BY MULTI-LEVEL METHODS by John R. Van Rosendale May 1980 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA- CHAMPAIGN URBANA, ILLINOIS 61801 Supported In part by the Air Force Office of Scientific Research Grant AFOSR 75-7854 and submitted in partial fulfillment of the requirements of the Graduate College for the degree of Doctor of Philosophy. Digitized by the Internet Archive in 2013 http://archive.org/details/rapidsolutionoff1021vanr iii ACKNOWLEDGMENTS Too many people contributed to the successful completion of this thesis for me to thank them all explicitly, but I will try to mention those who seemed to me to have contributed most. First I must thank my advisor, Bob Skeel, who supported my research and tried to make sense of it even as I sank into a quagmire of Hilbert spaces. Whatever rigor and professionalism is evidenced here is largely due to his help. Secondly, I would like to thank the many friends who supported me; first among these, Dennis Gannon, who shared many of his ideas about finite elements with me. I would also like to thank Paul Saylor, who first introduced me to multi-level methods. (At times I have not been sure whether to thank him for that.) Dave Ferguson, who helped on parts of the manuscript also deserves thanks. The other people who's support was more psychological must also be mentioned, my parents, Mitch and Diane, Bert, Rick, Qasim (who insisted his name be included), and many others. Their encouragement was very important. Finally I would like to thank Barbara Armstrong for her excellent typing in somewhat trying circumstances. She does not deserve to be last on the list. iv TABLE OF CONTENTS Page 1. INTRODUCTION 1 1.1. Scope of Thesis ■ 1 1.2. History of Multi-Level Methods 3 1.3. The Finite Element Approach 6 2. PRELIMINARIES 10 2.1. Elliptic Equations 10 2.2. Finite Element Spaces 14 2.3. Linear Equations 17 3. QUASI-UNIFORM GRIDS 21 3.1. Introduction 21 2 3.2. L Convergence 22 3.3. Computational Cost 30 4. LOCALLY REFINED GRIDS 36 4.1. Introduction 36 4.2. Notation 41 4.3. Algorithms 61 4.4. Complexity 77 4.5. Interpolation 89 4.6. Approximation 103 BIBLIOGRAPHY 116 VITA 118 1 . INTRODUCTION 1.1. Scope of Thesis Multi-level, or multi-grid, methods have existed for more than fifteen years, but have only recently become popular. Thus, there remain many unanswered theoretical questions. Several of these are addressed in this thesis, the major one being the computational complexity of multi-level methods for locally refined finite element grids. These are grids for which the ratio of diameters of the largest and smallest elements is very large, becoming unbounded on the family of discretizations of interest. The principal result of this thesis is that multi-level methods can approximately solve the linear algebraic systems for such finite element discretizations in 0(N) operations for a grid with N elements. Such asymptotically optimal complexity bounds have been given previously for multi-level methods on non-locally-refined grids, but not for locally refined grids. This result provides theoretical justification for the use of multi-level methods on locally refined grids, which is becoming increasingly popular. The proof of this result, which occupies almost half of this thesis, is also of interest for several reasons. Perhaps the most significant of these is that the rate of convergence shown for the multi-level iteration used depends only on local properties of the finite element spaces involved. Global properties, including elliptic regularity, are never invoked. As a consequence this result provides an independent verification of the 0(N) complexity of multi-level methods for non-convex domains, established previously by Bank and Dupont (1978). It also opens the door to easily computable bounds on the convergence rate of multi-level methods on irregular finite element grids. In all, three multi-level algorithms are considered in detail in this thesis. The first, analyzed in chapter three, is similar to methods considered by Nicolaides (1977) and by Bank and Dupont (1978) , but is simpler than these others. It is shown to converge in 0(N) operations for N-point quasi-uniform finite element grids on convex domains, through a convergence proof similar to those given by the above authors. This proof is given here since it is less complex than an earlier one based on Fourier analysis of the truncation error. The algorithm for quasi-uniform grids is treated here partly as background for the analysis of algorithms for locally refined grids in chapter four. It is also of independent interest that an 0(N) bound can be shown for this algorithm even though it is simpler than those in the literature. Two algorithms are considered in chapter four, both applicable to locally refined grids on arbitrary polygonal domains. The first of these is quite similar to the method used in the adaptive finite element code of Bank and Sherman (1978). It is the closest algorithm to theirs that could be analyzed by the theoretical techniques here. This method is shown to converge in 0(N) operations subject to a condition on the rate of growth of the dimensions of the finite element spaces used. The necessary growth condition requires that the finer finite element spaces used by the multi-level algorithm have far more elements than the coarser spaces used. Such a condition is also imposed in the Bank and Sherman code though it is less stringent there due to their technique of "level compression," Bank and Sherman (1978). Since these growth conditions can be quite restrictive, a second algorithm is also considered in chapter four. It is a modification of the first, and it will be shown to converge in 0(N) operations under a considerably weaker growth condition on the dimensions of the finite element spaces used. Though the theoretical interest in this modification is apparent, the author is at. present unsure of the practical utility of this second algorithm. One pays a large penalty for the weaker growth conditions through an increase in the hidden constant in the 0(N) work bound. 1.2. History of Multi-Level Methods Multi-level or multi-grid methods are usually thought of as originating in the work of Fedorenko (1961, 1964), although related ideas had been put forth earlier. The idea of using several discretizations in the solution of a PDE (partial differential equation) is actually quite old. Southwell (1940) suggested using the solution of an elliptic problem on a coarse finite difference grid to obtain a starting value for relaxation on a finer grid. He also recommended a technique called "block relaxation" which is closely related to multi-level methods. In this technique the person doing the relaxation calculation (in the era before computers) would find a region of the finite difference grid where the residuals were all of the same sign. A constant would then be added to the value of the trial solution throughout this region to reduce the average residual on this region as much as possible. As a consequence , after this correction the residuals would oscillate in sign on this region. Southwell observed that as long as there were no large areas where the residual did not alternate in sign, relaxation converged rapidly. Of course, though quite effective, such a heuristic approach was not of much help in the early days of computers, so simpler methods like SOR took over. The idea of making block corrections essentially died out until Fedorenko (1961) observed that rather than making corrections in isolated blocks as suggested by Southwell, one could make a global correction by solving a so-called coarse grid residual problem. For example, suppose that in iteratively solving Poisson's equation on a grid with mesh size h one has at the n-th iteration a trial vector u £ G, and residual — n h r £ G, where G, is the vector space of mesh functions on the grid. — n h h If A is the discrete Laplacian being used, h _h . h r = f - A, u — n — h — n where f^ is the data vector. Now let G~, be the vector space of mesh functions on a coextensive coarser grid with mesh size 2h and suppose we have an interpolation mapping X 2h : G 2h * G h • One can ask for a correction vector v £ G_, such that the improved — zn solution h h x h 2h U J.1 = u + T ou v — n+1 — n 2h — would have smaller residual. This idea has been considered independently by others, Wachspress (1974). It is attractive since v lies in the lower dimensional space G__ and so should be relatively easy to compute. zh 9Vi Fedorenko' s fundamental contribution was to note that v could be chosen as /i n -, \ a 2h T 2h h (1.2.1) A 2h v = I h r^ 9Vi where I : G, -> G_ is an interpolation mapping and A ? , is a discrete Poisson operator on G„, . Since (1.2.1) is just a discretization of Poisson's equation with data I, r , he realized that (1.2.1) could be h — n approximately solved by relaxation and the recursive use of corrections on a still coarser grid G,, . Fedorenko's theoretical work showed that for a model problem multi-level iteration reduced the convergence error by an amount per iteration independent of the mesh size. Thus, reducing the convergence error to a magnitude comparable to the truncation error required 0(N log N) operations on an N point grid, the same as ADI. Shortly thereafter Bakhvalov (1966) considered beginning the solution process on a very coarse grid, using the solution there as a starting value for multi-level iteration on a slightly finer grid, and so on, gradually descending down to the level of refinement desired. This idea, which was probably old when Southwell was young, changes the work bound from 0(N log N) to 0(N). No other current method, either direct or iterative, achieves this asymptotically optimal work bound. Besides improving the asymptotic work bound on multi-level iteration, Bakhvalov extended the method to a large class of PDEs on a square finite difference grid. After his work, almost no work seems to have been done on multi-level methods until this decade. Then several people became interested in this method including Brandt, Hackbusch and Nicolaides. Brandt in particular is largely responsible for popularizing this method and developing it into a practical computational tool. The other two people mentioned have been working more on the theoretical questions involved in this method and on generalization to finite element grids. In the last few years there has been a rapid growth of interest in multi-level methods, and there is considerable ongoing research. This research is directed both at the theoretical problems, Bank and Dupont (1978), Hackbusch (1978), Nicolaides (1977), and at developing algorithms for the many applications areas. Representative work of the latter type may be found in Bank and Sherman (1978), Brandt (1977), Brandt, Dendy and Ruppel (1978), Nicolaides (1975), and Poling (1977). If the author is not mistaken, we should see continued research in this area for many years and may eventually see multi-level methods as the method of choice for most parabolic and elliptic problems. 1.3. The Finite Element Approach Design and analysis of multi-level methods can be done from either a finite difference or finite element point of view. In one sense these two approaches are almost the same; finite element methods are a special class of finite difference method. Nevertheless, the way people think about finite element methods and the kinds of grids they use for these methods are so different from the usual way people understand and use finite difference methods that we may think of these classes as completely separate. This thesis is devoted to analysis of multi-level methods from the finite element point of view. Several of the many reasons the author had for adopting this point of view will be described in this section. While the multi-level method originated as a means of solving finite difference equations and was, until recently, solidly wedded to the finite difference approach, much of the current work on multi-level methods is being done in the finite element context. Starting with Hackbusch (1977) and Nicolaides (1977) , who seem to have been the first to consider multi-level methods as a means of solving finite element equations, there has been increasing interest in this approach. The reasons for this are not complicated. Finite element theory is in many respects superior to finite difference theory. A great many problems in the finite difference approach either do not arise in the finite element approach or are easily handled there. Questions of stability, the use of irregular grids, treatment of corner singularities, treatment of curved boundaries and construction of error estimates are all much simpler in the finite element context. Aside from the obvious (to the author) superiority of the finite element method, there are special reasons unique to multi-level algorithms for doing their theory in the finite element setting. The first of these is that the finite element error estimates are based on approximation results rather than on asymptotic expansions. What this means for multi-level methods is that, unlike some of the early finite difference arguments where rapid convergence of the iteration could be shown only for sufficiently small size, no such proviso is ever needed in the finite element context. A second, more important reason is that multi-level algorithms always involve more than one discretization. Thus one requires mappings from the space of grid functions or finite element functions of each discretization to those of 8 other discretizations. In the finite difference context such mappings are done by interpolations which require separate error analyses. In the finite element context the grids can be designed so that the finite element space on a coarser grid will be a subspace of the finite element spaces on finer grids. Then the mapping from a coarser grid to a finer grid will be just the natural injection. To go the other way, some kind of projection is needed, but many natural ones are available. The fact that different finite element discretizations can be related in this way yields a great simplification of the theory of multi-level methods. Unfortunately only the theory is simplified, not the programming. Rather than an interpolation formula, one has a similar formula relating the basis functions of one finite element space to those of another. It is interesting that in many cases these change of basis formulas arising in the finite element setting are exactly the same as the interpolation formulas already in use in the finite difference setting. There is another reason why studying multi-level convergence is important for this thesis. In chapter four we treat multi-level methods for locally refined grids. While it may be possible to analyze multi-level methods for locally refined finite difference grids in special cases, it seems unlikely that an analysis as general as that given here would be possible outside finite element theory. For locally refined grids, finite element theory also helps in coding. One never needs to choose special formulas for points where the mesh size changes abruptly or worry about truncation error. There are also beginning to be reliable a posteriori error estimates for finite elements enabling one to design robust adaptive algorithms, Babushka and Rheinboldt (1978). In short, multi-level methods, adaptive refinement, and finite element theory mesh together quite nicely. The main limitation of the finite element approach is that it is far more natural and powerful for self-adjoint problems than for the general case. Neither finite difference nor finite element methods work very well for non-self-adjoint problems, but it is frequently easier to see what to do on regular finite difference grids than on general finite element grids. While it would have been possible to extend many of the results of this thesis to the non-self-adjoint case by adding various restrictions, that was a more ambitious project than the author was prepared to attempt. If one takes the results here as indicative of the kind of results that could be proved in the general non-self-adjoint case, one should not go far wrong. 10 2. PRELIMINARIES In this chapter we collect notations and basic results from finite element theory. Since this material is completely standard, the reader familiar with finite elements can probably skip this chapter entirely, referring back to it later as required. The reader less familiar with the finite element approach to partial differential equations may find this chapter helpful. 2.1. Elliptic Equations Throughout this thesis we will be concerned with solving self-adjoint second order boundary-value problems. The problem of most interest will be the Dirichlet problem (2.1.1a) Lu = f on Q (2.1.1b) u = on 8ft , 2 where ft is a bounded open domain in the plane IR with piecewise smooth boundary 9ft. L may be any self-adjoint second order operator, for example the Laplacian, or more generally any operator in divergence form Lu = -V • (a V u) + bu where a, b are C functions on ft. In order for this operator to be elliptic on ft, we require that there be constants a., a, b such that 3. — a — a < b < b hold throughout ft. We will also be considering the Neumann problem (2.1.2a) Lu = f on ft (2.1.2b) u = on 9ft , n where u is the derivative of u in the direction normal to the boundary. n 11 In this case ellipticity requires that b be positive on some neighborhood in ft. Bank and Dupont consider also the singular Neumann problem, b = 0, but we exclude it here. The finite element approach to boundary-value problems is based in part on the functional analytic theory of partial differential equations, in which the boundedness or invertibility of a differential operator is expressed in terms of Sobolev norms. These norms are a 2 natural generalization of the L norm, defined for non-negative integers I by ||u|| = ( Z ||D a u|| 2 2 ) 1/2 , 0<|a|<£ L where a = (a , a ) is a multi-index with a. , a > and j ot | = a + a~ D a = 9 a l a 2 8x 9y This definition makes sense for all functions on ft whose distribution 2 °o derivatives up to order I exist and lie in L . The completion of C (ft) in || • || is denoted H (ft) and is a Hilbert space relative to the inner product (u,v) = / ( Z (D a u)(D°v)) . ft 0<|a|<£ If & is a negative integer || • |L is defined by duality, ||u|| £ = sup - ( ^- -i l|v 'U veH (ft) Here and throughout, (•, •) is the L inner product while || • || with no 2 subscript denotes the L norm. Sobolev spaces can also be defined for arbitrary real Z either through the interpolation theory of Hilbert 12 spaces or by Fourier transforms. Details may be found in Oden and Reddy (1976). Equations (2.1.1) and (2.1.2) implicitly require the solution u to be twice dif ferentiable, which is overly restrictive for many purposes Consider the problem (2.1.3a) Lu = f on Q (2.1.3b) u = on T (2.1.3c) u = on I\ , n 2 where 9S7 = V U r , which includes both (2.1.1) and (2.1.2) as special cases. Let H (fi) be the subspace of H (Q) satisfying the (essential) Dirichlet boundary condition (2.1.3b) but otherwise unrestricted. If (2.1.3) is satisfied then (Lu,v) = (f,v), v G Hg(fi) . The notation here means that equality holds for all test functions v in H (fi) . Defining the bilinear form a(u,v) = (Lu,v) , we can ask for the element u in H (fi) such that E (2.1.4) a(u,v) = (f,v), v6Hj(i]) . It is only necessary that u lie in H since Green's theorem can be used E to shift one derivative from u to v. A function u satisfying (2.1.4) is called a weak solution of (2.1.3), and (2.1.4) is called the weak form of the equation. The advantage here is that the weak solution exists quite generally and agrees with the classical solution satisfying equation (2.1.3) pointwise when the latter exists. 13 The condition for existence and uniqueness of the solution of (2.1.4) is that the problem be elliptic. This means that the energy norm ||| • ||| defined by |||u||| 2 = a(u,u) for u E H (ft) is equivalent to the first Sobolev norm || • || . That is, Jj 1 ffjM't < INII < ^H . u ^ H^ft) , for fixed , G~ > independent of u. If the problem (2.1.3) is elliptic, a simple calculation shows that for & a non-negative integer l|u|l £+2 > c ||Lu|| £ , u e H £+2 (ft) , for some c > 0. The principal results in elliptic theory go the other 3 I way. For I > - -~ and any f e H (ft) the weak solutions of (2.1.1) or (2.1.2) satisfy the regularity estimate (2.1.5) NI £+2 < C H f ll*, for some c > 0, assuming that the domain is smooth and the coefficients 00 of L are C . For the technical definition of smooth domains see Babushka and Aziz (1972) or Oden and Reddy (1976). The regularity estimate (2.1.5) implies existence, uniqueness, and that the solution has two more derivatives in the mean square sense than its data. In particular, if f is infinitely dif f erentiable, u will be as well. These results hold only for smooth domains. For domains that are only piecewise smooth the corners introduce singularities that limit the regularity of the problem. On a convex polygonal domain the 2 I solution will lie in H but not necessarily in H for £ > 2. In this case we have the regularity estimate 14 (2-1.6) IMU^IIf"* 3 for - "z- < £ _< 0. For non-convex polygonal domains (2.1.6) will not hold for I = because the reentrant corners, that is corners with interior angles exceeding tt, introduce severe singularities. In the worst case of 3 1 a slit domain, £ must be restricted to (- -~, - ■=■). The solution has less than one and a half derivatives in the mean square sense. 2.2. Finite Element Spaces The finite element discretization of (2.1.3) begins with the selection of a finite dimensional subspace M of H (ft). The discrete E solution to (2.1.3) is taken as that function u E M satisfying the finite element equation (2.2.1) a(u,) = (f,), J . The approximation properties of each of these spaces M. depend on the size of the elements {e} of M. and the element geometry, as well as the degree of the finite element space. Element geometry matters only to the extent that it must not degenerate as j increases. For standard finite elements it suffices that all elements In M. be convex and have 3 interior angles bounded away from and IT uniformly in j . The size of an element e is measured by its diameter, d , defined as d e = sup ||x-y|| 2 , x,yee I where ||»|| „ denotes the Euclidean norm. The family of spaces {M.} is I J said to be quasi-uniform if there is a corresponding family of parameters {h.} and a > 1 such that (2.2.2) - h. < d < ah. O 3 - e - j for all elements e of M. and all j. If not, that is if the constant O 16 must depend on j, the family of spaces {M } is said to be locally refined, While strictly speaking these concepts apply only to an infinite family of finite element spaces, we will talk about an individual finite element space being quasi-uniform or locally refined. This usage is extremely convenient though it should be born in mind that this really makes sense only with respect to membership in an infinite family of spaces {M. }. This definition of the mesh parameter h. associated with a quasi-uniform finite element space M. is not completely standard. Often h. is taken as the maximum of the element diameters. However, the J extra freedom allowed by (2.2.2) makes no difference in the error estimates since a is a fixed constant, and is quite convenient. In particular this freedom makes it easy to have a quasi-uniform family ■f^ili /n 1N parameterized by h rather than j. In this case the diameters tl n£ (,U, L) of elements in M, are required to satisfy — h < d < Gh . a — e — The standard finite element error estimates give the rate of convergence in any Sobolev norm of the finite element approximations in a quasi-uniform space M in terms of the mesh parameter h and the degree of M . For smooth solutions the convergence rate may be made arbitrarily high by taking finite element spaces {M, } of sufficiently high degree. But on a polygonal domain the corner singularities limit the rate of convergence attainable by any family of piecewise polynomial finite element spaces. The most important error estimate in the next chapter will be the L error estimate 17 (2.2.3) ||u-u h || 2 < Ch 2 ||f|| 2 , L L which holds for quasi-uniform spaces {M, } of degree at least one when h 2 the partial differential equation is H regular, that is, when (2.1.6) holds with I = 0. Thus, this bound applies to the Dirichlet problem (2.1.1) and the Neumann problem (2.1.2) on any convex polygonal domain. The convergence in energy for these two problems is of first order on any convex polygonal domain. On non-convex domains it is of less than first order. For a domain with a slit, the worst case ordinarily arising, it is of order one-half. Much of the utility of the finite element method comes from the fact that these meager convergence rates attained on polygonal domains are easily improved through the addition of singular functions or by doing local refinement to cope with the corner singularities. In chapter four the use of multi-level methods on locally refined finite element grids will be considered. A simple error estimate showing the quality of approximation of finite element methods on such grids will be given there. 2.3. Linear Equations By selecting a basis for the finite element space M, the finite element equation (2.2.1) is translated into a linear system. Let M be a finite element space for problem (2.1.3) and let l

. , i=l - 1 X 18 where u^ is the coordinate vector of u relative to this basis. The convention of using underbars for vectors will be adhered to throughout. Given the basis {<}>.}, we can define the mass and stiffness 1 matrices M and K with elements m ij = C* ± . Oj) k ij = a( V V respectively. M is just the Gram matrix of the basis relative to the 2 L inner product while K is the Gram matrix in the energy inner product. To form the linear system corresponding to the Ritz equation (2.2.1), let f e L (fi) and let F be given by F ± = (f, (JO, i = 1,..., N . Then (2.2.1) is equivalent to the linear system (2.3.1) Ku = F where jl is the coordinate vector of the finite element solution. Numerical solution of the linear system (2.3.1) is aided by several properties the matrices K and M have. First, M will always be symmetric positive definite and since L is self-adjoint K will be symmetric positive definite as well. Second, for the usual finite element bases M and K will be quite sparse. Finally, both M and K are usually well conditioned. The common finite element bases are largely motivated by these last two considerations which are of great importance for the numerical inversion of (2.3.1) or (2.3.2) whether this is done directly or iteratively. To make these considerations precise, consider a family of finite element spaces {M . } -.of increasing dimension. We require that J C) N - each space M., 1 < i < °°, have a basis {(J). }. , that is local in the j — i i=l 19 sense that the number of non-zero elements in each row of the corresponding mass and stiffness matrices is bounded uniformly in j . This will be true if the number of basis functions in {<£>. }._, whose supports intersect the supports of any given basis element <$. is bounded uniformly in j . Usually, under these conditions the condition number k(M.) of the mass matrix M. , 1 j< j < °°, will be uniformly bounded. For this to hold the proper scaling must be chosen and the condition numbers of the element mass matrices must be uniformly bounded. This is true for the common finite element bases as long as the element geometry does not degenerate. Then if the basis elements are scaled as ||^ j) || = 1, 1 < i < N. , 1 < j < °° , which we assume throughout, it is immediate that the eigenvalues of M. will lie in an interval [— , o] for O > 1 independent of j. Since the stiffness matrices {K. } are exactly analogous to the matrices arising in finite difference discretizations, their condition numbers cannot be uniformly bounded. Just as in the finite difference case the maximum eigenvalue of K. satisfies (2.3.2) A (K.) < ch" 2 max j — mm under ordinary circumstances, where h . is the diameter of the mm smallest element of M. . The minimum eigenvalue of K. can be bounded 3 3 below by a simple argument. We have T T T x K.x_ x K.x. x K.x A (K.) = min — ^— > min J min — ^~ J xx x M.x x x _ -,_ = Am: 1 K.) XAK.) 1 3 J 1 3 > A. (MT 1 K.)/a . ~ 1 J 3 20 The operator M K is the linear operator corresponding to the finite element equation (2.2.1). Since the finite element eigenvalues approximate the eigenvalues of the PDE from above, A (K.) must be bounded below uniformly in j . Thus one ordinarily has (2.3.3) k(K.) < ch" 2 2 - mm or for a quasi-uniform spaces JVL . (2.3.4) k(K.) < ch" 2 . The difficulty in iteratively inverting the linear system (2.3.1) arises directly from this growth of the condition of the stiffness matrix as h decreases. These bounds on the condition of M and K are proved in Fried (1971) and may be found also in Strang and Fix (1973). 21 3. QUASI-UNIFORM GRIDS In this chapter a multi-level algorithm simpler than those considered in previous theoretical work will be given. It will be shown to be of optimal order, producing a good solution on a quasi-uniform finite element grid in 0(N) operations on an N point grid. 3.1. Introduction In this chapter we look at a multi-level algorithm for self-adjoint elliptic boundary value problems on quasi-uniform finite element grids. This algorithm is related to algorithms considered by Nicolaides, Bank and Dupont , Hackbusch and others, but is less complex than these other algorithms. Like many of those considered by others, the algorithm here is of optimal order, producing a second order accurate solution in 0(N) operations on an N point grid. The optimal order convergence of the algorithm in this chapter was originally shown for finite difference discretizations by a cumbersome Fourier technique, and then for finite element discretizations in the same way. Here an easier proof, following the ideas of the above authors, will be given. The proof is based particularly on the work of Bank and Dupont but employs several minor simplifications. Throughout we follow, for the most part, their notation. It is convenient here since little reference to algebraic quantities such as vectors or matrices will be made. 22 2 3.2. L Convergence In this section and the next we examine the convergence of a simple multi-level algorithm. This algorithm can be used to obtain finite element approximations to the solutions of the elliptic problems (2.1.1) and (2.1.2) on a convex polygonal domain ft. The restriction to convex domains is necessary since we need to use the 2 2 L error estimate (2.2.3) which requires H regularity. Using curved elements as in Bank and Dupont [1978] the convergence results given here extend immediately to piecewise smooth domains free of reentrant corners. For details on the treatment of curved elements, the interested reader is referred to their paper. The algorithm considered here uses a nested family of quasi-uniform finite element spaces iM.}.-, nested in the sense that M. C M. ,. , 1 < j < °o . 3 3+1 To approximately solve the finite element equations for one of these spaces M , k _> 2, this algorithm uses a simple simultaneous displacement smoothing iteration applied to the equations for M, plus approximate solution of related problems in M, - . If k - 1 > 2 these approximate solutions are obtained recursively by applying the same algorithm to the equations for M, ,. The recursion continues, finally requiring approximate solutions to the equations for M . These can be found directly or by some convenient iteration. The multi-level method considered here is related to algorithms for finite elements by Nicolaides, Hackbusch, and Bank and Dupont and to methods for finite difference equations considered by Brandt, Hackbusch, Nicolaides, Bakhvalov, Fedorenko and others. The algorithm here differs 23 from algorithms described by others in that the solution to the equations for M, requires only two approximate solutions to the equations for M, _, • In other multi-level algorithms, to solve the equations for M, approximate solution of related problems in M, , is used in two ways. First it is used as a means of obtaining a good starting value for iteration on M . Subsequently it is used as an acceleration procedure for iteration on M, . k The first approximate solution of the equations for M, - is used to obtain a starting value for iteration on M, , and considerable accuracy is required. It is necessary to reduce the convergence error of the problem on M so it is comparable to the truncation error on M, n . In subsequent solutions of related problems in M, , used to k-1 r k-1 accelerate convergence of the iteration in M less accuracy is required. K. For example, in one of the proofs given by Bank and Dupont and in the related proof of theorem 4.1 in chapter four of this thesis it suffices to solve the related problem with a relative accuracy of 33%. The point of this chapter is that if all of the related problems in M are solved to the same accuracy, namely the error in the approximate solutions must be comparable to the truncation error in M, . , then only two approximate solutions in M are required. To approximately solve the elliptic problem (2.1.1) or (2.1.2) in M, , k _> 2, the algorithm begins by obtaining a good approximate solution to the finite element equations for M . This approximate solution is used as a starting value for a smoothing iteration on M . K. Though this iteration could be continued until adequate convergence was attained this would be inefficient. Instead, after a fixed number of iterations, m, independent of k, we terminate the iteration. 24 The resulting approximate solution is then corrected by subtracting from it the projection of the convergence error onto M . This corrected approximate solution in M is then taken as the accepted solution. The reason this is effective goes back to the work of Southwell. He observed that at least initially relaxation rapidly reduces the residual, while the error itself may be only slowly reduced. At the end of our algorithm's smoothing iteration the residual will be small relative to the convergence error. This means that the convergence error is quite smooth. In this case the convergence error can be well approximated in the coarser space M and the final correction will remove most of it. Brandt describes this same phenomenon in terms of the reduction of Fourier components of the convergence error. Most convergence proofs, including ours, follow Fedorenko who considered the reduction of the error components in the direction of eigenfunctions of the discrete partial differential operator. Fedorenko refers to "good" meaning smooth eigenfunctions and "bad" meaning oscillatory eigenfunctions, For practical purposes many possible smoothing iterations would be acceptable. Jacobi, Gauss-Seidel, and symmetric SOR are most often used. The author expects that improvement could be made through the use of Chebyshev acceleration, although in numerical tests so far changing to more subtle iterative schemes has had little effect for the well behaved self-adjoint problems considered here (Bank and Sherman [1979], Brandt [1977]). In this section the smoothing iteration will be taken as a simple but computationally unattractive simultaneous displacement iteration considered also by Bank and Dupont. The advantage of this 25 iteration is that the resulting convergence is easy to analyze and is basis independent. In the next section we will show that this smoothing iteration can be replaced by a more computationally attractive iteration without altering the convergence result. In this case the algorithm will be shown to be of optimal order, producing a second order accurate solution in 0(N) operations on a finite element grid with N elements. Let (M . } . , be a nested family of quasi-uniform finite J 3 = 1 element spaces, and suppose we have associated mesh parameters ih,} satisfying (3.2.1) h. = p 1 " 3 !^ , j > 1 , for p > 1 a fixed constant. Let u G M,, j > 1, be corresponding finite element solutions of problem (2.1.1) or (2.1.2). That is, u G M. are functions satisfying (3.2.2) a(u (j) ,(b) = (f, e M . From the error estimate (2.2.3) (3.2.3) ||u-u (j) ||< chJllfH for fixed c > where u is the true solution. In the space M. associated 3 with the finite element equation (3.2.2) we have eigenf unctions (') N * {i> . } ^ , where N. is the dimension of M., and associated eigenvalues i i^l 3 J U^}. 3 with i i=l < A< j) < A< j) <...< X ( T j) . 1 — 2 — — N. 3 That is 26 a(^ j) ,<}>) ■ A< j) (^ j) ,4>) , 1 independent of j, we have the bound (3.2.4) A X7 < c hT 2 . To solve the finite element equation in M, , k >^ 2, the algorithm to be considered produces a sequence of iterates u , u 9 ,..., u in M beginning with an approximate solution u_ to the finite nrrJL k. U element equations in M, - . For some fixed positive integer m, the approximate solution u in will be taken as the final answer. The m+1 iterates u. , 1 < I _< m, are produced by a simple smoothing iteration while the last one, which will be the accepted solution, is the result of a correction involving the approximate solution of a related 2 problem in M, , . For a fixed constant c_ > and data f in L (Q) the k-1 algorithm is as follows: 1. Let u G M be an approximate solution to (3.2.2) satisfying (k-1) I, . u -u < (3.2.5) HV U Hi c (h k-l } ll f ll (k-1) where u is the solution of (3.2.2). 2. Let u„ , 1 <_ £ <_ m be the element of M satisfying (3.2.6) (u £ -u £ _ lf ) , G M . m k-1 27 That is , v is given by a(v,cf>) = -(r,(f>) , £ M k _ 1 where, r £ M, is defined by k (r,<|>) = (f ,<))) - a(u ,({)), f£M, . in k. Let v be an approximation to v satisfying (3.2.7) ||v-v|| < c Q h^ ||r|| and let u , .. be given by m+I (3.2.8) u .. = u - v m+1 m The corrected approximate solution u , . is taken as the final answer. m+1 If we can show that the approximate solution, u . , produced is accurate enough uniformly in k >^ 2, the use of this algorithm recursively to solve the subproblems in steps one and three will be immediately justified. This is the content of the following theorem. Theorem 3.1. For any c n > there is a positive integer m, independent of k, such that the approximate solution u ., produced by m+1 algorithm (3. 2. 5)-(3. 2. 8) satisfies (3.2.9) H« 1|H -. ll. i=l N It remains to estimate the result of the correction in step 3, We have N A. r = E a. A. (1 - -r— ) \b , .-IX A T l 1=1 N and by the error estimate (2.2.3) llv-eJIich^Hrll; so by (3.2.7) Vill-HVlli »•.-'« + »«i -2-2- ( c + c_) (c(l + p 2 ) + c. p 2 ) , - c Q and the result follows. □ In this proof we have carried the constants p and c explicitly to show the dependence of m on these quantities. From (3.2.11) m 4 apparently goes as — . This sharp dependence of the number of iterations C m on the mesh ratio is typical of multi-level algorithms, but the dependence on c is not. For more common multi-level algorithms, such as that considered by Nicolaides [1977] and that considered in chapter 4 for locally refined grids, m depends on c only as |log(c n )|. Thus if one wishes to reduce the convergence error well below the truncation error, the algorithm here cannot be recommended. However, in the usual case where it suffices to have convergence error comparable to the truncation error this algorithm should perform about as well as other multi-level methods. 30 3.3. Computational Cost While iteration (3.2.6) used in step two of algorithm (3.2.5)- (3.2.8) was convenient in the last section since it made the analysis there basis independent, it is impractical. Each of the iterations, (3.2.6), required the solution of a linear system involving the mass matrix because of the implicit definition on the left side. To avoid this, we will show in this section that the smoothing iteration (3.2.6) can be replaced by a simple simultaneous displacement iteration without impairing the convergence result. This analysis follows Bank and Dupont [1978]. With this replacement the overall algorithm is of optimal order, producing a second order accurate solution in 0(N) operations for a finite element grid with N elements. Given a basis {(J).}. , for M, , for each u £ M, there corresponds l i=l k k N a coordinate vector u G |R . As in the last section, we have temporarily dropped the subscript k on the dimension N of M and on eigenvectors and eigenvalues. Let b(», •) denote the bilinear form T b (u,v) = u. v The norm induced by this bilinear form is the same as the Euclidean norm on corresponding vectors, 11-11? - b(u,u) = ||u|| 2 . I Now let {_$.} be the eigenvectors of the stiffness matrix K with corresponding eigenvalues 0< X < x «...« x N The eigenvectors {^.} are coordinate vectors for eigenfunctions {i|>.} satisfying 31 (3.3.1) a(ifi., $) = A. b(^., ) , G M XXX K. and the orthogonality relations a($ , $.) = bC^, J.) - for i ^ j . We wish to consider algorithm (3. 2.5)-(3. 2, 8) of the last section but with the smoothing iteration (3.2.6) replaced by the simultaneous displacement iteration (3.3.2) b(u £ - u^_ 1 , ((.) = -±- (a(u £ _ 1 , cf>) - (f, $)) , there is a positive integer (k) m independent of k > 2 such that the approximate solution u * — m+1 produced by algorithms (3. 2. 5 (- (3. 2. 8) with (3.2.6) replaced by (3.2.2) satisfies ll Vl -u«||) = - a(e ,<}>), f £ M, m k b(r,<}>) = - a(e m , ) , ♦ G M • r = Mr ; 2 T — 1 T 1 rll = r Mr = (M r) M(M f) ~T -1. = r M r ««llllP 2 -«ll«l6 As in the proof of Theorem 3.1 Now Vi'l- (c(1 + p2) + c o p2) h k l|r|1 < (C(l + p 2 ) + C Q p 2 )a 1/2 h 2 Hr^ 1=1 A N <||eJL max (X(l - ^-) m ) 0 - N-qII 2 = ^ He l! b A/ A/ and since the maximum eigenvalue L, satisfies the bound (2.3.2) the rest N of the proof goes through exactly as before. □ This theorem shows that by twice approximately solving the lower dimensional linear system for M .. and doing a fixed number of simultaneous displacement iterations on M. one can obtain a second order accurate k solution in M . Since the convergence shown by Theorem 3.2 is independent K. of k >_ 2, the two related problems in M can be solved by recursively applying the same algorithm to their smaller linear systems. Continuing the recursion one ends up repeatedly having to solve the linear system for JVL . This can be done either directly or by any convenient iterative method. To estimate the computational cost of this algorithm observe that from (3.2.1) the dimensions of the finite element spaces {M.} must satisfy C, p ** 1 N 1 c 2 p J , 1 £ j < °° , for c , c_ > independent of j. This assumes that h , the mesh size of M , is fixed and the degree of the highest degree polynomial on any element is bounded uniformly in j. Then, if s_ is the storage required for the direct or iterative inversion of the linear system for M , the total storage must be proportional to k k s + Z N. < s n + c E p 3 . 3 — 2 . - 3=2 J J=2 2k < s„ + c, - ° 2 i-p- 2 34 for p > 1. Since s is fixed, assuming the mesh ratio p is greater than one, the total storage required depends only linearly on the dimension of the solution space M. . k To bound the computation time, observe that the time required to do one of the simultaneous displacement iterations (3.3.2) for M. is proportional to N.. Let t be the computation time for the approximate solution of the linear system for M . Then taking into account that we must recursively apply the algorithm twice to M, - , four times to M ,_ 9 » and so on the computation time T is bounded by V k lr ' (3.3.4) T < 2 t A + cm I 2 J N. " ° 3=2 3 < 2 k t + cc m E 2 k " j p 2k j-2 ^k 2k „ , 2 vk--j £ 2 t + cc 2 mp E (-j) J . j=2 P We have neglected the computational cost of carrying out the mappings of functions on one grid to another in steps one and three of algorithm (3. 2. 5)- (3. 2. 8) . Despite the fact that the spaces {M.} are nested, the mappings required in steps one and three are not free since in general each space M. uses a different basis. It is however easy to see that the above bound on the computation time is unaltered by taking into account this additional cost. Now from (3.3.4) we have the bounds 0(N, ) for p > /2 k T < 0(N log(N )) for p = /2 log 2 . N/ XOg p for 1 < p < /2 . ^ k 35 The natural choice p = 2 gives an optimal order algorithm and is convenient in programming. 36 4. LOCALLY REFINED GRIDS This chapter, which contains the principal results of this thesis, presents two multi-level algorithms achieving 0(N) complexity bounds on locally refined grids. One of these attains this optimal complexity bound under weaker than expected restrictions on the dimensions of the finite element spaces used by the method. These optimal order convergence results are a consequence of an approximation theorem proved in this chapter. Aside from its use in extending multi- level convergence theory to locally refined grids, this approximation result is of interest since it relies only on local properties of the finite element spaces used and is not based on elliptic regularity. As a consequence it provides an independent verification of the 0(N) convergence of multi-level methods on non-convex domains shown previously by Bank and Dupont . 4.1. Introduction There are a variety of reasons why one might wish to use locally refined grids. They can be used to maintain high order accuracy in the face of singularities caused by discontinuous coefficients or the corners of the domain, or to resolve boundary layers. They may simply result from an adaptive discretization. They can also be used to give good resolution of the part of the solution of greatest interest without incurring the computational cost of a fine global grid. This situation occurs, for example, on unbounded domains. As long as the ratio of the diameter of the largest element to that of the smallest remains uniformly bounded on the family of finite element spaces of interest, multi-level convergence results for quasi-uniform grids, such as those in the last 37 chapter, remain valid. However, in many cases of interest this ratio of element diameters becomes unbounded on the family of finite element spaces being used. In this case, while results for quasi-uniform spaces, such as those in the last chapter, yield an 0(N) work bound for each finite element space, the constant hidden in the 0(N) work bound will depend on the ratio of element diameters. Such bounds are clearly unsatisfactory for locally refined finite element spaces. While nearly all of the convergence results in the literature apply only to the case of quasi-uniform finite difference or finite element grids, multi-level methods are being used increasingly for locally refined grids, Brandt (1977), Bank and Sherman (1978). In fact, one of their principal advantages is that their rate of convergence remains very good on locally refined grids. By contrast other iterative methods, whose rates of convergence depend on the condition number of the stiffness matrix, do poorly on locally refined grids. For example, the number of iterations required with conjugate gradient iteration goes as /K(K), where K(K) is the spectral condition number of the stiffness matrix, Axelsson (1977). This condition number satisfies K(K) > c h" 2 — mm where h . is proportional to the diameter of the smallest element in the mm grid. Thus, K(K) may grow much faster as a function of the number of unknowns on a locally refined grid. This presents a very severe difficulty for nearly all iterative methods. As will be shown in this chapter, multi-level methods do not suffer from this difficulty. The algorithms in this chapter differ considerably from that in the last chapter, aside from their application to locally refined grids. In the algorithm of the last chapter, to approximately solve the elliptic 38 problem on the finest grid, a related problem on the next coarser grid was solved twice. It was solved once at the beginning to obtain a starting value for the simultaneous displacement iteration and once at the end as a final correction. In the algorithms of this chapter and in other algorithms in the literature, solution of a related problem on the next coarser grid is used more frequently, typically after every few smoothing iterations on the finest grid. Thus, the basic operation in such an algorithm is best thought of as a multi-level iteration consisting of a simple relaxation iteration periodically accelerated by the solution of a related problem on the next coarser grid. Algorithms like those considered in this chapter have been described for quasi-uniform grids by Bank and Dupont, Nicolaides, Hackbusch, and others. They are also related to methods for locally refined grids suggested by Brandt and very closely to the method used in the adaptive finite element code of Bank and Sherman. An algorithm similar to that of chapter three could have been used, but probably would not have been as practical as the algorithms here. One could argue that in any case it is best to consider different modifications of the standard approach separately rather than in combination. To our knowledge there is only one convergence result in the literature for multi-level methods on locally refined grids. This is given in Bank and Dupont (1978) where a two-level multi-level iteration is shown to converge at a rate independent of the mesh size. Since this result relies only on local properties of the finite element space and not on global approximation properties, it applies automatically to locally refined grids. The limitation to this result is that it deals only with their two-level scheme. With such a scheme, as one decreases the mesh size, the dimensions of both the finer and coarser finite 39 element spaces grow. In their algorithm the linear system for the coarser finite element space must be approximately solved once per multi-level iteration. Whether this is done directly or by one of the common iterative methods, the cost of approximately solving this linear system grows faster than linearly in the dimension of the linear system. As a consequence, such a two-level scheme can never be of optimal order. It may, however, be quite attractive as a practical algorithm since inversion of the lower dimensional linear system for the coarser space may be much cheaper than the cost of inverting the linear system for the finer space in which the solution is sought. The use of direct methods to solve these linear systems for the coarser space is particularly attractive since after the first multi-level iteration only back solves are required. One might expect that a simple generalization of their convergence proof for the two-level case could be used to show optimal order convergence of a multi-level method using many levels on locally refined grids. In some cases it can, for example, for Poisson's equation on certain finite element grids. In general, however, the error reduction per iteration, y E (0, 1), of their two-level scheme, while independent of the mesh size, depends strongly on the PDE and on certain properties of the discretization. The inductive argument used to show the optimal order convergence of multi-level methods using many levels requires that we be able to choose a sufficiently small value of y. For example, for the first algorithm in this chapter, y must be chosen in (0, — ) . Since there is no easy way of controlling y in their two-level scheme, it seems to be quite hard to get general optimal order convergence 40 results from a direct generalization of their two-level result. This is unfortunate, since the argument showing the convergence of their two-level algorithm is far simpler than that given here. The results here and the technique of their proof are actually much more closely related to the second kind of algorithm considered in their paper. This second algorithm is a recursive algorithm using many grid levels just as the algorithm in the last chapter did. Their algorithm is shown to converge in 0(N) operations on quasi-uniform finite element spaces. Unlike the results in chapter three, their results apply also to non-convex domains, on which the elliptic problem 2 is less than H regular. The analysis of algorithms for locally refined grids here, which also applies to non-convex domains, is largely modeled on their analysis of the second algorithm in their paper. The notations here have also been kept almost the same as theirs. Analysis of multi-level methods is harder than that of most iterative methods since it involves both algebraic and approximation questions. The algebraic and complexity questions for locally refined grids are slightly harder, since the dimensions of the finite element spaces used in the multi-level algorithm must be taken into account more carefully. By contrast, the approximation questions for locally refined grids seem to be a great deal harder. In the first half of this chapter we consider the algebraic and complexity questions involved in treating multi-level methods on locally refined grids. The analysis here largely follows that in the paper by Bank and Dupont, while the algorithms themselves are motivated by the adaptive finite element code of Bank and Sherman. Two algorithms are considered; the first is similar to that in the Bank and Sherman code. The second is a modification yielding an improved complexity bound. This modification, which appears to be 41 new, allows greater freedom in the choice of finite element spaces used and may have significant practical applications. The second half of this chapter deals with the approximation questions involved in treating locally refined grids. The approximation result proved there is the main result of this thesis. In previous work on multi-level methods, the approximation questions were answered by using the elliptic regularity of the problem to show smoothness of the functions being approximated. Then the approximation error could be estimated using standard finite element techniques including duality arguments and the Bramble-Hilbert lemma. The approximation result here uses a completely different kind of argument making no use of the regularity of the problem. This is why it applies to non-convex 2 domains where the regularity is weaker than H . Perhaps of greater interest from a computational point of view is the fact that this approx- imation result relies only on local properties of the finite element spaces and the PDE operator. In consequence, it would be quite easy to compute rigorous upper bounds on the rate of convergence of the algorithms in this chapter. Up until now there were good heuristic bounds on the rate of convergence of multi-level methods on finite difference grids, Brandt (1977), and in a few cases rigorous bounds. No previous bounds for the rate of convergence of multi-level methods on irregular finite element grids appear to be available. 4.2. Notation To develop the algorithms considered in this chapter, it is necessary to go into greater detail concerning the finite element spaces 42 used than was required in the last chapter. Additional notation must also be developed. These are the tasks of this section. To keep the mathematical issues in this chapter clear, our analysis will be restricted to C triangular and rectangular, elements. Though the analysis in this chapter seems to be extendable to all or nearly all families of finite elements, the details of the proofs here can be considerably harder for some families of elements. Using the C square and triangular elements considered here permits considerable simplifi- cation of the analysis and may also cover the cases of greatest interest, since the majority of finite element computation seems to be done with these elements. Linear C triangular elements, called Turner triangles, are used in the Bank and Sherman adaptive multi-level code, so at least for polygonal domains, the solution space used by their code is covered by the theory here. However, as will be seen, their algorithm violates the assumptions of the theory developed here in a number of other ways. There are really two reasons why more detailed specification of the finite element spaces used in the multi-level algorithm is needed here than was required in the last chapter. The first is that locally refined spaces are intrinsically more complex than quasi-uniform spaces since they are not characterized by a single mesh parameter. The second is that, unlike the analysis in the last chapter, which was based on standard finite element error estimates, the analysis here is based on an approximation result developed here. The proof of this result depends on detailed analysis of the finite element spaces used. While this result could have been proven more abstractly, the approach here is certainly more straightforward. The reader with a background 43 in finite element theory should have little trouble in seeing the generalizations to elements with higher order continuity. The two elements considered here are the standard C triangular and square elements. The C triangular elements of degree k, k _> 1, use as element trial space the space P of polynomials in x and y K. of degree at most k. The C square elements of degree k, k >_ 1, are tensor product elements based on the trial space (L consisting of poly- nomials whose terms are of degree k or less in x and y separately. Both types of elements may be taken as nodal finite elements with node placement as in figures 1 and 2. K= K=2 K = 3 K=* Figure 1. C Lagrangian Triangular Elements • • (i • • « • • • • K= K=2 K=3 Figure 2. C Tensor Product Elements 44 The nodal parameters used are just the function values at these points, providing a substantial simplification of the analysis in this chapter over that which would be required for the more complex Hermite elements. Equally important from our point of view, for any k _> 1 , either type of element constitutes an affine family, Ciarlet (1979). What this means is that we may fix a reference element in the plane, on which we have a reference trial space L, which will be either V or () for triangular K. K. or square elements, respectively. For any other element of the same type anywhere in the plane we can find an affine transformation from the reference element to the given element. This affine transformation carries the region occupied by the given element, the nodes, and the element trial space from the reference element to any other element of the same type anywhere in the plane. Affine transformations are mappings of the form x 1 = Ax + b where b is a fixed coordinate vector and A is a nonsingular 2><2 matrix. For triangular elements all affine transformations may be used, carrying the reference element into the entire six parameter family of triangles in the plane. For square elements, A may be taken as a multiple of the identity. That is, there is no need for rotation and distortion in this case. The point of all this is that affine transformations and the use of reference elements provide an elegant simplification of the theory. The reference element and its element trial space basically contain all properties intrinsic to this type of finite element. On the other hand, the affine transformations contain all information about the geometry, mesh size, and distortion of any particular element in the 45 finite element grid. Thus, analysis of a finite element space is conveniently factored into two parts, a property we will use to advantage in section 4.5. Aside from the choice of element, two types of grids need to be considered. These may be conveniently labeled "aligned" and "unaligned" grids. Aligned grids are those which are standard in finite element theory, where each edge of an element is the edge of another, or part of the boundary. Unaligned grids are those where the edge of an element may be only part of the edge of another. See figures 3, 4, and 5 for examples of these two types of grids. In general, triangular grids may always be chosen to be aligned, without sacrificing efficiency, but it is often necessary to choose grids of square or rectangular elements to be unaligned, or much of the benefit of local refinement will be lost. To be admissible, the finite element spaces must be C and satisfy the essential, that is Dirichlet, boundary conditions. On aligned grids there is no problem obtaining C continuity. One simply requires agreement of the nodal parameter on the edge between adjacent elements. This does not work for unaligned grids. To simplify this problem we restrict the unaligned grids as follows. Assume that every edge of an element is either part of the boundary, the edge of another element, or the union of the edges of two other elements. This last situation is shown in figure 6. In this last situation we consider the single element to dominate the two elements it adjoins across this edge. For convenience call the single element the "larger" since it almost always is larger than the other two. Then we may say that the "larger" dominates the 46 Figure 3. Aligned Triangular Grid El Figure 4. Unaligned Triangular Grid Figure 5. Unaligned Grid with Square Elements P' ii • i P< n • i Figure 6. Adjacent, Unaligned Elements 47 other two in the sense that its nodes on this edge are the real nodes of the finite element space there. The nodes of the "smaller" elements may be called parasitic since their values will be taken from the larger element. Parasitic nodes of the "smaller" elements in figure 6 not coinciding with nodes of the "larger" element are labeled with a "P." It is easy to see that this system preserves the required C continuity and is fairly straightforward to code. One more item needs to be discussed before we consider construction of the locally refined spaces needed in this chapter. This is the angle condition. It was mentioned in section 2.3 that the mass matrix will usually be well conditioned uniformly in the level of refinement of the grid as long as the element geometry does not degenerate. While this is not always true for locally refined grids, it is still important that the element geometry not degenerate. One way of stating the required condition is to say that the minimum interior angle of any triangle in the grid must be bounded below uniformly in the level of refinement. More formally, if {T.j. , is the family of 3 3=1 triangulations being considered, we want (4.2.1) inf { min } > j>l T^T. J where 6 is the smallest interior angle of the triangle T. There are several alternate formulations of this condition in the finite element literature, Oden and Reddy (1976), but this one is adequate for our purposes. Obviously, no such condition is needed on square grids. We next turn to consideration of the class of finite element grids on which the multi-level algorithms of this chapter can be shown to converge rapidly. This class, which includes both locally refined and 48 quasi-uniform spaces, is quite large and may be thought of as an indexed family {M } _ of spaces. It consists of all spaces which can be a oea generated by the refinement process described below and which satisfy an additional condition given in section 4.4. An 0(N) complexity bound holds uniformly on the class of spaces in {M )__. satisfying this Ot OfcA condition uniformly. In general, this class of spaces with a uniform 0(N) complexity bound is large enough to show the improved convergence to singular solutions obtainable by discretization with locally refined grids. A partial exception to this is the case of point singularities, since the best discretizations for these problems do not satisfy the conditions of section 4.4 uniformly. For these problems 0(N) bounds can be shown by the theory here only for suboptimal, but still very good locally refined grids. Fixing a coarse space M , the refinement process described below will generate sequences {M „ . } . , , 3 e B of finite element spaces, 3,1 3=1 where the index 3 ranges over the uncountable family of choices allowed in this refinement process. The family of spaces {M } is just the union of all spaces which can be produced by this refinement process, The condition mentioned above, which is described more fully in section 4.4, is that for each fixed 3 G B, the dimensions of the spaces {M ./. . 3,3 3=1 in this sequence of refinements grow geometrically. The 0(N) complexity bound given here depends only on the problem, the coarsest finite element space M , and this geometric growth rate. 49 The multi-level algorithms here actually use only the finite family of spaces {M .}. ,, for fixed 3 e B and fixed k > 2, in their 3,j 3=1' - computation. Since 3 is fixed, it is convenient to drop it, although we oo retain j, and in fact consider the entire sequence of refinements (M.}._ since that simplifies the statements of some of the theorems here. One word of warning is needed, however. Though the spaces in the sequence {M.}. n become infinite dimensional, they do not usually become dense in H' and thus convergence to the true solution will not ordinarily be £ oo obtained in {M.}._ . In order to get a sequence of solutions converging to the true solution, one must go back to consideration of the family a a£A OO Description of the families (M . } of locally refined grids is essentially the same for grids based on square and triangular elements. We consider first the triangular case since it is the case of greatest interest in this chapter. In this case, the family {M.}._, of finite r i 00 element spaces will be based on a nested family of triangulations lT.j. ... By nested we mean that each triangle T of T., j > 2, is contained in a r "> °° triangle of T. ,. Once we have chosen the triangulations It.;. .. , and 3-1 J J =1 the particular element to be used, the family of finite element spaces r i °° W.i. . is completely determined. J 3=1 The construction of the nested family of triangulations It j. , j J =1 may be done by a process Bank and Sherman call regular subdivision. In regular subdivision of a triangle T, it is divided into four smaller triangles by pairwise connecting the midpoints of the edges together as shown in figure 7. Let T be a fixed, aligned triangulation of Si. Now inductively suppose T. has been chosen for some i > 1. Select a subset T, of x.. We 3 ~ J J 50 may take T as the (unaligned) triangulation of fi formed by regular subdivision of the triangles in T. while leaving the remaining triangles of T unaltered. Proceeding is this way, the entire family of oo triangulations {t,}. , may be constructed, J 3=1 Figure 7. Regular Subdivision of a Triangular Element In our approach the subsets T. C t , j > 1, of triangles to be refined cannot be chosen arbitrarily, but most satisfy certain conditions necessary for the operation and analysis of the algorithms here. Let £2.,, C fi be the subdomain convered by triangles in T., for j+l J j _> 1. For simplicity of notation it is convenient to set (L = fi. The first condition required is (4.2.2) CD, , j > 1 . J+l J - This condition means that for j > 1, only triangles of T. which were formed by regular subdivision of triangles in T. - may be subdivided in creating T . _ . Equivalently, for j _> 1» each triangle of T. - is contained in a triangle of T,. Now assume that the diameters d™ of the triangles T in the coarsest triangulation T satisfy (A. 2. 3) - h n < d < a h, , a 1 — T — 1 for some fixed a > 1 and mesh parameter h . As before set 51 (4.2.4) h. = p 1 ' 2 h ± , j > 1 , where p is a mesh ratio, which for the regular subdivision process here is two. Conditions (4. 2. 2)- (4. 2. 4) imply that the diameters of triangles in T . , j >_ 1, must satisfy (4.2.5a) - h. < d < a h. for T C ft. o J - T - j j (4.2.5b) ~ h. < d_ < a) h. for T C ft. \ n ., a i — T — i i i+I with 1 < i < j - 1 . The overbar here denotes closure, necessary since triangles are closed, and our domains {ft.}. , are not. Intuitively, the triangles in T are all of about the same size. By (4.2.3) or (4.2.5a) their diameters are all comparable to h . The triangles in T come in two sizes: fine ones in ft~ with diameter comparable to h_, produced by regular subdivisions, and coarse ones outside ft„ with diameters about h . In general, T. may be thought of as having j sizes of triangles, the finest lying in ft.. Of course, this intuitive point of view should not be pushed too far. Because of the constant a in (4.2.5) these sizes may overlap a great deal. Only one other condition is really essential for triangular grids. That is, that the number of elements adjoining a given element along one edge be at most two. The "two" here is arbitrary, and any other finite bound could be chosen, but in practice two is probably all one wants. The coding is difficult enough as it is. This restriction, that at most two "smaller" elements adjoin one "larger" element across an edge, together with the angle condition (4.2.1) limits the rate of change of the mesh size as one can easily show. 52 Locally refined grids of square elements may be constructed along the same lines as the triangular grids just considered. One begins with an aligned square grid G on ft and then generates a family of finer .00 (unaligned) grids (G.} , using an analogous regular subdivision process. This process is shown in figure 8. The subdomains (ft./. , can be defined analogously, and the analogs of (4. 2. 3)-(4. 2. 5) will hold, now with the notation "d " instead of "d ", using "e" for "element." Since such e T ° square grids are only being considered in this section and will thereafter be dropped for simplicity, there is no need to go into these details more deeply. Figure 8. Regular Subdivision of a Square Element A problem arises on locally refined grids not present on quasi-uniform grids. Unless special care is taken, the condition numbers of the mass matrices for a family of locally refined finite element spaces may not be uniformly bounded. For aligned quasi-uniform grids and the usual types of elements, including those considered here, the uniform boundedness of the condition number of the mass matrices is equivalent to the angle condition (4.2.1). Unaligned quasi-uniform grids are rarely considered, for obvious reasons, but the same thing can be shown for such grids. 53 The problem with locally refined grids is that the number of basis functions which overlap with a given basis function may be very large, becoming unbounded as the grid is refined. This violates the hypotheses of the result of Fried described in section 2.3, and in fact mass matrices with condition numbers becoming unbounded do occur. A simple finite element grid where this problem occurs is shown in figure 9. Here bilinear elements are used with the standard "hat function" or "pagoda" basis functions. Notice that the conditioning problem in this example occurs despite the fact that the mesh size changes slowly in the sense that the diameters of adjacent elements differ at most by a factor of two. Proper scaling of the basis functions mitigates the problem but does not eliminate it. sei; i-+ ++ f+ i f f+ +i + 4 Figure 9. Finite Element Space with Poorly Conditioned Mass Matrix The only real solution is to arrange things so that each basis function overlaps only a bounded number of elements, with the bound independent of the level of refinement. This can be done either by modifying the basis elements to prevent excessive overlap, or by selecting 54 only grids where such overlap does not occur. The second approach is probably preferable in most cases since it tends to be simpler. Consider the approach of modifying the basis functions for the simple bilinear elements. The situation is similar for the- other square and triangular elements here. Ordinarily, the nodal parameter at the intersection of the four elements shown in figure 10 would be associated with the "hat function" whose support is the union of the four elements. When one or more of these four elements is refined, the basis function associated with its center node can be modified to have support on a smaller region. Up to rotation, there are four cases to consider, as shown in Figure 11. Clearly the case where all four elements are refined need not be dealt with since then the same considerations apply on the resulting smaller elements. The dark lines in figure 11 show the boundaries of the supports of the modified basis functions. The squares marked with asterisks may be subject to further refinement without violating the requirement that the mesh size change slowly. The point is that all squares so marked lie outside the support of these basis functions, so the number of elements these basis functions overlap, and hence the number of other basis functions they overlap, will be bounded by a fixed constant independent of the level of refinement. The only difficulty with this approach is that computations of the element integrals is somewhat more complicated since these basis functions are more complex. This approach was suggested to the author by D. Gannon, who has done considerable work with rectangular tensor product elements on locally refined grids, Gannon (1980). 55 <» Figure 10. Support of a Basis Function #- 1 * * * * * X ^ )K * * * * * Figure 11. Modified Basis Functions 56 In the second approach, one retains the simple Lagrangian basis functions, "hat functions," and their higher order analogues, but restricts or modifies the grids so that the problem of one basis function overlapping an increasingly large number of others never arises. With triangular grids modification is not difficult. This is the approach used in the Bank and Sherman code. Beginning with an unaligned triangulation T . , 1 _< j < °°, belonging to a nested family of triangulations IT.;. , satisfying the conditions given so far, their algorithm constructs an aligned triangulation x., refining T., on which this condition number problem does not arise. This is done in two steps. In the first step each triangle with the property that there are vertices of other triangles in the middle of more than one of its edges is regularly refined, as in figure 12. Once all such triangles are refined, including those produced in carrying out this process, the only triangles in the grid are those with either one vertex in the middle of an edge or no vertices in the middle of its edges. The former are refined by a different subdivision process, called "green" refinement, shown in figure 13. The resulting aligned triangulation T. is then used to generate the finite element space M. used by the multi-level algorithm. Notice that green refinement may decrease the smallest interior angle present on the grid. This would be a problem if successive green refinements were applied. In Bank and Sherman's approach this is avoided r i °° by first generating the unaligned triangulations Ix.j and then r ~ 1°° forming the corresponding aligned triangulation It . J- . _ , rather than taking T , j > 1, as a refinement of x.. This ensures that the angle condition (4.2.1) is satisfied. Even so, the triangulations IT./.* used by their code are not generally admissible in the theory here, even when 57 A Figure 12. Regular Refinement of Triangles with Edges Containing Vertices Figure 13. Green Refinement of a Triangle the underlying unaligned triangulations It./., satisfy all of our assumptions. The problem is that in producing an aligned triangulation T . , j >_ 1, a triangle may be green refined while in a finer unaligned triangulation T., i > j, this same triangle may be regularly refined. r~ 1°° If so, the triangulations It./ will not be nested, nor will the resulting r 1°° finite element spaces {M.i. ,. While this does not seem to be a problem J J=l for their code, it is certainly a problem for the theory here. The cure for this seems to be to require the unaligned triangulations It.} to contain only triangles with a vertex of another triangle in at most 58 one of their edges. If so, in going from {x } to the corresponding r ~ -.00 aligned triangulations It) only green refinement will be done. Since the triangles of T , j > 1, green refined in producing T. must lie outside fi. (the portion of T. lying in fi. is aligned), there is no danger of subsequently regularly refining them. Further, it is easy to see that those triangles of T. produced by green refinement must be green refined in producing T , i > j , so the finite element spaces based r~-|°° on the family of triangulations it)., will be nested and the theory of this chapter will go through for them. The advantage of this approach is that the condition number of the mass matrix for a finite element space based on an aligned triangulation will be uniformly bounded, as long as the angle condition is satisfied. This is so for linear elements since each pyramid basis function overlaps only those elements meeting at the vertex with which its nodal parameter is associated. The support is thus a polygon formed from these elements, as in figure 14. The other kinds of basis functions required for higher order elements, those associated with nodes on the edges or interiors of triangles, are not a problem either, since their supports are smaller than those associated with the vertices. Figure 14. Support of a Basis Function on an Aligned Triangulation 59 Though the above approach is quite attractive for triangular grids, since aligned triangulations are convenient in programming, it does not generalize directly to square elements, since there is no analog of green refinement for these elements. Instead, we must either use the modified basis functions described above, or restrict the rate of change of the mesh size so that the conditioning problem does not occur when using the standard basis functions. The simplest way to describe the required restriction on the rate of change of the mesh size is in terms of the process for constructing the grids iG . J- ._ of square elements we wish to consider. Say that two elements, e , e„ £ G., 1 ■< j < oo } touch if they have any point in common; that is, if they share a vertex or part of an edge as in figure 15. Let G\ , j _> 1, be the set {e G G.: e C ft.}. That is, G\ is the set consisting of the J J J most refined elements of G.. Define the interior of G\ , denoted 3 3 int(G|) as the set of elements of G\ not touching any elements in G.\ G\ . Then a sufficient condition for the mass matrices for the locally 3 3 r i °° r r> 1°° refined spaces {M.}._ , based on the grids \G.}._ , to be uniformly well conditioned is G. C int(G|) , j > 1 . j J ~ This condition suffices for the mass matrices to be uniformly bounded since no basis functions centered at nodes outside of G\ have supports 3 extending into int(G'. ). It follows that they only overlap elements outside ft., and so they only overlap a bounded number of elements. 60 II I ' < 1 Figure 15. Touching Elements With the two types of elements considered here and the variety of ways of ensuring the uniform boundedness of the condition numbers of the mass matrices, a number of options are available. However, it would be unwieldy to treat the theory in the rest of this chapter in the same generality, particularly the work in section 4.5, which relies on detailed analysis of the finite element spaces used. For this reason, we restrict the analysis that follows to C spaces on aligned triangular grids. This choice is motivated partly by convenience, and partly by the use of such spaces in the Bank and Sherman code. There would be no difficulty in treating the other types of finite element spaces considered here, however, the interpo- lation mappings of section 4.5 are slightly more complicated to analyze on unaligned grids. 61 4.3. Algorithms In this section we describe the multi-level algorithms to be considered in this chapter and begin the analysis of their convergence. The first of these algorithms is identical to one of those described by Bank and Dupont for quasi-uniform grids and is similar to algorithms described by Nicolaides, Hackbusch, Brandt, Bokhvalov, and others. It is also similar to the algorithm used in the adaptive finite element code of Bank and Sherman though there is a significant difference here which will be discussed. The second algorithm to be considered here appears to be new. It can be used to show the 0(N) complexity bound given in the next section under weaker assumptions than otherwise necessary, Both of the algorithms described here differ substantially from that described in the last chapter. There, to approximately solve the finite element equations on a member M., j ^2, of a family {M.} J ] j = of finite element spaces, approximate solution of a related problem in M _ was required twice; once at the beginning to obtain a starting value for smoothing iteration on M., once at the end of the iteration as a final correction. The algorithms in this chapter and other multi-level algorithms in the literature use approximate solution of the equations for M more often. After obtaining a starting value from M._ , the accepted solution is obtained through a sequence of multi-level (outer) iterations. Each of these outer iterations consists of a number of smoothing (inner) iterations on M., followed by approximate solution of a related problem in M. . Thus, approximate solution of the related problem in M is required an indeterminate number of times depending on how far the outer iteration on M. is carried. J 62 This kind of algorithm is somewhat more complicated than that of chapter three. However it made sense to analyze the application of this kind of algorithm to locally refined grids rather than the application of algorithms like those in the last chapter for several reasons. The main one was the desire to establish the optimal order convergence of an algorithm as close as possible to those being used in practice. This was more fully accomplished than it would have been if we had examined algorithms more like those in chapter three. It might also be pointed out that relatively little of the complexity of this chapter is caused by the choice of algorithm. Description of the algorithms here and analysis of their convergence properties will involve the eigenf unctions and eigenvalues. Following the notation in chapter three let N. be the dimensions of M., 1 <. j < °°. Then for the problem (2.1.3) we have eigenf unctions (i) N i (1)i^i U>) } . -, in M. and eigenvalues {at/. - satisfying i i=l j 1 1= 1 aO^ J) , ), e M. . Since more than one finite element space may be involved in our considerations, we will retain the superscript j. As in chapter three we may replace the L inner product by the bilinear form b.( # , •) defined there, where we now use the subscript j to keep track of the finite element space in {M.}._ giving rise to this bilinear form. With this -eplacemi {Xf^}.^ satisfying i i=l replacement we have as before eigenfunctions {i|j. } -Z.-, an d eigenvalues a(^ j) , cjO = A< j) b.(^ j) , ♦), ^ 2, of a family (M.}., of locally refined spaces, multi-level iteration will need to be applied to all of the spaces {M.}. .; thus, the iteration will be described for any space M. in J J=2 j 3 3=1 Beginning with a trial solution u_ in M. the outer iteration 0,n j r (i ) V 00 produces a sequence of iterates in j converging to the true discrete solution u in M.. Each of these outer iterations is J similar to the algorithm of the last chapter. Beginning with a trial solution u_ £ M. an outer iteration produces an improved solution 0,n j U 0^n+1 given by u (j) = u (j) 0,n+l m+l,n ' where {u n J }„ , is a sequence of inner iterates. All but the last of J6,n ic=l these inner iterates are generated by a simultaneous displacement smoothing iteration on M.. The last uses a correction based on the 3 approximate solution of a related problem in the coarser finite element space M. . . More precisely, beginning with a trial solution u n £ M., j-1 0,n j an outer iteration consists of the steps: 64 1. For I = 1,..., m define u„ £ M by (a(u , , «J>) - (f, W),46H J( j 2. Letting v £ M._ 1 be given by (4.3.2) a(v, (j)) = a(u , $) - (f , ) , € U , m,n j-1 let v be an approximation to v and set (4.3.3) ii = u - v . m+1 , n m , n The quality required in the approximation of v by v will be specified in theorems 4.1 and 4.3. The intuitive meaning of these two steps is that in the first the more oscillatory components of the convergence error are rapidly reduced by the smoothing iteration. In the second step the smoother components of the error, for which the simultaneous displacement iteration is ineffective, are reduced through approximate solution of the related problem (4.3.2). Consequently convergence will be rapid for all components of the error. To make these considerations precise, let . £ (0, 1) be an arbitrary parameter. Let J (6.) be the integer such that x (j) < e _ A (j) for i < j(e ) i — 3 N. — 2 3 x (j) > e> x (j) for ± > j( j . l j N. 3 Similarly let J(0.) be the integer such that 65 A? j) < 9. X^ j) for i < J(0.) , i N. L (j) > 6. \[p for i > J(6.) . J Then we can define the subspace S. 0. S. V, of M. by 3 3 J J J S. = span {\p^\ i < J(6.)} , 0. = span ii>^\ i > J(6.)} , S. = span {$f j) , i < J(8.)} , 3 i - j 0. = span {ipf j) , i > J(6.)} . J i J The spaces S. and S. consist of relatively smooth functions while 0. 3 3 3 and 0. contain more oscillatory functions. For any 9. £ (0, 1) one has: M. = 5. © 0. = S. © 0. 3 3 ^ 3 3 w J In fact, 0. is the orthogonal complement of 5. in M. in both the energy and L inner products, while 0. is the orthogonal complement of S. in 3 3 energy and the inner product b^ (•, •). For the time being only the subspaces 5. and V. will be important in analyzing the convergence of the multi-level iteration (4. 3.1)-(4. 3.5) , since the approximation question involved will be deferred until sections 4.5-4.7. There S. and 0. will play an important role. To establish the rapid convergence of iteration (4. 3.1)-(4. 3. 3) we will show that error components in the oscillatory space 0. are rapidly reduced in step 1, while components in the smooth space S. are reduced in step 2. The multi-level iteration (4. 3.1)-(4. 3. 3) differs from that in the Bank and Sherman adaptive finite element code primarily in the smoothing inner iteration used. Their code uses a symmetric Gauss-Seidel 66 smoothing iteration rather than the simultaneous displacement iteration used here. While comparisons for quasi-uniform grids have shown that the choice of smoothing iteration is not very important, Brandt (1977), Bank and Sherman (1979), for locally refined grids the choice is significant, In Fourier terms, the iteration (A. 3.1) rapidly annihilates Fourier error components whose wave length is comparable to h , the finest mesh spacing of M.. Away from ft., the most refined region, this iteration has little effect. The reason for this is that the residual on the right hand side of 4.3.1 is multiplied by the constant factor - Z7T\ • By contrast y in Jacobi or Gauss-Seidel iteration each component of the residual is multiplied by the inverse of the corresponding diagonal entry of the stiffness matrix. In other words, in these latter cases the scaling 2 2 factor looks like -ch, , rather than -ch., where h, - is a function local j local which takes on a different value for each component of the residual vector, depending on the size of the corresponding elements. As a result of this difference, the smoothing iteration used in the Bank and Sherman code is effective everywhere on the grid, rather than just on the most refined parts. This is important in their code since they allow any element of M j _> 2, to be refined in creating a finer space M.. We do not, and the simultaneous displacement iteration considered here is perfectly adequate in this simpler case. The choice of smoothing iteration is intimately connected with the decomposition of M. into smooth and oscillatory components. The decomposition here, M. = S. (+) 0.. is natural only for simultaneous displacement iteration. The hard part of the analysis here is the approximation question involved in demonstrating the effectiveness of 67 step 2 of the multi-level iteration (4. 3.1)-(4. 3.3) . For the decomposition M. = S. (•+} V. and the case considered here where only the finest elements of M._ 1 are refined in creating M., this approximation question can be answered. The author was unable to handle the more general case where any element of M. , may be refined. j-l Though this approximation result will not be shown until section 4.6, we now state it since it is needed for the analysis in this section. Let {M.}. , be a family of locally refined spaces satisfying the assumptions of section 4.2. Let j >_ 2 and 6. S (0, 1) be arbitrary, and for any T) € S . let T) e M be given by a(n - n, ) = o , <(> e M , that is, rj is the elliptic projection of n. onto M._ . Then (4.3.4) llln-n|||ic el /4 !||n||, , where c > is independent of j , 0., and r\. This inequality, which is the main result of this thesis, shows that eigenf unctions ty . in M whose eigenvalues are relatively small can be well approximated in the coarser finite element space M. ,. J-l With these preliminaries behind we are ready to consider the first convergence theorem for the multi-level iteration (4. 3. l)-(4. 3. 3) . This theorem is the same as one found in Bank and Dupont for the case of quasi-uniform grids and is similar to theorems proved by Nicolaides, Hackbusch, and others. Our notation largely follows Bank and Dupont although there are some differences. 68 Theorem 4.1 . There exists y G (0, 1) and a positive integer m, both independent of j ^2, such that if the approximate solution v of (4.3.2) satisfies (4.3.5) ll!v-v!|| lllv H>n -" 6 M , m j-1 a(ii - n , ) = , G M . m j-1 Since £ is the projection relative to the energy inner product mi - gii i wgii < a - <> >■ nie in and by (4.3.4) n - njll < c e 1/4 Hln.HI 1 c e 1/4 |||. | Since v = r\ + E, and e = n. + E, , m m m e m - v||| < c^ 174 |||e ||| + (1 - 6.) m |||e : Then since e in = e - v, m+1 m llle^lll 1 ll|e m - v||| + [||v - v||| _ 2 using iteration (4. 3.1)- (4. 3. 3) with j = j causes the energy norm of the convergence error to be reduced by the factor y £ (0, 1) according to this theorem. Also suppose we have an arbitrary iteration on M. . which reduces the energy norm of the J convergence error on M. _, by the same factor y €= (0, 1) at each iteration, J Then this iteration on M. can be used to approximately solve the J related problem (4.3.2) on M. - arising in using the multi-level J iteration (4. 3.1)- (4. 3. 3) on M. . Beginning with the trial solution J v Q = . oo in M. the iteration on M. .. produces a sequence of iterates {v }„ _ j -i j -i * *=0 satisfying flK-v||| 2. It follows that if j -1 > 2 the iteration on M. ., can be taken as the multi-level iteration (4. 3.1)-(4. 3. 3) with j = j n ~l- In this case we must solve in turn related problems on M. ~. Again this can be done recursively with iteration (4.3. l)-(4. 3. 3) if j -2 > 2. To make these ideas clear, we may consider the following informal procedure. Calling this recursive procedure with j = j causes n of the multi-level iterations (4. 3.1)-(4. 3. 3) to be carried out on M. . J Procedure Outer_l_on_M. (n) If (j = 1) {solve the equations on M. directly; then return } else {repeat n times{ 1. do m smoothing inner iterations on M. according to (4.3.1). 2. approximately solve the related problem (4.3.2) on M by calling Outer_l_on_M (2); then make the correction (4.3.3) }} return The computational cost of this multi-level iteration, which will be analyzed in the next section, depends on the dimensions of the spaces {M.}._ . In order to have a cost per iteration proportional to the number of unknowns, the cost of the smoothing inner iteration on 72 M must dominate the cost of recursively solving the related problems on j the coarser spaces. This will only be true if the dimensions of the spaces (M } , grow fast enough. Since these dimensions will grow quite slowly for many families {M ./._-. of locally refined grids this presents a problem. To ameliorate this difficulty, we consider now a second approach to solving the related problem (4.3.2) in the multi-level iteration (4. 3. l)-(4. 3. 3) . This second approach is based on a slightly sharper version of theorem 4.1. First we require a simple algebraic lemma. Lemma 4.2 . Let x,y be constants in (0, 1). If m is given by [ 1 - y m = *- I xy then Proof. We have (1 - x) : y . so Then so .. vin+1 _ (1 " X) . - 1 = 1+ (1 -x) +...+ (1 -x) m (1 - x) - 1 > (m + 1)(1 - x) m , (1 - x) m+1 - 1 < -x(m + 1)(1 - x) m 1 > (1 - x)"* 1 + x(m + 1)(1 - x) m = (1 + mx)(l - x) m , (1 - x) m < (1 + mx) - 1 < (1 + Z ) 1 = y . D 73 Using this lemma, we can show the following theorem, which for our purposes improves on theorem 4.1. All notation is as before. Theorem 4.3 . Let r G (0, 1) be fixed and let y G (0, 1) and j > 2 be arbitrary. Then there exists a positive integer m independent of j such that if the approximate solution v of (4.3.2) satisfies (4.3.8) |||v - v||| < r y |||v||| , then (4.3.9) III V!-"^ Ill is a constant independent of j and y. Proof . The proof is the same as that of theorem 4.1, except for the choice of the parameters m, 9., y at the end of the proof. Instead of inequality (4.3.7) we have (4.3.11) III e^ III 1 (c^ 174 + (1 - e.) m + ry) |||e ||| , 2 because the factor ry in (4.3.8) replaces the factor y in (4.3.5). Since we wish to show c 9^ /4 + (1 - B.) m + ry < y , it is enough to choose G . and m such that c n e 1/A + (i - e.) m < (i - r)y . For convenience assume c_ _> 1, without loss of generality, and define y by _ (1 - r)y y = 2c Then we may choose as and m as 74 6 j =y m = i - y V Using lemma 4.2, c.e}' 4 + a - e.) m < o^>l + (l_^)l< (1 . r)Y , establishing the first part of the theorem. For the second part we have m < e.y -5 y 2c 1 - r 2c ,5 -3 < 1 + (rT 2-) 5 Y" 5 < (i + (y^) 5 ) y" 5 , completing the proof. D Using this theorem, we can now analyze a modified version of the multi-level iteration given. This modified iteration will be the same as (4.3.1)-(4.3. 3) except that m, the number of smoothing inner iterations, will now be made to depend on the level, j. In applying this modified iteration to M. , j_ >_ 3, the related problem in M, _ will be solved by doing only one of the modified outer iterations on M. _. However, to J adequately solve this related problem on M, - using only one outer J iteration, rather than two as before, the number of smoothing inner iterations must be increased. 75 An informal statement of this modified outer iteration follows. Calling this procedure with j = j causes n multi-level iterations to be performed on M. . Here m( # , •) is an integer valued function to be given following the description of this procedure. Procedure 0uter_2_on_M. (n, j ) If (j = 1) {solve the equations on M. directly; then return } else {repeat n times{ 1. do m(j , j ) smoothing inner iterations on M. according to (4.3.1). 2. approximately solve the related problem (4.3.2) on M by calling 0uter_2_on_M. (1, j ), then make the correction (4.3.3). }} return It remains only to choose the function m(* , •) required by this procedure. Let r €E (0, 1) be arbitrary and let m > be the corresponding constant from inequality (4.3.10) of theorem 4.3. Also, let Y n e (0, 1) be arbitrary. Then we may set (4.3.12) m(j, j Q ) = -5 5 «-V m o y o r Since the argument needed to show that the multi-level iteration in 0uter_2 converges rapidly is somewhat more complex than that for r 1°° Outer 1, we state this as a theorem. Let iu^ j ., be successive — 0,n n=l iterates produced by the multi-level iteration in 0uter_2. Then we have 76 Theorem 4.4 . Let J > 2 be arbitrary. Then with m(*, •) given by (4.3.12) the iterates {u. } , satisfy 0,n n=l (j n ) (j n ) I^W ° "I £T UI-o.n " " I"' - for all n >_ 1. Proof. In order for each outer iteration on M. , i. > 2, to reduce the energy norm of the convergence error by the factor y , according to theorem 4.3, it suffices to do m of the inner iterations (4.3.1) on M. for some m < m Q Y q 5 , and then approximately solve the related problem (4.3.2) with relative accuracy ry . Examination of (4.3.11) shows that doing extra inner iterations never hurts so the choice -5 m = m (J > 3 Q ) = Lm Y is fine. Thus, for the case where j_ = 2 and the related problem on M. . is solved directly the proof is done. Otherwise, we solve the 3 related problem on M. with one outer iteration. Choosing Y = rY„ 3 in theorem 4.3 it suffices to do m inner iterations on M. n for some Jo" 1 m < m Q (rY ) and then solve the related problem on M. with relative accuracy r Y Q . Doing /• -, n I -5-5 m = m(j -l, j ) = Lm Q Y Q r inner iterations is certainly sufficient, so the proof is done for j n = 3. Continuing the induction in the obvious way completes the proof. D 77 The computational cost of the multi-level iterations described in this section and construction of an attractive overall algorithm for obtaining solutions to the elliptic problem (2.1.3) will be considered in the following section. 4.4. Complexity In this section an estimate of the computational cost of the multi-level iterations described in the last section is given, and these iterations are then used as building blocks for computationally attractive algorithms for solving finite element equations on locally refined grids. Under certain restrictions these algorithms will be shown to be of optimal order, producing a good solution to the finite element equations for M k 2l 1, in 0(N, ) operations, where N is the dimension of the finite element space M, . This complexity analysis is somewhat harder than the corresponding analysis in chapter three, since the complexity of the multi-level iterations in the last section depends on the dimensions of the finite element spaces {M.}. .. used. This complication 3 j=l was absent from the last chapter, since for the family {M.}. , of J J=l quasi-uniform spaces there, the dimensions of the spaces were essentially determined by the dimension of the coarsest space M and the mesh ratio p . A second minor difficulty also arises in proving the optimal order convergence of the algorithms here. A multi-level algorithm is considered to be of optimal order if, when applied to any finite element space M, in {M.}._,, it produces a good solution u in $(N,) operations, k J j— 1 k To define what "a good solution" means, an error bound applicable to locally refined grids is needed. We now proceed to derive such a bound. It is important here that the error bound derived not be overly 78 pessimistic. If it were, the results here would be correspondingly weakened. To have shown a multi-level method produces a solution containing a convergence error smaller than the expected truncation error means very little if one expects far more truncation error than actually results. Only a simple error estimate is required here. Such an estimate can be derived directly from the Bramble-Hilbert lemma. First we need the notion of a Sobolev seminorm. For I a non-negative integer, the seminorm | • | .is the same as the Sobolev norm || • || . except that it IT if contains only derivatives of order £, not those of lower order. That is, for u £ H (fi) , |u| = ( I l^ull 2 /' 2 . h*(B) |o|-£ l Just as for the Sobolev norm, these seminorms can be defined for non- integral I in a more complicated way. Now let M be a finite element space of C Lagrangian elements, based on a triangulation T of Q and satisfying conditions (4. 2. 1)- (4. 2. 5) Suppose the element trial functions are polynomials of degree k - 1. Then M is said to be a finite element space of degree k. For u G H(J2), 1 < I < k, the nodal interpolant u G M of u is well defined E 1 and we have Lemma 4.5 . For any triangle T £ T and for s = 0, 1 U " U T I s - a s d T l U l I 1 H S (T) S L H^CT) where d is the diameter of T. 79 The variable i here need not be an integer. This lemma is a special case of the Bramble-Hilbert lemma, which may be found in, for example, Strong and Fix or Oden and Reddy. r i °° Applying this lemma to the spaces in a family {M.}. , of 3 J =1 locally refined finite element spaces satisfying the assumptions of section 4.2 yields the required error bound. Assume that the 1+a differential equation (2.1.3) is H regular for a^ (0, 1), that 1+a is, that its solutions lie in H . Then the following error estimate holds on M. , i > 1. J - Theorem 4.6 . Let u be the finite element solution of (2.1.3) in M . . Then J (4.4.1) |||u-u^||| is independent of u, j, and the family of locally refined 00 spaces {M . } . i i=l Proof . Let u be the interpolant of u in M.. Then III- - u (j) ||| = min |||u- v||| < |||u-u< j) vSM. Since the energy and H norms are equivalent, (4.4.2) l||u-u (j) ||| < c Mu-u^M . H Applying the above lemma on each triangle in T . , M (j) M 2 .2 v ,2 (1+a) , .2 I u " V 2 - a n Z d T u l+« 1 L 2 (fi) °^T. T H 1+a (T) Thus, 80 (j)l 2 ^ J- v ^ 2a I I 2 U " U T I 1 < °1 ^ d U 1_L~ (j)ii 2 . 2 _ . ,2a , .2 u " u t J II ! < c E d | u | H X (fi) T^T. H ±W (T) 2 2 2 for c = max{a n , a } , assuming d < 1 for all triangles T £ T. Now J combining this with (A. 4. 2) the result follows. D Since this bound contains d , in effect the local mesh parameter, it is much tighter for locally refined grids than bounds in terms of the maximum element diameter. Moreover, the exponent of d here is optimal. For sufficiently smooth u the reverse inequality, bounding the error from below, may be expected to hold although this is difficult to show. Establishing that the multi-level methods of this section produce solutions whose convergence error is comparable to the truncation error expected from this error bound may be regarded as more than adequate for most purposes. In particular, for problems with singular solutions the improved convergence of finite element spaces using the proper locally refined grids may be readily shown using this bound. The multi-level algorithms here will inherit this improved convergence to singular solutions on properly refined grids. With this error bound established we are ready to turn to the design of multi-level algorithms for solving the elliptic problem (2.1.3) 00 on locally refined grids. Let IM.}._ 1 be a family of locally refined finite element spaces satisfying the assumptions of section 4.2. The ~(k) algorithms we wish to consider will produce a good solution u in any number M of this family, k >_ 1. The first step in these algorithms is 81 to solve the finite element equations for M, directly to obtain an approximate solution u to the problem (2.1.3). After this, solutions u £ M,, 2 < j < k, will be produced, each a good approximate solution of (2.1.3). Each solution u , 2 <_ j <_ k, is obtained by beginning with the previous solution u G M , considering it instead as a function in M . , and then carrying out a sequence of multi-level iterations to obtain the improved solution u G M.. This multi-level iteration may be taken as either of the algorithms of the last section. Two questions arise immediately: what conditions must the above process satisfy in order to obtain an accurate solution u , and what determines the computational complexity of this process? With the error bound (4.4.1) already proven, the first question is easier. ~ (k) We wish to obtain a solution u which contains only convergence error comparable to the expected truncation error. That is, we want the convergence error to satisfy III 00 ( k )|ll • t v a 2ql I I 2 \ 1/2 WW - « III 1 c ( E d |u| ) u T^T, l H ™(T) k ~(k) In order to have an algorithm whose solution u satisfies this inequality, it is convenient to ask that the intermediate solutions u , 1 j< j £ k - 1, also satisfy this type of inequaltiy. That is, for 1 < j < k, (4.4.3) |||5 ( J ) -u^HI ||| i a + 2 1 ^)- 1 iBstt-" - .«>« . Then (4.4.3) holds for 1 < j < k. Proof . By assumption IH~(1) (Din . f - ,2a | ,2 ,1/2 ll u " u III 1 c ( Z d u ) re Tl T H 1+a (T) Now inductively suppose |||~(j) (j)ii| < f v ,2a | ,2 .1/2 III 11 " u III 1 c n ( Z d T u 1-Kv ' TET T H 1+a (T) for some j , 1 <_ j <^ k. Then |||5 (j) -.||| < |||5 (3) -u (3) ||| + |||u (j) -.||| ; so |||5«> - ulH < 2c ( E df |u| 2 ^ ) 1/2 rer. T h 14< *(t) using the error bound (4.4.1). Then, since III (J +1 ) III * t v A 2a I I 2 ^/ 2 ll u - u lll 1 c n ( E d T u 14a ) ° TGT..- T H 14a (T) 3+1 we have (4.4.5) |||u^. U ^ +1 >||| - u||! <2c.( I df \u\ 2 ) l > 2 • a much stronger inequality for large a. Thus, even though one must work harder per element for smooth problems, it will still be cheaper to produce a good answer. Now we are ready to consider the design of multi-level algorithms based on theorem 4.7 and the multi- level iterations of section 4.3. The multi-level iteration used may be that described in either procedure Outer_l or 0uter_2. In either case the dimensions of the finite element space {M . s . -, will be of central importance. If N. is the dimension of M., j _> 1, optimal complexity bounds can be shown under the assumption that the dimensions {N.}. . grow at least 3 1=1 geometrically. That is, (4.4.6) N... ' <5N. , j > 1 , for fixed 6 > 1. Actually, this can be weakened slightly. It is enough oo that we have a sequence of numbers {x,}._ growing at least geometrically, * j+1 1 6x. , j > 1 , for 6 > 1, and that N. and x. are comparable, J J — x. < N . < c x. c j — 3 - J for fixed c > 0. For simplicity, assume that (4.4.6) holds. Using either of the multi-level iterations of section 4.3, the storage required to 85 approximately solve the finite element equations for M , k _> 1, can be estimated as follows. Let c be the storage required for direct solutions of the equations for M . Then the total storage required is proportional to k k I N. + c, < N. Z 6 J + c, 3=2 3 ^ k j = 2 1 < \ t^t + C l • Since c is fixed, for any 6 > 1 the storage required is proportional to N, . k Next consider the use of the multi-level iteration described by procedure Outer_l. According to theorem 4.1, when this multi-level iteration is applied to M., j > 1, the convergence error is reduced by a factor Y at every iteration, where y < 1 is independent of j. Actually, the proof of theorem 4.1 shows Y can be chosen arbitrarily, though this is unimportant, since for any y < 1 ~we can find V ^> 1 such that Y V < (1 + 2 1+ V 1 . Let y < 1 an d an integer v satisfying this inequality be fixed. Then ~(k) according to theorem 4.7, a solution u in M , satisfying (4.4.5) can K. be produced by carrying out V of the multi-level iterations described by Outer_l on each space M. , 2 _< j <^ k. For each M. , 2 j< j j< k, the approximate solution produced in the next coarser space M is used as a starting value for this iteration. To estimate the complexity of this procedure, note that in applying Outer_l to M., j >_ 2, 0(N.) operations are required plus the cost of recursively solving the related 86 problem in M , . This related problem is solved in Outer_l applied to M., by recursively applying Outer_l twice to the equations for M , , if j - 1 > 1. Continuing the recursion one sees that the total cost in doing V multi-level iterations of this type on the equations for M. is proportional to v( Z 2 j_i N + 2 j_1 cj < VN.( Z (|) j_i + (|) j_1 cj , i=2 1 l ~ J i=2 6 6 2 where c„ is the cost of solving the linear system for M . Since c„ is a fixed constant, the computation time T(M.) to do v of these multi- J level iterations on M. must satisfy f 0(N.) , 6 > 2 T(M.) < J 0(N. log (N.)) , 6=2 log 2 L0(N. log 6 ) , 1 < 6 < 2 . Then following the procedure described above where the approximate solution ~(k) u in M is obtained by successively obtaining good solutions in M., 1 <^ j <_ k, the total computation time T to obtain u must satisfy (\) T = Z T(M.) < I (N, log (N, )) ,6 = 2 j=l J \ k log 2 [(N k l0g 6 ) , 1 < 6 < 2 . For many problems a bound like this is quite adequate. For example, in doing local refinement along a plane singularity or boundery in a three dimensional problem, the dimensions {N.}. , would at least J J=l quadruple at each level. Since the two dimensional results here readily extend to three dimensional problems, this example is important. In other problems, for example, refining along a line sigularity in two 87 or three dimensions, the dimensions typically only double at each level, yielding the N log (N) bound. Still other problems arise for which the growth of these dimensions is even slower. In refining around point singularities, these dimensions will usually not grow geometrically, and the results of this section are of little use. This does not necessarily mean that multi-level methods should be avoided in these cases, only that theoretical results showing their effectiveness have not yet been given. When the dimensions of the spaces M. grow geometrically according to (A. 4. 7), but the rate of growth, 5, is not greater than 2, the complexity results above can be improved by replacing the multi-level iteration given by Outer_l with that given by 0uter_2. In equation (4.3.12), governing the number of inner iterations used on each level in 0uter_2, we are free to choose both y n an d r. For convenience, y n may be chosen as With this choice in applying 0uter_2 to the equation for M^ , J (j ) 2 <_ j <_ k, only one iteration of 0uter_2 is needed to get u satisfying (4.4.4) from u . The other free parameter, r, in (4.3.12) will be chosen as r = 6 u . Then for any 6 > 1, r 6 (0, 1), so theorem 4.4 shows that applying the multi-level iteration described by Outer 2 to M. , 2 < i. < k, reduces V °" the convergence error by Y n P er iteration as required. From equation (4.3.12), governing the number of inner iterations on each level, we can 88 ~ ( V ~ ( V 1} estimate the computation time to produce u from u , j > 2. This time, T(M. ) is proportional to j -5(j n -j) j -5 j _j E r N + c < N E (f ) + c. j=2 J Z J j-2 5 2 j -I . . < N. E (6 6 ) j j + c J j-2 X < N. (1-6 6 ) 1 + c 9 . Thus, the amount of computation required to produce u from u j > 2, is 0(N. ). Then, since J k -1 -1 E N. < N (1 - 6 X ) j =2 ^0 for any 6 > 1, the total computation time is 0(N ). This is an optimal bound, showing 0(N) complexity under the weaker assumption that the dimensions IN.}., grow geometrically, with any growth rate 6 > 1. Unfortunately, the dependence of this 0(N) bound on 6 is quite severe, though this could be improved by choosing the various constants differently, and by using a more effective inner iteration. The hidden constant in the 0(N) bound just given contains not only the factors (1-6 ) and (1-6 ) seen above, but also the constant m of (4.3.12), which depends indirectly on 6. By the proof of theorem 4.3 it can be seen that m n goes as 89 (1 - r) 5 = (1 - 6 6 ) 5 Thus, including the dependence on 6, the complexity bound is 1 0((1 - 6 V 6 - (1 - 6" 1 )" 1 N) . _1 /£ _C. _-| _ 1 C. While the factor (1-6 ) (1 - 6 ) already exceeds 10 for 6=2, in practice one would expect the dependence on 6 to be far less severe than this, although the appropriate numerical experiments have not yet been carried out. This completes the analysis of the algebraic aspects of the multi-level algorithms of this chapter. It remains to prove the approximation result (4.3.4), on which the complexity results just given are based. 4.5. Interpolation This section begins the analysis of the approximation properties of locally refined finite element spaces needed to justify the convergence and complexity results of the last two sections. For convenience this analysis is split into two parts. The first of these, given in this section, deals with properties of the interpolation mapping from M. to 00 M. ,, i > 2, for M., M. , members of the family {M.}. , of locally J-l - 3 3-1 3 3=1 refined spaces. The second part of this analysis, given in the following section will use the results of this section to obtain the required approximation theorem. While the analysis in the next section is completely general, the analysis here is specific to the C triangular elements and grids based on aligned triangulations being considered. However, there is no reason to believe that similar results could not be shown for more general 90 grids and for most common elements including triangular elements with higher continuity, square or quadrilateral elements, and elements with curved edges. The results here were not proven in greater generality primarily because it did not seem to aid understanding of the multi-level algorithms here and because the notations and arguments required would have been burdensome to both author and reader. Proof of results analogous to those in this section for most reasonable finite element spaces is not difficult. The work in this section is closely related to the work Bank and Dupont did in analyzing their two-level scheme discussed briefly in section 4.1. Much of the notation here and one of the lemmas, lemma 4.10, are taken directly from their paper. The difference here is mainly in emphasis. While lemma 4.10 plays a fairly small role here, it is central to the convergence proof for the two-level scheme in their paper. Lemma 4.11, which is not from their paper, is really much more central here. The reason for considering first the interpolation mapping is the usual one: while projection with respect to the energy inner product is global, interpolation gives a local projection more easily analyzed. Yet the two mappings are sufficiently alike that properties of one can be derived from the other. This is a common point of view, but the way in which properties of the elliptic projection are derived here from the interpolation projection is unusual. Instead of bounding the approxi- mation error in the interpolate and using this as a bound on the approximation error in the elliptic projection, we begin by establishing properties of the null space of the interpolation mapping. From these 91 properties, in the next section we will derive properties of the null space of the energy projection. Finally these lead to the required approximation result. This proof may seem circuitous, but the author was unable to establish an approximation result such as (4.3.4) in any other way. The problem encountered with the usual finite element approach is not hard to see. Since elliptic projection is global, an error estimate, such as theorem 4.6, for a locally refined space M, must involve the maximum element diameter h . However, the eigenvalue bounds involve max h . . Thus, one will be able to prove inequalities such as (4.3.4) but mm . h the constant c occuring will contain some power of the factor — . min This factor becomes unbounded for the locally refined grids being considered, so such a result would be too weak to prove 0(N) complexity bounds like those in the last section. Since the finite element spaces considered in this chapter are based on C Lagrangian elements, the interpolation mapping from M. to M. , , j > 2, with M., M. , members of {M.}. ,, is automatically well J-l - J 3-1 3 3=1 defined. For any u e M . , u is evaluated at the points corresponding to the nodal parameters of M. ... Then there is a unique u T €= M. . attaining 3-1 I 3-1 the same values at these points. One of the problems encountered when using elements with higher continuity is that this interpolation may not be well defined. There are, for example, C elements using first derivatives as nodal parameters and almost all C elements use second derivatives. In these cases care must be taken that the nodes of M. .. J-l lie at points where the corresponding values of functions in M. are well defined. In those rare cases where nodal interpolation is not well 92 defined, more complex interpolation projections based on various local averaging techniques must be used. Now let M , j _> 2, denote the null space of the interpolation from M, to M, .. M. is the subspace of M, consisting of functions vanishing at the nodes of M. .. Much of this section and part of the j-l next are devoted to establishing properties of these null spaces. For now we note that M. yields a direct sum decomposition of M.: J J w j-l This is easy to show, for example, by dimensionality arguments. Besides the direct method of defining the interpolate u T G M . of a function u G M., there is also a more abstract way which I J-l J we now consider. The interpolation mapping from M. to M,_ considered here is completely local in the sense that the value of the interpolate u,. £ M. , on a triangle T 6 T. , depends only on the value of u £ M. I J-l J-l J on T. For this reason the interpolation mapping from M. to M. , J J-l decomposes into separate mappings from M. | to Ml for each T G T -_-i* Rather than considering separately each of the function spaces M.| and M | T as T ranges over M. we can consider only a few canonical function spaces on T_ and the interpolation mappings between these canonical spaces. These canonical spaces on T^ can then be transformed into the spaces M. | and M | through the affine transformation F . As in section 4.2 let L be the element trial space on the reference triangle T . We must consider three different refinements is. of T_: regular refinement, green refinement, and the null refinement, 93 in which T_ is not altered. There are three possible green refinements R but we need only one, say the one with the new edge leading to the top vertex. Regular Green Null Figure 16. Refinements of T R Since each of these refinements of T D is a triangulation, finite element spaces on each of them can be generated in the usual way from the reference element trial space L. Let K be the finite element space on T generated in this R R way from the regular refinement of T and let K r be the finite element R Vj space generated by this green refinement. Null refinement just reproduces L. All three of these spaces generated by refinement, K , R K c , and L contain L as a subspace. Moreover, there are interpolation mappings I_ , I„ , and I. from K , K n , and L, respectively to L. The R U L R (jt last is of course just the identity on L. Just as for the interpolation mapping from M. to M. , we denote the kernels of I., and I„ by K. and j j-1 R b R K r , respectively. The kernel of I. is the trivial function space {0}. Now let T G T. , be arbitrary. In going from T. n to T, the J-1 .1-1 J triangle T is regular, green, or null refined. Consider the same type of refinement of T . Then there is an affine transformation R F : T -> T carrying the subtriangles of T generated by this refinement 1 R R to the subtriangles created in refining T. When T is regular or null refined, any of the six affine transformations have this property. When T is green refined two of the six affine transformations carry the green 94 subtriangles of T_ to those of T. The other four affine transformations K correspond to the other possible green refinements of T. To see the importance of all this, let I be the interpolation mapping from ML to M . | and let F' be the dual of F , carrying functions on T to functions on T , R by F^,: H^T) ■* H 1 ^) F.J,(u) = u o F T . For each of the three types of refinement of T we get a commutative diagram. For regular refinement we have F' L ' T M. ,1 J-1'T t X R t X T T R "j'T K «— ^ — M For green refinement we have F' L h i M f< 1-1 'T t 1 ! F' K G —^ W ll T And finally for null refinement we have 95 f; L •* — M. ,1 j-l'T ■ F T L ^-^-M.| T That these diagrams are in fact commutative is a consequence of the fact that F carries not only T and its subtriangles into T and its subtriangles, but carries nodes of T and its subtriangles into R corresponding nodes of T and its subtriangles. We omit the straight- forward but messy verification. The lemmas which follow just elucidate the properties of these three commutative diagrams and the various function spaces involved, The first result here shows that functions in L cannot well approximate •\ ^\ nonzero functions in K. or K , the kernels of I and I . The kernel, R G R G {0}, of I. contains no nonzero functions so the question does not arise /\ ^ there. The functions in K_ and K n vanish at the nodes of T- while R G R nonzero functions in L cannot vanish at all nodes of T_. The results R follow from this and the finite dimensionality of all spaces involved. For convenience define the interpolation n L : C°(T R ) -> L . The mappings I_ , I n , and I. described above are just restrictions of R G L 7T. to the subspaces K , K , and L of C (T ) , respectively. L R G R Lemma 4.8 . If w £ K or w 6 K , for any v £ L R 96 (4.5.1a) l|w-v|| >c||w|| L Z (T R ) l/(T R ) (A. 5.1b) |w - v| . >_ c | w | .. , E l (T R ) ' ti L CT R ) for c > independent of v and w. AAA, /\ Proof . Let K = K + K„. Then K is contained in the null space of tt. . Since tt. is the identity mapping on L, (4.5.2) K n L = 3 . Now consider the functional f(w, v) = || w - v || „ L 2 (T R ) defined on K x L. Let C be the compact subset {(w, v): ||w|| = ||v|| =1} L Z (T R ) L^(T R ) of K x i m Since f is continuous it attains its intimum of some (w_, v_) G c, but by (4.5.2) we have or |W Q - v || > c > , U L Z (T R ) |w Q - v || 2 > c 1 1 w 1 1 U U L Z (T R ) ° L^(T R ) Since f(w, v) is minimal at (w , v_) , (4.5.1a) holds on all of C, and then on K x L by linearity. For (4.5.1b) note that since the functions in K. are continuous 1 ~ i i and piecewise C , the only functions u £ K. for which |u| , =0 are the H constant functions. The only constant function in K is the function identically zero, since K is contained in the null space of TT. and the function values at the vertices of T^ are nodal parameters for L. Thus, K 97 • | - is norm on K. Because of this, the verification of H (4.5.1b) goes through exactly the same as for (4.5.1a). □ 2 For any triangle T C IR , each affine transformation F : T ->■ T induces a dual mapping F' carrying functions on T back to IK T functions on T^ : R F^: hV) -> H 1 ^) by Fj(u) = u o F T . One can then ask for the norm of this dual mapping F' either as a mapping from H (T) to H (T„) or from L (T) to L (T ). For the H seminorm the answer is that |u| , and |u ° F | are essentially IT(T) T T 1 ^) the same since the change in area over which | • | .. is computed and H change in the size of the derivatives are opposite effects and cancel perfectly (in two dimensions). Only triangles T €E T are considered in the following lemma both because they are the only ones of interest here and because the requirements of section 4.2 force all triangles in T . . to have diameter comparable to h. , , simplifying the statement J-l J-l of the lemma. Lemma 4.9. Let T G T , . j > 2 and let F be an affine transformation J-l - from T_ onto T. If u e H (T) , then u* = u ° F e H (T D ) and R R (4.5.3a) - |u| < |u*| < c |u| . C H X (T) H X (T R ) H X (T) 98 (4.5.3b) £ HI 2 iNll^Hz 1 C HI 2 l/(T) J l/(T R ) L Z (T) for fixed c > 1 independent of u, j, and T. 2 Proof . F = Ax + b for some 2x2 matrix A and b e R . Then as shown in Oden and Reddy (1976, eqs. (6.167) and (6.168)), for any s > (A. 5.4a) |u*| < ||a|| S Idet A|~ 1/2 lul * I I q — Hill I I I q ■ H S (T R ) H S (T) (4.5.4b) lul < Ha" 1 !! 8 Idet a| 1/2 |u*| * ' I I o — II III I lie ' H S (T) H S (T R ) where ||A|| is the spectral norm. Now let d T and d be the diameters of R T and T , respectively and let p_ and p_ be the diameters of the R T T R inscribed circles of T and T . Then as given in Oden and Reddy, lemma 6.2, the following simple bounds on the norms of A and A " apply llA||< P T R IK x |li We also have P T , ar ea (T) det A = area (T R ) The triangle T satisfies the angle condition (4.2.1) uniformly for T S t. - and j >_ 2. Using this condition it can be shown by simple geometric arguments that 99 — d < p c T - P T - d < p m c T - H T R R ^ d^ < area(T) < c d 2 \&\ < area(T R ) < c d^ , R R for c > 1 independent of T £ T, and j >_ 2 . Then since d = 1 by R assumption it follows that I|a||< c d T liA" 1 !! < e d" 1 |det A| 1/2 _ 2. Proof . The proof will be done locally over the triangles of T . _ . Let T S t. . be arbitrary. Then there exists an affine transformation J-1 F : T -*■ T. Letting w € $. be arbitrary we can define w* w ° F . Then T R 2 *■ w* £ K„ or w* £ K even when T £ T. , and w* = 0. Also we have R 6 j-1 {v* = v ° F for v £ M. } c L where the inclusion may be strict because of boundary conditions or interelement continuity conditions. Thus, inf |v ° F - w*| > inf |v* - w*| 1 > c |w*| 1 v£M._ 1 T H X (T R ) v*GL H X (T R ) " H X (T R ) inf ||v o F - w*|| ? > inf ||v* - w*|| ? > C ||w*|| , vSM._ 1 i L Z (T R ) v*GL L^(T R ) L Z (T R ) where the inequalities on the right are from lemma 4.8. By lemma 4.9 these inequalities can be converted into 101 inf |v - w| 1 ii c | w | i v€M._ H (T) H (T) inf v - w II „ > c ||w|| „ , ^J-l L 2 (T) " L 2 (T) where the constants here are Independent of d T « Combining these relations yields 2 2 2 inf || v - w|| . >_ inf |v - w| . + inf ||v - w|| _ veM._ 1 h (T) "" v£M. h (T) vGM L (T) > c(|w| 2 +||w|| 2 ) = c ||w|| H X (T) L (T) H (T) Then summing over the triangles of x 2 2 inf || v - w 1 1 = inf E ||v - w|| , v€M. , H v£M. , T^T. , H (T) 3-1 j-1 j-1 > E inf ||v - w|| 2 " T^t. , vSM. . H (T) j-1 j-1 >_ c E ||w|| = c ||w|| , . T^T H^T) H So using the equivalence of the H and energy norms inf |||v - w||| > c MI w Ml . J-1 This is equivalent to (4.5.5). To see this, assume |||w|| = 1. Then inf |||v - w||| > inf |||v - w||| > c |||w||| = c , |||v|||=l vSM vSM. , J_i J-1 so if |||v||| = |||w||| = 1 a(v, w) = l-\ |||v-w||| 2 < 1- c- (1- c) |||v||| |||w||| . The general case of (4.5.5) follows by linearity. □ 102 One more lemma is needed, which will be, for our purposes, the most important. It states that functions in M , the kernel of the interpolation mapping from M to M,_ 1 , are all quite oscillatory. Intuitively, the relatively smooth functions in M. agree quite closely with their interpolate in M. ,. Only the most oscillatory functions in M. can have the zero function as their interpolate. This observation J can be made precise quite easily through the use of lemma 4.9. For our purposes a function is "oscillatory" when its H norm or seminorm is 2 ~ much larger than its L norm. That functions in M, are oscillatory in this sense is shown by the following lemma. Lemma 4.11. If w €E M . , i > 2 then j' ~ (4.5.6) -llwll , < h. Iwl - < c M ii 2 - 3 i II- c ||w|| ? , L" J L with c > 1 independent of w and j To Proof . It is enough to show this locally on the triangles of T . - . see this, suppose that for all w S M. and T £ T ._ (A. 5. 7) i||w|| < h |w| < c || w|| L (T) J H (T) L (T) ' with c > 1 independent of w, T, and j. Then for any w £ M. \ Z ll w H 2 2 1 h i ^ I w| < c 2 E ||w|| 2 , c tet l (t) 3 tgt._ 1 h (T) ^-i L (T) which is the same as (4.5.6). To show (4.5.7), let w G M. and let T be an arbitrary triangle in T . .. Then there is an affine transformation F: T_. -*- T such that J-l R w* = w o F is in K or K„. Since • , is a norm on K_, and on K and R G tt 1 K b n 103 u 6 K or u £ L R 6 ^ Hull 2 1 independent of u. Then 1 i|W*|| 2 £ |w*| < c ||w*|| 2 l/(T R ) H X (T R ) l/(T R ) so using lemma 4.9 - II 1 h , l w l i 1 c ll w ll 2 J H (T) L L (T) J H (T) L (T) completing the proof. □ 4.6. Approximat ion In this section the approximation inequality (4.3.4) needed to justify the convergence and complexity results of sections 4.3 and 4.4 will be proven. This inequality generalizes a result of Bank and Dupont for quasi-uniform grids. Assume the elliptic equation is H regular for a ^ (0, 1). Then their result is that if (M.}._ is a family of quasi-uniform spaces and S. and 0., j >_ 2, are defined as before, the following inequality holds for T) €= S . , |||n-n|||£ce«/ 2 ||| n ||| , where r\ is the elliptic projection of n. on M. ... A domain with a slit, the worst case ordinarily arising, yields a = -r- and their result becomes 104 llln-n|||