14

Data Structures for

Unstructured Mesh
Generation

14.1 Introduction
14.2 Some Basic Data Structures
Linear Lists • A Simple Hash Table
14.3 Tree Structures
Binary Trees • Heaps • Binary Search Tree • Digital Trees
14.4 Multidimensional Search
Searching Point Data • Quadtrees • Binary Trees for
Multidimensional Search • Intersection Problems
Luca Formaggia 14.5 Final Remarks

14.1 Introduction
The term data structures, or information structures, signifies the structural relationships between the
data items used by a computer program. An algorithm needs to perform a variety of operations on the
data stored in computer memory and disk; consequently, the way the data is organized may greatly
influence the overall code efficiency.
For example, in mesh generation there is often the necessity of answering queries of the following
kind: give the list of mesh sides connected to a given node, or find all the mesh nodes laying inside a
certain portion of physical space, for instance, a sphere in 3D. The latter is an example of a range search
operation, and an inefficient organization of the node coordinate data will cause the looping over all
mesh nodes to arrive at the answer. The time for this search operation would then be proportional to
the number of nodes n, and this situation is usually referred to by saying that the algorithm is of order
n, or more simply O(n). We will see later in this chapter that a better data organization may reduce the
number of operations for that type of query to O(log2 n), with considerable CPU time savings when n
is large.
The final decision to use a certain organization of data structure may depend on many factors; the
most relevant are the type of operations we wish to perform on the data and the amount of computer
memory available. Moreover, the best data organization for a certain type of operation, for instance
searching if an item is present in a table, is not necessarily the most efficient one for other operations,
such as deleting that item from the table. As a consequence, the final choice is often a compromise. The
fact that an efficient data organization strongly depends on the kind of problem at hand is probably the
major reason that a large number of information structures are described in the literature. In this chapter,
we will describe only a few of them: the ones that, in the author’s opinion, are most relevant to
unstructured mesh generation. The reader interested in a more ample surveys may consult specialized

©1999 CRC Press LLC

texts, among which we mention [10, 2] for a general introduction to data structures and related algo-
rithms, and [11, 20, 17] for a more specific illustration of range searching and data structures relevant
to computational geometry.
It is a commonly held opinion that writing sophisticated data structures is made simpler by adopting
programming languages that allow for recursion, dynamic memory allocation, pointer and structure data
types. This is probably true, and languages such as C/C++ are surely among the best candidates for the
purpose. However, all the data structures presented in this work may be (and indeed they have been)
implemented in Fortran, and a Fortran implementation is often more cumbersome but normally not
less efficient than the best C implementation. I am not advocating the use of Fortran for this type of
problem — quite the contrary — but I wish to make the point that also Fortran programs may well
benefit from the use of appropriate information structures.
This chapter is addressed to people with a mathematical or engineering background, and only a limited
knowledge of computer science, who would like to understand how a more effective use of data structures
may help them in developing or improving a mesh generation/adaption algorithms. Readers with a strong
background in computer science will find this chapter rather trivial, apart from possibly the last section
on multidimensional searching.

14.2 Some Basic Data Structures

In this section we review some basic data structures. It is outside the scope of this work to give detailed
algorithmic descriptions and analysis. We have preferred to provide the reader with an overview of some
of the information structures that may be profitably used in mesh generation/adaption procedures,
together with some practical examples, rather than to delve into theoretical results.
First, some nomenclature will be given. A record R is a set of one or more consecutive memory locations
where the basic pieces of information are kept, in separate fields. Many authors use the term node instead
of record. We have chosen the latter to avoid possible confusion with mesh nodes. The location &R of
record R is the pointer to the position where the record is stored, while *p indicates the record whose
location is p. In the algorithm descriptions, I will use a C-like syntax, for example i++ is equivalent to i =
i+1. Moreover, with the expression A.b we will indicate the attribute b, which can be either a variable or
a function, associated to item A. For instance R.f may indicate the field f of the record R.

14.2.1 Linear Lists

Quoting from Knuth [10] “A linear list is a set of records whose structural properties essentially involve
only the linear (one-dimensional) relative position of the records.” In a list with n records, the record
positions can be put in a one-to-one correspondence with the set of the first n integer numbers, so that
we may speak of the ith element in the list (with 1 < i < n), which we will indicate with L.i. The type of
operations we normally want to perform on linear lists are listed in the following.
Record Access RA(k). This operation allows the retrieval of the content of a record at position k.
Record Insert RI(r,k). After this operation the list has grown by one record and the inserted record will
be at position k. The record previously at position k will be at k + 1 and the relative position of
the other records remains unchanged.
Record Delete RD(k). The kth record is eliminated, all the other records remain in the same relative
position.
14.2.1.1 Stacks and Queues
Linear lists where insertion and deletion are made only “at the end of the list” are quite common, so
they have been given the special name of deques, or double-ended queues. Two special types of deques
are of particular importance: stacks and queues.

©1999 CRC Press LLC

 With stack, or LIFO list, it is indicated a linear list where insertion, deletion, and accesses are made
