1. Introduction

There are some very powerful techniques for obtaining solutions to partial differential equations that govern fluid flow problems, i.e., finite difference, finite volume, and finite element methods. In order to utilize one of those techniques, it is necessary to replace the spatial domain of the interest into a structured or unstructured discrete domain, depending on the solution method employed to solve the problem. This spatial domain conversion process is very important for many reasons. Depending on how the domain is discretized and grid points are controlled, the solutions obtained using an appropriate solution scheme may be not so accurate, due to many errors, such as discretization, and interpolation errors. Even though considerable progress has been made in grid generation, the generation of good quality structured or unstructured grid systems in geometrically complex, three-dimensional spatial domains remains a difficult problem.

The grid generation techniques for computational fluid dynamics has reached such a level of maturity that it can be considered as a solved problem. Although this is the case, each technique has its particular strengths and weaknesses. A well-defined standard of an optimal grid, due to numerical properties of the solution scheme is still not clear. Thus, further investigation is needed in this field.

In this paper, current structured and unstructured grid generation techniques are surveyed and presented. The strengths and weaknesses, especially for unstructured grid generation, are also presented and discussed. Since the purpose of this paper is to present various techniques to generate grids, which are often involved with complicated mathematics, the detailed and theoretical derivations for those techniques are avoided as much as possible. The main goal of this paper is to show currently well developed and advanced structured and unstructured grid generation techniques.
 

2. Structured Grid Generation

Even though there are some advantages of employing structured grids due to their simplicity, availability of code, and suitability for multi-grid and finite difference methods, it is very difficult to generate a structured grid for complex configurations, such as a complete aircraft. In case of structured grid based computations, grid generation is considered the most significant part of finite-difference and finite-volume methods because the grid points used for the computations strongly affect the accuracy, efficiency and ease with which these methods generate solutions. Depending on how a structured grid is generated, finite-difference or finite-volume method may not even be used to compute the solutions to a problem.
 

Types of Grid systems

All possible grid systems that can be used by employing finite-difference or finite-volume methods can be classified as follows: structured, unstructured, or mixed. This classification is based on how the grid points are connected to each other. A structured grid system in turn can be classified as either a single or composite grid. A composite grid can be classified as completely discontinuous, partially discontinuous, partially continuous or completely continuous. In this section of the paper, the main required conditions and techniques used to generate structured grids will be explored.

If a finite-difference or finite-volume method is used with a structured grid to obtain solutions to fluid flow problems, there are certain conditions to be satisfied when generating a structured grid. From [1], they are:

1. The total number of grid points in the grid should be kept to the minimum needed to obtain solutions of the desired accuracy. This condition should be met for both structured and unstructured grids and very significant for computational efficiency. This condition can be achieved by clustering grid points in regions where they are most needed, i.e., the places where physics are rapidly changing, e.g., a shock wave in a supersonic flow, and scattering them elsewhere. If the grid points are not clustered in the regions where needed, solutions obtained may not have meaningful physics due to a low accuracy. The classical example is a boundary layer computation. If not enough grid points are used where a boundary layer is expected to occur, the boundary layer sometimes cannot be even seen. Currently, based on my experience in the field of computational fluid dynamics, the grid points are often clustered near solid wall boundaries. If too many grid points used overall, even currently well advanced computers may not able to handle the computation at all due to a huge memory requirement. Thus, this condition is very important.
2. The grid should be boundary conforming. That is, one set of grid lines should always coincide with the physical boundary of the spatial domain regardless of the geometric complexity. This rule is often met, although it is very difficult to generate a boundary conforming grid for a highly curved surface.
3. Grid lines that intersect a boundary should intersect that boundary perpendicularly so that derivative boundary conditions can be implemented more easily and accurately. In case where an adiabatic wall condition is used, the Neumann boundary condition is desired, which is the temperature gradient normal to the wall surface is equal to zero. This boundary condition can be easily implemented with a low error if the grid lines intersect a boundary orthogonally. At the interior of the spatial domain, the angle of intersection between grid lines only needs to be nearly orthogonal, but must be somewhere between 45 and 135 degrees.
4. The spacings between grid points should change slowly from a region where grid points are concentrated to a region where grid points are sparsely distributed, especially in regions where gradients of the flow are large. This condition is important because Fourier components, which make up the solution reflect and refract at interfaces where grid spacings change.
5. One set of grid lines should align with the flow direction. This condition is important for convection-dominated flows when the aspect ratio of the control volume about each grid point is very high and/or when the thin layer Navier-Stokes equations are used to study such flows.

Completely discontinuous composite (CDC) grids

