Parallel Algorithms for Adaptive Mesh Refinement.
Computational methods based on the use of adaptively constructed nonuniform meshes reduce the amount of computation and storage necessary to perform many scientific calculations. The adaptive construction of such nonuniform meshes is an important part of these methods. In this paper, we present a parallel algorithm for adaptive mesh refinement that is suitable for implementation on distributed-memory parallel computers. Experimental results obtained on the Intel are presented to demonstrate that for scientific computations involving the finite element method, the algorithm exhibits scalable performance and has a small run time in comparison with other aspects of the scientific computations examined. It is also shown that the algorithm has a fast expected running time under the parallel random access machine (PRAM) computation model.
Introduction
. Adaptive mesh refinement techniques have been shown to
be very successful in reducing the computational and storage requirements for
solving many partial differential equations [10]. Rather than use a uniform mesh
with grid points evenly spaced on a domain, adaptive mesh refinement techniques
place more grid points in areas where the local error in the solution is large. The
mesh is adaptively refined and/or unrefined during the computation according to
local error estimates on the domain. This technique is much more efficient than
the use of uniform meshes when the solution is changing much more rapidly in
some areas than in others.
The adaptive construction of these nonuniform meshes is a crucial part of
adaptive mesh solution methods and has been examined by many researchers, for
example, [3], [10], [11], [12], [13], [14], [15], [16], and [18]. Typically, one begins
with an initial mesh conforming to a particular geometry. This mesh is selectively
refined, based on local error estimates, to construct a mesh that satisfies a certain
error tolerance. Most research has focused on meshes composed of simplicial ele-
ments: line segments in one dimension, triangles in two dimensions, or tetrahedra
in three dimensions. This paper focuses primarily on two-dimensional simplicial
meshes. However, the algorithms and analyses presented here are applicable to
other dimensions and to nonsimplicial meshes.
In this paper, we present a new parallel algorithm for the adaptive construction
of nonuniform meshes. This algorithm is well suited for implementation on
medium-grained distributed-memory parallel computers such as the Intel DELTA.
The algorithm is based on the simplicial bisection algorithm given by Rivara [15].
Our algorithm is scalable in that it has an expected run time that is a very slowly
growing function of the triangles in the mesh.
This work was supported in part by the Office of Scientific Computing, U.S. Department of
Energy, under Contract W-31-109-Eng-38.
y The address of the first author is Computer Science Department, University of Tennessee,
Knoxville, TN 37996. The address of the second author is Mathematics and Computer Science
Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439.
To demonstrate the performance of this algorithm, we present experimental
results obtained on the Intel DELTA. The results demonstrate that, for practical
scientific calculations, the algorithm exhibits scalable performance and a run time
that is much smaller than the other computations necessary for the entire solution
method.
This paper is organized as follows. In x2 we review methods for adaptive mesh
refinement. In x3 we present our algorithm and analyze its expected run time under
the P-RAM computation model. A medium-grained distributed-memory version
of the algorithm is detailed in x4. We discuss our experimental results from the
Intel in x5. Finally, in x6 we summarize this research and discuss possible
future work.
2. Adaptive Refinement Methods. The finite element method has proven
to be extremely effective in the computation of approximate solutions to partial
differential equations (PDEs). Our focus is on adaptive, or local, refinement strategies
for generating finite element meshes. This approach can obtain much more
accurate solutions to these problems than a uniform mesh with the same number
of elements.
The computation of an approximate solution to a PDE consists of three main
tasks: (1) the construction of the finite element mesh, (2) the assembly of a sparse
linear system, and (3) the solution of this linear system. Although we are not
explicitly discussing the last two tasks, they must be kept in mind. In particular,
any method of local error estimation requires the approximate solution on a given
element mesh.
With a parallel implementation, we must remember that it is essential that
any adaptive mesh algorithm be integrated with parallel algorithms for the matrix
assembly and the solution of the resulting linear system. In addition, with
an adaptive strategy the assignment (or partitioning) of elements and vertices to
processors must be updated with each modification of the mesh to ensure the
continued efficient execution of the matrix assembly and linear system solution.
Solve the PDE on T k
Estimate the error on each triangle
while the maximum error on a triangle is larger than the given tolerance do
Based on error estimates, determine a set of triangles, S k , to refine
? Divide the triangles in S k , and any other
triangles necessary to form T k+1
Solve the PDE on T k+1
Estimate the error on each triangle
endwhile
Fig. 1. A framework for the adaptive solution of PDEs
Consider the general adaptive mesh algorithm [10] shown in Figure 1. We
begin by assuming that we have an initial element mesh given by the triangulation
consistent with the geometry of the problem domain. Our attention is focused
on the step in the algorithm where the current mesh T k is adaptively refined (the
step denoted by the ? in Figure 1).
Suppose that some arbitrary subset of triangles, S k , of T k is marked for refine-
ment. We have developed and implemented parallel algorithms for constructing
a new mesh T k+1 that satisfies the required changes in the mesh. To keep our
presentation clear and brief, we assume that the set S k contains triangles marked
only for refinement, not for unrefinement. However, our software is able to unrefine
triangles that have been previously refined.
The refinement of the mesh must maintain several important properties, given
that finite element approximations are used. First, we require that each mesh T k
be conforming (or compatible). That is, the intersection of any two triangles in T k
should be a single vertex, a line segment connecting two vertices, or the empty set.
A side of a triangle is called
-nonconforming if it has s ? 0 vertices between
any two endpoints. A triangle is called compatible if none of its sides are 1
nonconforming. Examples of conforming and nonconforming meshes are given in