only at one end. For example, a list where the operations allowed are RA(n), RI(n+1), and RD(n), i.e.,
all the operations made on the last list position, is a stack. The insert operation is often called a “push,”
while the combination of RA(n) and RD(n) is referred to as a “pop” operation.
In a queue, also called FIFO list, the elements are inserted at one end, and accessed and deleted at the
other end. For instance, a linear list where only RD(1), RA(1), and RI(n+1) operations are allowed is a
queue.
The stack is a very common data structure. It occurs every time we wish to “accumulate items” one
by one and then retrieve them in the inverted order. For instance, when in a triangulation process we
are searching the nodes that lie inside a sphere, every time a new node is found it may be pushed onto
a stack. At the end of the search, we may “pop” the nodes from the stack one by one. We have so far
identified a linear list by its properties and the set of operations that may be performed on it. Now, we
will investigate how a linear list could be actually implemented, looking in some detail at the implemen-
tations based on sequential and linked allocation.
14.2.1.2 Sequential Allocation
The method of sequential allocation is probably the most natural way of keeping a linear list. It consists
in storing the records one after the other in computer memory, so that there is a linear mapping between
the position of the record in the list and the memory location where that record is actually stored. With
sequential allocation, direct addressing is, therefore, straightforward. A sequentially allocated list broadly
corresponds to the ARRAY data structure, present in all high-level computer languages. In the following,
we use the C convention that the first element in array A is A[0].
As an example, let us consider how to implement a stack using sequential allocation. One possibility
is to store the stack S in a structure formed by two integers. S.max and S.n indicating the maximum and
the actual number of records on the stack, respectively, and an array S.A[max] containing the records.
Unless we know beforehand that the program will never try to store more than S.max records on the
stack, we need to consider the possibility of stack overflowing. When such condition occurs, we could
simply set an error indicator and exit the push function. A more sophisticated approach would consider
the possibility of increasing the stack size. In that case, we will probably store an additional variable sgrow
indicating how much the stack should increase if overflow occurs. In that situation, we could then allocate
memory for an array sized S.max + sgrow, adjourn S.max to the new value, and move the old array on the
new memory location. We must remember to verify that there are enough computer resources available for
the new array. If not, we have a “hard” overflow and we can only exit the function with an error condition.
We have just considered the possibility of letting the stack grow dynamically. What about shrinking it
when there is a lot of unused space in S? We should first decide on a strategy, in order to avoid growing
and shrinking the stack too often, since these may be costly operations. For instance, we could shrink
only when S.max – S.n > 1.5sgrow. The value of sgrow may itself be a result of a compromise between memory
requirement and efficiency. A too small value could mean performing too many memory allocation/deallo-
cations and array copying operations. Too large a value will imply a waste of memory resources. We will not
continue this discussion further. We wanted only to show how, even when dealing with a very simple
information structure such as stack, there are subtle details that could be important for certain applications.
The sequential implementation just described may be readily modified to be used also for a general
double-ended queue Q. Figure 14.1 shows how this may be done. We use an array Q.A[max], plus the
integer quantities Q.n, Q.max, Q.start, and Q.end, respectively, indicating the actual and maximum
number of records in the deque and the position of the initial and final record in the array. In Table 14.1,
we illustrate the algorithms for the four basic operations, RI(1), RI(n+1), RD(1), and RD(n). As a matter
of fact, Q.n is not strictly necessary, yet it makes the algorithms simpler. When an overflow occurs we
may decide to grow the structure by a given amount, and the same considerations previously made for
stacks apply here.

©1999 CRC Press LLC

 FIGURE 14.1 How to implement a double-ended queue (deque) using sequential allocation.

TABLE 14.1 An Example of Algorithms for Inserting and Removing

Records from the Sequentially Allocated Deque Illustrated in Figure 14.1.
Deque::RI (R,1) Deque::RI (R,n+1)

1. (n > max) a OVERFLOW; 1. (n > max) a OVERFLOW;

2. start = (start + max – 1) mod max; 2. end = (end + 1) mod max;
3. A [start] = R; 3. A [end] = R;
4. n + +. 4. n + +.

Deque::RD (1) Deque::RD(n)

1. n = 0 a UNDERFLOW; 1. n = 0 a UNDERFLOW;
2. start = (start + 1) mod max; 2. end = (end + max – 1) mod max;
3. n – –. 3. n – –.

FIGURE 14.2 Adding a record in a sequentially stored list.

14.2.1.3 Linked Allocation

What happens if we have to add a record at a random location inside a sequentially stored list? Figure 14.2
graphically shows that we should move “in place” a slice of the array. This procedure requires, in general,
O(n) operations and therefore it should be avoided. In order to increase the flexibility of a linear list by
allowing an efficient implementation of random record insertion and deletion, we need to change the
way the structure is implemented. This may be done by adding to each record the link to the next record
in the list. For instance, we could add to R a field R.next, containing the location of the successive record.
A list which uses this type of layout is called a linked list, or, more precisely a singly linked list. There are,
in fact, many types of linked lists. If we have also the link to the previous record R.prev, we have a doubly
linked list that permits enhanced flexibility, as it allows one to sequentially traverse the list in both
directions and to perform record insertions in O(1) operations. Figure 14.3 illustrates an example of a
singly and of a doubly linked list.
A disadvantage of linked allocation is that direct addressing operations are costly, since they require
traversing the list until the correct position is reached. Moreover, a linked list uses more memory space
per record than the corresponding sequential lists. However, in many practical applications direct address-
ing is not really needed. Furthermore, with a linked list it is normally easier to manage the memory
requirements dynamically and to organize some sharing of resources among different lists.

©1999 CRC Press LLC

 FIGURE 14.3 An illustration of a singly and of a doubly linked list.

TABLE 14.2 Algorithmic Implementation of the Addition and Deletion

of a Record from a Doubly Linked Circular List
Dcllist::RI (R, Q). Dcllist::RD (R).
Insert record R in list, after record Q Delete record R from list

1. p = Q.next; 1. r = R.prev;
2. R.next = p; 2. p = P.next;
3. (*p).prev = &R; 3. (*r).next = p;
4. Q.next = &R; 4. (*p).prev = r.
5. R.prev = &Q.

It is often convenient to use a variant of the linked list, called a circular linked list. In a circular (singly
or doubly) linked list every record has a successor and a predecessor and the basic addition/deletion
operation has a simpler implementation. There is also usually a special record called header that contains
the link to the first record in the list, and it is pointed to by the last one. Table (14.2) shows a possible
algorithm for the implementation of the basic addition/deletion operations on a circular doubly linked
list L. The memory location for a new record could be dynamically allocated from the operating system,
where we would also free the ones deleted from the list. However, this type of memory management
could be not efficient if we expect to have frequent insertions and deletions, as the operations of allocating
and deallocating dynamic memory have a computational overhead. Moreover, it cannot be implemented
with programming languages that do not support dynamic memory management. It is then often
preferable to keep an auxiliary list, called list of available space (LAS), or free list, which acts as a pool
where records could be dumped and retrieved. At start-up the LAS will contain all the initial memory
resources available for the linked list(s). The LAS is used as a stack and is often singly linked. Here, for
sake of simplicity, we assume that also the LAS is stored as a doubly linked circular list. Figure 14.4 shows
graphically an example of a doubly linked circular list and the corresponding LAS, plus the operation
required for the addition of a record. In the implementation shown in the table we have two attributes
associated with a list L, namely L.head, and L.n, which gives the location of the header and the number
of records currently stored in the list, respectively. Consequently, LAS.n indicates the number of free
records currently available for the linked list(s). In Table (14.3) we illustrate the use of the LAS for the
insert and delete operation. We have indicated with R.cont the field where the actual data associated with
R is kept.
It remains to decide what to do when an overflow occurs. Letting the list grow dynamically is easy:
we need to allocate memory for a certain number of records and join them to the LAS. The details are
left to the reader. If we want to shrink a linked list we can always eliminate some records from the LAS
by releasing them to the operating system. Again, we should take into account that many fine grain
allocations/deallocations could cause a considerable computational overhead, and a compromise should
be found between memory usage and efficiency. We have mentioned the possibility that the list of available
storage could be shared among many lists. The only limitation is that the records in all those lists should
be of equal size. Linked lists may be implemented in languages, such as Fortran, that do not provide
pointer data type. Pointers would be substituted by array indices, and both the linked list and the LAS
could be stored on the same array. The interested reader may consult [2] for some implementation details.

