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
At
A
A
A
A
A A A
A
At
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