SOME REFINEMENT ALGORITHMS AND DATA STRUCTURES FOR REGULAR LOCAL MESH REFINEMENT  RANDOLPH E. BANKy , ANDREW H. SHERMANz , AND ALAN WEISERx Abstract. In this paper, we consider local mesh re nement algorithms and data structures for nite element methods for linear elliptic partial di erential equations in the plane. Quadrilateral and triangular are treated in a uni ed fashion. Because we restrict the local re nement to be regular, the resulting nite element systems are always sparse, and the re nement algorithms can be implemented eciently, in time proportional to the number of elements. 1. Introduction. The nite element method for computing approximate solutions to linear elliptic partial di erential equations in the plane [12] is usually considered to have three parts: subspace selection, matrix assembly, and matrix solution. If a particular approximate solution is deemed too inaccurate, improvement may be obtained by selecting a larger subspace, using either mesh re nement or higher-order basis functions. In this paper, we consider various aspects of mesh re nement, particularly in the context of accuracy improvement in the presence of singularities in the solution to the partial di erential equation. A nite element mesh is a subdivision of the domain of the di erential equation into a number of polygons, called elements. In this paper, we deal only with meshes formed of shape-regular quadrilateral or triangular elements, and we assume that is, in fact, speci ed as the union of such elements. Extensions to allow for curved edges in the boundary of are straightforward. Throughout this paper, we assume that the nite element basis functions are piecewise linear for triangular meshes or bilinear (isoparametric) for quadrilateral meshes. In part, the reason for this is simplicity, since many of our remarks apply as well to higher order bases. Moreover, since we are concerned with mesh re nement in the presence of singularities, it may often be the case that little bene t can be gained through the use of higher order, more costly nite elements. For discussion of accuracy improvement through the use of higher order elements for problems with smooth solutions or point singularities, see [1, 14]. Mesh re nement strategies are of three types. In global mesh re nement, every element in the mesh is re ned to obtain a ner mesh. In most cases, this is the simplest strategy to implement, but it can be very wasteful in the sense that many elements are generated away from the areas of interest. A modi cation of the global strategy is semi-local mesh re nement, in which elements in one or more selected cross-sections of the mesh are re ned. In certain cases this strategy may be implemented as easily as global re nement and may be less wasteful. However, complications occur when the desired re nement does not correspond to a natural cross-section of the mesh. The last re nement strategy is (adaptive) local mesh re nement, in which only selected elements are re ned. This is an attractive strategy for problems involving either sharp fronts or point singularities, because the re nement can be restricted to those portions of where it is required. However, for adaptive local mesh re nement to be ecient, it is necessary that: (i) the best (or nearly best) elements to re ne can be selected cheaply;  Scienti c Computing, R. Stepleman et al (eds.), IMACS/North-Holland, 1983. y Department of Mathematics, University of California at San Diego, La Jolla, California z Exxon Production Research Company, Houston Texas 77001. x Exxon Production Research Company, Houston Texas 77001. 1 92093. (ii) the resulting systems of linear equations retain their sparseness as the mesh is re ned; and (iii) the local re nement process itself can be implemented eciently. For our purposes, we will assume that (i) holds; there is developing an extensive literature on this selection problem, and the reader is referred to such work as [14, 4, 8] for details. In this work we focus on (ii) and (iii). For example, if the system of equations becomes increasingly dense as the mesh is re ned, then the advantage gained by re ning only a small number of elements is partially lost. (In global re nement, all matrices retain the original sparseness.) And if the local re nement process cannot be implemented eciently, then the actual cost of computing an approximate solution to the di erential equation may actually be higher with local re nement than with the other two strategies, even in cases where the other strategies are extremely wasteful. An outline of this paper is as follows. Section 2 introduces notation and other basic information common to all the meshes we consider. Sections 3 and 4 develop formal criteria for regular meshes and re nement strategies that guarantee that (ii) holds for quadrilateral and triangular meshes, respectively. Section 5 describes a sample re nement algorithm that, we will show, generates re ned meshes satisfying the criteria developed in Sections 3 and 4. Section 6 discusses the implementation of the re nement algorithm, concentrating mainly on the required data structures, and demonstrating that the algorithm, in some sense, optimally satis es (iii). An appendix outlines the proof of several statements made in Sections 3 and 4. We note that some of the information in this paper has appeared previously, in less uni ed form, in [14, 7, 6, 5]. 2. De nitions and Notation. Let be a bounded, polygonal region of IR2 . The domain is assumed to be decomposed into a small number of macro elements ti , 1  i  nt0 which form a tessellation T0 of . The macro elements are either triangles or quadrilaterals, and we require that is ti ; tj 2 T0 , ti 6= tj , then ti \ tj is either empty, a common vertex, or a shared edge.1 We assume all macro elements are shape regular; that is, each element's side lengths are bounded above and below by xed multiples of that element's diameter, and all interior angles are uniformly bounded above 0 and below . In particular, quadrilateral elements must be strictly convex. Shape regularity does not require global quasi-uniformity of the mesh. Globally, we denote the vertices of T0 by vk , 1  k  nv0 and the edges by ek , 1  k  ne0 , As with the usual types of nite element computations (e.g., matrix assembly), it is important in our algorithms to have also have local notation to describe the edges and vertices of a given element. Thus we say that element ti contains ke vertices ij , 1  j  ke and ke edges "ji , 1  j  ke , where ke = 3 for triangles and ke = 4 for quadrilaterals (cf. Figure 2.1). Thus, for 1  i  nt0 , and 1  j  ke , ij = vk for some k , 1  k  nv0 and "ji = e` for some `, 1  `  ne0 . Throughout this paper, we will view local designations (e.g., ij ) and global designations (e.g., vk ) as interchangeable names for a unique entity, and use whichever designation makes more sense in context. To ensure the shape regularity of the elements in T0 is inherited by elements created during the re nement process, we restrict attention to bisection-type mesh re nement. A triangular element ti is subdivided into four geometrically similar triangles, called the sons of ti , by pairwise connecting the midpoints of the three edges 1 t denotes the closure of t . i i 2 i3 S  S "1 2 " S  t S  S  S S 2 1 3 i i i i i "i "4i i1 "1i i4 "3i ti i2 i3 "2i . Local Notation. Fig. 2.1 of t . Element t is referred to as the father of its sons. In the case of quadrilaterals, the sons are created by joining midpoints of opposing sides of t (see Figure 2.2). In this case, the sons are not geometrically similar to their father (except for the important special case in which t is a parallelogram), but they remain shape regular. Indeed, they become \nicer" in the sense that the bilinear mapping which maps the unit square (or reference element) to a son of t is more nearly ane than the corresponding mapping for t itself. In Figure 2.2, the numbers in the corners of the quartet of sons are to be interpreted as the superscript j in the local designations  + of the vertices of t + , 0  k  3. They de ne the geometric orientation of each member of the quartet relative to the father. In particular, note that vertices of each son in Figure 2.2 are locally ordered in a fashion consistent with Figure 2.1. It is this consistency, and the fact that the sons of an element bear known, xed relationships to their father, and not the speci c conventions we adapt, that are important in the algorithms we discuss. In the case of triangles, we will use a second type of re nement called \green" re nement, in which a vertex of t is connected to the midpoint of the opposite edge (Figure 2.3). While green re nement creates elements which can be less shape regular than their father, we will restrict its use so that no element created by green re nement is ever re ned itself. An admissible mesh [4, 3] is a collection of re ned and unre ned elements, which is recursively de ned as follows: (i) T0 is an admissible mesh; and (ii) if T is admissible and t 2 T is re ned, then T [ sons of t is admissible.2 Given an admissible mesh T , and vertex v which is a corner of each unre ned element it touches is called regular, and all other vertices are called irregular (Figure 2.4). Because admissible meshes allow irregular vertices, they are not necessarily tessellations of . The irregularity index [4] of a mesh T is the maximum number of irregular vertices on a side of any unre ned element. A k-irregular mesh has irregularity index no greater than k. For example, the mesh in Figure 2.4 is a 3-irregular mesh. Vertices lying on @ (which are necessarily regular) are called boundary vertices; all other vertices are called interior vertices. Each edge of an element t is either a boundary edge of or part of the perimeter of one or more other elements in the mesh. We de ne the neighbor of t across edge " to be the smallest element with one edge which completely overlaps " . For convenience, we denote both this neighbor and its global element number by  . It is i i i i i i j si si k k i i j i 2 Here sons of t are created only by regular re nement. 3 i j i j i i3 3 T  T  T  T  T tsi  T  T  T  T T T  T  T  T  T tsi  T  T  T  T tsi tsi  T  T  T  T  T  T T  T +3 1 11 2 22 +1 i1 +2 33 3 1 i4 4              tsi      tsi +3 3 tsi +2 21 34 1 Fig. 2.2.                  i2 3 E i 3 4 1 4 i1 2 E E E E E E E E E E E E E E E E E 2 3 tsi +1 2 1 2 2 i Bisection-Type Mesh Re nement. important to note that the neighbor relation is nonsymmetric and time-dependent. The level `i of element ti is de ned inductively as follows `i =  `fi 1 if ti 2 T0 + 1 if ti 2= T0 ; where tf is the father of ti . If we regard as a pseudo-element whose \sons" are the macro-elements of T0 , the process of creating an admissible mesh can be viewed as the creation of an element tree, in which the root is , other nodes correspond to macro-elements or elements created during the re nement process, and the branches lead downward from an element to its sons. All internal nodes in the tree have exactly four sons, except for , which has nt0 sons. Leaves of the tree correspond to unre ned elements. The level `i of ti then i 4 i3 T 1 1T  T  T  T  T  T  T  T  T  T  T  T  T tsi +1 T tsi   T  T 2 2 3T 2 3 1  i Fig. 2.3. i Green Re nement (One of Three Possibilities). t t t Fig. 2.4. t t t t Irregular Vertices. corresponds to the distance from t to in the tree; that is, the length of the shortest part from t to the root. In general, it is advantageous to \regularize" an admissible mesh T by restricting the number of irregular vertices on each edge. There are several reasons for this: simplifying computations such as matrix assembly and mesh re nement, increasing approximation power by insuring that neighboring elements are not of vastly di ering sizes, and guaranteeing that each element is in the support of a bounded number of basis functions. (cf. Sections 3 and 4). There are a number of ways to accomplish this regularization, but we shall mainly consider the following 1-irregular rule [14, 6, 5] and several of it variants. i i 1-Irregular Rule: Re ne any unre ned element for which any of the sides contains more than one irregular vertex. We shall see that applying such a rule as often as possible to an admissible mesh leads to a regularized mesh having several important properties described in the next two sections. The algorithm described in Section 5 incorporates both this rule T 5 and the ke ? 1 neighbor rule: ke ? 1 Neighbor Rule: Re ne any element with ke ? 1 neighbors that have been regularly re ned. 3. Quadrilateral Meshes. In this section, we consider bisection-type local mesh re nement for quadrilateral meshes (e.g, [14, 4, 2]). Let T0 be the given initial mesh and T be an admissible mesh generated from T0 . For convenience, in this section we assume T0 consists of a single square macro element. Thus we can use bilinear, rather than isoparametric basis functions in our discussions, and the initial adjacency information about elements in T0 is not required. These simpli cations do not have a signi cant impact on the issues addressed in this section. Let S = S (T ) denote the space of C 0 piecewise bilinear functions associated with T. S = f j  is continuous and jt2T is bilinearg We de ne a basis B = fbig for S as the unique set of elements of S satisfying the Lagrange conditions bi (vj ) = ij for all regular vertices vj in T . Here ij is the Kronecker delta. The basis B is well de ned [14, 4]. This choice of basis gives a natural correspondence between basis functions and the regular vertices in the mesh. Figure 3.1 shows a typical case in which the elements forming the support of bi are denoted by dots. t t vi t t t t Fig. 3.1. Support of bi . In a mesh with no irregular vertices, the support of a given basis function bi is limited to those unre ned elements which have vi as a vertex. There are thus at most four nonzero, linearly independent basis functions in any given element. These two properties imply that the nite element sti ness matrix will be sparse and the elementwise assembly procedure will be similar for all elements. The following example [14] shows that these desirable properties cannot be guaranteed for general admissible meshes. Let T0 be (0; 1)  (0; 1), and let T` be the (` + 1)-level mesh resulting from the regular re nement of exactly those elements containing the point p = (2=3; 2=3). T4 is shown in Figure 3.2. Since 2=3 = 1 ? 1=2 + 1=4 ? 1=8 + 1=16 ? : : :, T` is formed by alternately re ning the upper right or lower left son of the of the smallest currently re ned element containing p. Each mesh T` , `  3 is 2-irregular. 6 Fig. 3.2. T4 . The regular vertices in T` are the ten boundary vertices and the interior vertices ` vi =( ) p i ; pi ; pi = i  j X 1 j =0 2 1  ; i `: Then the smallest unre ned element in T` containing is in the support of the basis functions corresponding to i , 1   , and the resulting sti ness matrix will be essentially dense for large values of . The element assembly process will also be complicated due to the large numbers of nonzero basis functions in many elements. This provides the motivation for considering subsets of the class of admissible meshes for which certain desirable properties can be guaranteed. Our basic strategy is to take an arbitrary admissible mesh T and to generate a more re ned mesh T having the desired properties by applying one or more re nement rules. Let T have unre ned elements f i g, regular vertices f i g, and associated nite element space S with basis B = f i g. For example, a mesh T generated using the 1-irregular rule will satisfy the following properties (see [14]): Q1. T is 1-irregular; Q2. T uniquely contains the fewest elements of any 1-irregular mesh containing T; Q3. There are at most four basis functions j having support in any element i ; Q4. The restrictions to any element i of those basis functions functions nonzero in i are linearly independent; Q5. The support of any basis function i intersects the support of at most twelve other basis functions j ; Q6. jf i gj  9jf i gj; Q7. f i g can be partitioned into mutually disjoint sets 1 2 12 such that for any distinct i and j in the same set k , supp( i ) \ supp( j ) contains no elements. The e ect of properties Q1-Q7 is to ensure the eciency of the overall solution process. Q1, Q3, and Q4 combine to allow the assembly procedure for the sti ness matrix that is much simpler than what might be required for T . Q1 implies that the levels of neighboring elements di er by at most one, and allows the use of mesh data structures that conserve storage without a severe penalty in computational e ort (see Sections 5 and 6). Q2 shows that no mesh with these two properties can have fewer p v i ` ` 0 0 t 0 b 0 v 0 0 0 0 0 0 b t t b v 0 t 0 b t 0 0 0 0 t 0 V ;V ;:::;V v 0 v 0 V 7 b 0 b 0 0 elements than T . Q6 ensures that the sti ness matrix for T is not too much larger than the one for T . and Q5 shows that it is not signi cantly denser. (in fact Q5 and Q6 are usually very pessimistic; normally T has only a few extra elements, and rows of the sti ness matrix have nine nonzeroes rather than the thirteen suggested by Q5.) In the terminology of [3], Q6 shows that the overlap index of S is bounded by twelve. Q7 is of interest in applying iterative methods to the solution of the system for T . For instance, Q7 implies that an SOR iteration for the system can be broken up into separate vectorizable iterations for each V , since the new coecient for a basis function in V does not depend on any other new coecient for a basis function in V . In the Appendix, we outline a proof that for an arbitrary admissible mesh T , Q4 holds if and only if 0 0 0 0 0 k k k Q8. There exist no unre ned elements t1 , t2 and t3 with `1 < `2 < `3 , all containing a common irregular vertex on their boundaries. For quadrilaterals, the 1-irregular rule is equivalent to the following big element rule: Big-Element Rule: Whenever Q8 is violated, re ne t1 . We do not know of other rules that always require no more mesh re nement than the 1-irregular rule but are still sucient to guarantee properties Q3-Q6. For example, we might use the following middle element rule: Middle-Element Rule: Whenever Q8 is violated, re ne t2. The resulting mesh would satisfy Q3 and Q4 and might contain fewer elements than T , but no property analogous to Q5 could be insured. A 2-irregular rule analogous to the 1-irregular rule cannot guarantee Q3, since the mesh illustrated in Figure 3.2 would satisfy such a rule. In the case of quadrilateral elements, the 1-irregular rule by itself is sucient to guarantee Q1-Q7. Nonetheless, it may be desirable to also impose the 3-neighbor rule. By re ning an element with irregular vertices on three sides, one adds only two additional vertices but obtains four additional basis functions in S . Meshes generated using a combination of the two rules satisfy properties analogous to Q1-Q7, with, in some cases, larger constants. 0 0 4. Triangular Meshes. In this section, we consider bisection-type local mesh re nement for triangular meshes (e.g, [6, 5]). As in the case of quadrilaterals, our goal is to develop a procedure which is ecient for mesh generation and use in matrix assembly and solution. Let T0 be a given initial mesh and T be an admissible mesh generated from T0 . For convenience, in this section we assume T0 consists of a single triangle. Let S = S (T ) denote the space of C 0 piecewise linear functions associated with T . S = f j  is continuous and jt 2T is linearg We de ne the Lagrange basis B = fb g for S as in the case of quadrilateral elements. Given an admissible mesh T , we can obtain a more re ned mesh T , with unre ned elements ft g, regular vertices fv g, subspace S and basis B = fb g, by application of the 1-irregular rule. i 0 0 i 0 0 i 8 0 0 i The properties of the 1-irregular rule when applied to triangular meshes di er somewhat from the case of quadrilateral elements. First, unlike quadrilaterals, re nement of triangular elements does not always generate new regular vertices. Let the ` + 1 level mesh T` be generated by successive re nement of the center element of the smallest quartet of currently unre ned elements. T3 is illustrated in Figure 4.1. All the T` are 1-irregular. However, the only regular nodes for ` > 2 are the six boundary nodes.
A
A
A
A
A
A
A A
A