©1999 CRC Press LLC

FIGURE 14.4 An example of a doubly linked circular list and of the associate list of available storage. The operations
involved in the addition of record D after position Q are graphically illustrated.

TABLE 14.3 Record Addition and Deletion from a Doubly Linked Circular
List, Using a List of Available Space for Record Memory Management
Insert data x in list L in a record placed
after record Q Delete R from list L

1. LAS.n = 0 → OVERFLOW; 1. n = 0 → UNDERFLOW;

2. p = (LAS.head).next ; 2. RD(R);
3. R = *p; 3. n– –;
4. LAS.RD(R); 4. Q = *LAS.head ;
5. LAS.n – – ; 5. LAS.RI(R,Q);
6. R.cont = x; 6. LAS.n ++.
7. RI(R,Q)
8. n ++.

14.2.2 A Simple Hash Table

It may be noted that searching was not present among the set of operations to be performed on a linear
list. This is because linear lists are not well suited for this type of application. We will now introduce a
data structure used in unstructured grid generation and grid adaption procedures and that is better
designed for simple search queries.
Let’s first state in a general form the problem we wish to address. Let us assume that we need to keep
in a table H some records that are uniquely identified by a set of keys K = {k1, k2, K, kl} and let us indicate
with R.ki the ith key associated to R, respectively. The type of operations we want to perform on H are
as follows:

©1999 CRC Press LLC

 1. Search if a record with given keys is present in the table;
2. Add a new entry to the table;
3. Delete an entry from the table.
A possible implementation that may allow efficiently solving the problem is to consider one of the keys
as the principal key. Without loss of generality, we assume that the first key k1 is the principal key, and
in the following it will be simply referred to as k. Let U be a set of keys and k ΠU a generic key of the
set. We now build a function h(k), called hashing function,* h(k): U → { 0, …, m – 1 } , that assigns to
each key of U an integer number between 0 and m – 1. We have various ways of building the hashing
table H, depending whether h is a one-to-one mapping or not. However, before proceeding further, let
us consider a practical example.
Assume that we want to keep track of the triangular faces of a 3D tetrahedral mesh, when the mesh
layout is constantly changing, for example during a mesh generation or adaption process. A face F could
be identified by its node numbering {i1, i2, i3}, and to make the identification unique we could impose that

k ≡ k 1 = min (i 1, i 2, i 3), k 3 = max (i 1, i 2, i 3), k 2 = {i 1, i 2, i 3} – {k 1, k 2}

Since k is an integer number, and we expect that the k’s will be almost uniformly distributed, a simple
and effective choice for h(k) is the identity function h(k) = k. A hash table H may then be formed by an
array H.A[m], where m is the maximum node numbering allowed in the mesh. The array will be addressed
directly using the key k. Each table entry H.A[k] will either contain a null pointer, indicating that the
corresponding key is not present in the table, or a pointer to the beginning of a linked list whose records
contain the remaining keys for each face with principal key k, plus possible ancillary information. If we
use doubly linked lists, each entry in the table may store two pointers, pointing to the head and the tail
of the corresponding list, respectively. In practice, each array element acts as a header for the list.
Figure 14.5 illustrates this information structure, where, for simplicity a 2D mesh has been considered.
If we use doubly linked list, add and delete operations are O(1) while simple search is a O(1 + n/m)
operation, where n indicates the number of records actually stored in the table. Since in many practical
situations, such as the one in the example, the average value of n/m is relatively small (≈ 6 for a 2D
mesh), the searching is quite efficient. As for memory usage, if we assume that no ancillary data is stored,
we need approximately 2mP + nmax[2P + (l – 1)I] memory locations, where P and I are the storage
required for a pointer and an integer value, respectively, while l is the number of keys (3 in our case),
and nmax is the maximum number of records that could be stored at the same time in the linked lists. In
this case, nmax is at most the maximum number of faces in the mesh. All chained lists have records of the
same size, therefore a common LAS is normally used to store the available records. Some memory savings
may be obtained by storing the first record of each chained list directly on the corresponding array
element, at the expense of a more complicated bookkeeping.
The structure just presented is an example of a hash table with collision resolved by chaining. The
term collision means the event caused by two or more records that hash to the same slot in the table.
Here, the event is resolved by storing all records with the same hash function in a linked list. This is not
the only possibility, however, and many other hashing techniques are present in the literature, whose
description is here omitted. In the previous example we have assumed that we know the maximum
number of different keys. What can be done if we do not know this beforehand, or if we would like to
save some memory by having a smaller table? We need to use a different hash function. There are many
choices for hash functions that may be found in the literature. However, for the type of problems just
described, the division method, i.e.,

h(k) = k mod m

*In general, h may be a function of all keys, i.e., h = h(K). For sake of simplicity, we neglect the general case.

©1999 CRC Press LLC

 FIGURE 14.5 A simple hash table to keep track of triangular faces.

is simple and effective. Going back to our example, if we choose m = 103, then faces {104, 505, 670} and
{342, 207, 849} have the same hash value h = 1, even if their principal key is different (104 and 207,
respectively). In order to distinguish them, we need to store also the principal key in the chained linked
list records, changing the memory requirement to approximately 2mmax P + nmax [2P + lI]. Comparing
with the previous expression, it is clear that this method is convenient when nmax < < mmax. In which
particular situations would a hash table like the one presented in the example be useful? Let us assume
that somebody has given you a tetrahedral grid, without any indication of the boundary faces. How do
you find the boundary faces? You may exploit the fact that each mesh face, apart from the ones at the
boundary, belongs to two tetrahedra, set up a hash table H of the type just described, and run the
following algorithm.
1. Loop over the elements e of the mesh
1.1. Loop over element faces
1.1.1. Compute the keys K for the face and the principal key k
1.1.2. Search K in H
• If K is present then delete the corresponding record
• Otherwise add to H the record containing the face keys
2. Traverse all items in the hash table and push them onto stack F
The stack F will now obtain all boundary faces. A similar structure may be used also to dynamically store
the list of nodes surrounding each mesh node, or the list of all mesh sides and many other grid data. We
have found this hash table structure very useful and rather easy to program.
The implementation just described is useful in a dynamic setting, when add and delete operations are
required. In a static problem, when the grid is not changing, we may devise more compact representations
based on sequential storage and direct addressing. Again, let’s consider a practical problem, such as storing
a table with the nodes surrounding each given mesh node, when the mesh, formed by n nodes and ne
elements, is not changing. We use a structure K with the following attributes:
K.n = n, number of entries in the table;
K.IA[n + 1], the array containing the pointer to array JA;
K.JA[3ne], the array containing the list of nodes.