Figure

2. If the mesh is conforming, then only one basic type of finite element is
necessary. Otherwise, several special element types are required, and/or a more
complicated matrix assembly. Note, however, that the use of triangles does not
restrict one to linear finite elements; one can use higher-order basis functions in a
triangulation.
Fig. 2. On the left, a conforming mesh; on the right, a nonconforming mesh
A second requirement is that the mesh T k be graded (or smooth). That is,
adjacent triangles should not differ dramatically in area. A nonsmooth mesh could
result in the finite element approximation being very far from the continuous solution

A final requirement is that all angles in the mesh be bounded away from 0
and -. The latter condition is necessary because the discretization error in a finite
element approximation has been shown to grow as the maximum angle approaches
- [1]. We would like to avoid small angles because the condition number of the
matrices arising from mesh elements has been shown to grow as O( 1
' min is the smallest angle in the mesh [4].
2.1. Related Work. A number of mesh refinement algorithms have been
shown to maintain the mesh properties given above. In this section we briefly
review the three most widely used of these refinement methods. To begin, we note
that there are two methods used to subdivide a triangle: bisection and regular
refinement. In bisection, a vertex of the triangle is connected to the midpoint of
the opposite side of the triangle, as in Figure 3, forming two triangles of equal area.
In regular refinement, the midpoints of the sides of the triangle are connected, as
in

Figure

3, to form four similar triangles.
Fig. 3. On the left, a triangle divided with bisection; on the right, a triangle divided by using
regular refinement
The regular refinement algorithm of Bank, Sherman, and Weiser [3] has been
used very successfully in the software package PLTMG [2]. Triangles are divided
by using regular refinement and temporary bisections of selected triangles to make
the mesh conforming. The bisected triangles are merged before the mesh is refined
again. By merging the bisected triangles at each level, the method guarantees
that each triangle in T k+1 either is similar to a triangle in T 0 or is a bisection of a
triangle similar to a triangle in T 0 . Clearly, the angles in T k+1 are bounded away
from 0 and -.
The mesh refinement algorithm 5.6 of Rivara [15] uses bisections of triangles
across the largest edge (dividing the largest angle) and selective divisions across
smaller edges. This approach has been shown to yield triangulations, T k , whose
smallest angle is bounded by at worst one-half the smallest angle in T 0 [17]. A
detailed discussion of this algorithm is given in the following subsection.
The newest-node algorithm of Sewell is also based on bisection, but without
the restriction on bisecting the longest edge [10]. In this algorithm, a triangle
is always bisected by using its newest node. The propagation inherent in the
bisection and regular refinement algorithms is avoided by refining triangles only in
pairs. However, because of the pair restriction, it is possible that a triangle may
never be able to be refined. In the experiments run by Mitchell [10], this difficulty
did not arise.
Mitchell compared these three methods in a series of numerical experiments
and found that it was difficult to choose a consistently superior algorithm [10].
In addition, he found that all three methods were superior to using uniform refinement
except on smooth problems. Given the similar performance of the three
methods, we choose to discuss the bisection algorithm in detail in this paper for
three reasons: (1) it is simpler from an implementation standpoint than the regular
refinement algorithm; (2) it manifests the propagation inherent in both the bisection
and regular refinement algorithms and, therefore, demonstrates the ability of
our algorithm to handle such propagations; and (3) it does not have the potential
for having "unrefinable" nodes as in the newest-node algorithm. We note, how-
ever, that our algorithms are applicable to all three algorithms. In addition, only
a simple modification to our parallel implementation is required to implement the
newest-node algorithm.
We note that Williams [19] has developed a voxel database approach to the
parallel mesh refinement problem. Our approach, which we present in the following
section, differs from his approach in that we have explicit parallel runtime bounds.
In addition, we have designed our approach to yield data structures that are more
suitable to the assembly of the sparse linear systems that arise from these meshes
as well as the solution of these sparse linear systems by sophisticated iterative and
direct sparse factorization methods.
denotes triangles not yet refinedg
children of refined trianglesg
while do
Bisect each triangle in Q i across its longest edge
Bisect each triangle in R i across a nonconforming edge
All incompatible triangles embedded in
are placed in R i+1
All other incompatible triangles are placed in Q
endwhile
Fig. 4. The bisection algorithm
2.2. The Bisection Algorithm. In Figure 4 we present the bisection algo-
rithm. This algorithm is slightly altered, for ease of presentation, from Algorithm
5.6 as presented by Rivara in [15]. However, this modified algorithm yields the
same final mesh as the original algorithm presented by Rivara.
Fig. 5. The process of the bisection algorithm is shown from left to right. In the initial
mesh, the shaded triangles are refined; subsequently the shaded triangles are refined because they
are not compatible.
To illustrate the bisection algorithm, we give an example of the propagation
in

Figure

5. Note that the refinement could propagate through unmarked triangles
not adjacent to marked triangles before finishing. Rivara, however, has shown that
this loop will terminate in a finite number of iterations. We denote this number
of iterations by L P . In general, L P depends on the characteristics of the mesh
being refined. Rivara also has shown that each triangle in T k embeds 1, 2, 3, or 4
triangles of the resulting compatible mesh, T k+1 . We show the possible 2, 3, or 4
resulting triangles in Figure 6. We formalize the following useful result from [15].
Theorem 2.1. During the execution of the bisection algorithm, no side of a
triangle may be divided more than once.
Fig. 6. The possible divisions of a single triangle in the bisection algorithm
Proof: Vertices may be created only when a triangle, as member of Q i , bisects
its longest edge-and then only if a vertex does not already exist in that location.
Once a triangle is in Q i , its children are excluded from Qm for steps m ? i and
thus may not create vertices.
Therefore, if a triangle t a containing the edge e creates a vertex on e at step
no further vertices will be created by that triangle or its children. Given the
creation of the vertex, a triangle t b adjacent to t a and sharing e would not create
another vertex on e when it is a member of some Q In addition, the children
of t b cannot create any vertices. Thus, at most one vertex may be inserted on an
edge during refinement. 2
3. Parallel Adaptive Refinement. In this section, we present a parallel
algorithm for adaptive refinement that correctly implements the bisection method.
We illustrate the key aspect of this algorithm: the synchronization necessary for the
correct parallel execution of the bisection algorithm. Finally, we give an analysis
of the algorithm under the P-RAM computation model.
First, we need the following definitions. Let ng be the
set of vertices in the mesh and mg be the set of polygons.
We assume that the final mesh consists only of triangles (i.e., a conforming mesh).
However, intermediate meshes can be nonconforming, hence we allow for these
nonconforming elements in our definition. Let E) be the graph associated
with the mesh, with edges g. Let be the dual
graph associated with the mesh, where g.
The refinement algorithm will be formulated within the context of the dual
graph. To begin the P-RAM analysis, we assume that at any given time we have
as many processors as we have triangles and that triangle t a is assigned to the
processor p a . For the analysis that follows, the specific P-RAM computational
variant does not make a difference; one may assume that the CREW P-RAM
model is used. Some synchronization must be managed during the execution of
the algorithm to maintain the correct neighbor information in both the graph G
and the dual graph D as they are modified. Thus, each processor p a must keep
track of the current neighbors of t a in D. We note that the correct neighbor
information for G can be constructed in a straightforward way from D.
To illustrate the synchronization required for the correct execution of the parallel
algorithm, we note the two ways that neighbor information can be corrupted.
U
Fig. 7. On the left, two processors creating a vertex at the same location; on the right, a
possible corruption of neighbor information
First, two different processors must not create vertices at the same location when
bisecting their triangles. If two vertices are created at the same location, then a
postprocessing step must be included to merge these vertices. We eliminate the
need for postprocessing by proper synchronization. In Figure 7 we see an example
of two processors, P 1 and P 2 , creating two vertices at the same location. Second,
we must ensure that outdated neighbor information is not propagated. For exam-
ple, in the same figure we see that triangle U 1 may believe that triangle W is its
neighbor, rather than triangle W 1 , if triangles U and W are refined simultaneously.
The key observation is that both of these synchronization problems can be
avoided if only triangles from an independent set in D are refined simultaneously.
An independent set, I, is a subset of triangles of T such that no two triangles from
I are adjacent in D. Once these triangles are refined and neighboring triangles
are notified, another independent set can be chosen for refinement based on the
correctly updated neighbor information. In the following subsection, we consider
two possible approaches for computing these independent sets in parallel.
3.1. Two Methods for Computing Independent Sets in Parallel. For
the purpose of the running time analysis for the refinement algorithm, we review
two approaches for computing the independent sets. Both of these approaches
require that the graph be of bounded degree-which is true for the problem we
consider. The first approach uses an assignment of random numbers to the vertices
of a graph to obtain a sequence of independent sets that is a slowly growing function
of the size of the graph. The second approach is to compute a graph coloring and
use this coloring to generate the independent sets. The advantage of the coloring
approach is that we can guarantee that the number of colors, and thus the number
of independent sets, is independent of the size of the graph for a bounded degree
graph. However, the computation of this coloring requires the use of the first
random number approach; therefore, the coloring is useful only if it is used enough
times to justify the initial expense.
First we consider the use of independent random numbers to generate the
independent sets. Suppose we wish to compute the sequence of independent sets
for the set of triangles T 0
, a subset of T , in the corresponding subgraph D 0
For each triangle t a in D 0
we assign a distinct, independent random number ae(t a ).
We choose an independent set I from T 0
according to the following rule: t a 2 I if for
each of its neighbors t b in D, we have that either (a) t
or (b) ae(t a ) ? ae(t b ). We
then update the set of triangles under consideration by deleting the independent
I. We are now free to generate the next independent set in the
sequence using the same rule. This process continues until T 0
is the empty set.
The expected number of independent sets is given by the following lemma.
Lemma 3.1. Let D 0
be a bounded degree, undirected graph with n vertices.
Suppose each vertex t in D 0
is assigned a unique independent random number ae(t).
Consider the sequence of independent sets generated by the above rule. The expected
number of these independent sets is bounded by EO(log n=log log n).
Proof: This bound is a consequence of Corollary 3.5 in [7]. 2
We now consider a second approach for obtaining the sequence of independent
sets. First, note that any valid coloring of the graph can be used to
generate the required independent sets. Recall that the function oe
is an s-coloring of D, if oe(t a ) . Thus, a sequence
of s independent sets can be generated from an s-coloring of D by assigning all
triangles of the same color to one of the sets.
To efficiently compute this coloring, we use the parallel greedy heuristic presented
in [7]. An outline of this heuristic is presented in Figure 8. The independent
sets required for this heuristic can be generated by using the random number
method described above. The greedy step in the heuristic is the color assignment;
the smallest consistent color for t is the smallest color not assigned to a neighbor
of t.
While T 0
do
Choose an independent set I from T 0
Color I in parallel by choosing the smallest
consistent color oe(t) for each t 2 I
I
Fig. 8. Outline of a parallel greedy coloring heuristic
The advantage of using a coloring to generate the independent sets is that for
a bounded degree graph the maximum number of colors is independent of the size
of the graph. We include this well-known result as the following lemma.
Lemma 3.2. Consider a bounded degree graph, D, of maximum degree \Delta.
The parallel greedy coloring heuristic computes an s-coloring of D with s -
Proof: Every vertex t is colored in the greedy heuristic by assigning it the smallest
consistent color. Since, at worst, every neighbor of t is a different color, the
maximum color assigned t required is the degree of t plus one. Thus, the maximum
color assigned by the greedy heuristic to any vertex in D is \Delta
In sum, we have available two methods for generating the sequence of independent
sets required for the parallel refinement algorithm. For the following P-RAM
running time analysis, it turns out that the best running time bound is obtained
by maintaining a coloring of the dual graph comprising the triangles to be re-
fined. However, in practice, the overhead associated with maintaining the coloring
is not advantageous. Hence, the first approach is used in the practical algorithm
presented in x4.
3.2. A P-RAM Adaptive Refinement Algorithm. In Figure 9 we present
a P-RAM algorithm that avoids the synchronization problems discussed above, by
simultaneously refining triangles from independent sets in D. Note that the independent
sets used for refinement are also used to update the coloring. This update
is required because the dual graph is modified after the bisection of a triangle. In
the remainder of this section, we show that this algorithm avoids the two possible
synchronization problems and has a fast run time.
Based on local error estimates, a set of triangles, Q 0 , is marked for refinement
Each triangle, t j , in Q 0 is assigned a random number, ae(t j )
The subgraph D(Q 0 ) is colored by the parallel greedy coloring heuristic
While do
While do finner loopg
Choose an independent set in D, I, from (Q i
Simultaneously bisect each of the triangles in I
embedded in across its longest edge
Simultaneously bisect each of the triangles in I
embedded in R i across a nonconforming edge
Each new triangle, t j , is assigned the smallest consistent
new processor
Each processor owning a bisected triangle updates this
information on processors owning adjacent triangles
Endwhile
incompatible triangles embedded in
incompatible triangles
Endwhile
Fig. 9. Parallel algorithm for refinement
We assume that the initial dual graph, D, is of bounded degree. In fact,
because the triangulation of a surface is of primary interest, we assume that each
triangle edge is shared by at most two triangles in the initial triangulation. In this
case we have that maximum degree of the initial conforming mesh is three. The fact
that D has bounded degree not only is useful in the following runtime proof, but
also is useful in practice. Design of data structures and software is simplified if the
maximum number of neighbors in the graph is bounded by a small constant. We
now show that degree of any intermediate, nonconforming dual graph is bounded
by at most twice the initial maximum degree.
Lemma 3.3. The dual graph, D, is of bounded degree at all times during
the execution of the algorithm. In fact, the degree of a vertex in D never exceeds
six. As a result, the maximum number of colors required at all times during the
execution of the algorithm to color D is seven or less.
Proof: From Theorem 2.1 each triangle edge is divided at most once. Therefore,
a triangle can at most double the number of its neighbors. By Lemma 3.2, if the
maximum degree of the dual is six, at most seven colors will be required during
the execution of the refinement algorithm. 2
Because the degree of the dual graph remains bounded, we note that the work
assigned to each processor during one pass through the inner loop of this algorithm
can be done in constant time under the P-RAM computational model. We now
show that the two possible corruption problems discussed above cannot occur.
Lemma 3.4. Neighbor information in the dual graph is correctly updated
during the execution of the refinement algorithm.
Proof: The proof is by induction. We assume that the initial neighbor information
is correct and that the neighbor information is correct following step a is
being refined at step i, by the properties of the independent set none of its neighbors
in D are being refined. The triangles, t a 1
and t a 2
, resulting from bisection of t a
have correct information about their neighbors. The former neighbors of t a can,
therefore, be notified of the refinement of t a and be given the correct information
about their new neighbors. Thus, following step i of the refinement algorithm the
modified neighborhood information for D is correct. 2
vertices will be created at the same position during the
execution of the refinement algorithm.
Proof: Again, the proof is by induction. We assume that all vertices are unique
initially and following step 1. For a vertex to be created at the same position
at step i by two different processors, one of two situations must occur: (1) two
processors must simultaneously refine the same edge, or (2) a processor must refine
a previously refined edge because it has not been notified that a vertex has been
created on that edge. The first condition is prevented by the definition of the
independent set-no adjacent triangles are refined simultaneously. The second
condition is prevented by the correct notification of neighbor information in D
ensured by Lemma 3.4. 2
Finally, we give a bound on the expected running time of the refinement algorithm

Theorem 3.6. Recall that L P is the number of loop iterations in the serial
bisection algorithm in Figure 4. The algorithm given in Figure 9 terminates in a
finite number of steps and has an expected run time under the P-RAM computational
model of EO( logjQ 0 j
log logjQ 0 j
Proof: First, we consider the expected running time to compute the initial coloring
of D(Q 0 ). By Lemma 3.1 this time is EO( logjQ 0 j
log logjQ 0 j
Next, we consider the running time for the inner loop of the refinement algorithm
at step i. Define the graph, D is the set of triangles
to be refined at this step. The set F i is the subset of edges (t a ; t b ) from F
with t a ; t b 2 S i . By Lemma 3.3, we know that D i is always a bounded degree graph.
Also by Lemma 3.3 the number colors required, and thus number if independent
sets, is bounded by a constant. Hence the work assigned to any processor in the
inner loop (the bisection of its triangle, updating the coloring, and the neighbor
notification) takes time bounded by a constant independent of the mesh size.
Finally, we must show that it takes L P iterations of the outer loop to form a
conforming mesh. Clearly, every triangle that becomes incompatible at step i is
refined at step i just as in the sequential algorithm in Figure 4. Thus, the
number of iterations of the outer loop in each algorithm is identical, L P .
Hence, the total expected running time for the entire algorithm is bounded by
log logjQ 0 j
We close this section with several notes about this running time analysis.
First, since typically in this context the initial mesh to be refined was obtained
from a previous level of refinement, the initial coloring step would not be required.
Instead, a coloring of the mesh could be maintained between levels of refinement.
Using this information, the P-RAM running time of the algorithm would be O(L P ).
Fig. 10. On the left, the shaded triangle is marked for refinement; on the right, the resulting
conforming mesh after refinement. Note that the refinement has propagated through every triangle
in the mesh but one.
Finally, we close the analysis with some comments about L P , the length of
propagation. We point out that it is easy to construct a worst case example where
j). For example, in Figure 10 we give an example that
can be generalized to illustrate this worst case behavior. However, another way of
looking at this example is to assume that this particular mesh was generated from
previous refinements starting from a single triangle. In this case, the average length
of propagation, L P , over all levels of refinement, is actually constant. Furthermore,
as we note with the experimental results presented in x5, the average number of
independent sets required to obtain a conforming mesh appears to bounded by a
small constant and independent of the size of the mesh. Thus, we believe that this
possible worst case behavior of L P is not the ominous problem that it appears it
could be in a practical implementation.
4. Distributed-Memory Implementation. For use on a practical parallel
computer, we must modify the P-RAM algorithm analyzed in the preceding section.
Rather than assigning a single triangle or vertex to each processor, we assign a set of
vertices and triangles to each processor. The vertices V are partitioned into disjoint
subsets owns the subset V j , and we have that
We choose to partition the vertices rather than the triangles because we have found
that it makes the finite element evaluation, mesh refinement, and sparse matrix
assembly and solution (if necessary) more straightforward and efficient. Based on
the partitioning of V , we determine a partitioning of
disjoint
subsets where processor j owns the subset T j . In practice, one can assume that at
least one vertex of triangle in T j is in the set V j .
For communication purposes, each processor, j, stores the set of triangles
is the set of triangles adjacent to a
triangle in T j in the dual graph D, and T is the set of triangles containing a
vertex in V j . In addition, processor j stores the set of vertices -
is the set of vertices contained by all triangles in -
Given the sets -
processor j has all the information necessary to
evaluate all finite elements that have vertices in V j , assemble complete rows and/or
columns of a sparse matrix associated with each vertex in V j , and perform the
parallel refinement algorithm (yet to be specified) on the triangles in T j . We
illustrate these sets for some processor j in Figure 11. In this figure we have
partitioned the vertices by the geometric cuts represented by the orthogonal dashed
lines. The vertices in the interior of the four dashed lines have been assigned to
processor j-the set V j -and are shown as filled vertices. The set -
is the set of
unfilled and filled vertices, T j is the set of shaded triangles, and -
is the set of
unshaded and shaded triangles.
Fig. 11. An illustration of the sets -
and -
maintained on processor j. The set of vertices
assigned to processor j, V j
, is shown as the set of filled vertices. The set of shaded triangles is
the set T j
. The union of the set of shaded and unshaded triangles is -
; the union of the set of
filled and unfilled vertices is -
In

Figure

12 we present a practical version of the P-RAM algorithm given in

Figure

9. The algorithm ensures that vertices are not created at the same location
and that the sets -
T j on each processor j are correct. Note that, in this
modified algorithm, if a triangle or vertex is created on processor j, processor
j is its owner. The inner and outer loops in the P-RAM algorithm have been
combined into a single loop for greater efficiency; the separate loops allowed for
a clearer presentation of the runtime bounds, but that is not necessary in this
section.
Based on local error estimates an initial set of triangles, Q,
is marked for refinement
Each triangle, t a 2 Q, is assigned a random number, ae(t a )
While (Q [ R) 6= ; do
Choose an independent set in D, I
, from the
triangles in (Q [ R), where I
Each processor, j, bisects the triangles in I j
embedded in Q across its longest edge
Each processor, j, bisects the triangles in I j
embedded in R across a nonconforming edge
For each new triangle, t b , a new random number, ae(t b ), is chosen
Each new triangle, t b , created on processor j is added to T j
Each new vertex, v k , created on processor j is added to V j
For each triangle, t a 2 I j notification of bisection is sent to each processor
l for which ((adj D (t a
Each processor receives notification and updates its sets -
incompatible triangles embedded in W
incompatible triangles not in R
Endwhile
Fig. 12. A practical parallel algorithm for refinement
For this algorithm, independent sets are chosen according to a slightly different
rule from the rule used in the P-RAM algorithm. The triangle t a is in I j if, for each
of its neighbors t b in D, one of the following hold: (a) t
or (c) ae(t a ) ? ae(t b ). This modification allows two triangles on the same processor
to be refined on the same step. The computation of the independent sets requires
no communication because each processor has all the necessary information in -
for this computation. Communication of the random numbers is not necessary if
the seed given the pseudo-random number generator used to determine ae(t a ) is
based solely on a. Thus, the only communication necessary in the algorithm is the
notification of bisections and the global reduction required to determine whether
;. For further efficiency, the notification messages can be packed so
that each processor receives at most one message from another processor during
each time through the while loop.
Because the modified algorithm in Figure 12 uses essentially the same synchronization
scheme presented in x3, collisions are avoided, and neighbors in D on
separate processors are not simultaneously bisected. Thus, we have the following
theorem.
Theorem 4.1. All changes made by other processors to the triangles/vertices
in the sets -
V j on each processor j are received so that these sets are kept
updated throughout the algorithm in Figure 12.
Proof: The sets T j and V j are updated correctly because only processor j can
bisect triangles in this set or create new vertices in this set. Any changes to
triangles or vertices in -
attributable to changes in triangles or vertices in
are directly communicated in the algorithm.
The remaining portion of -
attributable to T (V j ), is accounted for because,
is bisected on another processor, then notification
of this bisection will be sent and received.
Finally, to show that the vertex neighbor information is correct, we note that
the neighbor information is correct on the subgraph of D induced by -
the neighbor information contained in the subgraph of G induced by -
must also
be correct because -
5. Experimental Results. In this section computational results are presented
that demonstrate that the parallel refinement algorithm is scalable and
that its execution time is negligible compared with that of other computations
required to solve a PDE.
The parallel refinement algorithm is implemented as a subroutine library that
can be called by an application program. Chameleon [5] is used to achieve portability
across several architectures, including the Intel DELTA, which is the focus
of this section. Note that in addition to the refinement algorithm, the subroutine
library also includes a similarly constructed, parallel unrefinement algorithm.
Because the unrefinement algorithm is necessary in many applications, including
one of those used here, and its performance is similar to the refinement algorithm,
results from it are included here as well. Results are presented for the parallel
refinement algorithm for two different two-dimensional PDEs: Poisson's equation
and the equations for linear elasticity. These problems are solved on two different
geometries.
5.1. Test Problems. Our first set of test problems models Poisson's equa-
tion
where S is a square domain and a linear finite element approximation is used. The
function f(x) is a Gaussian charge distribution centered at a point (S x
the domain. The mesh is selectively refined according to the energy norm [10] until
the estimate of the local error on each triangle is less than a specified tolerance.
Further, the point (S x ; S y ) is moved several times, and a new solution/mesh is
found from the old solution/mesh. This movement requires significant mesh refinement
around the new charge position and definement around the old position
while the remainder of the mesh remains relatively constant. The parallel conjugate
gradient method preconditioned by an incomplete factorization is used to
solve the sparse linear systems that arise [6].
The second problem considered is the linear elasticity equations for the plane
stress problem, given (without inclusion of a load) as
@x@y
@x@y
where u and v are the x and y displacements, respectively. These equations are
solved on a rectangular region with a central hole. One side of the region is constrained
to have zero displacement, and a constant traction is applied to the opposite
side. Again, linear finite elements are used to approximate these equations.
The mesh is selectively refined according to the energy norm until the local error
estimate for each triangle is less than a specified tolerance. The linear systems are
solved by using the same code used for the Poisson problem.
In each problem set, the initial coarse mesh has approximately 200 nodes
except when running on 128 and 256 nodes; in these cases the initial coarse meshes
are approximately 2 and 4 times larger, respectively. For each of these problems,
by carefully choosing the maximum tolerance for the local error estimator, one can
determine the maximum number of vertices in the solution meshes. The following
two problem sets have been constructed such that the final solution mesh for each
successive problem has roughly twice as many vertices/triangles as in the previous
problem. Information about the two problem sequences is given in Tables 1 and 2.

Table
A sequence of test problems based on the Poisson problem
Maximum Number Maximum Number Ratio of Area of
of Vertices in the of Triangles in the Largest Triangle to
Name Adaptive Mesh Adaptive Mesh Smallest Triangle
POISSON1 2,673 5,268 256
POISSON2 5,176 10,260 512
POISSON3 10,238 20,330 512
POISSON4 20,296 40,412 1,024
POISSON6 80,116 159,872 2,048
POISSON7 159,758 318,948 2,048
A sequence of test problems based on the linear elasticity equations for the plane stress problem
Maximum Number Maximum Number Ratio of Area of
of Vertices in the of Triangles in the Largest Triangle to
Name Adaptive Mesh Adaptive Mesh Smallest Triangle

Tables

3 and 4 give the number of refinement steps required for each problem
during the solution process. A refinement step consists of finding an approximate
solution to the PDE on the current mesh, T k , by solving the sparse linear system
arising from the finite element model, computing estimates for the local error
at each triangle, and then refining T k according to these estimates to obtain the
conforming mesh, T k+1 . One observes that, not unexpectedly, it takes more mesh
refinement steps to construct the larger meshes. In addition, the number of iterations
through the loop in the algorithm in Figure 12 is given. The number of
iterations should be at least L P and perhaps a slowly growing function of the mesh
size because we use the random number rule to generate the independent sets.
One notes that the number of loop iterations needed is a slowly growing function
of the number of processors and problem size. This result indicates that one can,
in general, achieve scalable performance, as may be expected from Theorem 3.6.
For the POISSON problem set, the charge location was moved twice; this
movement meant that at two solution steps the mesh was not only refined around
the charge, but also unrefined around the old charge position. However, there
were still more refinement operations/steps than unrefinement operations. No
unrefinement was necessary in the ELASTIC problem set; the load function was
unchanged.
A good partitioning of the vertices for each of these problems is necessary
for the new algorithm to perform efficiently. Many good partitioning methods
are available; a geometric partitioning algorithm [8] was chosen for this work.

Figure

13 shows the average number of partitions that are adjacent to a given
partition. This information gives some sense of the number of processors each
processor shares triangles with and must, therefore, exchange information with.

Figure

14 shows the percentage of the total triangle edges that have endpoints on
two different processors. This data gives some sense of the number of triangles each
Number of refinement steps and loop iterations for the sequence of Poisson test problems
Number of Number of Average Number
Name Processors Refinement Steps of Loop Iterations
POISSON5
POISSON6
POISSON8 128 24 2.46

Table
Number of refinement steps and loop iterations for the sequence of linear elasticity test problem

Number of Number of Average Number
Name Processors Refinement Steps of Loop Iterations
processor has that must be coordinated with another processor. Note that these
values initially rise rapidly, as one would expect, until approximately processors
are in use. For larger numbers of processors, these values increase very slowly.
5.2. Experiments. The experiments were run on up to 256 nodes of the
Intel DELTA. The parallel computer is a mesh-connected, 16 \Theta 32 array
of Intel i860 microprocessors. 1 In all of the experiments, the reported times are
given in seconds. The operations rates indicate the number of bisections and vertex
deletions (note that vertex deletions correspond to unrefinement and constitute a
small percentage of the total) per second.
1 Note that because of constraints on the amount of time available to us on the DELTA, the
512-processor case was not run. We believe, however, that the results convincingly demonstrate
the effectiveness of the parallel refinement algorithm.
POISSON
Number
of
Neighbors6Number of Processors
Fig. 13. The average number of partitions each partition is adjacent to in the final mesh
POISSON
100 2003Number of Processors
Percentage of
Cross-Edges
Fig. 14. The percentage of cross-edges in the final mesh
To demonstrate the scalability of the new algorithm and implementation, we
designed the problems from each test set to have nearly equal numbers of vertices
from the final mesh assigned to processors. This fact can be seen in Tables 1
and 2, which show how many processors each problem was run on and the size of
the final meshes. Each of the test problems is refined in localized regions of the
mesh; therefore, some processors have more refinement work than others. This
load imbalance is reflected in Tables 5 and 6, which give the average number
of operations per processor per step and the average of the maximum number
of operations on a single processor per step. The average number of operations
falls as the number of processors increases; this decrease results because more
refinement steps are taken to achieve the same number of vertices per processor in
the final mesh. The average maximum number of operations increases because, as
the mesh size increases, more refinement is concentrated in the same size area in
which a limited number of processors are working. Recall that the entire mesh is
repartitioned after each refinement step. Thus, this concentration of new elements
is continuously redistributed to processors with fewer elements.
However, even given these handicaps, the results demonstrate that the algorithm
performs quite well. Figure 15 shows the average number of refinement
operations per second per processor as a function of the number of processors. If
The number of refinement operations for the Poisson problem sequence
Average Number of Average of the
Number of Operations per Maximum Number of
Name Processors Processor Operations per Processor
POISSON5
POISSON6
POISSON8 128 111 573
POISSON
100 200
Number of Processors300700
Number of
triangles
per second
Fig. 15. The average number triangles refined per processor per second
refinement were occurring uniformly on all processors, one could expect this rate to
be nearly constant; however, in the test problems, as in most practical problems,
this is not the case. Figure 16 shows a more interesting rate, the maximum number
of refinement operations per second on an individual processor. One would expect
this rate to remain constant, or nearly so, if the algorithm is perfectly scalable.
With the POISSON problem set, one sees very little degradation in the maximum
refinement rate. One might expect some degradation resulting from the
increasing number of neighbors each processor must exchange information with
as the number of processors increases. However, with the POSSION problem set
the increase in the number of neighbors is offset by the rapidly increasing maximum
number of operations per processor (given in Table 5). With the ELASTIC
problem set one observes this expected degradation because the maximum number
of operations per processor is increasing only moderately. Prior to reaching
processors, the maximum rate of refinement is rapidly changing because of the
The number of refinement operations for the linear elasticity problem sequence
Average Number of Average of the
Number of Operations per Maximum Number of
Name Processors Processor Operations per Processor
POISSON
100 200
Number of Processors300700
Number of
triangles
per second
Fig. 16. The maximum number of triangles refined per second on an individual processor
increasing communication requirements as the the number of processor neighbors
and the percentage of cross-edges increases. After reaching 16 processors, the
number of processor neighbors and the percentage of cross-edges stabilizes, and
one sees approximately a 20% degradation in the rate of refinement from 16 to 256
processors.
Results given in Figure 17 demonstrate that, for a reasonably complex set of
problems, the time to solve the linear systems dominates the time to refine the
mesh for any number of processors. In fact, the total refinement time is always
less than 4 percent of the total execution time. Note that a linear system is solved
after each level of refinement. So, for example, the total execution time shown
for 128 processors in Figure 17 includes the assembly and solution of 13 sparse
linear systems. The time represented by the white region in the bar graph is
composed almost entirely of the sum of the times required for the repartitioning of
the mesh after each level of refinement. This partitioning time includes the time
Number of Processors
Time (sec)
Execution
Total
Matrix Solution Time
Refinement Time
Other (Partitioning, etc.)
Fig. 17. A comparison of the refinement times with the times required for all other aspects
of the problem solution for the linear elasticity problem sequence
to move vertices, triangles, and the data associated with them between processors.
We note that the implementation of our partitioning heuristic is preliminary; we
believe that these times can be significantly reduced.
To examine the total running time in more detail, we consider one problem,
ELASTIC9, run on 256 processors. In Figures we show the time required
to solve the linear system and the number of nonzeros in the assemblied matrix
as a function of the refinement level. Initially the matrix size is doubling after
every level of refinement, since most triangles are bisected at each refinement step.
However, for the last several refinement levels only small areas of the mesh are
being refined. As a result, the interpolated solution from the previous mesh is an
excellent initial guess to the solution on the refined mesh, and only a small number
of conjugate gradient iterations are required to obtain a solution that satisfies the
specified tolerance for the relative residual.
For the same problem, ELASTIC9 run on 256 processors, we show in Figure 19
the time required to refine the mesh as a function of the refinement level. In
the figure we show the number of vertices that have been added to the mesh at
that level of refinement. Note that refinement time continues to increase after the
number of vertices reaches a maximum. This effect can be explained by noting that
the areas of the computational domain on which refinement is occurring become
confined to fewer processors as the mesh is refined. Recall from Figure 16 that it
is the maximum rate of refinement on a processor that is constant. Thus, because
the refinement is occurring on a smaller number of processors, the average rate of
refinement is worse at the higher levels of refinement This behavior explains the
Linear
System
Solution
Time
SOLUTION TIME
Matrix
(millions)
Refinement Level
Fig. 18. The time required to solve the linear systems and the number of nonzeros in the
assembled linear system at each level of refinement for problem ELASTIC9 run on 256 processors
of the Intel
VERTICES ADDED(sec)
Total
(thousands)
Created
Vertices
Refinement
Refinement Level
Fig. 19. The refinement times and the number of vertices added at each level of refinement
for problem ELASTIC9 run on 256 processors of the Intel
decrease in the average rate of refinement as a function of the number of processors
as shown in Figure 15.
6. Concluding Remarks. We have described a parallel algorithm for the
adaptive refinement of meshes. This algorithm was shown to run in provably fast
time under a P-RAM model of computation. In addition, we described an efficient
method of implementation for this algorithm on a practical, distributed-memory
parallel computer. We then gave results for two problems that demonstrate the
scalable nature of this algorithm.
The results given in this paper are for a two-dimensional triangular mesh.
The use of independent sets for parallel synchronization, however, generalizes to
the three-dimensional case as well as other refinement algorithms. The next logical
step in this work is to develop theoretical results for three-dimensional tetrahedralizations
as well as a practical, parallel implementation for three dimensions. In
addition, we note that the use of higher-order basis functions is straightforward
in this methodology; in fact, we include this functionality in the current parallel
implementation [9].



--R


PLTMG: A Software Package for Solving Elliptic Partial Differential Equa- tions
Refinement algorithms and data structures for regular local mesh refinement
Condition of finite element matrices generated from nonuniform meshes
Users Manual for Chameleon Parallel Programming Tools




A comparison of adaptive refinement techniques for elliptic problems
An algorithm for adaptive refinement of triangular element meshes
On a data structure for adaptive finite element mesh refinements
Algorithms for refining triangular grids suitable for adaptive and multigrid techniques



A lower bound on the angles of triangles constructed by bisecting the longest side
A finite element program with automatic user-controlled mesh grading
A dynamic solution-adaptive unstructured parallel solver
--TR

--CTR
J. P. Surez , P. Abad , A. Plaza , M. A. Padrn, Computational aspects of the refinement of 3D complex meshes, Proceedings of the international conference on Computational methods in sciences and engineering, p.615-618, September 12-16, 2003, Kastoria, Greece
Jorn Behrens, Adaptive Atmospheric Modeling: Scientific Computing at Its Best, Computing in Science and Engineering, v.7 n.4, p.76-83, July 2005
Miguel A. Padrn , Jos P. Surez , ngel Plaza, Adaptive techniques for unstructured nested meshes, Applied Numerical Mathematics, v.51 n.4, p.565-579, December 2004
parallel algorithms;distributed memory computers;unstructured mesh computation;sparse matrices;adaptive mesh refinement