AA

A
A
A
A
A A
A
A A


A A

A A

A A
Fig. 4.1. T3 . Second, the 1-irregular rule for triangular meshes does not guarantee exactly three basis functions will be nonzero in each element. Consider the situation in Figure 4.2. The mesh satis es the 1-irregular rule, but the four basis functions corresponding to the vertices marked by dots are nonzero in t.
A
A
A
A
A
A
A A
A t

A
At A
A
A

A
A
t A
A
A A
t A

A
A
A
A
A

A
A A
t A
Fig. 4.2. Four Nonzero Basis Functions in t. Given a 1-irregular mesh T , we can generate a more re ned mesh T by applying, wherever possible, the following green rule: 0 00 Green Rule: With as few elements as possible, triangulate any unre ned element with an irregular vertex on one or more of it sides. The three situations in which the green rule can occur are shown in Figure 4.3. The following analogues of Q1-Q7 hold: T1. T2. T0 T0 T; is 1-irregular; uniquely contains the fewest elements of any 1-irregular mesh containing 9 
A
A
A
t A
A
t At
A A

A
t At
A
t A
A
A
A A

A
A
A  A

A
A
A
A
A
A (a) (b) (c) Fig. 4.3. The Green Rule. T3. There are at most three basis functions having support in any element 2T . T4. The restrictions to any element of those basis functions functions nonzero in are linearly independent; T5. The support of any basis function intersects the support of at most twelve other basis functions ; T6. jf gj  13jf gj and jf gj  2jf gj; T7. f g can be partitioned into mutually disjoint sets 1 2 9 such that for any distinct and in the same set , supp( ) \ supp( ) contains no elements. T5 and T6 are usually pessimistic. The most common number of nonzeroes in a row of the sti ness matrix is seven, and for most meshes encountered in practice T contains fewer than twice as many elements as T . The big element rule and the 1-irregular rule are not equivalent for triangular meshes; the big element rule may force more re nement than the 1-irregular rule. While the green rules solves the problem of irregular vertices by regularizing them at the conclusion of the re nement process, other strategies are possible. For example, a rule like Re ne the neighbors of any element whose re nement is dictated by accuracy considerations. used in conjunction with, say, the big element rule, forces the irregular vertices to lie in regions where re nement is not required for reasons of accuracy. In that case, regularizing may not be necessary. The 2-neighbor rule, used in conjunction with the 1-irregular rule, leads to 1irregular meshes in which each remaining irregular vertex is located at the midpoint of an edge of a unique unre ned element. When the green rule is applied to such a mesh, the only situation which occurs is that of Figure 4.3a. Analogues of T1-T6 hold for this process, but the constants are somewhat larger. Note that the green rule could also be used with quadrilateral elements, re ning quadrilaterals with irregular vertices on their sides into quadrilaterals and triangles. The three special cases for quadrilateral meshes satisfying the 1-irregular and 3-neighbor rules are shown in Figure 4.4. One can also use the green rule in conjunction with the -irregular rule for any xed 1. Analogues of T1-T6 would hold for such meshes, with di ering constants. The number of special cases requiring irregular re nement as in Figure 4.3 increases with increasing . The amount by which the shape regularity of the result elements 00 00 bj 00 ti 00 ti 00 ti 00 bi 00 bj 0 ti ti vi 00 ti 0 0 ti V ;V ;:::;V 0 vi 0 vj Vk 0 0 bi bj 00 k k > k 10 t t t
A
A
A A
(a) t (b) t
@
@
 