©1999 CRC Press LLC

 FIGURE 14.6 The structure used for searching the point surrounding each mesh node in the static case.

Figure 14.6 graphically shows how the structure works. The indices {IA[i], K, IA[i + 1] – 1} are used to
directly address the entries in array JA that contain the numbering of the nodes surrounding node i
(here, we have assumed that the smallest node number is 0). The use of the structure is straightforward.
The problem remains of how to build it in the first place. A possible technique consists of a two-sweep
algorithm. We assume that we have the list of mesh sides.* In the first sweep we loop over the sides and
we count the number of nodes surrounding each node, preparing the layout for the second pass:
1. For ( i = 0, i ≤ n ; i + + ) IA[i] = 0
2. Begin sweep 1: loop over the mesh sides i1, i2
2.1. For i Π{ i 1, i 2 } ) IA[i] ++
3. For (i = 1, i ≤ n, i ++) IA[i]+ = IA[i – 1]
4. For (i = n – 1, i ≥ 1, i – –) IA[i] = IA[i – 1]
5. IA[0] = 0
6. Begin sweep 2: loop over the mesh side i1, i2
6.1. For ( i Π{ i 1, i 2 } )
6.1.1. JA[IA[i]] = i1 + i2 – i
6.1.2. IA[i] ++
7. For (i = n, i ≥ 1, i – –) IA[i] = IA[i – 1]
8. IA[0] = 0
It is worth mentioning that this structure is also the basis of the compressed sparse row format, used to
efficiently store sparse matrices.

*The algorithm that works using the element connectivity list, i.e., the list of the nodes on each element, is only
a bit more complicated, and it is left to the reader.

©1999 CRC Press LLC

 FIGURE 14.7 An example of the representation of a generic tree.

14.3 Tree Structures

The data structures seen so far are not optimal when the following operations are to be made on the data:
1. Simple search: search if a record is present in the structure;
2. Range search: find all data which are within a certain range of values;
3. Tracking the maximum (or minimum) value: return the record which has the maximum value of
a given field.
Those queries are normally better answered with the adoption of a tree-type structure. Before explaining
why this is so, let us give some nomenclature. Formally, a tree T is a finite set whose elements are called
nodes (here we cannot avoid the ambiguity with mesh nodes), such that
1. There is a special node called root, which will be normally indicated by T.root;
2. The other nodes may be partitioned into n disjoint sets, {T1,K, Tn}, each of which is itself a tree,
called subtrees of the root.
If the order of the subtrees is important, then the tree is called an ordered tree. The one just given is
a recursive definition, and the simplest tree has only one node. In that case the node is called a leaf node.
A tree takes its name by the way it is usually represented. Actually, the most used representation, such
as the one shown in Figure 14.7, is an upside-down tree, with the root at the top. The tree root is linked
to the root of its subtrees, and the branching continues until we reach a leaf node. Sometimes null links
are also indicated, as in Figure 14.7. In all graphical representations of a tree contained in the remaining
part of this work, we have omitted the null links.
The number of subtrees of a node is called the degree of the node. A leaf has degree 0. A tree whose
nodes have at most degree n is termed an N-ary tree. The root of a tree T is said to be at level l = 0 with
respect to T and to be the parent of the roots of its subtree, which are called its children and are at level
l = 1, and so on. If a node is at level k it takes at least k steps to reach it starting from the tree root, and
the highest level reached by the nodes of given tree is called the height of that tree. A tree is called complete
if it has the minimum height possible for the given number of nodes.
An important case is that of binary trees. In a binary tree each node is linked to at most two children,
normally called the left and right child, respectively. Binary trees are important because many information
structures are based on them and also because all trees may in fact be represented by means of a binary
tree [10]. Tree structure is applied to grid generation in Chapter 15.

©1999 CRC Press LLC

 TABLE 14.4 Inorder and Preorder Traversal of a Binary Tree
inorder_walk (T ) preorder_walk (T )

1. T.root = NIL a RETURN; 1. T.root = NIL a RETURN;

2. inorder_walk (T.left); 2. VISIT (T.root);
3. VISIT (T.root); 3. preorder_walk (T.left);
4. inorder_walk (T.right); 4. preorder_walk (T,right);

Note: VISIT indicated whatever operation we wish to do on the tree node,

for example printing its content.

14.3.1 Binary Trees

In a binary tree each node may be considered as a record N whose two fields N.left and N.right contain
a pointer to the left and right subtree root, respectively. Therefore, a binary tree may be thought of as a
different layout of the same record used in a doubly linked list. Indeed, all considerations about using a
LAS for records management apply straight away to binary trees. While there is an unequivocal meaning
to the term “traversing a linear list,” i.e., examine the records in their list ordering, this is no longer true
for a tree structure. So, if we want to list all the nodes in a binary tree, we must first decide which path
to follow. There are basically three ways of traversing a binary tree, depending whether the root is visited
before traversing the subtrees (preorder traversal), between the traversing of the left and the right subtree
(inorder, or symmetric, traversal), or after both subtrees have been traversed (postorder traversal).
Table 14.4 shows two possible recursive algorithms for inorder and preorder traversal, respectively, the
extension to postorder being obvious. We have used recursion, since it allows to write concise algorithms.
We should warn, however, that recursion causes some computational overhead. Therefore, nonrecursive
algorithms should be preferred when speed is an issue.
14.3.1.1 How to Implement Binary Trees
We have already seen that a binary tree implementation is similar to that of a doubly linked list. Each
node is represented by a record which has two special fields containing pointers that may be either null
or pointing to the node subtrees. If both pointers are null, the node is a leaf node.
However, if the tree is a complete binary tree, a more compact representation could be adopted, which
uses an array and no pointers, but rather integer operations on the array addresses. We will use a data
organization of this type in the next section dedicated to a special complete binary tree: the heap. A more
general discussion may be found in [10].

14.3.2 Heaps
Often, there is the necessity to keep track of the record in a certain set that contains the maximum (or
minimum) value of a key. For example, in a 2D mesh generation procedure based on the advancing front
method [14] we need to keep track of the front side with the minimum length, while the front is changing.
An information structure that answers this type of query is a priority queue, and a particular data
organization which could be used for this purpose is the heap. A heap is normally used for sorting
purposes and indeed the heap-sort algorithm exploits the properties of a heap to sort a set of N keys in
O(Nlog2N) time, with no need of additional storage. We will illustrate how a heap may also be useful as
a priority queue.
We will indicate in the following with > the ordering relation. A heap is formally defined as a binary
tree with the following characteristics: If k is the key associated with a heap node, and kl and kr are the
keys associated to the not-empty left and right subtree root, respectively, the following relation holds:

k >= k1 and k >= kr

©1999 CRC Press LLC
 FIGURE 14.8 Tree representation of a heap and an example of its sequential allocation.

As a consequence, the key associated to each node is not “smaller” than the key of any node of its subtrees,
and the heap root is associated to the “largest” key in the set. We have placed in quotes the words “largest”
and “smaller” because the ordering relation > may in fact be arbitrary (as long as it satisfies the definition
of an ordering relation), and it does not necessarily correspond to the usual meaning of “greater than.”
An interesting feature of the heap is that, by employing the correct insertion and addition algorithms,
the heap can be kept complete, and the addition, deletion, and simple query operations are, even in the
worst case, of O(log2 n), while accessing the “largest” node is clearly O(1).
A heap T may be stored using sequential allocation. We will indicate by T.n and T.max the current
and maximum number of records stored in T, respectively, while T.H[max] is the array that will hold
the records. The left and right subtrees of the heap node stored at H[i] are rooted at H[2i + 1] and H[2i
+ 2], respectively, as illustrated in Figure 14.8. Therefore, by using the integer division operation, the
parent of the node stored in H[j] is H[(j – 1)/2]. The heap definition may then be rewritten as

H [ j ].key < = H [ ( j – 1 ) ⁄ 2 ].key ( 0 < j < n ) (14.1)

When inserting a new node, we provisionally place it in the next available position in the array and
we then climb up the heap until the appropriate location is found. Deletion could be done employing a
top-down procedure, as shown in Figure 14.9. We consider the heap rooted at the node to be deleted,
and we recursively move the “greatest” subtree root to the parent location, until we reach a leaf where
we move the node stored on the last array location. Finally, a bottom-up procedure analogous to that
used for node insertion is performed.*
Since the number of operations is clearly proportional to the height of the tree, we can deduce that,
even in the worst case, insertion and deletion are O(log2 n). Simple and range searches could be easily
implemented with a heap as well. However, a heap is not optimal for operations of this type.

*The terms top-down and bottom-up refer to the way a tree is normally drawn. So, by climbing up a tree we
reach the root!

©1999 CRC Press LLC

FIGURE 14.9 Deletion of the root of a heap. (a) The “highest” subtree root is recursively “promoted” until we reach
a leaf. (b) The last node is placed in the empty leaf and (c) it is sifted-up to the right place by a succession of exchanges
with its parent, until the final position (d) is reached.

14.3.3 Binary Search Tree

Better techniques for simple and range searching make use of a binary search tree. Indicating again with
k, kl, and kr the keys associated with a node, its left and right subtree roots, respectively, a binary search
tree is an oriented binary tree where the following expression is satisfied:

k1 < = k and kr > k. (14.2)

As before, > indicates an ordering relation. It should be noted that we must disambiguate the case of
equal keys, so that the comparison may be used to discriminate the records that would follow the left
branch from the ones that would go to the right. Inorder traversal of a binary search tree returns the
records in “ascending” order.
The simple search operation is obvious. We recursively compare the given key with the one stored in
the root, and we choose the right or left branch according to the comparison, until we reach either the
desired record or a leaf node. In the latter case, the search ends unsuccessfully. In the worst case, the
number of operations for a simple search is proportional to the height of the tree. For a complete tree
the search is then O(log2 n). However, the shape of a binary tree depends on the order in which the
records are inserted and, in the worst case, (which, for example, happens when a set of ordered records
is inserted) the tree degenerates and the search becomes O(n). Fortunately, if the keys are inserted in
random order, it may be proved the search is, on average, still O(log2 n) [10].
Node addition is trivial and follows the same lines of the simple search algorithm. We continue the
procedure until we reach a leaf node, of which the newly inserted node will become a left or right child,
according to the value of the key comparison. Node deletion is only slightly more complicated, unless
the deleted node has fewer than two children. In that case the deletion is indeed straightforward if the
node is a leaf, while if it has a single child, we can slice it out by connecting its parent with its child.

©1999 CRC Press LLC

FIGURE 14.10 An example of binary search tree and a graphical illustration of the operations necessary to delete
a “single parent” node (0.25) and a node with two children (0.70). The tree has been created inserting the keys in
the following order: 0.37, 0.25, 0.50, 0.10, 0.70, 0.15, 0.20, 0.55, 0.75, 0.52, 0.74.

In the case that the deleted node has two children, we have to find its successor S in the inorder
traversal of the tree, which has necessarily at most one child. We then slice S out of the tree and put it
at the deleted node location. The resulting structure is still a binary search tree, Figure 14.10 illustrates
the procedure. It can be proved that both insert and delete operations are O(log2 n) for complete binary
search trees. Other algorithmic details may be found, for example, in [2].
A binary search tree may be kept almost complete by various techniques, for example, by adopting
the red-black tree data structure [7] or the AVL tree, whose description may be found in [22].

14.3.4 Digital Trees

In a binary search tree the discrimination between left and right branching is made by a comparison
between keys. What happens if we make the comparison with fixed values instead? Let us suppose that
the keys are floating point numbers within the range [A, B). We will say that a tree node N is associated
to the interval [a, b) if all keys stored on the tree rooted at N fall within that range. Then, we can put
[a, b) = [A, B) for the root of tree T, and we may recursively build the intervals associated to the subtrees
as follows: given a tree associated to the interval [a, b), the left and right root subtrees will be associated
to the intervals [a, r) and [r, b), respectively, where r Π[ a, b ) is a discriminating value, which is usually
taken as r = (a + b)/2. This data structure is called digital search tree. Adding a node to a digital binary
tree is simple, and resembles the algorithms used for binary search trees. We start from the root and
follow the left or right path according to the result of the test k > r? The difference with binary search
trees is that the discriminant r now has a value that does not depend on the record currently stored at
the node, but on the node position in the tree structure. We want again to point out that the result of
the discriminating test must be unique: that is why the intervals are open at one end. In this way, the
case k = r will not be ambiguous, by leading to follow the path on the right.
Deleting a node is even more trivial, since every node of a digital search tree can be placed at the root
position! Therefore, when deleting a node N, we just have to substitute it with a convenient leaf node of
the subtree rooted at N (of course, we do not move the nodes, we just reset the links). For example, if
N is not a leaf node (in which case we could just release it), we would traverse in postorder the subtree
rooted at N, as the first “visited” node will certainly be a leaf, and substitute N with it. In order to keep
a better balanced tree in a highly dynamic setting, when many insertion and deletion operations are
expected, we could keep at each node N a link to a leaf node of the tree rooted at N with the highest
level l. That node will be used to substitute N when it has to be deleted. With this technique the algorithms
for insertion/deletion are just a little more involved, and it could be useful to store at each node N also
the link N.parent to the parent node. With a digital tree we cannot slice out of the tree a single-child
node as we have seen in binary trees, since this operation will cause a change of node level in the subtree
rooted at the sliced-out location, and the discriminant value r, which is function of the node level, will