There are some major advantages using completely discontinuous composite grids. The current major CDC grid widely used is chimera grid. Due to trivial patching, the grid generation process is very easy and straightforward. Also, since the structure of each single grid can be different from each other, optimization of each single grid is possible and this is the easiest grid to refine when a local grid refinement is desired. Furthermore, since composite grids are composed of a series of single grids, it is possible to do computations on one single grid at a time, thus reduces memory requirements.

Although many advantages of using CDC grids as listed above, they are the most difficult to use in obtaining solutions when compared to other types of composite grids because interpolation schemes are needed to transfer information from one grid to another when two or more single grids overlap. Also, properties such as conservation and monotonicity must be ensured in regions where two or more single grids overlap, which is very difficult to do.
 

Partially discontinuous composite (PDC) grids

The grid generation steps for partially discontinuous composite grids are very similar to completely discontinuous composite grids, except the fact that partially discontinuous composite grids do not overlap. Since this is the case, PDC and CDC grids share same advantages and disadvantages. Only difference is that PDC grids are more difficult to generate but are easier to use than CDC grids. The main difficulty in using PDC grids is implementing boundary conditions at interfaces where different single grids meet.
 

Partially and completely continuous composite grids

In completely continuous composite grid systems, all grid lines and all of their derivatives of every order are continuous at all interfaces where different single grids meet. In practice, these grid systems are not employed because the continuity of all derivatives must be ensured when generating the grids. Instead, partially continuous composite grids are often used because of their efficiency. For the partially continuous composite grids, only up to second derivatives are continuous in this grid system.

The major advantage of using partially continuous composite grids is the ease of implementing boundary conditions at the interfaces where different grids meet. In fact, boundary conditions are not even necessary because computations can be carried across them. This is possible because all of the grid lines are continuous up to the second derivatives at the interfaces.

On the other hand, these grids are the most difficult ones to generate since the structure and number of grid points in each single grid must satisfy certain compatibility conditions in order to ensure continuity at the interfaces. And with these conditions, it is impossible to optimize each grid for different location and to do local grid refinement.
 

2.1 Grid Generation Methods

There are two major ways to generate grids. They are differential equation methods and algebraic methods. Differential equation methods are difficult to use because a system of partial differential equations must be solved in order to distribute grid points within the spatial domain. Because these methods require a significant amount of computational time to solve those complex partial differential equations, algebraic methods are often employed to generate grids. The algebraic methods generate the grids by interpolating between boundaries of the spatial domain and therefore computationally much more efficient.

Depending upon the complexity and dimensionality of the spatial domain, one of the following three methods can be used to generate grids. They are the two-boundary method, the four-boundary method, and the six-boundary methods. These methods are very similar to each other and can generate grid lines that intersect boundaries orthogonally by using algebraic methods mentioned above. They are also highly efficient and give precise controls over the distribution of grid points in the spatial domain when used in conjunction with stretching functions.

The two-boundary method is intended for problems in which it is only necessary to map correctly two arbitrary-shaped boundaries of the spatial domain, e.g., a two-dimensional converging diverging nozzle. If the remaining boundaries are straight lines in the two-dimensional case or flat surfaces in the three-dimensional case, the all boundaries can be mapped correctly. Similarly, the four-boundary method is used when four boundaries must be conformed correctly and the six-boundary method is used when all of the six boundaries must be mapped correctly. The four-boundary method is an extension of the two-boundary method and the six-boundary method is an extension of both the two and four-boundary methods. These grid generation methods are highly efficient and they are able to provide very precise controls over the distribution of grid points when stretching functions are used.

There are some problems when using those three methods described above. First, slope discontinuities present in the boundaries of spatial domains propagate into the interior of the grid systems generated by using these methods. These discontinuities in the slopes of the grid lines are undesirable since they can lead to errors in the solution. To correct this problem, some technique should be employed to smooth the grid. One way is to apply a Laplacian operator to the region near the discontinuity. Other major problem is that when connecting curves are discretized to form grid points, the orthogonality, which was forced at the boundaries, may be lost. The detailed explanation can be found in the reference [1].
 

3. Unstructured Grid Generation

Although unstructured grid generation is considered to be well-advanced over the last fifteen years, further improvements in efficiency and robustness are desired. Many unstructured grid generation techniques are purely based in the field of computational geometry. The main classical example is the Delaunay triangulation, which we are very familiar from the class. Many grid construction algorithms utilize this popular and well known technique. But unfortunately, the concept of an optimal triangulation or algorithm from a computational geometry point of view does not always coincide with the view from a computational fluid dynamics point of view, and thus many of the computational geometry results have found little use in the area of unstructured grid generation. In this section, a preliminary discussion of several computational geometry algorithms, which are most relevant to unstructured grid generation techniques, is presented first. As mentioned previously, some theoretical details will not be presented.
 