(c) Fig. 4.4. The Green Rule for Quadrilaterals. is degraded also increases. 5. A Local Mesh Re nement Procedure. In this section we present an algorithm for computing locally re ned meshes of the types described in Sections 3 and 4 [6]. For quadrilateral elements, our algorithm implements the 1-irregular and 3-neighbor rules, and for triangles, it implements the 1-irregular and 2-neighbor rules. Our scheme assumes that a logical-valued function DV T EST is available which indicates, for a given element in the mesh, whether the element is to be re ned. DV T EST can either be a user speci cation of a xed re nement pattern (say, by using element level numbers to control re nement) or be the output of a self-adaptive mechanism within the code which uses local error indicators. An element in the mesh may be re ned either because DV T EST indicates it should be re ned, or because it violates the ke ? 1-neighbor rule or the 1-irregular rule. Note that a given element may satisfy both rules when it it is processed, but violate a rule later in the re nement process due to the re nement of one or more of its neighbors. Thus it is clear that we must examine certain elements more than once. Procedure REFINE f If this is the rst call to REF IN E , THEN [nt nt0 ]; i 1; While (i  nt) Do [For j 1 to ke Do [If ij is unre ned Then [If ij has ke ? 1 or ke re ned neighbors Then DIV IDE (ij ); Else If `i > ` + 1 Then DIV IDE (ij )]]; If DV T EST (ti) Then DIV IDE (ti ); i i + 1]g; Procedure DIVIDE f si nt + 1; nt nt + 4; create ts +j , 0  j  3, along with associated verticesg; Logical Function DVTEST f. . . g j i i Fig. 5.1. Procedure The re nement procedure REF IN E 11 REFINE . is presented in Figure 5.1. Notice that is a one-pass algorithm in which elements are processed in the order they are created; suppose the re nement of tk is forced by the ke ? 1-neighbor or 1-irregular rule because of the re nement of ti (for whatever reason). REF IN E learns this fact (and acts on it) by testing the neighbors of an element against the ke ? 1-neighbor and 1-irregular criteria before the element itself is tested for re nement in DV T EST . REF IN E tk PP PP P ti  PP aa   aa   a  a  HH  HH  H Fig. 5.2. tk PP P j
PP ata PP
a  aa 
 aa   a  a  HH  HH  H Re nement of ti . For example, in the situation depicted in Figure 5.2, when ti is re ned, its four sons are added to the end of the element list. When son tj is processed, its neighbor tk violates the 1-irregular rule, so tk is re ned and it's sons are added to the end of the element list. If the re nement of tk causes new violations of the 1-irregular or ke ? 1-neighbor rules, this will be discovered and remedied when the sons of tk are processed. Since a given element has at most ke neighbors, we test (and possibly re ne) at most ke + 1 elements for any reason at any step of the process. We will show that generating the data to carry out any given step of REF IN E requires constant time. Thus, the complexity of REF IN E is linear in the number of elements. Since the 2-neighbor and 1-irregular rules are satis ed by triangular meshes generated in REF IN E , each remaining irregular vertex is the sole edge midpoint in some unre ned element. Each such element is re ned by GREEN into two triangles. Note that procedure REF IN E can be reentered using the mesh generated by a previous call to REF IN E . (This requires setting nt nt0 outside REF IN E .) It is thus easy to incorporate REF IN E into adaptive re nement schemes in which mesh generation is alternated with equation solution. Information regarding which elements are to be re ned at a given stage is passed to REF IN E via DV T EST . In the case of triangular elements, green triangles are logically removed from the mesh before invoking REF IN E , in order to avoid the possibility of creating new triangles with excessively small angles. 6. Data Structures and Implementation. In order to carry out procedure REF IN E and the other procedures normally associated with nite element calculations, we must be able to generate and retrieve certain data about elements and vertices in the mesh. We describe this data below in the form of functions. For a given element ti , we must be able to compute: F1. knots: knots(j; i) = k where ij = vk , 1  j  ke F2. neighbors:  j  = i if "ji is an interior edge ; 1  j  k n(j; i) e < 0 if "ji is a boundary edge 12 F3. father: f (i) = k where tk is the father of ti F4. sons: s(i) = F5.  si if ti is re ned 0 if ti is not re ned level: `(i) = `i For a given vertex vk we must be able to compute: F6. vertex type: 8 1 if vk is a user boundary vertex > > > > interior vertex > < 23 ifif vvkk isis aa user regular vertex vty(k) = > 4 if v is a regular boundary interior edge vertex k > > 5 if v is a regular interior center vertex > k > : 6 if vk is irregular User vertices are those associated with T0 (1  k  nv0 ). Regular interior edge vertices are those regular interior vertices which were created as the midpoint of an element edge during the re nement process. Regular interior center vertices were created as the center vertex in the re nement of a quadrilateral element. F7. vertex fathers: vf (j; k) = mj if vk was obtained as the midpoint of the edge with endpoints vm , 1  j  2 j Vertex fathers are well de ned for vertices of types 3, 4, and 6. If vk is a center vertex, it was created during the re nement of a particular quadrilateral element, and we can require its global vertex number k to be the largest of the nine vertices involved in that re nement step. The vertex fathers of vk are then de ned as an arbitrary pair of opposing edge midpoints in the re ned element. F8. coordinates: vx(k) = x-coordinate of vk vy(k) = y-coordinate of vk Since algorithm REFINE relies on integer triangle adjacency information, it does not require the x and y coordinates of the vertices; they are required, however, for the nite element matrix assembly process. The data structure associated with and created by REFINE corresponds to a re ned element tree,,and consists of four integer arrays: two short arrays IUSR and 13 IBOUND, and two larger arrays IRE and IV ERT . IBOUND is an array of length 2  nb0 , where nb0 is the number of boundary edges in T0 , and contains information about the boundary conditions. IUSR is a 2 ke  nt0 array whose i-th column contains knots(j; i) and n(j; i), 1  j  ke , for macro element ti . Let maxv (respectively maxt) be the total number of vertices (respectively elements) contained in the nest mesh generated. A re ned triangular mesh with maxv vertices will have maxt  8=3 maxv re ned and unre ned elements, while a quadrilateral mesh will have maxt  4=3 maxv total elements. Then IRE will be of size ke  maxt and IV ERT will be of size 2  maxv. 2  maxv real words of storage will be required for the x and y coordinates in the matrix assembly process. i 1 ? i 2 ? 3 ? 4 s(i) Fig. 6.1. 1 ? 2 ? 3 s(i) Macro Element Nodes in IRE ; Quadrilaterals (left) and Triangles (right). Suppose nt0 = 3 mod 4. (If not, renumber the macro elements starting at 2, 3, or 4 rather than 1 to accomplish this.) Columns 1 ? nt0 of IRE contain information about macro elements, as indicated in Figure 6.1. The remaining storage in IRE is devoted to re ned elements. If ti+j , 0  j  3 form a quartet of sons of a regularly re ned element, then columns i ? (i +3) of IRE contain information about the quartet as a whole. Because nt0 = 3 mod 4, given tk with k > nt0 , the decomposition (6.1) i = k ? k mod 4 j = k mod 4 both gives the starting column index i of the block of storage allocated to the quartet, and establishes that tk is the j -th son of tf (k) . If ti+j , 0  j  3, form a quartet of quadrilaterals, columns i ? (i + 3) of IRE contain the information indicated in Figure 6.2. If the four elements form a quartet of triangles, the the corresponding columns of IRE contain the information indicated in Figure 6.3. The information contained in both IRE blocks is the same, except that nine global vertex numbers are stored for quadrilaterals, while only six are needed for triangles. In the case of triangular elements, green triangles can be appended to the end of the IRE array. The son eld of the father should be set to point at the storage block allocated to the pair, and given a negative value to distinguish green triangle sons from regular triangle sons. For each pair of green triangles, we must store four global vertex numbers and the father of the pair, so two columns of IRE are sucient to contain all the required data for the pair. Information about vertices is stored in IV ERT , with column k corresponding to vertex vk . The type of vertex information varies according to vertex type, as illustrated in Figure 6.4. Here IBC is simply a pointer into IBOUND that allows determination of the appropriate boundary conditions. The IV F 1, IV F 2 and F elds require some explanation, since they contain the key adjacency information about the mesh. When a new vertex vk is created in the interior of , the situation is typically one of those shown in Figure 6.5 for triangular elements. In both cases, the creation 14 vknots(4;f (i)) vknots(4;i+2) ti+3 vknots(3;f (i)) ti+2 vknots(3;i) vknots(1;i+3) ti vknots(1;f (i)) vknots(3;i+1) ti+1 vknots(2;i) vknots(2;f (i)) i i+1 i+2 i+3 1 f (i) `(i) ? knots(3; i) 2 knots(2; i) knots(3; i + 1) knots(4; i + 2) knots(1; i + 3) 3 knots(1; f (i)) knots(2; f (i)) knots(3; f (i)) knots(4; f (i)) 4 s(i) s(i + 1) s(i + 2) s(i + 3) Fig. 6.2. Re ned Quadrilateral Node in IRE . of vk is caused by the creation of the quartet of regular elements ti+j , 0  j  3, i = 0 mod 4. Since vk is irregular F is set to m, and IV F 1 is set to i + j , where vk is the midpoint of edge "jf (i) ("jf+1 for quadrilateral elements) of element tf (i) . (i) A decomposition of IV F 1 as in (6.1) give the address of the storage block for the quartet in IRE , and establishes the geometric relationship of vk to the quartet. If the situation is as in Figure 6.5a, then F is updated when tm is re ned, but vk remains irregular. If the situation is as in Figure 6.5b, vk will become regular if element tm is re ned. In this case the IV F 2 eld in the same way as IV F 1, but for the second quartet. If vk is a center node in the re nement of a quadrilateral into elements ti+j , 0  j  3, then the IV F 1 eld is set to i, since the geometric relationship of vk to ti+j is established by its vertex type. The actual re nement process consists of lling in entries in IRE and IV ERT . When an element is re ned, all the information required for IRE is trivially known, except for the global vertex numbers for the edge midpoints. (The center node in the quadrilateral case must be created, so its global vertex number is known.) For each edge, we must either create new irregular or (regular) boundary vertices, or change current irregular vertices into regular ones. To do this, we must have the ability to compute knots(j; i), n(j; i), 1  j  ke , f (i), s(i), and `(i) for any element currently in the mesh. When a new vertex is created or changes status, all the necessary information is present, having been generated in the process of determining whether or not an irregular vertex already existed on the edge. 15 vknots(3;f (i)) T  T  T  T  T  T  T t +3  T  T  T v v (2 ) (1 ) T T    TT  TT   T T  t  T T   T T   T T  t +1 t +2  T T   T T   T T   T T i knots ;i knots ;i i i i vknots(1;f (i)) vknots(3;i) +1 (1 ) (1 ( )) ( + 1) i 1 2 3 knots ` i knots s i s i Fig. 6.3. ? =1 0 +2 (2 ) (2 ( )) ( + 2) i () () () f i vty vknots(2;f (i)) vty knots ;f i knots ;i ;f i Re ned Triangle Node in =2 =3 0 1 0 ? vty vty =4 1 2 knots Nodes in s i =5 =6 1 ? 1 0 vty IV F IV F ;i ;f i IRE . vty IV F I BC Fig. 6.4. i knots s i IV F I BC +3 (3 ) (3 ( )) ( + 3) i ;i IV F F IV ERT . We now discuss brie y the computation of functions F1-F8. If 1   0 , ( ), ( ) and ( ) are found by table look-up in and . () and ( ) are trivially known. For 0 , one makes the decomposition (6.1); then ( ), ( ), ( ) and ( ) can be found by table look-up. Finding neighbors is more complicated. For triangles, if = 0 mod 4 (the center element) then ( ) are trivially known. If 6= 0 mod 4, the one neighbor (the center element) is known. The other two are found by checking the status of the two edge midpoint vertices of ( ) which are corners of . If such a vertex is a regular boundary vertex (type 3) or an irregular vertex (type 6), then ( ) is found in the or elds, respectively, of . If the vertex is a regular interior vertex (type 4), then both 1 and 2 are analyzed as in (6.1). One of 1 or 2 points to the quartet containing , and the other points to the neighboring quartet. Since we know how this vertex is geometrically related to each quartet, we can easily determine how the two quartets are related. For quadrilateral elements, two neighbors are known, and the other two can be determined by a procedure analogous to that for triangular elements. For a given vertex , its vertex type can be determined by checking the signs of i knots j; i n j; i s i ` i knots j; i I U SR I RE nt f i i > nt s i f i ` i i n j; i i tf i ti n j; i I BC F I V E RT IV F IV F IV F vk 16 IV F ti 