©1999 CRC Press LLC

FIGURE 14.11 An example of a binary search tree (a), a digital search tree (b), and binary trie structure (c), for a
given set of data.

change. As a consequence the structure resulting from this operation will, in general, not be a digital
search tree anymore.
The name of this data structure derives from the fact that it is in principle possible to transform a key
into an ordered set of binary digits d0, d1, K, dm so that, at tree level l, the decision to follow the right
or left path could be made by examining the lth digit. In particular, for the example just shown, the
decision could be made by considering the lth significant digit of the binary representation of the number
(k – a)/(b – a), and following the left path if dl = 0.
Directly related to the digital search tree is the trie structure [11], which differs from the digital search
tree mainly because the actual data is only stored at the leaf nodes. The trie shape is completely inde-
pendent from the order in which the nodes are inserted. The shape of a digital tree, instead, still depends
on the insertion order. On the other hand, a trie uses up more memory, and the algorithms for adding
and deleting nodes are a little more complex. Contrary to a binary search tree, both structures require
the knowledge of the maximum and minimum value allowed for the key.
Figure 14.11 shows an example of the search structures seen so far for a given set of keys k Π[ 0, 1 ). For
the digital and binary trie we have put beside each node the indication of the associated interval.

14.4 Multidimensional Search

In mesh generation procedures there is often the necessity of finding the set of nodes or elements lying
within a specified range, for instance, finding the points that are inside a given sphere. There also arises
the need to solve geometric intersection problems, such as finding the triangular faces that may intersect
a tetrahedron.
Needless to say, these are fundamental problems in computational geometry (cf. Part III of this
handbook) and there has been a great deal of research work in the last years aimed at devising optimal
data structures for this purpose. There is not, however, a definite answer. Therefore, we will again
concentrate only on those structures suited for mesh generation procedures, where we usually have
dynamic data and where memory occupancy is a critical issue.

©1999 CRC Press LLC

14.4.1 Searching Point Data
Given a set of points P in Rd, where d is either 2 or 3, we will consider the following queries:
• Point search: is the point P present in the P?
• Range search: that are the points of P that lie inside a given interval I ⊂ R d ?

In addition to those operations, we may want to be able to efficiently add and delete nodes to the set.
For the case d = 1, it was shown in the previous section that a binary search tree could efficiently answer
these queries. It would be natural to ask whether it can be used also for multidimensional searches.
Unfortunately, binary search is based on the existence of an ordering relation between the stored data,
and there is normally no way of determining an ordering relation between multidimensional points. In
fact, we will see that in principle an ordering relation may be found, for instance using a technique called
bit interleaving, but in practice this procedure is not feasible, as it would require costly operations, both
in terms of computation and memory.
The most popular procedures for handling multidimensional searches are either based on hierarchical
data structure or on grid methods [20]. We will illustrate some of the former, and in particular data
structures based either on binary trees quadtrees, or octrees. For sake of simplicity, we will consider only
a Cartesian coordinate system and the two-dimensional case, the extension to 3D being obvious.

14.4.2 Quadtrees
The quadtree is “4-ary” tree whose construction is based on the recursive decomposition of the Cartesian
plane. Its three-dimensional counterpart is the octree. There are basically two types of quadtrees, depend-
ing whether the space decomposition is driven by the stored point data (point-based quadtrees) or it is
determined a priori (region-based quadtrees). Broadly speaking, this subdivision is analogous to the one
existing between a binary and a digital search tree. We will in the following indicate with B the domain
bounding box defined as the smallest interval in R2 enclosing the portion of space where all points in P
will lie. Normally, it can be determined because we usually know beforehand the extension of the domain
that has to meshed. For sake of simplicity, we will often assume in the following that B is unitary, that
is B ≡ [ 0, 1 ) × [ 0, 1 ) . There is no loss of generality when using this assumption, as an affine transfor-
mation can always be found that maps our point data set into a domain enclosed by the unitary interval.
14.4.2.1 Region-Based Quadtrees
A region-based quadtree is based on the recursive partitioning of B into four equally sized parts, along
lines parallel to the coordinate axis. We can associate to each quadtree node N an interval N.I = [a, b)
× [a, b) where all the points stored in the tree rooted at N will lie. Each node N has four links, often
denoted by SW, SE, NW, NE, that point to the root of its subtrees, which have associated the intervals
obtained by the partitioning. Point data are usually stored only at leaf nodes, though it is also possible
to create variants where point data can be stored on any node. Figure (14.12) illustrates an example of
a region quadtree. The particular implementation shown is usually called PR-quadtree [20,19].
Point searching is done by starting from the root and recursively following the path to the subtree
root whose associated interval encloses the point, until we reach a leaf. Then the comparison is made
between the given node and the one stored in the leaf. Range searching could be performed by examining
only the points stored in the subtrees whose associated interval has a non-empty intersection with the
given range. Details for point addition/deletion procedures may be found in the cited reference. The
shape of the quadtree here presented, and consequently both search algorithm efficiency and memory
requirement, is independent of the point data insertion order, but it depends on the current set of points
stored. If the points are clustered, as often happens in mesh generation, this quadtree can use a great
deal of memory because of many empty nodes. Compression techniques have been developed to overcome
this problem: details may be found in [15]. In unstructured mesh generation the region quadtree, and
the octree in 3D, is often used not just for search purposes [12], but also as a region decomposition tool
(see also Chapter 22). To illustrate the idea behind this, let us consider the example shown in Figure 14.13,

©1999 CRC Press LLC

 FIGURE 14.12 An example of a region-based quadtree.

FIGURE 14.13 Domain partitioning by a region quadtree. The quadtree contains the boundary points. The parti-
tions associated with the quadtree nodes are shown with dotted lines.

