- Email: [email protected]
A divide-and-conquer algorithm for grid generation
APPLIED NUMERICAL MATHEMATICS Applied Numerical Mathematics 14 (1994) 125-134
ELSEVIER
A divide-and-conquer Randall L. Dougherty
algorithm for grid generation * and James M. Hyman **
Theoretical Division, MS B284, Center for Nonlinear Studies, Los Alamos, NM 87545, USA
Abstract We present a new algorithm for generating a logically rectangular curvilinear mesh for a region in the plane bounded by four smooth curves. The curvilinear grid is generated using a divide-and-conquer algorithm where at each step a single grid curve or line is generated using position and derivative information derived from previously generated curves. The grid curves are produced by shape-preserving cubic Hermite interpolation that allow for corners or kinks along the boundary. After all grid curves have been generated, the mesh is modified by a rezoning algorithm to improve smoothness, orthogonality, resolution of a particular function defined in the region, or other desirable qualities. Key words: Grid generation
1. Introduction
To solve numerically partial differential equations in an irregular region, one usually starts by covering the region with a discrete grid; the differential operators and solution are then approximated by difference operators and discrete solution values solved for at the grid points. A good grid reduces the errors in the numerical approximation of the differential operators and often allows the equations to be accurately solved with fewer mesh points than would otherwise be required. In this paper we will describe a new method based on a divide-and-conquer (DAC) algorithm for generating a logically rectangular curvilinear mesh for solving two-dimensional problems in the plane bounded by four smooth boundary curves. In this algorithm, at each step a single grid curve or line is generated using position and derivative information derived from previously generated curves. After all grid curves have been generated, the mesh is modified by a rezoning
* Current address: Department of Mathematics, [email protected]. * * Corresponding author. E-mail: [email protected]. 0168.9274/94/$07.00
Ohio
State
0 1994 Elsevier Science B.V. All rights reserved
SSDI 0168-9274(93)E0057-9
University,
Columbus,
OH
43210, USA.
E-mail:
126
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
algorithm to improve smoothness, orthogonality, resolution of a particular function defined in the region, or other desirable qualities. The DAC algorithm is simpler to implement than most of the existing methods, such as the ones surveyed in [5,6], and offers a quick efficient methodology to generate a grid for very distorted regions. For a simple rectangular region, a tensor-product grid (i.e., one with points (xi, yj), where the indices i and j vary independently) is often suitable. For complicated regions, a grid with an arbitrary interconnection structure has many advantages, but the algorithms and software for solving equations on these unstructured grids are more complicated and difficult to develop. An intermediate possibility is the logically rectangular grid, where the grid points are indexed by i and j and the neighbors of grid point pij = (xij, yij> are pi f I,j and pi,j + 1. Such a grid (or a combination of such grids) can be used on a wide variety of regions. Software for working with functions defined on such grids (PDE solvers, graphics packages, differentiation and integration routines, etc.) is readily available Because of the simplicity of the data structure and the versatility of these logically rectangular grids, our initial analysis and implementation of the DAC algorithm is for this case. In the next section we give a brief description of the DAC algorithm, leaving out many of the details, to give a broad overview of how the algorithm functions. We then present a more detailed mathematical description of the method and describe some numerical examples.
2. Algorithm: overview The DAC algorithm sequentially builds a grid by starting with four smooth boundary curves and then successively generates new grid curves (the points pij for one fixed i or for one fixed j) one at a time; each curve provides information used to generate succeeding curves. The grid curves are determined by shape-preserving piecewise-cubic Hermite interpolation [l] through the known points where the curve intersects preceding curves, including the initially specified boundary curves. Once all grid curves have been generated, a rezoning algorithm [3] can be applied to the resulting grid to improve smoothness, orthogonality, resolution of a particular function in the region, or other desirable properties. The DAC algorithm generates grid points pij = (xij, yij) for 1 < i GM, 1 G j G N, that cover a region with specified boundary curves; originally the points pij for i = 1, i = M, j = 1, or j = N are given. At each stage of the algorithm the points pij have been already determined for certain i’s and for certain j’s, and a new grid curve in the i-direction (i.e., all points pij for some fixed j) or in the j-direction (all pij for some fixed i) will be generated. Fig. 1 shows an example of how the DAC algorithm starts with the initial boundary curves (Fig. l(a)> and fi11s in the interior grid, one curve at a time. At each stage, the new grid curve is added that is, in some sense, most needed and will most reduce the distance between the existing grid curves. Note that this approach is similar to what many people do instinctively when asked to generate a grid when given only the boundary curves. Many of the heuristic decisions made in designing the algorithm were chosen to mimic this instinctive process. The first step in each stage is to decide which new grid curve to add. In order to allow the greatest flexibility in choosing the new curve, while providing the most guidance in the choice of succeeding curves, we choose to add the grid curve which will most reduce the distance
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (19941 125-134
127
Fig. 1. Starting with the initial boundary curves, new grid curves are added one at a time. Each new grid curve is chosen to reduce the maximum distance between the existing curves and is orthogonal to the boundary curves. After the initial grid is completed, it is improved by rezoning. (a) The four initially specified boundary curves. (b) The first new grid line is added in the j-direction to bisect the longer dimension at the position where the distance between the upper and lower boundary is minimized. (c) More new curves are added in the j-direction until the distance between them is comparable to the distance between the i-direction curves. (d) The first new curve in i-direction that bisects the curves in the j-direction when the maximum distance between the i-curves is greater than the distance between the j-curves. (e) Additional curves are added until the desired resolution is reached. (f) The grid in Fig. l(e) is rezoned to improve the smoothness of the grid.
between the current existing grid curves. That is, we will add a grid curve that is farthest from the existing grid curves. The sense in which we mean farthest here is based on defining the distance between curve number j and curve number j + 1 in the i-direction to be the maximum of the Euclidean distances dist(~,~, P,,~+ 1> and dist(p,, P~,~+ 1> (and similarly for the curves in the j-direction). For example, in the case of an O-grid (a grid with two opposite boundary curves forced to coincide, such a grid is used for an annulus or similar region), this strategy causes several radial curves to be generated before any curves in the other direction are generated. To add grid curve number j in the i-direction, when the points pij for certain values of i are already known, a shape-preserving parametric cubic Hermite interpolant [l] through the known points is defined and then evaluated at the remaining values of i. At the knots, the tangent vectors (derivatives) along known grid curves in the same direction as the curve being generated and are easily approximated by finite difference methods [4]. The resulting interpolants have continuous first derivatives, and because we have the freedom to choose a tangent vector at
128
R.L. Dougherfy, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
Fig. 2. The derivative vectors used by the parametric interpolant must be scaled by arc length to insure that the interpolated grid points are evenly spaced. (a) The unit tangent and normal vectors to a grid curve are used to define a local coordinate system to eliminate any bias toward a given coordinate system. (b) The tangent and normal vectors used by the parametric interpolant are scaled by arc length in proportion to the lengths of the desired grid curves.
each of the known points, we can require extra properties of the grid curve, such as orthogonality to previous curves or having a specific angle of intersection with the boundary curves. These tangent vectors could also be specified in advance along boundary or reference grid lines to incorporate some desired properties into the grid. When interpolating a tangent vector along a grid curve to approximate the tangent at other points along the curve, instead of using the standard Cartesian coordinates of the vector, we use the coordinate system determined by the unit tangent and normal to the grid curve to eliminate bias toward a specific coordinate system. The latter point is illustrated in Fig. 2(a), because the vector being interpolated is almost orthogonal to the grid curve at the two points where it is known, it will end up being almost orthogonal to the grid curve at the nearby intermediate points. The initial boundary curves for Fig. 1 (Fig. 2(b)) is another such example; interpolating the known vectors in the standard coordinate system would result in zero vectors at the intermediate points. Fig. 2(b) also illustrates another potential problem: in this case, if we interpolated the derivative vectors in the curved coordinate system without normalizing them to the appropriate length scale, then the interpolated tangent vectors would have the same length as the original vectors. The resulting grid curve would have unevenly spaced grid points (too few near the ends and too many in the middle). Therefore, when interpolating these vectors, we also rescale them in proportion to the lengths of the desired grid curves; this will be specified more precisely in the next section. To prevent the interpolant from oscillating between the knots, we use a monotonicity-preserving parametric cubic Hermite interpolant [l]. The monotonicity preservation insures that between the existing grid points (knots) the interpolant preserves the sign of the tangents at the knots. This interpolation is used both for generating the new tangent vectors and for producing the grid points on the new curve. We use a piecewise monotonicity constraint that requires the interpolants to be monotone between each pair of consecutive given data points; as a result, the interpolation will also preserve positivity and negativity of the given data. The sign preservation is important for the interpolation of derivative values; the normal component of these values should have a constant sign along a grid line, and the interpolation must preserve this.
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
129
This completes our overview description of the DAC algorithms for obtaining one new grid curve. The process is repeated until the desired grid resolution is reached. The grid generated by the above process is excellent for moderately distorted regions and can be used directly in solving PDEs. However, in extreme cases it may have rough areas or even singularities caused by the grid curves crossing. When this occurs, we apply the rezoning techniques of [3] to make the grid smoother and, optionally, improve its orthogonality or its resolution of a given function. The rezoner may optionally change the angles at which the grid curves intersect the boundary curves or the positions of the boundary points on the boundary curve. We now recap the procedure for the grid generated in Fig. 1 for the DAC algorithm applied to the boundary curves given in Fig. 2(b). Fig. l(a) shows the grid points on the boundary. Fig. l(b) shows the first grid curve added, in the j-direction; the even spacing of grid points on the curve is the result of the resealing during the derivative interpolation. As shown in Fig. l(c), because the length of the top boundary curve is much longer than the length of the side boundaries, six more curves in the j-direction are added before one is added in the i-direction (Fig. l(d)). Note that, since the original boundary curves were orthogonal to each other, the interpolated tangent vectors will also be orthogonal to each other; however unless the interpolation algorithm is designed to retain this property, the orthogonal&y will be lost in the process, as can be seen in Fig. l(e), which shows the completed DAC grid. Finally, Fig. l(f) shows the results of rezoning, where orthogonality was not specified and therefore not preserved. (The grid curves would not have been so curved at the ends if we had allowed the rezoner to move the boundary points along the boundary curves.>
3. Algorithm: details In this section we again describe the DAC algorithm for generating the initial grid, this time giving details that were left out of the preceding section and describing more options. Initially the boundary points of the grid must be specified. One can also specify one or more of the interior grid curves; this is useful when one wants the grid to be a tensor-product grid on some subregion, or when one wants to force the grid to be finer in a subregion. We refer to the boundary or interior curves specified in advance as reference curves. The derivatives at the grid points on the reference curves (in the direction of the grid curves not specified) may be specified in advance (in which case the shape-preserving interpolation routines will not change them); if not, the program will approximate them by interpolation [4]. If the reference curves are orthogonal, then the algorithm may be required to insure that all computed tangent vectors are orthogonal to the given or previously computed grid curves, with magnitudes determined by interpolation. The derivatives for boundary curves may be specified while the interior derivatives are left unspecified. When using local Lagrange polynomial to approximate tangent derivatives, we have the freedom to choose the independent variable with respect to that which we are interpolating. Interpolating with respect to i or j does not work well for uneven grid spacing. To handle uneven grid spacing, we define an independent variable r for interpolation in the i-direction and an independent variable s for the j-direction; that is, we define values r(1) < r(2) < . * * < r(M) and s(1)
130
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
lative maximum distance between the grid lines worked extremely well. Here we set r(1) = 0 and r(i + 1) = r(i) + maxdist( pi;, Pi+l,j) + EL, i where j runs over those values for which grid curve number j in the i-direction is a reference curve. The EL regularization term ensures that r(i + 1) > r(i) even if some grid points coincide; L is the total length of the boundary, and E is a small constant (but no smaller than the square root of machine roundoff). We define s analogously. Other definitions of r and s are possible (and allowed by the DAC code). Using the same distance parameters Y and s for interpolation along all grid curves may not be appropriate if the grid point spacing varies widely from curve to curve. Therefore, we also define an independent variable rj for interpolation along curve number j in the i-direction, and an independent variable si for curve number i in the j-direction. If grid curve number j in the i-direction is a reference curve, then define Pj by F;(l) = 0 and ?;(i + 1) = F;(i) + dist(Pij, ~~+i,~) + E[ r(i + 1) - r(i)], and rescale to [O, 11 by letting r;(i) = Fj(i)/Fj(M). (Note that approximating the arc length may result in inappropriate or at least noninvariant results if the two dimensions in the plane are incomparable; however, a curvilinear grid is rarely appropriate for such regions.) We obtain rj for other j’s by yet another application of monotonicity-preserving piecewise-cubic Hermite interpolation. Note that the monotonicity being preserved is not monotonic&y of each rj (i.e., rj(i + 1) > rj(i)), so these rj values must be checked and, if necessary, adjusted (we adjust a nonmonotone rj by taking a weighted average of it with what one would have obtained by piecewise-linear instead of piecewise-cubic interpolation). The values sJj> are defined in the same way. We now describe the algorithm for adding one new grid curve in more detail than in the preceding section. First we decide which curve to add, by finding the grid curve whose distance from the previously generated grid curves (in terms of the distances given by Y or s> is largest. If this distance is less than a certain threshold value (which we set to 0.15 min{r(M> - r(l), s(N) - s(l)) for simple domains), then the divide-and-conquer approach is stopped and the remaining curves in the i-direction are added successively; otherwise, the curve most distant from existing grid curves is added. Suppose we have decided to add grid curve number j in the i-direction (the other case is analogous). For each i for which we have already determined grid point pii, we need a tangent vector vii for the new curve in the i-direction. By using finite difference approximations, we compute tangent vectors uiViJfor those j’ for which we have curve number j’ in the i-direction, as well as tangent vectors witi, along the curves in the j-direction for those i’ for which we have these curves. The finite difference formulas used should give tangent vectors which have approximately the same direction as those of the piecewise-linear interpolant; the second-order centered formulas for nonuniform grids [4] have this property. Some of the vectors Vi; may have been specified in advance; the rest must be produced by interpolating the known values Vi;‘. We now describe the algorithm for computing vij for a fixed i. (If Uij is required to be
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
131
orthogonal to wij but its magnitude is not given, then this magnitude is computed by the following algorithm.) For each j’, let ej, be the (signed) angle from the vector (1, 0) to wij,. Define i, and i, to be, respectively, the greatest i’ < i and the least i’ > i such that grid curve number i’ in the j-direction is known. (Let i, = 1 if i = 1, and let i, = M if i = M.) For each i’ for which grid curve number j’ in the i-direction is known, define Bij, to be the result of rotating ujjVby -dj,, and let Gij,= qij,/dist(pilj., piZjO. Interpolate the values cij, with respect to the independent variable si, using monotonicity-preserving piecewise-cubic Hermite interpolation, to get CI,,; then define vii to be fiij dist(Pi,j, piZj> rotated by yj. Special actions must be taken when wijJ or dist(Pilj,, piZj,) is zero (or nearly zero). For example, if the first grid curve to be constructed for an O-grid is in the nonradial direction, then we will have i, = 1, i, = M, and pij, =pMjMj’ for all j’. In this case we can omit the scaling and obtain uij by interpolating the values cij, and then rotating by ej. The case where wijV= 0 is somewhat trickier. If this occurs for isolated values of j’, then we define ej, by setting it equal to oj,, for some j” near j’, or by using, say, the vector from pi,jl_k to pi j,+k for some small k instead of wijf. If all of the values wjjf are zero (as happens when generating a polar grid as an O-grid by mapping one of the four boundary curves to the origin), then the angles 4, may be obtained from vii, instead of Wij~where vii, is known, and 8j can be determined by mterpolation and carefully adding the right multiples of 27r to the angles Oj, to make the interpolation work out correctly (we do so by examining winding numbers of the boundary curves). Once the vectors uij are determined, the new grid curve is constructed by shape-preserving (or partially shape-preserving, since some of the vectors vii may be unmodifiable) cubic Hermite interpolation, this time with respect to rj. This process is repeated until the initial grid has been constructed, and then, if desired, the grid can be rezoned to improve smoothness, orthogonality, etc.
4. Examples In this section we describe three examples of grids generated by the DAC algorithm and the effect of rezoning on these grids. The first example (Figs. 3(a) and 3(b)) shows that the initial algorithm produces the expected result when generating a polar grid in Cartesian coordinates. This example is trivial when the grid is generated in polar coordinates, but is difficult when posed in Cartesian coordinates with the singular boundary. The boundary curves used are from (0, 0) to (1, 01, from (1, 0) around the unit circle back to (1, 01, from (1, 0) to (0, 01, and constant at (0, 01, respectively. (The edge mapped entirely to the origin is one of the special cases mentioned in the previous section; winding numbers were used for the interpolation of tangent vectors across this edge.) Because the grid is defined in Cartesian coordinates, the rezoning algorithm does not recognize or preserve radial symmetry. The second example (Figs. 3(c) and 3(d)) is a more difficult case where we found it useful to change the constant used in the monotonicity algorithm (described in [l]> from 2 to 0.6 and require the grid curves to be orthogonal to the boundaries. Without these changes, the resulting initial grid was more irregular; however, this was mild enough that the rezoning algorithm produced a grid indistinguishable from Fig. 3(d).
132
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) X25-134
The third example (Figs. 3(e) and 3(f)) is a case where an orthogonal grid was desired, and the rezoning parameters were specified accordingly. Note that, since the reference curves were orthogonal, the divide-and-conquer algorithm attempted to make all curves orthogonal, but this
(d)
Fig. 3. Two example grids generated by the algorithm. Initial DAC grid and (d) DAC grid after rezoning.
(a) Initial
DAC grid and (b) DAC grid after rezoning.
CC>
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
Fig. 3 (continued). rezoning.
Another
example
grid generated
by the algorithm.
133
(e) Initial DAC grid and (f) DAC grid after
became increasingly more difficult as more curves were added. However, the rezoner had no trouble restoring the orthogonal&y that was lost in the initial grid generation.
5. Conclusion
The DAC algorithm described in this paper, together with the rezoning described in [3], produces highly acceptable grids for a wide variety of situations. The numerous parameters allow direct control of the grid generation process; the default values for all of these parameters are robust enough to allow easy generation of a grid that will be acceptable in most cases. The DAC algorithm is appropriate for most simply-connected regions (and even some that are not simply connected if branch cuts are introduced) for which it is reasonable to use a single logically rectangular grid. More complicated or irregular regions can be divided into subregions, with a logically rectangular grid on each subregion. We have found it useful to combine the initial grid generator and rezoner with an interactive graphical display package. This allows the user to interact with the programs, displaying the results of using the current parameters and then changing one or more parameters if necessary to obtain an improved grid. Once the user has produced a suitable grid (note that there is a high correlation between visual acceptability and suitability for PDE computations), the results can be saved and passed on to other software. The methods of this paper can be generalized to three dimensions with moderate effort. Just as the problem of generating a two-dimensional grid was reduced to that of generating a single grid curve in the plane, so the problem of generating a three-dimensional grid reduces to the problem of generating a surface in space with some curves already specified. This is little different from generating a grid in the plane, except that the problem of generating suitable grid curves in space will be more fragile than that in the plane; probably several variants of
134
R.L. Dougherty,J.M. Hyman /Applied Numerical Mathematics14 (1994) 125-134
shape-preserving interpolation, such as convexity preserving interpolation [l], will have to be tried before a suitable one is found. The rezoning step probably should be applied to each surface as it is generated (without moving the given curves) as well as to the full three-dimensional grid.
6. Acknowledgement This work was performed under the auspices of the U.S. Department of Energy under contract W-7405-ENG-36 and the Office of Basic Energy Sciences, Applied Mathematical Sciences Program KC-07-01-01.
7. References [l] R.L. Dougherty, A.S. Edelman and J.M. Hyman, Nonnegativity-, monotonicity-, or convexity-preserving cubic and quintic Hermite interpolation, Math. Camp. 52 (1989) 471-494. [2] R.L. Dougherty and J.M. Hyman, DAC: a grid generation program using the divide-and-conquer algorithm, in preparation. [3] W.D. Henshaw and J.M. Hyman, Static rezone methods on logically rectangular grids, Los Alamos National Laboratory Report (1986). [4] J.M. Hyman and B. Larrouturou, The numerical differentiation of discrete functions using polynomial interpolation methods, in: J.F. Thompson, ed., Numerical Grid Generation for Numerical Solution of Partial Differential Equations (Elsevier North-Holland, New York, 1982) 487-506. [5] J.F. Thompson, ed., Numerical Grid Generation, for Numerical Solution of Partial Differential Equations (Elsevier North-Holland, New York, 1982). [6] J.F. Thompson, Z.U.A. Warsi and C.W. Mastin, Numerical Grid Generation: Foundations and Applications (Elsevier North-Holland, New York, 1985).
ELSEVIER
A divide-and-conquer Randall L. Dougherty
algorithm for grid generation * and James M. Hyman **
Theoretical Division, MS B284, Center for Nonlinear Studies, Los Alamos, NM 87545, USA
Abstract We present a new algorithm for generating a logically rectangular curvilinear mesh for a region in the plane bounded by four smooth curves. The curvilinear grid is generated using a divide-and-conquer algorithm where at each step a single grid curve or line is generated using position and derivative information derived from previously generated curves. The grid curves are produced by shape-preserving cubic Hermite interpolation that allow for corners or kinks along the boundary. After all grid curves have been generated, the mesh is modified by a rezoning algorithm to improve smoothness, orthogonality, resolution of a particular function defined in the region, or other desirable qualities. Key words: Grid generation
1. Introduction
To solve numerically partial differential equations in an irregular region, one usually starts by covering the region with a discrete grid; the differential operators and solution are then approximated by difference operators and discrete solution values solved for at the grid points. A good grid reduces the errors in the numerical approximation of the differential operators and often allows the equations to be accurately solved with fewer mesh points than would otherwise be required. In this paper we will describe a new method based on a divide-and-conquer (DAC) algorithm for generating a logically rectangular curvilinear mesh for solving two-dimensional problems in the plane bounded by four smooth boundary curves. In this algorithm, at each step a single grid curve or line is generated using position and derivative information derived from previously generated curves. After all grid curves have been generated, the mesh is modified by a rezoning
* Current address: Department of Mathematics, [email protected]. * * Corresponding author. E-mail: [email protected]. 0168.9274/94/$07.00
Ohio
State
0 1994 Elsevier Science B.V. All rights reserved
SSDI 0168-9274(93)E0057-9
University,
Columbus,
OH
43210, USA.
E-mail:
126
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
algorithm to improve smoothness, orthogonality, resolution of a particular function defined in the region, or other desirable qualities. The DAC algorithm is simpler to implement than most of the existing methods, such as the ones surveyed in [5,6], and offers a quick efficient methodology to generate a grid for very distorted regions. For a simple rectangular region, a tensor-product grid (i.e., one with points (xi, yj), where the indices i and j vary independently) is often suitable. For complicated regions, a grid with an arbitrary interconnection structure has many advantages, but the algorithms and software for solving equations on these unstructured grids are more complicated and difficult to develop. An intermediate possibility is the logically rectangular grid, where the grid points are indexed by i and j and the neighbors of grid point pij = (xij, yij> are pi f I,j and pi,j + 1. Such a grid (or a combination of such grids) can be used on a wide variety of regions. Software for working with functions defined on such grids (PDE solvers, graphics packages, differentiation and integration routines, etc.) is readily available Because of the simplicity of the data structure and the versatility of these logically rectangular grids, our initial analysis and implementation of the DAC algorithm is for this case. In the next section we give a brief description of the DAC algorithm, leaving out many of the details, to give a broad overview of how the algorithm functions. We then present a more detailed mathematical description of the method and describe some numerical examples.
2. Algorithm: overview The DAC algorithm sequentially builds a grid by starting with four smooth boundary curves and then successively generates new grid curves (the points pij for one fixed i or for one fixed j) one at a time; each curve provides information used to generate succeeding curves. The grid curves are determined by shape-preserving piecewise-cubic Hermite interpolation [l] through the known points where the curve intersects preceding curves, including the initially specified boundary curves. Once all grid curves have been generated, a rezoning algorithm [3] can be applied to the resulting grid to improve smoothness, orthogonality, resolution of a particular function in the region, or other desirable properties. The DAC algorithm generates grid points pij = (xij, yij) for 1 < i GM, 1 G j G N, that cover a region with specified boundary curves; originally the points pij for i = 1, i = M, j = 1, or j = N are given. At each stage of the algorithm the points pij have been already determined for certain i’s and for certain j’s, and a new grid curve in the i-direction (i.e., all points pij for some fixed j) or in the j-direction (all pij for some fixed i) will be generated. Fig. 1 shows an example of how the DAC algorithm starts with the initial boundary curves (Fig. l(a)> and fi11s in the interior grid, one curve at a time. At each stage, the new grid curve is added that is, in some sense, most needed and will most reduce the distance between the existing grid curves. Note that this approach is similar to what many people do instinctively when asked to generate a grid when given only the boundary curves. Many of the heuristic decisions made in designing the algorithm were chosen to mimic this instinctive process. The first step in each stage is to decide which new grid curve to add. In order to allow the greatest flexibility in choosing the new curve, while providing the most guidance in the choice of succeeding curves, we choose to add the grid curve which will most reduce the distance
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (19941 125-134
127
Fig. 1. Starting with the initial boundary curves, new grid curves are added one at a time. Each new grid curve is chosen to reduce the maximum distance between the existing curves and is orthogonal to the boundary curves. After the initial grid is completed, it is improved by rezoning. (a) The four initially specified boundary curves. (b) The first new grid line is added in the j-direction to bisect the longer dimension at the position where the distance between the upper and lower boundary is minimized. (c) More new curves are added in the j-direction until the distance between them is comparable to the distance between the i-direction curves. (d) The first new curve in i-direction that bisects the curves in the j-direction when the maximum distance between the i-curves is greater than the distance between the j-curves. (e) Additional curves are added until the desired resolution is reached. (f) The grid in Fig. l(e) is rezoned to improve the smoothness of the grid.
between the current existing grid curves. That is, we will add a grid curve that is farthest from the existing grid curves. The sense in which we mean farthest here is based on defining the distance between curve number j and curve number j + 1 in the i-direction to be the maximum of the Euclidean distances dist(~,~, P,,~+ 1> and dist(p,, P~,~+ 1> (and similarly for the curves in the j-direction). For example, in the case of an O-grid (a grid with two opposite boundary curves forced to coincide, such a grid is used for an annulus or similar region), this strategy causes several radial curves to be generated before any curves in the other direction are generated. To add grid curve number j in the i-direction, when the points pij for certain values of i are already known, a shape-preserving parametric cubic Hermite interpolant [l] through the known points is defined and then evaluated at the remaining values of i. At the knots, the tangent vectors (derivatives) along known grid curves in the same direction as the curve being generated and are easily approximated by finite difference methods [4]. The resulting interpolants have continuous first derivatives, and because we have the freedom to choose a tangent vector at
128
R.L. Dougherfy, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
Fig. 2. The derivative vectors used by the parametric interpolant must be scaled by arc length to insure that the interpolated grid points are evenly spaced. (a) The unit tangent and normal vectors to a grid curve are used to define a local coordinate system to eliminate any bias toward a given coordinate system. (b) The tangent and normal vectors used by the parametric interpolant are scaled by arc length in proportion to the lengths of the desired grid curves.
each of the known points, we can require extra properties of the grid curve, such as orthogonality to previous curves or having a specific angle of intersection with the boundary curves. These tangent vectors could also be specified in advance along boundary or reference grid lines to incorporate some desired properties into the grid. When interpolating a tangent vector along a grid curve to approximate the tangent at other points along the curve, instead of using the standard Cartesian coordinates of the vector, we use the coordinate system determined by the unit tangent and normal to the grid curve to eliminate bias toward a specific coordinate system. The latter point is illustrated in Fig. 2(a), because the vector being interpolated is almost orthogonal to the grid curve at the two points where it is known, it will end up being almost orthogonal to the grid curve at the nearby intermediate points. The initial boundary curves for Fig. 1 (Fig. 2(b)) is another such example; interpolating the known vectors in the standard coordinate system would result in zero vectors at the intermediate points. Fig. 2(b) also illustrates another potential problem: in this case, if we interpolated the derivative vectors in the curved coordinate system without normalizing them to the appropriate length scale, then the interpolated tangent vectors would have the same length as the original vectors. The resulting grid curve would have unevenly spaced grid points (too few near the ends and too many in the middle). Therefore, when interpolating these vectors, we also rescale them in proportion to the lengths of the desired grid curves; this will be specified more precisely in the next section. To prevent the interpolant from oscillating between the knots, we use a monotonicity-preserving parametric cubic Hermite interpolant [l]. The monotonicity preservation insures that between the existing grid points (knots) the interpolant preserves the sign of the tangents at the knots. This interpolation is used both for generating the new tangent vectors and for producing the grid points on the new curve. We use a piecewise monotonicity constraint that requires the interpolants to be monotone between each pair of consecutive given data points; as a result, the interpolation will also preserve positivity and negativity of the given data. The sign preservation is important for the interpolation of derivative values; the normal component of these values should have a constant sign along a grid line, and the interpolation must preserve this.
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
129
This completes our overview description of the DAC algorithms for obtaining one new grid curve. The process is repeated until the desired grid resolution is reached. The grid generated by the above process is excellent for moderately distorted regions and can be used directly in solving PDEs. However, in extreme cases it may have rough areas or even singularities caused by the grid curves crossing. When this occurs, we apply the rezoning techniques of [3] to make the grid smoother and, optionally, improve its orthogonality or its resolution of a given function. The rezoner may optionally change the angles at which the grid curves intersect the boundary curves or the positions of the boundary points on the boundary curve. We now recap the procedure for the grid generated in Fig. 1 for the DAC algorithm applied to the boundary curves given in Fig. 2(b). Fig. l(a) shows the grid points on the boundary. Fig. l(b) shows the first grid curve added, in the j-direction; the even spacing of grid points on the curve is the result of the resealing during the derivative interpolation. As shown in Fig. l(c), because the length of the top boundary curve is much longer than the length of the side boundaries, six more curves in the j-direction are added before one is added in the i-direction (Fig. l(d)). Note that, since the original boundary curves were orthogonal to each other, the interpolated tangent vectors will also be orthogonal to each other; however unless the interpolation algorithm is designed to retain this property, the orthogonal&y will be lost in the process, as can be seen in Fig. l(e), which shows the completed DAC grid. Finally, Fig. l(f) shows the results of rezoning, where orthogonality was not specified and therefore not preserved. (The grid curves would not have been so curved at the ends if we had allowed the rezoner to move the boundary points along the boundary curves.>
3. Algorithm: details In this section we again describe the DAC algorithm for generating the initial grid, this time giving details that were left out of the preceding section and describing more options. Initially the boundary points of the grid must be specified. One can also specify one or more of the interior grid curves; this is useful when one wants the grid to be a tensor-product grid on some subregion, or when one wants to force the grid to be finer in a subregion. We refer to the boundary or interior curves specified in advance as reference curves. The derivatives at the grid points on the reference curves (in the direction of the grid curves not specified) may be specified in advance (in which case the shape-preserving interpolation routines will not change them); if not, the program will approximate them by interpolation [4]. If the reference curves are orthogonal, then the algorithm may be required to insure that all computed tangent vectors are orthogonal to the given or previously computed grid curves, with magnitudes determined by interpolation. The derivatives for boundary curves may be specified while the interior derivatives are left unspecified. When using local Lagrange polynomial to approximate tangent derivatives, we have the freedom to choose the independent variable with respect to that which we are interpolating. Interpolating with respect to i or j does not work well for uneven grid spacing. To handle uneven grid spacing, we define an independent variable r for interpolation in the i-direction and an independent variable s for the j-direction; that is, we define values r(1) < r(2) < . * * < r(M) and s(1)
130
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
lative maximum distance between the grid lines worked extremely well. Here we set r(1) = 0 and r(i + 1) = r(i) + maxdist( pi;, Pi+l,j) + EL, i where j runs over those values for which grid curve number j in the i-direction is a reference curve. The EL regularization term ensures that r(i + 1) > r(i) even if some grid points coincide; L is the total length of the boundary, and E is a small constant (but no smaller than the square root of machine roundoff). We define s analogously. Other definitions of r and s are possible (and allowed by the DAC code). Using the same distance parameters Y and s for interpolation along all grid curves may not be appropriate if the grid point spacing varies widely from curve to curve. Therefore, we also define an independent variable rj for interpolation along curve number j in the i-direction, and an independent variable si for curve number i in the j-direction. If grid curve number j in the i-direction is a reference curve, then define Pj by F;(l) = 0 and ?;(i + 1) = F;(i) + dist(Pij, ~~+i,~) + E[ r(i + 1) - r(i)], and rescale to [O, 11 by letting r;(i) = Fj(i)/Fj(M). (Note that approximating the arc length may result in inappropriate or at least noninvariant results if the two dimensions in the plane are incomparable; however, a curvilinear grid is rarely appropriate for such regions.) We obtain rj for other j’s by yet another application of monotonicity-preserving piecewise-cubic Hermite interpolation. Note that the monotonicity being preserved is not monotonic&y of each rj (i.e., rj(i + 1) > rj(i)), so these rj values must be checked and, if necessary, adjusted (we adjust a nonmonotone rj by taking a weighted average of it with what one would have obtained by piecewise-linear instead of piecewise-cubic interpolation). The values sJj> are defined in the same way. We now describe the algorithm for adding one new grid curve in more detail than in the preceding section. First we decide which curve to add, by finding the grid curve whose distance from the previously generated grid curves (in terms of the distances given by Y or s> is largest. If this distance is less than a certain threshold value (which we set to 0.15 min{r(M> - r(l), s(N) - s(l)) for simple domains), then the divide-and-conquer approach is stopped and the remaining curves in the i-direction are added successively; otherwise, the curve most distant from existing grid curves is added. Suppose we have decided to add grid curve number j in the i-direction (the other case is analogous). For each i for which we have already determined grid point pii, we need a tangent vector vii for the new curve in the i-direction. By using finite difference approximations, we compute tangent vectors uiViJfor those j’ for which we have curve number j’ in the i-direction, as well as tangent vectors witi, along the curves in the j-direction for those i’ for which we have these curves. The finite difference formulas used should give tangent vectors which have approximately the same direction as those of the piecewise-linear interpolant; the second-order centered formulas for nonuniform grids [4] have this property. Some of the vectors Vi; may have been specified in advance; the rest must be produced by interpolating the known values Vi;‘. We now describe the algorithm for computing vij for a fixed i. (If Uij is required to be
R.L. Dougherty, J.M. Hyman /Applied
Numerical Mathematics
14 (1994) 125-134
131
orthogonal to wij but its magnitude is not given, then this magnitude is computed by the following algorithm.) For each j’, let ej, be the (signed) angle from the vector (1, 0) to wij,. Define i, and i, to be, respectively, the greatest i’ < i and the least i’ > i such that grid curve number i’ in the j-direction is known. (Let i, = 1 if i = 1, and let i, = M if i = M.) For each i’ for which grid curve number j’ in the i-direction is known, define Bij, to be the result of rotating ujjVby -dj,, and let Gij,= qij,/dist(pilj., piZjO. Interpolate the values cij, with respect to the independent variable si, using monotonicity-preserving piecewise-cubic Hermite interpolation, to get CI,,; then define vii to be fiij dist(Pi,j, piZj> rotated by yj. Special actions must be taken when wijJ or dist(Pilj,, piZj,) is zero (or nearly zero). For example, if the first grid curve to be constructed for an O-grid is in the nonradial direction, then we will have i, = 1, i, = M, and pij, =pMjMj’ for all j’. In this case we can omit the scaling and obtain uij by interpolating the values cij, and then rotating by ej. The case where wijV= 0 is somewhat trickier. If this occurs for isolated values of j’, then we define ej, by setting it equal to oj,, for some j” near j’, or by using, say, the vector from pi,jl_k to pi j,+k for some small k instead of wijf. If all of the values wjjf are zero (as happens when generating a polar grid as an O-grid by mapping one of the four boundary curves to the origin), then the angles 4, may be obtained from vii, instead of Wij~where vii, is known, and 8j can be determined by mterpolation and carefully adding the right multiples of 27r to the angles Oj, to make the interpolation work out correctly (we do so by examining winding numbers of the boundary curves). Once the vectors uij are determined, the new grid curve is constructed by shape-preserving (or partially shape-preserving, since some of the vectors vii may be unmodifiable) cubic Hermite interpolation, this time with respect to rj. This process is repeated until the initial grid has been constructed, and then, if desired, the grid can be rezoned to improve smoothness, orthogonality, etc.
4. Examples In this section we describe three examples of grids generated by the DAC algorithm and the effect of rezoning on these grids. The first example (Figs. 3(a) and 3(b)) shows that the initial algorithm produces the expected result when generating a polar grid in Cartesian coordinates. This example is trivial when the grid is generated in polar coordinates, but is difficult when posed in Cartesian coordinates with the singular boundary. The boundary curves used are from (0, 0) to (1, 01, from (1, 0) around the unit circle back to (1, 01, from (1, 0) to (0, 01, and constant at (0, 01, respectively. (The edge mapped entirely to the origin is one of the special cases mentioned in the previous section; winding numbers were used for the interpolation of tangent vectors across this edge.) Because the grid is defined in Cartesian coordinates, the rezoning algorithm does not recognize or preserve radial symmetry. The second example (Figs. 3(c) and 3(d)) is a more difficult case where we found it useful to change the constant used in the monotonicity algorithm (described in [l]> from 2 to 0.6 and require the grid curves to be orthogonal to the boundaries. Without these changes, the resulting initial grid was more irregular; however, this was mild enough that the rezoning algorithm produced a grid indistinguishable from Fig. 3(d).
132
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) X25-134
The third example (Figs. 3(e) and 3(f)) is a case where an orthogonal grid was desired, and the rezoning parameters were specified accordingly. Note that, since the reference curves were orthogonal, the divide-and-conquer algorithm attempted to make all curves orthogonal, but this
(d)
Fig. 3. Two example grids generated by the algorithm. Initial DAC grid and (d) DAC grid after rezoning.
(a) Initial
DAC grid and (b) DAC grid after rezoning.
CC>
R.L. Dougherty, J.M. Hyman /Applied Numerical Mathematics 14 (1994) 125-134
Fig. 3 (continued). rezoning.
Another
example
grid generated
by the algorithm.
133
(e) Initial DAC grid and (f) DAC grid after
became increasingly more difficult as more curves were added. However, the rezoner had no trouble restoring the orthogonal&y that was lost in the initial grid generation.
5. Conclusion
The DAC algorithm described in this paper, together with the rezoning described in [3], produces highly acceptable grids for a wide variety of situations. The numerous parameters allow direct control of the grid generation process; the default values for all of these parameters are robust enough to allow easy generation of a grid that will be acceptable in most cases. The DAC algorithm is appropriate for most simply-connected regions (and even some that are not simply connected if branch cuts are introduced) for which it is reasonable to use a single logically rectangular grid. More complicated or irregular regions can be divided into subregions, with a logically rectangular grid on each subregion. We have found it useful to combine the initial grid generator and rezoner with an interactive graphical display package. This allows the user to interact with the programs, displaying the results of using the current parameters and then changing one or more parameters if necessary to obtain an improved grid. Once the user has produced a suitable grid (note that there is a high correlation between visual acceptability and suitability for PDE computations), the results can be saved and passed on to other software. The methods of this paper can be generalized to three dimensions with moderate effort. Just as the problem of generating a two-dimensional grid was reduced to that of generating a single grid curve in the plane, so the problem of generating a three-dimensional grid reduces to the problem of generating a surface in space with some curves already specified. This is little different from generating a grid in the plane, except that the problem of generating suitable grid curves in space will be more fragile than that in the plane; probably several variants of
134
R.L. Dougherty,J.M. Hyman /Applied Numerical Mathematics14 (1994) 125-134
shape-preserving interpolation, such as convexity preserving interpolation [l], will have to be tried before a suitable one is found. The rezoning step probably should be applied to each surface as it is generated (without moving the given curves) as well as to the full three-dimensional grid.
6. Acknowledgement This work was performed under the auspices of the U.S. Department of Energy under contract W-7405-ENG-36 and the Office of Basic Energy Sciences, Applied Mathematical Sciences Program KC-07-01-01.
7. References [l] R.L. Dougherty, A.S. Edelman and J.M. Hyman, Nonnegativity-, monotonicity-, or convexity-preserving cubic and quintic Hermite interpolation, Math. Camp. 52 (1989) 471-494. [2] R.L. Dougherty and J.M. Hyman, DAC: a grid generation program using the divide-and-conquer algorithm, in preparation. [3] W.D. Henshaw and J.M. Hyman, Static rezone methods on logically rectangular grids, Los Alamos National Laboratory Report (1986). [4] J.M. Hyman and B. Larrouturou, The numerical differentiation of discrete functions using polynomial interpolation methods, in: J.F. Thompson, ed., Numerical Grid Generation for Numerical Solution of Partial Differential Equations (Elsevier North-Holland, New York, 1982) 487-506. [5] J.F. Thompson, ed., Numerical Grid Generation, for Numerical Solution of Partial Differential Equations (Elsevier North-Holland, New York, 1982). [6] J.F. Thompson, Z.U.A. Warsi and C.W. Mastin, Numerical Grid Generation: Foundations and Applications (Elsevier North-Holland, New York, 1985).