3.1 Computational Geometry Constructs and Algorithms
 

The Delaunay Triangulation

The Delaunay triangulation is the dual of the Voronoi Tessalation. A Voronoi Tessalation is the graph obtained by drawing the median line segments, which separate the plane into regions, which are closer to a given point than any other points. Another property of the Delaunay triangulation is known as the empty circumcircle property. This simply means that there should not be any points within a circumcircle of a triangle formed by three vertices, which make an "element". Another property of the Delaunay triangulation is a max-min angle triangulation. That is, out of all the possible triangulations of a given set of points, it is the triangulation, which gives the largest minimum angle for all triangular elements. Because of those properties listed above, the Delaunay Triangulation is expected to give well-shaped elements. Especially, it is very useful when used to generate three-dimensional grids because the empty circumcircle property can be easily extended to three-dimensional problems as the circumsphere associated with each tetrahedron.
 

Bowyer-Watson Algorithm

This algorithm basically is an incremental algorithm to generate unstructured grids by adding new points sequentially to an existing Delaunay triangulation. Instead of breaking all the triangles when a new point is added, the algorithm first searches for the main triangle in which the new point is inserted and its neighboring triangles. Then, all the possible triangles whose circumcircle is intersected by the new point are identified and removed from the current triangulation. Properties of the Delaunay triangulation guarantee that such a neighbor search is sufficient enough to locate all intersected triangles. Once all of the above steps are done, the algorithm re-triangulates by joining the new point to all vertices on the boundary of the polygonal cavity.

This algorithm was proven to be very efficient and useful, mainly due to its cheap computational cost. For unstructured mesh generation, the computational cost is near linear O(n) for two and three-dimensional cases.
 

Green-Sibson Algorithm

This algorithm is also based on the Delaunay triangulation and is very similar to the Bowyer-Watson algorithm, in the sense that a new point is inserted to an existing triangulation. But, there is a huge difference between two algorithms. That is, when the Green-Sibson Algorithm is utilized, only the triangle in which the new point is inserted is identified by search and the new triangulation is done simply by joining the new point to the three vertices of the triangle in which the new point is inserted, without removing any other edges of the neighboring triangles. The resulting triangulation is not necessarily Delaunay. Therefore, the rest of the process is to make the newly formed triangles

Delaunay by rearranging the grid lines in the vicinity of the new point. In case where one of these circumcircles associated with the new triangles contains a vertex, only the edge of the triangle which is shared by the neighboring triangles are removed and rearranged. Each time an edge is rearranged, two triangles are altered, and these must therefore be checked for the Delaunay criterion. The only edges, which need to be considered for swapping, are those, which border on a modified triangle and an outer triangles, which are not modified. The edge swapping procedure begins with the outer edges of the newly formed triangles, and propagates outwards until no further edges need to be swapped.
 

Tanemura-Merriam Algorithm

This algorithm approaches differently than the above two algorithms presented. The previous two algorithms represent a top-down approach, while this algorithm is a bottom-up approach by using the advancing-front Delaunay algorithm. The main idea is to construct the triangulation one triangle at a time, starting at the boundaries of the spatial domain, thus "advancing" or sweeping a front throughout the spatial domain. An arbitrary interior point is chosen, and the triangle formed by the two endpoints of the front edge and the interior point is constructed. If this triangle contains any other points within its circumcircle, an alternate point is chosen because it is not a Delaunay triangle. That is, the newly chosen point will be the point contained inside the newly formed circumcircle, which is closest to its center of the circumcircle.

By continuing this process until no points found within its circumcircle, the appropriate point will be eventually found. This new triangle formed by those two points on the front edge and this appropriate point will be a new triangle, and the front is advanced by removing the current edge from and adding the new edges to the front. This algorithm will stop once all the points in the spatial domain are used to form triangles.
 

Constrained Delaunay Triangulations

A constrained Delaunay triangulation is a triangulation which contains a set of prescribed edges, which are not necessary convex. The existence of constrained Delaunay triangulations ensures the validity of the Tanemura-Merriam algorithnm for arbitrary initial front. It also guarantees the possibility of modifying an existing Delaunay triangulation to include a set of edges, which define the boundaries of the domain to be triangulated. Unfortunately, an equivalent definition for constrained Delaunay triangulations in three-dimensions is not available.
 