where the line at the border with the shaded area represents a portion of the domain boundary. A region
quadtree of the boundary nodes has been built, and we are showing the hierarchy of partitions associated
to the tree nodes. It is evident that the size of the partitions is related to the distance between boundary
points, and that the partitioning is finer near the boundary. Therefore, structures of this type may be
used as the basis for algorithms for the generation of a mesh inside the domain, in various ways. For
instance, a grid may be generated by appropriately splitting the quad/octree partitions into triangle/tet-
rahedra [23]. Alternatively, the structure may be used to create points to be triangulated by a Delaunay
type procedure [21] (cf. Chapter 16). Finally, it can be adopted for defining the mesh spacing distribution
function in an advancing front type mesh generation algorithm [9] (cf. Chapter 17).
14.4.2.2 Point-Based Quadtrees
A point quadtree is a type of multidimensional extension of the binary search tree. Here the branching
is determined by the stored point, as shown in Figure 14.14. It has a more compact representation than
the region quadtree, since point data is stored also at non-leaf nodes. However, the point quadtree shape
strongly depends on the order in which the data is inserted, and node deletion is rather complex.
Therefore, it is not well suited for a dynamic setting. However, for a static situation, where the points

©1999 CRC Press LLC

 FIGURE 14.14 An example of a point-based quadtree. Nodes have been inserted in lexicographic order.

are known a priori, a simple technique has been devised [4] to generate optimized point quadtree, and
this fact makes this structure very interesting for static situations, since simple search operations become
O(log4 n) n being the total number of points stored. It should be mentioned that a procedure that allows
for dynamic quadtree optimization has also been devised [16]. Its description is beyond the scope of this
chapter.

14.4.3 Binary Trees for Multidimensional Search

A shortcoming of quad/octrees is that they are rather costly in terms of the memory required. It is
possible, however, to use a binary tree also for a multidimensional search. Indeed many of the ideas
illustrated for the one-dimensional case may be extended to more dimensions if we allow the discrimi-
nating function r at each tree level to alternate between the coordinates. We may use the same technique
adopted for the binary search tree, but we now discriminate according to the x coordinate at even-level
nodes, while nodes at odd level are used to discriminate according to the y coordinate. We have now a
two-dimensional search tree denoted k-d tree, which may be also considered as the binary tree counterpart
of a point-based quadtree. Figure (14.15) shows an example of a k-d tree. As the node where point D is
stored by a y-discriminator, its left subtree contains only points P that satisfy P.y < D.y. K-d trees are a
valid alternative to point-based quad/octrees. According to Samet [20]: “we can characterize the k-d tree
as a superior serial data structure and the point quadtree as a superior parallel data structure.” However,
they also share the same defects. Their shape strongly depends on the node insertion order, unless special
techniques are adopted [3], and node deletion is a quite complex operation also for k-d trees.
An alternative that has encountered great success in unstructured mesh generation is the alternating
digital tree (ADT) [1], which is a digital tree where, as in a k-d tree, the discrimination is alternated
between the coordinates. It differs from a k-d tree because here the discrimination is made against fixed
space locations. To each node N we associate an interval in [x1, x2) × [y1, y2), and all point data in the
subtree rooted at N will lie in that interval. The tree root is associated to the bounding interval B. If a
node is an x-discriminator, its left and right child will be associated to the intervals [x1, r) × [y1, y2) and
[r, x2) × [y1, y2), respectively, with r = (x1 + x2)/2. A y-discriminating node will act in a similar fashion
by subdividing the interval along the y axis. Figure (14.15) illustrates an example of an ADT tree.
The algorithms for node addition and deletion are analogous to the ones shown for the one-dimen-
sional digital tree. Simple searching is O(log2n) if the tree is complete. Unfortunately, the tree shape is
not independent of the order of node insertion, even if, in general, ADT trees are better balanced than
their k-d counterpart. For static data, while a special insertion order has been devised to get a balanced
k-d tree, no similar techniques are currently available for an ADT. Therefore, we may claim that ADTs
are better than k-d trees for dynamic data, while for a static situation, a k-d tree with optimal point

©1999 CRC Press LLC

FIGURE 14.15 An example of a two-dimensional k-d tree (top) and an ADT (bottom), on the same set of data.
Nodes at even levels discriminate the x coordinate, while odd nodes discriminate the y coordinate. In the k-d tree
the discrimination is made against the data stored at the node, while in the ADT structure we use fixed spatial
locations. Nodes have been inserted in lexicographic order in both cases.

insertion order is more efficient. Range searches in an ADT are made by traversing the subtrees associated
with intervals which intersect the given range.
A region decomposition based structure similar to ADT, where the data points are stored only at leaf
nodes is the bintree [20]: it has not been considered here because of its higher memory requirement
compared with ADT.
14.4.3.1 Bit Interleaving
For sake of completeness, we mention how, at least theoretically, a binary search tree may be used also
for multi-dimensional searching using a technique called bit interleaving. Let us assume that B is unitary.
Then, given a point P = (x, y) we may consider the binary representation of its coordinates, which we
will indicate as x0, x1, x2, K, xd and y0, y1, y2, K, yd. We may now build a code by interleaving the binary
digits, obtaining x0, y0, K, xd, yd and define an ordering relation by treating the code as the binary
representation of a number. The code is unique for each point in B, and we can use it as a key for the
construction of a binary search tree. This technique, however, is not practical because it would require
storing at each node a code that has a number of significant digits twice as large as the one required for
the normal representation of a float (three times as large for 3D cases!). It may be noted, however, that
the ADT may indeed be interpreted as a digital tree where the discrimination between left and right
branching at level l is made on the base of the lth digit of the code built by a bit interleaving procedure
(without actually constructing the code!).

14.4.4 Intersection Problems

Geometry intersection problems frequently arise in mesh generation procedures (see also Chapter 29).
In a front advancing algorithm, for instance, we have to test whether or not a new triangle intersects the
current front faces. General geometrical intersection problems may be simplified by adopting a two-step

©1999 CRC Press LLC

FIGURE 14.16 Intersection problem solved by means of an ADT structure. The subtree rooted at node B (shaded
nodes) does not need to be examined, since all subtree nodes correspond to rectangles that lie in the half-hyperspace
x1 ≥ 1/2, which cannot intersect R, since R.x2 ≤ 1/2.