A
A
A

vkA

tiAA
A
A

A A

A

A

A

tm A

A
A

A
(a) Fig. 6.5. Creation of Vertex
A

A
A

A
A

A t
A

i AA
v kA



tm A

A
A

A
(b) vk . its two entries in IV ERT . If k > nv0 , the vertex fathers are found by decomposing IV F 1 as in (6.1). Once the geometric relationship of the vertex to the quartet pointed at by IV F 1 is known, the vertex fathers can be looked up in IRE . Thus, each piece of data required by algorithm REFINE can be obtained in constant time. Another data structure which can be applied to arbitrary admissible meshes, but which does not obtain data in constant time, is presented in [10]. The nite element assembly procedure in addition requires the actual x and y coordinates of the vertices. These are easily generated from the vertex fathers when k > nv0 using vx(k) = vx(vf (1; k)) +2 vx(vf (2; k)) vy(k) = vy(vf (1; k)) +2 vy(vf (2; k)) Since vf (j; k) < k for k > nv0 , j = 1; 2, and the x and y coordinates for the vertices in T0 are given, the remaining coordinates can be generated in one pass through the vertices in the order in which they were created. We conclude with several remarks on storage requirements. For a given mesh with maxv vertices, a standard data structure for representing the nite element mesh is a list of pointers from unre ned elements to adjacent vertices and a list of pointers from vertices to adjacent unre ned elements [11]. For triangular elements, the number of unre ned elements is about 2 maxv, so the standard data structure requires about 12 maxv storage. IRE and IV ERT require about 10 maxv storage for the same mesh. For quadrilateral meshes the number of unre ned elements is about maxv, so the standard data structure requires about 8 maxv storage, as opposed to 7 1=3 maxv storage for IRE and IV ERT . Hence our data structures, which allow for local mesh re nement, actually require less storage than standard data structures for static meshes. 7. Appendix. Proofs of properties T6, T7 and the equivalence of Q4 and Q8. We now prove property T6: There are fewer than 13 times as many unre ned elements in the triangular mesh T as there are in T . Proof. Since T0 consists of one triangle, the number of unre ned elements in either T0 or T is three times the number of re ned elements plus one. Since T is admissible, each element t shares a corner with at most 13 di erent elements of the same level as t, so it suces to show that 0 0 17 A1. Each re ned element 0 2 T 0 shares a corner with a re ned element 2 T . with the level of equal to the the level 0 of 0 . A1 follows by induction on 0 . Let max be the largest level of an element in T 0 . 0 If = max ? 1 then 0 must be in T , so let = 0 . Otherwise if 0 2 T , 0 was forced to re ne by the 1-irregular rule, so there must be an element 0 sharing part of a side with 0 , with the level of 0 equal to 0 + 1. By induction, there is a re ned element 2 T sharing a corner with 0 with level 0 + 1. Then , the father of , must satisfy A1 (Figure 7.1). The proof is very similar to the proof of property Q6 given in [14]. t ` t t ` ` ` ` t ` t t t t = t s t s ` s s ` t s @@ @ @@ @@ @@ @@ @@ @ 0 @ @@ @@ @@ @@ 0 @@ @ @ s s t Fig. 7.1. Proof of Property T6. We now prove property T7: The regular triangular mesh vertices f i0 g for the mesh T 00 can be partitioned into mutually disjoint sets 1 2 9 such that for any distinct i0 and j0 in the same set k , supp( 0i ) \ supp( 0j ) contains no elements. Proof. For simplicity, we consider the case in which T0 consists of the triangle with vertices (0 0), (0 1) and (1 0). Let min( ) denote the minimum level of an unre ned element in T 0 for which vertex i is a corner. Let v V ;V ;:::;V v ; v ; V ; ` b b i v VL ` =f ij v ` min(i) = `g: Since T0 consists of one triangle, at most six edges in T 0 meet at i so i is not nonzero in any element in T 0 at level min( )+3 or greater. (Figure 7.2 depicts the worst case: min ( ) = k ? 2.) Recalling the coordinates of i are ( i i ), let v ` ` i b i ` v 0 1 V 2 V V x ;y = f i j ( i ? i )2` ?1 = 0 mod 3g = f i j ( i ? i )2` ?1 = 1 mod 3g = f i j ( i ? i )2` ?1 = 2 mod 3g v x y v x y v x y i i i : The vertices in sets 0 , 1 , and 2 when i = 3 are labeled 0, 1, and 2, respectively, in Figure 7.3. By construction, supp( i ) \ supp( j ) contains no elements for distinct k \ ` . Thus we can set i j 2 V V V ` b v ;v V b V L +3k V` =( V L ` [ V L +3 [ V L`+6 : : :) \ V k ` 18 @@ @@ @@ @@ @@ @@ @@ @ @@ @ @ @ @@ @@@@@@@ @@ @ @@ @@ @@ @@ @@ @@ @@ @ @ @ tj vi tk `min (i) = `(tj ). Fig. 7.2. 2 @@ 0 @1 @@ @@ 1 @2 @ 0 @@ @@ @@ 2 @ 0 @@1 @@2 @@ @@ @ @ @ 0 @@ 1 @@ 2 @@ 0 @ 1 Fig. 7.3. V 0, V 1, and V 2. for = 1 2 3, = 0 1 2. A similar construction shows property Q7. Note that can be constructed in a simple one pass algorithm processing vertices. Finding fewer sets which satisfy property T7 is possible (since T is planar, four sets suce, but algorithms for constructing such sets may be more expensive. We now outline a proof of equivalence of Q4 and Q8, i.e., ` ; ; k ; ; V1 ; : : : ; V9 00 A2. In each element in an admissible quadrilateral mesh T , the restrictions to of the which are nonzero in are linearly independent t bi t t holds if and only if A3. No irregular vertex in T touches unre ned elements 1 , 2 , and 3 with levels 1 2 3. v ` < ` t t t < ` Proof. First, suppose that A3 does not hold: there exists a vertex touching 1 , , and 3 of di erent levels. Suppose 1 is the largest element in T touching such a vertex. Each straight line segment (consisting of one or more edges) in T touches a regular vertex, so we can assume that 3 , 4 , and 5 are chosen to be the rst regular vertices away from on their straight line segments. Figure 7.4 depicts this situation, where 3 , 4 , and 5 may, in fact, be corners of , 3 , and 2 respectively. Since 2 1 and 3 1 , corners 1 and 2 of 1 must be regular as well. It follows that 1, 2 , 3 , 4 , and 5 must all be nonzero in element . Since their restrictions to are bilinear functions, A2 cannot hold. v t2 t t t v v v t v ` b v > ` b v t v v t t ` t b t 19 t > ` b b tv5 tv3 t t2 v1 t3 v t tv4 v2 t t1 Fig. 7.4. Condition A3 Violated. Now suppose A3 holds. Since fb g is a basis for S , to show A2, it suces to show that there are at most four b which are nonzero in any element t. Suppose that t has four regular corners v1 , v2 , v3 , and v4 . Then b1 , b2 , b3 , and b4 are nonzero in t, and all other b must be zero by the Lagrange conditions and continuity. Suppose t has exactly one irregular corner v4 . Then v4 must be on the side of a neighbor t2 of t, connecting corners v3 and v5 of t2 , with `2 less than the level of t. It follows that v3 and v5 are regular vertices, as are the other corners v1 and v2 of t. Then b1 , b2 , b3 , and b5 are nonzero in t, and all other b must be zero in t (Figure 7.5). The case when t has exactly two irregular corners is handled similarly. Note that t cannot have three irregular corners. Thus A2 holds in every case. i i i i v1 v2 t t t v t3 v4 t2 v5 t Fig. 7.5. t Has One Irregular Corner. Acknowledgments: We are grateful to the management of Exxon Production Research Company for permission to publish this paper, and to Marie Mason of EPR for eciently typing it. This research was supported in part by the Oce of Naval Research through contracts N00014-80-C-0645 (University of Texas at Austin), N00014-82-K-0197 (University of California at San Diego), and N00014-76-C-0277 (Yale University). REFERENCES [1] I. Babuska and M. R. Dorr, Error esitmates for the combined h and p versions of the nite element method, Tech. Rep. BN-951, Institute for Physical Science and Technology, University of Maryland, 1980. [2] I. Babuska and A. Miller, A posteriori error esitmates and adaptive techniques for the nite element method, Tech. Rep. BN-968, Institute for Physical Science and Technology, University of Maryland, 1981. [3] I. Babuska and W. C. Rheinboldt, Error estimates for aadptive nite element computations, SIAM J. Numer. Anal., 15 (1978), pp. 736{754. [4] , Reliable error estimation and mesh adaptation for the nite element method, in Computational Methods in Nonlinear Mechanics, North-Holland, New York, 1980, pp. 67{108. [5] R. E. Bank and A. H. Sherman, A multilevel iterative method for solving nite element 20 [6] [7] [8] [9] [10] [11] [12] [13] [14] equations, in Proc. Fifth Sympos. on Reservoir Engineering, Society of Petroleum Engineers of AIME, Dallas, 1979, pp. 117{126. , A re nement algorithm and dynamic data structure for nite element meshes, Tech. Rep. CNA-166, Center for Numerical Analysis, University of Texas at Austin, 1980. R. E. Bank, A. H. Sherman, and A. Weiser, On the regularity of local mesh re nement, in Proceedings of the IMACS Tenth World Conference, 1982, pp. 61{64. R. E. Bank and A. Weiser, Some a posteriori error estimators for elliptic partial di erential equations, Mathematics of Computation, 44 (1985), pp. 283{301. D. B. Gannon, Self Adaptive Methods for Parabolic Di erential Equations, PhD thesis, University of Illinois, 1980. W. C. Rheinboldt and C. K. Mesztenyi, On a data structure for adaptive nite element mesh re nements, ACM Trans. Math. Software, 6 (1980), pp. 166{187. R. B. Simpson, A survery of two dimensional nite element mesh generation, in Ninth Manitoba Conference on Numerical Mathematics and Computing, Manitoba, 1979, pp. 49{124. G. Strang and G. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cli s, New Jersey, 1973. J. R. Van Rosendale, Rapid Solution of Finite Eelement Equations on Locally Re ned Grids by Multi-Level Methods, PhD thesis, University of Illinois, 1980. A. Weiser, Local-mesh, local-order, adaptive nite element methods with a posteriori error estimators, Tech. Rep. 213, Computer Science Department, Yale University, 1981. 21