3.2 Unstructured Grid Generation Methods
 

Quadtree and Octree Based Grid Generations

One of the earliest and simplest methods for generating unstructured grids involves the use of quad and octrees in two and three dimensions respectively. An initial quad is formed which is large enough to cover the entire domain. Assuming a grid element-size distribution function exists, the quarter is recursively subdivided until all leaf quads are no larger than the local value of the element size distribution function. Subdivision in the domain interior can be proceed by ensuring that jumps in the sizes of neighboring leaf elements never exceeds two to one. At the physical boundaries, the quadtree must be made boundary conforming. This is done by moving the closest quadtree vertices to coincide with the boundaries. These methods are very simple and efficient to implement.
 

Advancing-Front Methods

The main idea of these methods is to construct grid element by element, adding new elements to previously generated elements, thus sweeping out a front across the entire domain. These methods divide the boundaries first, using an element size distribution function.

Once this is done, one edge is selected to form a new triangle by joining the two ends of the current edge either to a newly created point, or to an existing point on the front. The current edge is then removed from the front and new front is to one of the two edges of the new triangle, depending on their visibility.

One of the advantages of such an approach is the automatic point placement strategy, which generally results in high-quality elements throughout most of the flow field, due to the optimum positioning of these new points. Also, boundary integrity is guaranteed, since the boundary discretization are used as the initial condition.

Although there are many advantages using these methods, there are also disadvantages associated with them as well. One of the major disadvantages is related mainly to the efficiency. Since the advancing front techniques generates a grid one element at a time and there are about twice as many triangles as points in two-dimensional cases, some improvements are desired. Efficiency can be gained by determining all the potential triangles, which join this new point to the existing front each time a new point is generated. Much more efficiency can be gained in three-dimensional cases.
 

Delaunay Point-Insertion Methods

These Delaunay based methods use a set of single structured grid. Each individual grid represents a simple configuration, which can be patched locally to each other. By constructing a set of overlapping structured single grids this way and ignoring the grid lines within those grids, as well as the grid points which fall outside of the spatial domain, a set of grid points can be obtained which then can be used as the basis for a Delaunnay triangulation grid generation technique, e.g., Boyer-Watson or Green-Sibson algorithms, which are mentioned previously. Due to the overlapping mesh construction of the point set, close grid points may be produced, and grid smoothing is used as a post process in order to remove the local irregularities in the final grid spacing. The space and time requirements of these methods are typically O(nlogn). In practice, these methods have produced some of the most efficient grid generation codes available today, capable of generating one million three-dimensional elements per hour on present day workstations.
 

Other Unstructured Grid Generation Techniques based on Advancing Front and Delaunay Triangulation

Besides the unstructured grid generation techniques explained previously, there are a few more techniques which employ the same basic principles, i.e., advancing front and Delaunay Triangulation. Some well known and successful techniques include ‘Advancing Front Delaunay Triangulation’ and ‘Advancing Front Point-Insertion’ Methods. Those techniques are basically extensions and/or combinations of the above two major principal techniques. Computational cost of the advancing-front Delaunay triangulation is O(n²) at worst cases, but usually O(nlogn) is expected, which is about the same as the Delaunay Point-Insertion Methods.
 

4. Future Directions

I will end this paper with some promising future directions. There are several problems in the area of grid generation that should be addressed in the future. Currently, it usually takes two to six weeks to prepare the geometry and volume grid for a detailed Navier-Stokes calculation. This time should be reduced to hours, and at the same time, the grid generated should be flexible and robust enough to handle computations dealing with high complexities. The grid generating process also should be run in a more automatic fashion to reduce human errors and time.

In case of unstructured grids, additional work is needed in order to allow the unstructured grid generations for the viscous flow computations. Also, the current unstructured grids demand a huge memory requirement, due to the needs to keep track of the vertex and cell information and their connectivities. Therefore, ways to reduce this memory requirement should also be explored.
 
 

References

[1] Shih, T. I-P, et al, "Algebraic Grid Generation for Complex Geometries", International Journal for Numerical Methods in Fluids, Vol. 13, p. 1-31, 1991.

[2] Mavriplis, D. J., "Unstructured Mesh Generation and Adaptivity", ICASE Report No. 95-26.

[3] Pirzadeh, S., "Structured Background Grids for Generation of Unstructured Grids by Advancing-Front method", AIAA Journal, Vol. 31, No. 2, p. 257-265,      Feb. 1993.

[4] Bowyer, A., "Computing Dirchlet tessellations", The Computer Journal, Vol. 24, No. 2, 1981.