procedure. In the first step we associate with each geometrical entity of interest G its smallest enclosing
interval I G ≡ [ x 1G, x 2G ] × [ y 1G, y 2G ] , and we then build specialized data structure which enables one to
efficiently solve the following problem.
Given a set I of intervals (rectangles in 2D or hexahedra in 3D), find the subset H ⊂ I of all
elements of I which intersect a given arbitrary interval.
In the second phase, the actual geometrical intersection test will be made, restricted only to those
geometrical entities associated to the elements of H.
Data structures that enable solving efficiently this type of problem may be subdivided into two
categories: the ones that represent an interval in Rn as a point in Rn2, and those that directly store the
intervals. An example of the latter technique is the R-tree [20], which has been recently exploited in a
visualization procedure for three-dimensional unstructured grid for storing sub-volumes so that they
can be quickly retrieved from disk [13]. We will here concentrate on the first technique: i.e., how to
represent an interval as a point living in a greater dimensional space.
14.4.4.1 Representing an Interval as a Point
Let us consider the 2D case, where intervals are rectangles with edges parallel to the coordinate axis. A
rectangle R ≡ [ x 1, x 2 ] × [ y 1, y 2 ] may be represented by a point P Œ R 4 . There are different representa-
tions possible, two of which are listed in the following:
1. P ≡ [ x 1, y 1, x 2, y 2 ]
2. P ≡ [ xc, yc, dx, dy ] where xc = (x1 + x2)/2, dx = (x2 – xc), K
In the following we will consider the first representation. Once the rectangle has been converted into a
point, we can adopt either a k-d or an ADT data structure both for searching and geometrical intersection
problems. If we use an ADT tree, the problem of finding the possible intersections can be solved by
traversing the tree in preorder, excluding those subtrees whose associated interval in R4 cannot intersect
the given rectangle. Figure 14.16 shows a simple example of this technique.

14.5 Final Remarks

We have given an overview of some of the information structures that may be successfully adopted within
mesh generation schemes. This survey is, of course, not exhaustive. Many ingenious data structures have
been devised by people working on grid generation and related fields in order to solve particular problems
in the most effective way. We will in this final section, just mention a few of the efforts in this direction
that we have not had the possibility to describe in detail.

©1999 CRC Press LLC

 The List [10] structure (with uppercase L!) has been adopted [18] to control a hierarchy of grids for
multigrid computations. An edge-based structure [8] has been devised for storing mesh topology data,
which should be more efficient for Delaunay mesh generation algorithms. Doubly linked circular lists
are used for the implementation of a grid topology model that allows an efficient automatic block
detection for multiblock structured mesh generation procedures [5, 6].
Many other examples could be made. We hope that the reader now has and idea of how appropriate
data structures may help in devising efficient grid generation procedures.

References
1. Bonet, J. and Peraire, J., An alternating digital tree (ADT) algorithm for 3D geometric searching
and intersection problems, Int. J. Num. Meths. Eng., 31, pp. 1–17, 1991.
2. Cormen, T. H., Leiserson, C. E., and Rivest, R. L., Introduction to Algorithms, The MIT Electrical
Engineering and Computer Science Series, McGraw-Hill, 1990.
3. Fiedman, J. H., Bentley, J. L., and Finkel, R. A., An algorithm for finding best matches in loga-
rithmic expected time, ACM Transactions on Mathematical Software, 3(3), pp. 209–226, September
1977.
4. Finkel, R. A. and Bentley, J. L., Quad trees: a data structure for retrieval on composite keys, Acta
Inform. 1974, 4: pp 1–9.
5. Gaither, A., A topology model for numerical grid generation, Weatherill, N., Eiseman, P. R., Hauser,
J., Thompson, J. F., (Eds.), Proceedings of the 4th International Conference on Numerical Grid
Generation and Related Fields, Swansea, Pineridge Press, 1994.
6. Gaither, A., An efficient block detection algorithm for structured grid generation, Soni, B. K.,
Thompson, J. F., Hauser, J., Eiseman, P. R., (Eds.), Numerical Grid Generation in Computational
Field Simulations, Vol. 1, 1996.
7. Guibas, L. J. and Sedgewick, R., A diochromatatic framework for balanced trees, Proceedings of
the 19th Annual Symposium on Foundations of Computer Science, IEEE Computer Society, 1978,
pp. 8–21.
8. Guibas, L. J. and Stolfi, J., Primitives for the manipulation of general subdivisions and the com-
putation of Voronoï diagrams, ACM Transaction on Graphics, April 1985, 4(2).
9. Kallinderis, Y., Prismatic/tetrahedral grid generation for complex geometries, Computational Fluid
Dynamics, Lecture Series 1996-06. von Karman Institute for Fluid Dynamics, Belgium, March 1996.
10. Knuth, D. E., The Art of Computer Programming. Vol. 1, Fundamental Algorithms of Addison-Wesley
Series in Computer Science and Information Processing. Addison–Wesley, 2nd ed., 1973.
11. Knuth, D. E., The Art of Computer Programming, Vol. 3, Sorting and Searching of Addison–Wesley
Series in Computer Science and Information Processing. Addison-Wesley, 1973.
12. Lohner, R., Generation of three dimensional unstructured grids by advancing front method, AIAA
Paper 88-0515, 1988.
13. Ma, K. L., Leutenegger, S., and Mavriplis, D., Interactive exploration of large 3D unstructured-
grid data, Technical Report 96-63, ICASE, 1996.
14. Morgan, K., Peraire, J., and Peirò, J., Unstructured grid methods for compressible flows, Special
Course on Unstructured Grid Methods for Advection Dominated Flows, 1992, AGARD-R-787.
15. Ohsawa, Y. and Sakauchi, M., The BD-tree, a new n-dimensional data structure with highly
efficient dynamic characteristics, Mason, R.E.A., (Ed.), Information Processing 83, North-Holland,
Amsterdam, 1983, pp. 539–544.
16. Overmars, M. H. and van Leeuwev, J., Dynamic multi-dimensional data structures based on quad-
and k-d trees, Acta Informatica, 17(3), pp. 267–285, 1982.
17. Preparata F. P. and Shamos, M. I., Computational Geometry: An Introduction, Springer–Verlag, 1985.
18. Rivara, M. C., Design and data structures of fully adaptive, multigrid, finite-element software,
ACM Trans. Math. Soft., 10, 1984.

©1999 CRC Press LLC

 19. Samet, H., The quadtree and related hierarchical data structures, Computing Surveys, 1984, 16, pp.
188–260.
20. Samet, H., The Design and Analysis of Spatial Data Structures, Addison–Wesley, 1990.
21. Schroeder, W. J. and Shephard, M. S., A combined octree/Delaunay method for fully automatic
3D mesh generation, Int. Journal on Numerical Methods in Eng., 29, pp. 37–55, 1990.
22. Wirth, N., Algorithm + Data Structures = Programs, Prentice-Hall, Englewood Cliffs, NJ, 1976.
23. Yerry, M. A. and Shephard, M. S., A modified quadtree approach to finite element mesh generation,
IEEE Computer Graphics and Applications, 3, pp. 39–46, January/February 1983.

©1999 CRC Press LLC