Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
II. BOUNDARYCONFORMING COORDINATE SYSTEMS
1. Basic Concepts
To provide a familiar ground from which to view the general development to follow, consider first a twodimensional cylindrical coordinate system covering the annular region between two concentric circles:
Here the curvilinear coordinates (r,) vary on the intervals [r_{1},r_{2}] and [0,2], respectively. These curvilinear coordinates are related to the cartesian coordinates (x,y) by the transformation equations

(1) 
The inverse transformation is given by

(2) 
Note that one of the curvilinear coordinates, r, is constant on each of the phys1cal boundaries, while the other coordinate, , varies monotonically over the same range around each of the boundaries. Note also that the system can be represented as a rectangle on which the two physical boundaries correspond to the top and bottom sides:
The transformed region, i.e., where the curvilinear coordinates, r and the independent variables, thus can be thought of as being rectangular, and can be treated as such from a coding standpoint. These points will be central to what follows.
The curvilinear coordinates (r,) can be normalized to the interval [0,1] by introducing the new curvilinear coordinates (, ), where

(3) 
or

(4) 
The transformation then may be written

(5a) 

(5b) 
where now and both vary on the interval [0,1]. This is thus a mapping of the annular region between the two circles in the physical space onto the unit square in the transformed space, i.e., each point (x,y) on the annulus corresponds to one, and only one, point (,) on the unit square:
The bottom ( = 0) and top ( = 1) of the square correspond, respectively, to the inner and outer circles, r = r_{1}, and r = r_{2}. The sides of the square, = 0 and = 1 correspond to = 0 and = 2, respectively, and hence to the two coincident sides of a branch cut in the physical space. Therefore, boundary conditions are not to be specified on these sides of the unit square in the transformed space. Rather these sides are to be considered reentrant on each other with points adjacent to one, outside the square, being equivalent to points adjacent to the other, inside the square.
Conceptually, the physical region can be considered to have been opened at the cut = 0 and 2 and then deformed into a rectangle to form the transformed region:
Here, point correspondence across the reentrant boundaries (indicated by the dashed connecting line) in the transformed region is illustrated by the coincidence of the pair of circled points. This conceptual device and mode of illustration for the the point correspondence across reentrant boundaries will serve later for more general configurations.
These simple concepts extend to more complicated twodimensional configurations, the central feature being that one of the curvilinear coordinates is made to be constant on a boundary curve (as was r above), while the other varies monotonically along that boundary curve (as does ). The transformation to the rectangle is achieved by making the range and direction of variation of the varying coordinate the same on each of two opposing boundaries (as varies from 0 to 2 on each circle above).
The physical space thus transforms to the rectangle shown above regardless of the shape of the physical region. (It is not necessary to normalize the curvilinear coordinates to the interval [0,1], and in fact, any normalization can be used. In computational applications the normalization is more conveniently done to different intervals for each coordinate. The field in the transformed space is then rectangular, rather than square.) Familiar examples of this are elliptical coordinates for the region between two confocal ellipses, spherical coordinates for two spheres, parabolic coordinates for two parabolas, etc.
These, same concepts will be extended later to completely general configurations involving any number of boundary curves and branch outs. The extension to three dimensions follows directly, using boundary surfaces instead of curves, i.e., one curvilinear coordinate will be made constant on a boundary surface, with the other two forming a twodimensional coordinate system on the surface.
Returning to the concentric circles, if the functional dependence of on , and/or that of r on , had been made more general than the simple linear normalizations given by Eq. (4), the corresponding coordinate lines would have become unequally spaced in the physical space, while remaining as radial lines and concentric circles:
The transformation, from Eq. (1), is now given by

(6a) 

(6b) 
In this case the points on the inner and outer circular boundaries are not equally spaced around the circles in the physical space for equal increments of , although they remain equally spaced on the top and bottom of the unit square in the transformed space by construction. The spacing around these circles is determined by the functional dependence of on , and, since the points are located at equal increments of by construction, this functional relationship is defined by the placement of these points around the circles. This point, that the coordinate system in the field is determined from the boundary point distribution, will be central to the discussion of grid generation to follow. The distribution of circumferential lines is controlled here by the functional relationship between r and , which is not related to any boundary point distribution. Thus factors other than the boundary point distribution may be expected to be involved in grid generation, as well. That the point distribution on the boundaries may be controlled by direct placement of the points, while the coordinate line distribution in the field must be controlled by other means will also continue to appear in the developments to follow.
The onedimensional functional relationship between and in Eq. (6) requires that the relative distributions of boundary points around the inner and outer circles be the same. This restriction can be removed by making a function of , as well as of , while retaining the periodic nature of the dependence on . In this case the ooordinate lines of constant will no longer be straight radial lines, although they will continue to connect corresponding points on the inner and outer circular boundaries. Similarly the circumferential coordinate lines (lines of constant here) can be made to depart from circles by making r dependent on both and , but with the restriction that the dependence vanishes on the inner and outer circular boundaries (where = 0 and = 1, respectively, here).
Obviously certain constraints will have to be placed on the functions (,) and r(,) to keep the mapping onetoone. All of these considerations will reappear in the general developments that follow.
Finally, it should be realized that the intermediate use here of the cylindrical coordinates (r,) in defining the transformation between the curvilinear coordinates (,) and the cartesian coordinates (x,y) has been only in deference to the familiarity of the cylindrical coordinates, and such intermediary coordinates will not appear in general. The generalized statement for the simple configuration under consideration here is as follows: Find (x,y) and (x,y) in the annular region bounded by the curves x^{2} + y^{2} = and x^{2} + y^{2} = , subject to the boundary conditions
Specified monotonic variation of over [0,1] on x^{2} + y^{2} = and on x^{2} + y^{2} = with same sense of direction on each of these two curves.
It is the inverse problem that will be treated in fact, however, i.e., find x(,) and y(,) on the unit square in the transformed space (0 1, 0 1), subject to the boundary conditions
x(,0) and y(,0) specified on = 0 such that x^{2}( = 0) + y^{2},0) =
x(,1) and y(,1) specified on = 1 such that x^{2}( = 1) + y^{2}(,1) =
Periodicity in : x(1 + ,) = x(,)
y(1 +,) = y(,)
The simple form for the transformation given by Eq. (6) is made possible by choosing the same functional dependence of x and y on on the boundaries, = 0 and = 1. The familiar cylindrical coordinate system is thus a special case of the general grid generation problem for this simple configuration applicable to the region between two concentric circles, as is the elliptical coordinate system for two ellipses, etc.
2. Generalization
Generalizing from the above consideration of cylindrical coordinates, the basic idea of a boundaryconforming curvilinear coordinate system is to have some coordinate line (in 2D, surface in 3D) coincident with each boundary segment, analogous to the way in which lines of constant radial coordinate coincide with circles in the cylindrical coordinate system. The other curvilinear coordinate, analogous to the angular coordinate in the cylindrical system, will vary along the boundary segment and clearly must do so monotonically, else the same pair of values of the curvilinear coordinates will occur at two different physical points. (It should be clear that the curvilinear coordinate that varies along a boundary segment must have the same direction and range of variation over some opposing segment, e.g., as the angular variable varies from 0 to 2 over both of two concentric circles in cylindrical coordinates).
With the values of the curvilinear coordinates thus specified on the boundary, it then remains to generate values of these coordinates in the field from these boundary values. There must, or course, be a unique correspondence between the cartesian (or other basis system) and the curvilinear coordinates, i.e., the mapping of the physical region onto the transformed region must be onetoone, so that every point in the physical field corresponds to one, and only one, point in the transformed field, and vice versa. Coordinate lines of the same family must not cross, and lines of different families must not cross more than once.
In this chapter a twodimensional region will be considered in most of the discussions in the interest of economy of presentation. Generalization to three dimensions will be evident in most cases and will be mentioned specifically only when necessary. As noted above, the curvilinear coordinates may be normalized to any intervals, just as the radial and angular coordinates of the cylindrical coordinate system can be expressed in many different units. Since the interest of the present discussion is numerical application, it will be generally convenient to define the increments of all the curvilinear coordinates to be uniformly unity, and then to normalize these coordinates to the interval [1,N^{(i)}], where N^{(i)} is the total number of grid points to be used in the ^{i} direction. (The three curvilinear coordinates will be indicated as ^{i}, i = 1,2,3, in general. In two dimensions, however, the notation (, ) will often be used for the two coordinates ^{1} and ^{2} .) The computational field, i.e., the field in the transformed space, thus will have rectangular boundaries and will be covered by a square grid. (It will become clear later that the actual values of the increments in the curvilinear coordinates are immaterial since they do not appear in the final numer1cal expressions. Therefore no generality is lost in making the grid square and of unit increment in the transformed field.)
A. Boundaryvalue Problem  Physical Region
The generation of the curvilinear coordinate system may be treated as follows: with the curvilinear coordinates specified on the boundaries, e.g., (x,y) and (x,y) on a boundary curve (this specification amounting to a constant value for either or on each segment of , with a specified monotonic variation of the other over the segment), generate the values, (x,y) and (x,y), in the field bounded by . This is thus a boundary value problem on the physical field with the curvilinear coordinates (, ) as the dependent variables and the cartesian coordinates (x,y) as the independent variables, with boundary conditions specified on curved boundaries:
(In these discussions, the transformation is assumed to be from cartesian coordinates in the physical space. The transformation can, however, be from any system of coordinates in the physical space.)
B. Boundary value Problem  Transformed Region
The problem may be simplified for computation, however, by first transforming so that the physical cartesian coordinates (x,y) become the dependent variables, with the curvilinear coordinates (,) as the independent variables. Since a constant value of one curvilinear coordinate, with monotonic variation of the other, has been specified on each boundary segment, it follows that these boundary segments in the physical field will correspond to vertical or horizontal lines In the transformed field. Also, since the range of variation of the curvilinear coordinate varying along a boundary segment has been made the same over opposing segments, it follows that the transformed field will be composed of rectangular blocks.
The boundary value problem in the transformed field then involves generating the values of the physical cartesian coordinates, x(,) and y(,), in the transformed field from the specified boundary values of x(,) and y(,) on the rectangular boundary of the transformed field, the boundary being formed of segments of constant or , i.e., vertical or horizontal lines. With = constant on a boundary segment, and the increments in taken to be uniformly unity as discussed above, this boundary value specification is implemented numerically by distributing the points as desired along the boundary segment and then assigning the values of the cartesian coordinates of each successive point as boundary values at the equally spaced boundary points on the bottom (or top) of the transformed field in the following figure.
Boundary values are not specified on the left and right sides of the transformed field since these boundaries are reentrant on each other (analogous to the 0 and 2 lines in the cylindrical system), as discussed above, and as indicated by the connecting dotted line on the figure. Points outside one of these reentrant boundaries are coincident with points at the same distance inside the other. The problem is thus much more simple in the transformed field, since the boundaries there are all rectangular, and the computation in the transformed field thus is on a square grid regardless of the shape of the physical boundaries.
With values of the cartesian coordinates known in the field as functions of the curvilinear coordinates, the network of intersecting lines formed by contours (surfaces in 3D) on which a curvilinear coordinate is constant, i.e., the curvilinear coordinate system, provides the needed organization of the discretization with conformation to the physical boundary. It is also possible to specify intersection angles for the coordinate lines at the boundaries as well as the point locations.
3. Transformed Region Configurations
As noted above, the generation of the curvilinear coordinate system is done by devising a scheme for determination of the field values of the cartesian coordinates from specified values of these coordinates (and/or curvilinear coordinate line intersection angles) on portions of the boundary of the transformed region. Since the boundary of the transformed region is comprised of horizontal and vertical line segments, portions of which correspond to segments of the physical boundary on which a curvilinear coordinate is specified to be constant, it should be evident that the configuration of the resulting coordinate system depends on how the boundary correspondence is made, i.e., how the transformed region is configured.
Some examples of different configurations are given below, from which more complex configurations can be inferred. In these examples only a minimum number of coordinate lines are shown in the interest of clarity of presentation tation. In all of these examples, boundary values of the physical cartesian coordinates (and/or curvilinear coordinate line intersection angles) are understood to be specified on all boundaries, both external and internal, of the transformed region except for segments indicated by dotted lines. These latter segments correspond to branch cuts in the physical space, as is explained in the examples in which they appear. Such reentrant boundary segments always occur in pairs, the members of which are indicated by the dashed connecting lines on each of the configurations shown. Points outside the field across one segmentof such a pair are coincident with points inside the field across the other member of the pair. The conceptual device of opening the physical field at the cuts is used here to help clarify the correspondence between the physical and transformed fields. In many cases an example of an actual coordinate system is given as well. References to the use of various configurations may be found in the surveys given by Ref. [1] and [5], and a number of examples appear in Ref. [2].
A. Simplyconnected Regions
It is natural to define the same curvilinear coordinate to be constant on each member of a pair of generally opposing boundary segments in the physical plane. Thus, a simplyconnected region formed by four curves is logically treated by transforming to an empty rectangle:
Similarly, an Lshaped region could remain Lshaped in the transformed region:
Here, for instance, the cartesian coordinates of the desired points on the physical boundary segment 45 are specified as boundary conditions on the vertical line 45, in corresponding order, which forms a portion of the boundary of the transformed region.
The generalization of these ideas to more complicated regions is obvious, the transformed region being composed of contiguous rectangular blocks. An example follows:
The physical boundary segment on which a single curvilinear coordinate is constant can have slope discontinuities, however, so that the Lshaped region above could have been considered to be composed of four segments instead of six, so that the transformed region becomes a simple rectangle:
Here the cartesian coordinates of the desired points on the physical boundary 543 are the specified boundary values from left to right across the top of the transformed region. Whether or not the boundary slope discontinuity propagates into the field, so that the coordinate lines in the field exhibit a slope discontinuity as well, depends on how the coordinate system in the field is generated, as will be discussed later.
It is not necessary that corners on the boundary of the transformed region correspond to boundary slope discontinuities on the physical boundary and a counterexample follows next:
In this case, the segment 12 on the physical boundary is a line of constant , while the segment 14 is a line of constant . Thus at point 1 we have the following coordinate line configuration:
The lines through point 1 are as follows:
so that the angle between the two coordinate lines is at point 1, and consequently the Jacobian of the transformation (the cell area, cf. Chapter III) will vanish at this point. The coordinate species thus changes on the physical boundary at point 1. (Difference representations at such special points as this, and others to appear in the following examples, are discussed in Chapter IV.) Since the species of curvilinear coordinate necessarily changes at a corner on the transformed region boundary, the identification of a concave corner on the transformed region boundary with a point on a smooth physical boundary will always result in a special point of the type illustrated here. (A point of slope discontinuity on the physical boundary also requires special treatment in difference solutions, since no normal can be defined thereon. This, however, is inherent in the nature of the physical boundary and is not related to the construction of the transformed configuration.)
Some slightly more complicated examples of the alternatives introduced above now follow:
Still another alternative in this case would be to collapse the intrusion 2345 to a slit in the transformed region:
Here the physical cartesian coordinates are specified and are doublevalued on the vertical slit, 295, in the transformed region. The cartesian coordinates of the desired points on the physical boundary 29 are to be used on the slit in the generation of the grid to the left of the slit in the transformed region, while those on the physical boundary 59 are used for generation to the right of the slit. Solution values in a numerical solution on such a coordinate system would also be doublevalued on the slit, of course. This doublevaluedness requires extra bookkeeping in the code, since two values of each of the cartesian coordinates and of the physical solution must be available at the same point in the transformed region so that difference representations to the left of the slit use the slit values appropriate to the left side, etc. Difference representations near slits are discussed in Chapter IV. With the composite grid structure discussed in Section 4, however, this need for doublevaluedness, and the concomitant coding complexity, with the slit configuration can be avoided.
The point 9 here requires special treatment, since the coordinate line configuration there is as follows:
The coordinate lines through point 9 are as follows:
Here the slope of the coordinate line on which varies is discontinuous at point 9, and the line on which varies splits at this point. Such a special point will always occur at the slit ends with the slit configuration.
B. Multiplyconnected Regions
With obstacles in the interior of the field, i.e., with interior boundaries, there are still more alternative configurations of the transformed region. One possibility is to maintain the connectivity of the transformed region the same as that of the physical region, as in the following examples showing two variations of this approach using interior slabs and slits, respectively, in the transformed region. The slab configuration is as follows:
In coding, points inside the slab in the transformed region are simply skipped in all computations.
This configuration introduces a special point of the following form at each of the points corresponding to the slab corners in the transformed field:
The coordinate lines through 
point 7 are shown below: 

This type of special point, where the coordinate species changes on a smooth line, occurs when a convex corner in the transformed field is identified with a point on a smooth contour in the physical field. Both coordinate lines experience slope discontinuities at this point.
The slit configuration is as shown below:
(An obvious varition would be to have the slit vertical.) In this slit configuration, the point 5 and 6 are special points of the form shown on p. 26 characteristic of the slit configuration, and will require special treatment in difference solutions.
The transformed region could, however, be made simplyconnected by introducing a branch cut in the physical region as illustrated below:
Conceptually this can be viewed as an opening of the field at the out and then a deformation into a rectangle:
Here the coincident coordinate lines 12 and 43 form a branch cut, which becomes reentrant boundaries on the left and right sides of the transformed region. All derivatives are continuous across this cut, and points at a horizontal distance outside the rightside boundary in the transformed region are the same as corresponding points at the same horizontal distance on the same horizontal line inside the leftside boundary, and vice versa. (In all discussions of point correspondence across cuts, "distance" means distance in the transformed region). In coding, the use of a layer of points outside each member of a pair of reentrant boundaries in the transformed region holding values corresponding to the appropriate points inside the other boundary of the pair avoids the need for conditional choices in difference representations, as discussed in Section 6 of this chapter.
Boundary values are not specified on the cut. (This cut is, of course, analogous to the coincident 0 and 2 lines in the cylindrical coordinate system discussed above.) At the cut we have the following coordinate line configuration, as may be seen from the conceptional deformation to a rectangle:
so that the coordinate species and directions are both continuous across the cut.
This type of configuration is often called an Otype. Another possible configuration is as shown below, often called a Ctype:
Opening the field at the cut we have, conceptually,
with 1234 to flatten to the bottom of the rectangle. Here the two members of the pair of segments forming the branch cut are both on the same side of the transformed region, and consequently points located at a vertical distance below the segment 12, at a horizontal distance to the left of point 2, coincide with points at the same vertical distance above the segment 43, at the same horizontal distance to the right of point 3. The point 2(3) is a special point of the type shown on p. 26 for slit configurations.
The coordinate line configuration at the cut in this configuration is as follows:
where it is indicated that varies to the right on the upper side of the cut, but to the left on the lower side. The direction of variation of also reverses at the cut, so that although the species and slope of both lines are continuous across the cut, the direction of variation reverses there.
It is possible to pass onto a different sheet across a branch cut, and discontinuities in coordinate line species and/or direction occur only when passage is made onto a different sheet. It is also possible, however, to remain on the same (overlapping) sheet as the cut is crossed, in which case the species and direction are continuous, and this must be the interpretation when derivatives are evaluated across the cut, as is discussed in Section 5 to follow. These concepts are illustrated in the following figure, corresponding to the Ctype configuration given on p. 30:
In the present discussion of configurations, the behavior of the coordinate lines across the cut will always be described in regard to the passage onto a different sheet, since this is in fact the case in codes. It is to be understood that complete continuity can always be maintained by conceptually remaining on the same sheet as the cut is crossed. Much of this complexity can, however, be avoided with the use of an extra layer of points surrounding the transformed region as will be discussed in Section 6.
Although in principle any region can be transformed into an empty rectangular block through the use of branch cuts, the resulting grid point distribution may not necessarily be reasonable in all of the region. Furthermore, an unreasonable amount of effort may be required to properly segment the boundary surfaces and to devise an appropriate point distribution thereon for such a transformation. Some configurations are better treated with a computational field that has slits or rectangular slabs in it.
Regions of higher connectivity than those shown above are treated in a similar manner. The level of connectivity may be maintained as in the following illustration:
Here one slit is made horizontal and one vertical just for generality of illustration. Both could, of course, be of the same orientation. Slabs, rather than slits, could also have been used. The example has three bodies.
With the transformed region made simplyconnected we have, using two branch cuts, a configuration related to the 0type shown above for one internal boundary:
The conceptual opening here is as follows:
with segment 234567 opening to the bottom. Here the pairs of segments (12,87) and (34,65) are the branch cuts, which form reentrant boundaries in the transformed region as shown. In this case, points outside the right side of the transformed region coincide with points inside the left side, and vice versa. This cut is of the form described on p. 30, where both the coordinate species and direction are continuous across the cut. Points below the bottom segment 34, to the left of point 4, coincide with points above the bottom segment 65 to the right of point 5. This cut is of the form discussed on p. 31, for which the coordinate species is continuous across the cut but the direction changes there. There are a number of other possibilities for placement of the two cuts on the boundary of the transformed region, of course, some examples of which follow.
It is not necessary to reduce the connectivity of the region completely; rather, a slit or slab can be used for some of the interior boundaries, while others are placed on the exterior boundary of the transformed region.
One other possibility in two dimensions is the use of a preliminary analytical transformation of infinity to a point inside some interior boundary, with the coodinates resulting therefrom replacing the cartesian coordinates in the physical region. The grid generation then operates from these transformed coordinates rather than from the cartesian coordinates. This typically gives a fine grid near the bodies, but may give excessively large spacing away from the body.
Thus, for example, if points on the two physical boundaries shown below
are transformed according to the complex transformation
z' = 1/z
where z = x+iy and z' = x'+iy', infinity in the x,y system will transform to the origin in the x', y' system, as shown below.
Then with the grid generated numerically from the x', y' system the following configuration results:
References to the use of this approach are made in the survey of Ref. [1]. Somewhat related to this are various twodimensional configurations which arise directly from conformal mapping, cf. Ref. [6] and the survey of Ives on this subject, Ref. [7]. (Conformal mapping is discussed in Chapter X.)
C. Embedded Regions
In more complicated configurations, one type of coordinate system can be embedded in another. A simple example of this is shown below, where an 0type system surrounding an internal boundary is embedded in a system of a more rectangular form, using what amounts to a slit configuration.
The conceptual opening of this system is best understood in stages: First considering only the embedded 0type system surrounding the interior boundary, we have the region inside the contour 121369 opening as follows:
This then opens to the rectangular central portion of the transformed region shown above, with the inner boundary contour 878 collapsing to a slit. The rest of the physical region then opens as shown below:
These two regions then deform to rectangles and are fitted to the top and bottom of the rectangle corresponding to the inner system along the contours 1213 and 96 as shown.
Here points at a vertical distance below the segment 1112 are coincident with points at the same vertical distance below the segment 109 on the same vertical line, and vice versa, with similar correspondence for the pair of segments 1314 and 65. Points at a horizontal distance to the left of the segment B12, at a vertical distance above point 8, coincide with points at the same horizontal distance to the right of the segment 89, at the same vertical distance below point 8. Similar correspondence holds for the pair 713 and 76. Boundary values are specified on the slit 87.
The composite system shown on p. 40 can also be represented as a slit configuration in the transformed region:
with the inner system represented as 

and the lower side of the slit considered reentrant with the left half of the top boundary of the rectangle corresponding to the inner system, the upper side of the slit being reentrant with the right half of this top boundary of the inner region. Now the conceptual opening is as follows for the inner region:
Difference representations made above the slit thus would use points below the right half of the top of the inner region in the transformed region, etc. Similarly, representation made below the left half of the top of the inner region would use points below the slit. The slit is thus a "black hole" into which coordinate lines from the outer system disappear, to reappear as part of the inner system. The slit here, matched with the top of the inner system, is then clearly a branch cut, and passage through the slit onto the inner system is simply passage onto a different sheet.
Note that the embedded system has its own distinctive species and directions for the coordinate lines, entirely separate from the outer system. Thus for the inner region the directions are as follows:
while for the outer region they are as shown below:
Thus at a point on the upper interface, 1213, between the systems the lines are as follows:
while on the lower interface, 96 they are as follows:
Thus both coordinates reverse direction at the lower interface although the species is continuous, while both the species and directions are continuous across the upper interface. This again corresponds to passage onto a different sheet, for the interface between the inner and outer systems, i.e., the segments 1213 and 96, is actually a branch cut.
The points 9(12) and 6(13) here require special notice. For example, at point 9 the coordinate line configuration is as follows:
The lines through point 9 are as shown below
There are thus several changes in species and direction at this point. This type of special point embodies the form which always occurs with the slit configuration, shown on p. 26, and occurs here because the embedded region inside the contour 961312 is essentially contained inside a slit defined by the same set of numbers.
The above discussion refers to the slit configuration on p. 41. For the configuration on p. 40, the lines in the outer region are still as diagrammed on p. 43, but the lines in the inner region now are as follows:
The coordinate line species and direction given on p. 43 for the upper interface, 1213, thus applies here on the entire interface between the two regions.
An alternative treatment of the two special points is to place them inside cells as shown below:
This results in a sixsided cell surrounding each of these two points which requires special treatment as discussed in Chapter IV.
Embedded systems can also be constructed in the block configuration:
Here the top of the block, 78, in the outer system is reentrant with the corresponding segment, 78, on a portion of the top of the inner system. The left side of the block, 67, and the bottom of the block, 65, are similarly reentrant with single portions of the top of the inner system. Finally, the right side of the block, 5128, is reentrant with two portions, 512 and 128, of the top of the inner system. Points outside one of these segments in one system are thus located at corresponding positions inside the other segment of the reentrant pair in the other system. The slab sides, matched with the top of the inner system, are thus branch cuts between the inner and outer systems.
Here the coordinate lines proceed as follows for the outer system:
while those for the inner system are the same as before, as shown on p. 42. This means that on the left and right sides of the block, i.e., segments 67 and 58, the line directions are as follows:
and on the top and bottom, segments 78 and 65, the directions are as shown below:
There are thus changes in coordinate species and/or direction that are different on each side of the block.
The point 8 (and points 7,6 and 5) are special points of the following form:
The lines through the point are as shown below:
Here the special points occur in the field instead of on the boundary.
An example of a Ctype system embedded in another Ctype system is given next:
Here the conceptual opening is as follows: First, considering the system about the upper body, we have the following configuration:
which, with the body collapsed to a slit, opens to the rectangle in the center of the transformed region. Next consider the system about the other body:
This opens to a rectangle, with the body flattening to a portion of the bottom, which is fitted to the first rectangle along the segment 1113. Finally, the outermost portion opens as follows:
which opens to a rectangle which is fitted to the first one along the segment 1214.
Again the embedded region inside the contour 14121113 can be considered to lie inside a slit. This contour, which forms the interface between the inner and outer systems, is actually a branch cut between the two systems, across which there are discontinuities in coordinate species and directon in the same manner as was discussed above for the previous embedded system. Points below segment 1612 coincide with points below segment 1711 in this case. Points to the left of segment 1512, above point 15, are coincident with points to the right of segment 1511 below point 15. The slit here is formed of the segments 815 and 915. The coincident points 11 and 12 here must be taken as a point boundary in the physical region, i.e., fixed at a specified value. Several special points of the types discussed above are present here.
An alternative arrangement of the transformed region that corresponds to exactly the same coordinate system in the physical region is as follows:
Here points below segment 34, to the left of point 4, coincide with points above segment 65, to the right of point 5. When calculations are made on or above the segment 1214 on the larger block, points below this segment coincide with points below the corresponding segment on the smaller block. Similarly, when calculations are made on or below the segment 1311 on the larger block points above this segment coincide with points below the corresponding segment on the smaller block. Finally, points below the segment 78, to the left of point 8, on the smaller block are coincident with points above the segment 109, to the right of point 9.
This configuration displays explicitly the correspondence of the embedded region inside the contour 14121113 to a slit. Conceptually, coordinate lines from the main system disappear into the slit and emerge into the embedded system. These coordinate lines thus are continued from the main system onto another sheet representing the embedded system. This concept of embedded systems, with continuation onto another sheet through a slit adds considerable flexibility to the grid configurations and is of particular importance with multiple boundaries and in three dimensions. The composite structure discussed in Section 4 removes much of the coding complexity associated with systems of this type.
D. Other Configurations
Another arrangement of cuts, where the species of coordinate changes on a continuous line as the cut is crossed, is illustrated below. The transformed region in this case is composed of three blocks connected by the cuts.
Here points outside one section are coincident with corresponding points inside the adjacent section.
The coordinate line configuration on the interface on the right side of block A here is as follows:
This same type of configuration occurs, in different orientations, on each of the interfaces. These interfaces are branch cuts, so that passage onto the adjacent block amounts to passage onto another sheet in the same manner discussed above.
As a final configuration for consideration in two dimensions, the following example shows a case with fewer lines on one side of a slab than on the other. This does not necessitate the use of different increments of the curvilinear coordinates in the numerical expressions, because, as has been mentioned, these increments always cancel out anyway.
E. Threedimensional Regions
All the general concepts illustrated in these examples extend directly to three dimensions. Interior boundaries in the transformed region can become rectangular solids and plates, corresponding to the slabs and slits, respectively, illustrated above for two dimensions. Examples of threedimensional configurations can be found in the surveys given by Ref. [8] and [9].
It is also possible to use branch cuts, as illustrated above for two dimensions, to bring the interior boundaries in the physical region entirely to the exterior boundary of the transformed region:
Physical space
Computational space
The correspondence between the physical and transformed fields can, however, become much more complicated in three dimensions, and considerable ingenuity may be required to visualize this correspondence. For instance, the simple case of polar coordinates corresponds to a rectangular solid with two opposing sides having the radial coordinate constant thereon, and two reentrant sides on which the longitudinal coordinate is constant at 0 and 2, respectively (corresponding to the cut). The remaining two sides correspond to the north and south polar axes, so that an axis opens to cover an entire side. There is thus a line, i.e., the axis, in the physical region that corresponds to an entire side in the transformed region.
Threedimensional grids may be constructed in some cases by simply connecting corresponding points on twodimensional grids generated on stacks of planes or curved surfaces:
It should be noted, however, that this procedure provides no inherent smoothness in the third direction, except in cases where the stack is formed by an analytical transformation, such as rotation, translation or scaling, of the twodimensional systems. An example of such an analytical transformation of twodimensional systems is the construciton of a threedimensional grid for a curved pipe by rotating and translating (and scaling if the crosssectional area of the pipe varies) twodimensional grids generated for the pipe crosssection so as to place these transformed twodimensional grids normal to the pipe axis at successive locations along the axis:
Another example is the rotation of a twodimensional grid about an axis to produce an axisymmetric grid:
4. Composite Grids
All of the above concepts can be incorporated in a single framework, and the coding complexity can be greatly reduced, by considering the physical field to be segmented into subregions, bounded by four (six in 3D) generally curved sides, within each of which an individual coordinate system is generated. The overall coordinate system, covering the entire physical field, is then formed by joining the subsystems at the subregion boundaries. The degree of continuity with which this juncture is made is a design consideration in regard to the mode of application intended for the resulting grid.
This segmentation concept is illustrated in the figure below.
The locations of the interfaces between the subregions in the physical region are, of course, arbitrary since these interfaces are not actual boundaries. These interfaces might be fixed, i.e., the location completely specified just as in the case of actual boundaries, or might be left to be located by the grid generation procedure. Also the coordinate lines in adjacent subregions might be made to meet at the interface between with complete continuity:
with some lesser degree of continuity, e.g., continuous line slope only:
or with a discontinuity in slope:
or perhaps not to meet at all:
Naturally, progressively more special treatment at the interface will be required in numerical applications as more degrees of line continuity at the interface are lost. Procedures for generating segmented grids with various degrees of interface continuity are discussed later, and conservative interface conditions are given in Ref. [52], [53].
Now, with regard to placing these concepts in the framework of segmentation, the sides of an individual subregion (called a "block" hereafter) can be treated as boundaries on which the coordinate points, and/or the coordinate line intersection angles, are specified, just as is done for actual boundaries, or a side may be treated as one member of a pair of reentrant boundaries, i.e., one side of a branch cut in the physical region across which complete continuity is established. The other member of the pair may be another side (or portion thereof) of the same block or may be all (or part of) a side of an adjacent block in the physical field. Recall that it is not necessary for a coordinate to remain of the same species across a reentrant boundary, since the passage is onto a different sheet. This can introduce some coding complexity, but the treatment is straightforward, and in fact the coding can be greatly simplified by using an extra layer of points surrounding each block as is discussed in Section 6.
Some of the general concepts have been embodied in the twodimensional code discussed in Ref. [19] and in three recent threedimensional codes, Ref. [13] Ref. [14], and [51].
A. Simplyconnected Regions
The first Lshaped simplyconnected configuration on p. 21 can be interpreted as being composed of three blocks, with the sides of adjacent blocks forming pairs of reentrant boundaries:
or two blocks with a portion of a side of one block reentrant with an entire side of another block:
Here, and in the examples to follow, solid lines correspond to physical boundaries, while the dashed lines correspond to the interfaces between the blocks. The dashed arrows indicate the linkage between the interfaces. (Obviously, any single block can be broken into any number of blocks connected by reentrant boundaries across adjacent sides.) In contrast, the Lshaped configuration on p. 22 corresponds to the use of a single block. Similarly, the configuration on p. 24 can be formed with three blocks:
while the first configuration for the same boundary on p. 25 is formed with a single block.
The slit configuration on p. 25 can be formed of three blocks:
or two blocks with only a portion of the adjacent sides of two blocks forming a reentrant boundary:
B. Multiplyconnected Regions
The configuration with a single cut shown on p. 29 corresponds to the use of a single block with the left and right sides here being the members of a pair of reentrant boundaries:
The multiplyconnected slab configuration on p. 27 can be broken into four blocks:
Other decompositions should also be immediately conceivable. The slit configuration on p. 28 can be formed with two blocks, again with only portions of adjacent sides serving as reentrant boundaries:
or into four blocks, with entire sides as reentrant boundaries in all cases:
The doublebody region on p. 34 opens to a single block as shown there, with portions of sides as reentrant boundaries. A five block configuration would use only entire sides as reentrant boundaries, however:
There is no real advantage, however, to the fiveblock system here.
C. Embedded Regions
The segmentation concept is most useful in the construction of embedded coordinate systems. For instance, the system on p. 40 can be considered to be formed of three blocks as follows:
Here portions of adjacent sides of the two larger blocks are reentrant with each other, while each of the remaining portions of these sides is reentrant with half of one side of the smaller block. The left and right sides of the smaller block are reentrant with each other. This configuration could also have been constructed with eight blocks:
with only entire sides being involved in reentrant pairs as shown.
With embedded systems the coordinate species often changes as the reentrant boundary is crossed. These systems also show that the blocks need be physically adjacent only in the physical field, and it is in this sense that "adjacent" is always to be interpreted. The transformed (computational) field should always be viewed as only a bookkeeping structure. Various constructions are possible for the configurations on p. 48 and 50, and a two block structure was actually used on p. 50. A further example follows:
D. Threedimensional Regions
For general threedimensional configurations, it is usually very difficult to obtain a reasonable grid with the entire physical region transformed to a single rectangular block. A better approach in most cases is to segment the physical region into contiguous subregions, each bounded by six curved surfaces, with each subregion being transformed into a rectangular block. An individual grid is generated in each subregion:
These subregion grids are patched together to form the overall grads, as in the twodimensional cases discussed above. Examples of the use of this segmentation in three dimensions are found, in particular, in Ref. [11] and [12]. Others are noted in the survey given by Ref. [9].
As noted above, complete continuity can be achieved at the subregion interfaces by noting the correspondence of points exterior to one subregion with points interior to another. The necessary bookkeeping can be accomplished, and the coding complexity can be greatly reduced, by using an auxiliary layer of points just outside each of the six sides of the computational region, analogous to the procedure mentioned above for two dimensions. A correspondence is then established in the code between the auxiliary points and the appropriate points just inside other subregions. This approach has recently been incorporated in an internal region code, Ref. [13], and in two codes for general regions, Ref. [14] and [51]. This is discussed in more detail in Section 6.
General threedimensional regions can be built up using subregions as follows: First, point distributions are specified on the edges of a curved surface forming one boundary of a subregion:
and a twodimensional coordinate system is generated on the surface:
When this has been done for all surfaces bounding the subregion, the threedimensional system within the subregion is generated using the points on the surface grids as boundary conditions:
In three dimensions it is possible for a line, e.g., a polar axis, in the physical region to map to an entire side of the computational region as in the illistration below, where the axis corresponds to the entire left side of the block:
The system illustrated here could be one of several identical blocks joining together to form a complete system around the axis.
It is illustrated by an exercise that the occurence of a polar axis can be avoided, and this facilitates the construction of a block structure. Thus a surface grid, having eight "corners", analogous to the four "corners" on the circle in the 2D grid on p. 23, can be constructed on the surface of a sphere. This serves much better than a latitudelongitude type system for joining to adjacent regions. Similarly, the use of the four "corner" system, rather than a cylindrical system, in a circular pipe allows Tsections and bifurcations to be treated easily by a composite structure, c.f. Ref. [13].
Generally, grid configurations with polar axis should not be used in composite grid structures.
E. Overlaid Grids
Another approach to complicated configurations is to overlay coordinate systems of different types, or those generated for different subregions:
Here an appropriate grid is generated to fit each individual component of the configuration, such that each grid has several lines of overlap with an adjacent grid. Interpolation is then used in the region of overlap when solutions are done on the composite grid, with iteration among the various grids. This approach has the advantage of simplicity in the grid generation, in that the various subregion grids are only required to overlap, not to fit. However, there would appear to be problems if regions of strong gradients fall on the overlap regions. Also the interpolation may have to be constructed differently for different configurations, so that a general code may be hard to produce. Some applications of such overlaid grids are noted in Ref. [5].
5. Branch Cuts
As has been noted in the above discussion of transformed field configurations, it is possible for discontinuities in coordinate species and/or direction to occur at branch cuts, in the sense of passage onto another sheet. Continuity can be maintained, however, by conceptually remaining on the same overlapping sheet as the cut is crossed. All derivatives thus do exist at the cut, but careful attention to difference formulations is necessary to represent derivatives correctly across the cut. Although the correct representation can be accomplished directly by surrounding the computational region with an extra layer of points, as is discussed in Section 6, it is instructive to consider what is required of a correct representation further here.
A. Point Correspondence
Points on reentrant boundaries in the transformed region, i.e., on branch cuts in the physical region, are not special points in the sense used above. Points on reentrant boundaries, in fact, differ no more from the other field points than do the points on the 0 and 2 lines in a cylindrical coordinate system. Care must be taken, however, to identify the interior points coinciding with the extensions from such points beyond the field in the transformed space. This correspondence was noted above in each of the configurations shown above, being indicated by the dashed connecting lines joining the two members of a pair of reentrant boundaries. There are essentially four types of pairs of reentrant boundaries, as illustrated in the following discussion of derivative correspondence. In these illustrations one exterior point, and its corresponding interior point, are shown for each case. The converse of the correspondence should be evident in each configuration.
For the configurations involving a change in the coordinate species at the cut, not only must the coordinate directions be taken into account as the cut is crossed, but also the coordinate species may need to be interpreted differently from that established across the cut in order to remain on the same sheet as the cut is crossed. For example, points on an line belonging to section A in the figure on p. 52, but located outside the right side of this region, are coincident with points on a line of region B at a corresponding distance (in the transformed region) below the top of this region.
B. Derivative Correspondence
Care must be taken at branch cuts to represent derivatives correctly in relation to the particular side of the cut on which the derivative is to be used. The existence of branch cuts indicates that the transformed region is multisheeted, and computations must remain on the same sheet as the cut is crossed. Remaining on the same sheet means continuing the coordinate lines across the cut coincident with those of the adjacent region, but keeping the same interpretation of coordinate line species and directions as the cut is crossed, rather than adopting those of the adjacent region. As noted above, points outside a region across a cut the transformed space are coincident with points inside the region across the other member of the pair of reentrant boundary segments corresponding to the cut in the transformed space. The positive directions of the curvilinear coordinates to be used at these points inside the region across the other member of the pair in some cases are the same as the defined directions there, but in other cases are the opposite directions. As noted above, the coordinate species may change also.
For cuts located on opposing sides of the transformed region, the proper form is simply a continuation across the cut. Thus in the configuration on p. 29, with a computation site on the right side of the transformed region, i.e., on the upper side of the cut in the physical plane, we have points to the right of the site (below the cut in the physical plane) coinciding with points to the right of the left side of the transformed region (below the cut in the physical plane) as noted above. When derivatives and derivatives for use outside the right side of the transformed region are represented inside the left side, the positive directions of and to be used there are to the right and upward, respectively, as is illustrated below. (In this and the following figures of the section, the dotted arrows indicate the proper directions to be used at the interior points coincident with the required exterior points, i.e., on the same sheet across the cut, while solid arrows indicate the locally established directions for the coordinate lines, i.e., on a different sheet.)
With the two sides of the cut both located on the same coordinate line, i.e., on the same side of the transformed region as in the configuration on p. 30, however, the situation is not as simple as the above. In this case, when the computation site is on the left branch of the cut in the transformed region (on the lower branch in the physical region), the points below this boundary in the transformed region coincide with points located above the right branch of the cut (above the cut in the physical region) at mirrorimage positions, as has been noted earlier. The derivatives for use at such points below the left branch thus must be represented at these corresponding points above the right branch. The positive direction of for purposes of this calculation of derivatives above the right branch, for use below the left branch, must be taken as downward, not upward. There is a similar reversal in the interpretation of the positive direction of . This is in accordance with the discussion on p. 31. These interpretations are illustrated below:
In the configuration on p. 40, where two sides of a cut face each other across a void, there is really no problem of interpretation, since the directions in the configuration are treated simply as if the void did not exist. This correspondence is as shown below:
In all cases the interpretation of the positive directions of the curvilinear coordinates must be such as to preserve the direction in the physical region, i.e., on the same sheet, as the cut is crossed. In the cases where the coordinate species change at the cut, the situation is even more complicated. Thus on the left side, segment 67, of the slab interface between the inner and outer systems in the embedded configuration on p. 46, where the species changes across the cut, the correspondence is as follows:
Thus, when a derivative is needed outside the outer sytem, for use inside the left slab interface, the positive direction at the corresponding points inside the inner system must be taken to coincide with the negative direction of the inner system. Similarly, an derivative would be represented taking the positive direction to coincide with the positive direction of the inner system. In an analogous fashion, a derivative needed outside the inner system, for use inside the segment 67, would be represented at the corresponding point inside the outer system, i.e., to the left of the left slab side, but with the positive direction taken to be the positive direction of the outer system. An derivative would be represented similarly, taking the positive direction to be the negative direction of the outer system.
A derivative to the left of the right side of the slab in the outer system would be represented below segment 125 or 812, as the case may be, but with the positive direction taken to be the positive direction of the inner system. Similarly, an derivative would be represented taking the positive direction to be the negative direction of the inner system. For a derivative above the bottom of the slab in the outer system, the correspondence is to below the segment 56 inside the inner system, with the positive direction taken to be the negative direction of the inner system. The derivative is represented taking the positive direction to be the negative direction of the inner sytem. Finally, for derivatives below the top of the slab in the outer system, the correspondence is to below the segment 78 inside the inner system, with both the species and direction of the coordinates unchanged.
The proper interpretation of coordinate species and direction across branch cuts for all the other configurations discussed above can be inferred directly from these examples. A conceptual joining of the two members of a pair of reentrant boundaries in accordance with the dashed line notation used on the configurations given in this chapter will always show exactly how to interpret both the coordinate species and directions in order to remain on the same sheet and thus to maintain continuity in derivative representation across the cut. Examples of the proper difference representation are given in the following section. The complexities of this correspondence can be completely avoided, however, by using surrounding layers around each block in a segmented structure as discussed in the next section.
6. Implementation
As discussed above, the transformed region is always comprised of contiguous rectangular blocks by construction. This occurs because of the essential fact that one of the curvilinear coordinates is defined as constant on each segment of the physical boundary. Consequently, each segment of the physical boundary corresponds to a plane segment of the boundary of the transformed region that is parallel to a coordinate plane there. The complete boundary of the transformed region then is composed of plane segments, all intersecting at right angles. Although the transformed region may not be a simple sixsided rectangular solid, it can be broken up into a contiguous collection of such solids, here called blocks.
Now it is noted in Chapter III that the increments ^{i} cancel from all difference expressions, and that the actual values of the curvilinear coordinates ^{i} are immaterial. The coordinates in the transformed region can thus be considered simple counters identifying the points on the grid. This being the case, and the transformed region being comprised of a collection of rectangular blocks, it is convenient to identify the grid points with integer values of the curvilinear coordinates in each block, and thus to place the cartesian coordinates of a grid point in _{ijk}, where the subscripts (i,j,k) here indicate position (^{i},^{2}, ^{3}) in the transformed region. (In coding, a fourth index may be added to identify the block.) In each block, the curvilinear coordinates are then taken to vary as ^{i} = 1,2,...,I^{i} over the grid points, where I^{i} is the number of points in the ^{i}direction. Grid points on a boundary segment of the transformed region will be placed in _{ijk} with one index fixed.
Now each block has six exterior boundaries, and may also have any number of interior boundaries (cf. the slab and slit configurations of Section 3), all of which will always be plane segments intersecting at right angles, although the occur ence of interior boundaries can be avoided if desired by breaking the block up into a collection of smaller blocks as discussed in Section 4. The boundary segments in the transformed plane may correspond to actual segments of the physical boundary, or may correspond to cuts in the physical region. As discussed in Section 5, these cuts are not physical boundaries, but rather are interfaces across which the field is reentrant on itself. A boundary segment in the transformed region corresponding to such a cut then is an interface across which one block is connected with complete continuity to another block, or to another side of itself, several examples having been given above in this chapter.
Depending on the type of grid generation system used (cf., the later chapters), the cartesian coordinates of the grid points on a physical boundary segment may either be specified or may be free to move over the boundary in order to satisfy a condition, e.g., orthogonality, or the angle at which coordinate lines intersect the boundary.
To set up the configuration of the transformed region, a correspondence is established between each (exterior or interior) segment of the boundary of the transformed region and either a segment of the physical boundary or a segment of a cut in the physical region. This is best illustrated by a series of examples using the configurations of this chapter. The first step in general is to position points on the physical boundary, or on a cut, which are to correspond to corners of the transformed region (exterior or interior). As noted in Section 3, these points do not have to be located at actual corners (slope discontinuities) on the physical boundary.
For example, considering the twodimensional simplyconnected region on p. 23, four points on the physical boundary are selected to correspond to the four corners of the empty rectangle that forms the transformed region here:
Now, considering any one of these four points, one species of curvilinear coordinate will run from that point to one of the two neighboring corner points, while the other species will run to the other neighbor:
The corresponding species of coordinates will run to connect opposite pairs of corner points:
Since the curvilinear coordinates are to be assigned integer values at the grid points, ^{i} is to vary from 1 at one corner to a maximum value, I^{i}, at the next corner, where I^{i} is the number of grid points on the boundary segment between these two corners. Thus, proceeding clockwise from the lower left corner, the cartesian coordinates of the four corner points are placed in _{1,1}, _{1,J}, _{I,J}, and _{I,1}, where I^{1} = I and I^{2} = J.
The boundary specification is then completed by positioning I2 points on the lower and upper boundary segments of the physical region as desired, and J2 points on the left and right segments. The cartesian coordinates of these points on the lower and upper segments are placed in _{i,1} and _{i,J}, respectively, for i from 2 to I1, and those on the left and right segments are placed, respectively, in _{1,j} and _{I,j} for j from 2 to J1.
This process of boundary specification can be most easily understood by viewing the rectangular boundary of the transformed region, with I equallyspaced points along two opposite sides and J equallyspaced points along the other two sides, conceptually, as being deformed to fit on the physical boundary. The corners can be located anywhere on the physical boundary, of course. Here the point distribution on the sides can be conceptually stretched and compressed to position points as desired along the physical boundary. The cartesian coordinates of all the selected point locations on the physical boundary are then placed in as described above.
This conceptual deformation of the rectangular boundary of the transformed region to fit on the physical boundary serves to quickly illustrate the boundary specification for the doublyconnected physical field shown on p. 29, which involes a cut. Thus I points are positioned as desired clockwise around the inner boundary of the physical region from 2 to 3, and I points are positioned as desired, also in clockwise progression, around the outer boundary from 1 to 4. The cartesian coordinates of these points on the inner boundary are placed in _{i,1}, and those on the outer boundary in _{i,J}, with i from 1 to I. Note that here the first and last points must coincide on each boundary, i.e., _{I,1} = _{1,1} and _{I,J} = _{1,J}. The left and right sides of the transformed region (i=1 and i= I) are reentrant boundaries, corresponding to the cut, and hence values on these boundaries are not set but will be determined by the generation system. The system must provide that the same value appears on both of these sides, i.e., _{I,j} = _{1,j} for all j from 2 to J1.
The conceptual deformation of the rectangle for a Ctype configuration is illustrated on p. 31. Here, with I1 the number of points on the segments 12 and 34 (which must have the same number of points), I2 points are positioned as desired around the inner boundary in the physical region in a clockwise sense from 2 to 3, and the cartesian coordinates of these points are placed in _{i,1} for i from I1 to I1+I21.
The first and last of these points must be coincident, i.e. _{I1,1} = _{I1 + I21,1}. Now the top, and the left and right sides, of the rectangle are deformed here to fit on the outer boundary of the physical region. (In the illustration given, the two top corners are placed on the two corners that occur in the physical boundary, a selection that is logical but not mandatory.) The cartesian coordinates of the J points(positioned as desired on the segment 45 of the physical boundary) are placed in _{I,j}, proceeding upward on the physical boundary from 4 to 5 for j= 1 to J, and those on the segment 16 are placed in _{1,j}, but proceeding downward on the physical boundary from 1 to 6 for j1 to J. Finally, the cartesian coordinates of the I selected points on the physical boundary segment 65 are placed in _{i,J}, proceeding clockwise from 6 to 5 for i=1 to I. Since the same number of points must occur on the top and bottom of the rectangle, we must have I=2(I11)+I2. Here the portions of the lower side of the rectangle, i.e., i from 2 to I11, and from I1+I2 to I1 with j=1, are reentrant boundaries corresponding to the cut, and hence no values are to be specified on these segments. The generation system must make the correspondence _{i,1} = _{Ii+1,1} for i=2 to I11 on these segments.
The conceptual deformation of the boundary of the transformed regions also serves for the slab configuration on p. 27, where the interior rectangle deforms to fit the interior physical boundary, while the outer rectangle deforms to fit the outer physical boundary. On the inner boundary, the cartesian coordinates of J2J1+1 selected points on the segment 58 of the physical boundary are placed in _{I1,j} for j from J1 to J2, proceeding upward on the physical boundary from 5 to 8, where J1 and J2 are the jindices of the lower and upper sides, respectively, of the interior rectangle and I1 is the iindex of the left side of this rectangle. Similarly, J2J1+1 points are positioned as desired on the segment 67 of the physical boundary and are placed in _{I2,j}, where I2 is the iindex of the right side of the inner rectangle. Also I2I1+1 points on the segments 56 and 87 of the physical boundary are placed in _{i,J1} and _{i,J2}, respectively, for i from I1 to I2, proceeding to the right on each segment. The outer boundary is treated as has been described for an empty rectangle. Here there will be J11 coordinate lines running from left to right below the inner boundary, and JJ2 lines running above the inner boundary. Similarly, there will be I11 lines running upward to the left of the interior boundary and II2 lines to the right. Thus the specifications of the desired number of coordinate lines running on each side of the inner boundary serves to determine the indicies I1, I2, J1, and J2. Note that the points inside the slab, i.e., I1 < i < I2 and J1 < j < J2 are simply excluded from the calculation.
The slit configuration, illustrated on p. 28, can also be treated via the conceptual deformation, but now with a portion of a line inside the rectangle opening to fit the interior boundary of the physical region. This requires that provision be made in coding for two values of the cartesian coordinates to be stored on the slit. If the iindices of the slit ends, 5 and 6, are I1 and I2, respectively, then the cartesian coordinates of I2I1+1 points positioned as desired on the lower portion of the physical interior boundary, again proceeding from 5 to 6, are placed in a onedimensional array, while the coordinates of the same number of points selected on the upper portion of the physical interior boundary, again proceeding from 5 to 6, are placed in another onedimensional array. The first and last points in one of these arrays must, of course, coincide with those in the other. Then the generation system must read values into _{i,J1} for i from I1 to I2 (J1 being the jindex of the slit) from the former array for use below the slit, or values from the latter array for use above. (As has been noted, the use of a composite structure eliminates the need for these two auxiliary arrays.) Note that the index values I1 and I2 are determined by the number of lines desired to run upward to the left and right of the interior boundary, respectively, i.e., I11 lines on the left and II2 on the right. Similarly, there will be J11 lines below the interior boundary, and JJ1 above.
Configurations, such as those illustrated on pp. 2425, which involve slabs or slits that intersect the outer boundary are treated similarly, with points inside the slab again being simply excluded from the calculations. Also multiple slab or slit arrangements are treated by obvious extensions of the above procedures. Here the indices corresponding to each slab or slit will be determined by the number of points on the interior boundary segments and the number of coordinate lines specified to run between the various boundaries. For example, in the slit configuration shown on p. 33, the ends of the horizontal slit would be at iindices I1 and I2, where I11 lines run vertically to the left of the slit and there are I2I1+1 points on the slit. The vertical slit would be at i=I3 where there are I3I21 vertical lines between this slit and the horizontal slit (and II2 lines to the right). Similarly, if the jindices of the ends of the vertical slit or J1 and J2, there will be J11 horizontal lines below this slit and JJ2 lines above. With the jindex of the horizontal slit as J3, there will be J31 horizontal lines below this slit and JJ3 above. Provision will now have to be made in coding for two onedimensional arrays for each slit to hold the cartesian coordinates of the points on the segments of the physical interior boundaries corresponding to the two sides of each slit. Again this coding complexity is avoided in the composite structure.
The use of the conceptual deformation of the rectangle to setup the boundary configuration for the case with multiple interior boundaries on p. 34 should follow with little further explanation. Here there must be the same number of points on the pair of segments 23 and 67, which correspond to the two segments forming the interior boundary on the right. There must also be the same number of points on the pair, 34 and 56, corresponding to the cut connecting the two interior boundaries. Finally the number of points on the outer boundary must, of course, be the same as that on the bottom boundary. Note also that the values of the cartesian coordinates placed at 2 must be the same as are placed at 7; those at 3 must be the same as those at 6, and those at 4 the same as at 5. Values are not set on the cuts, of course, but the generation system must provide that values at points on the segment from 3 to 4 are the same as those on the segment 56, but proceeding from 6 to 5. Also values on the segments 21 and 78 must be the same, proceeding upward in each case.
Following the conceptual deformation of the rectangular boundaries of the transformed region and the indexing system illustrated above, it now should be possible to set up the more complicated configurations such as the embedded regions shown in Section 3C. As noted there, however, the most straightforward and general approach to such more complicated configurations is to divide the field into contiguous rectangular blocks, each of which has its own intrinsic set of curvilinear coordinates and hence its own (i,j,k) indexing system. The necessary correspondence between the individual coordinate systems across the block interfaces was discussed in some detail in Section 3C. This block structure greatly simplifies the setup of the configuration. For example, consider the 3block structure shown on p. 49 for the physical field shown on p. 48, for which the blocks are as follows:
Here the selected points on the right interior boundary (segment 8159) are placed in _{i,1} of the first block, for i from the iindex at 8 to that at 9,proceeding clockwise from 8 to 9 on the physical boundary. (The difference between these two iindices here is equal to the number of points on this interior boundary,less one.) Similary,the selected points on the left interior boundary (segment 45) are placed in _{i,1} of the second block for i from the iindex at 4 to that at 5, proceeding clockwise from 4 to 5 on the physical boundary. The selected points on the outer boundary of the physical region are placed in _{1,j} of the third block for j from 1 to J3, in _{i,J3} for i from 1 to I3, and in _{I3,j} for j from J3 to 1, proceeding from 16 to 1 to 2 to 14 on the physical boundary. Points on the remainder of the physical outer boundary are placed in _{1,j} of the second block for j from 1 to J2 and in _{I2,j} for j from J2 to 1, proceeding from 3 to 17 for the former and from 13 to 6 for the latter, and in _{1,j} of the first block for j from 1 to Jl and in _{I1,j} for j from Jl to 1, proceeding from 7 to 13 for the former and from 14 to 10 for the latter.
Since the three blocks must fit together we have I3=I2, (I1+1)/2 equal to the difference in iindices between 11 and 13 in the second block and to that between 12 and 14 of the third block. The quantities J1, J2, and J3 determine how many Ctype lines occur in each block, and can be chosen independently. Here the segment 1113 on the top of the first block interfaces with the corresponding segment on the top of the second block. The segment 1214, which forms the remainder of the top of the first block, interfaces with the corresponding segment on the bottom of the third block. Finally, the segment 1216, which forms the remainder of the bottom of the third block, interfaces with segment 1117, which forms the remainder of the top of the second block. The segments 34 and 65 on the bottom of the second block interface with each other in the order indicated, as do also the segments 78 and 109 on the bottom of the first block.
In coding, this block structure can be handled by using a fourth index to identify the block, placing an extra layer around each block, (i=0 and I+1, j=0 and J+1) and providing an imagepoint array by which any point of any block can be paired with any point of any other, or the same, block. Such pairs of points are coincident in the physical region, being on or across block interfaces, and consequently are to be given the same values of the cartesian coordinates by the generation system. This imaging extends to the extra layer surrounding each block, so that appropriate points Inside other blocks can be identified for use in difference representations on the block interfaces that require points outside the block, (cf. Section 5).
Interface correspondence then can be established by input by setting the imagepoint correspondence on the appropriate block sides, i.e., placing the (i,j,k) indices and block number of one member of a coincident pair of points in the imagepoint array at the indices and block number of the other member of the pair. This correspondence is indicated on the block diagram on pp. 8586 by the points enclosed in certain geometric symbols.
Thus, for the 3block configuration considered above, the indices (I1i+1, 1) and block number 1, corresponding to a point on the segment 910 of the first block, would be placed in the imagepoint array at the point (i,1) on the segment 78 of this block, and vice versa. A similar pairing occurs for points on the segments 34 and 56 of the second block. The indices (I2i+1, J2) and block number 2, corresponding to a point on the segment 1113 of the second block would be placed in the imagepoint array at the point (i,Jl) of the first block on the segment 1311 of that block, and vice versa. The indices (I3I1+i, 1) and block number 3 (a point on the segment 1214 of the third block) would be placed in the array at the point (i,J1) of the first block for a point on the segment 1214 of that block. Finally, the indices (i,1) and block number 3 (point on segment 1612 of the third block) would be placed in the array at the point (i,J2) of the second block for a point on the segment 1711 of that block. The remaining segments all correspond at portions of the physical boundary and hence do not have image points.
In the same manner the following image correspondence can be set between interior points and points on the surrounding layers in order to establish difference representations across the block interfaces: (This correspondence is indicated symbolically on the block diagram on pp. 8586 by geometric symbols.):
As noted, all of this information would be input into the imagepoint array. Then with values of the cartesian coordinates at the image points on the surrounding layer set equal to those at the corresponding object point inside one of the blocks, it is possible to use the same difference representations on the interfaces that are used in the interior.
The discussion given in this chapter should now allow the imagepoint input to be constructed for any configuration of interest. As noted, it is not necessary that the coordinate species remain the same as the interface is crossed. Thus, for instance, a point on the right side of one block could be paired with one on the bottom of another block. In such a case the image point of the point (I+1, j) the first block would be the point (j,2) inside the second block. Similarly the image of the point (i,0) below the second block would be the point (I1, i) inside the first block, The correct difference representation across interfaces is thus automatically established, eliminating the need for the concern with passage onto different sheets discussed in detail earlier in this chapter.
This greatly simplifies the coding, since with the surrounding layers and the use of the image points, all of the derivative correspondences are automatic and do not have to be specified for each configuration. It is only necessary to specify the point correspondence by input. This construction also allows codes for the numerical solution of partial differential equations on the grid to be written to operate on rectangular blocks. Then any configuration can be treated by sweeping over all the blocks. The surrounding layers of points and the image correspondence provide the proper linkage across the block interfaces. In an implicit solution the values on the interfaces would have to be updated iteratively in the course of the solution. The solution for the generation of the grid would similarly keep the interface and surrounding layer values updated during the course of the iterative solution.
This, of course, maintains completecontinuity across the block interfaces. If complete continuity is not required, then the surrounding layer is not required and the interfaces would be treated in the same manner as are physical boundaries. However, the surrounding layer and the point correspondence thereon discussed above might still be needed for the numerical solution to be done on the grid.
The extension of all of the above concepts and structures to three dimensions is direct, the illustrations having been given in two dimensions only for economy of presentation.
Exercises
1. Sketch the grid when the physical region shown below is transformed to an empty rectangle as indicated.
2. Locate the cuts on the grid in the physical region for the grids shown on pp. 3537.
3. Sketch the grid for a Ogrid, a Cgrid, and a slit configuration for a cascade arrangement (a periodic stack of bodies, e.g., turbine blades.)
4. For the configuration shown below, let body II transfrom form to a slit in a Ctype system about I. Sketch the grid lines and show the transformed region configuration.
5. Sketch the surface grid on a sphere with eight "corners" on the sphere. This is analogous to the 2D grid for the circle with four "corners" on p. 23. This configuration would be more appropriate for embedding a spherical region in a composite structure than would the usual polar system.
6. Diagram the transformed region configuration for a polar coordinate system between two concentric spheres. Here the polar axis will map to an entire side of the transformed region.
7. Sketch the surface grid on a circular cylinder with a hemispherical cap for two cases: (1) with a cylindricalcoordinate type grid on the hemisphere and (2) with four "corners" on the intersection of the hemisphere with the cylinder. The latter case would be more appropriate for a composite system.
8. For the twoslit configuration on p. 33, diagram a composite system composed of empty blocks. Show all cuts and the correspondence between pairs thereof
9. Sketch the grid and diagram the blocks for a composite twoblock system for a circular pipe Tsection. It is necessary to use a crosssectional grid of the type shown on p. 23, having four "corners" on the circle, since cylindrical grids would not join with line continuity.
10. Diagram the block structure and grid for a sixblock composite system for the region between two concentric spheres, based on the surface grid of Exercise 10. Note that no polar axis occurs with this configuration.
11. Consider a region between two boundaries, both of which are formed of cylinders with hemispherical caps, these being coaxial with one inside the other. Sketch the grid and diagram the blocks for a threeblock system, with one block corresponding to the annular region between the caps, for the following two configurations: (1) with the polar axis connecting the caps and (2) with no polar axis. In the latter case each of four sides of one block will correspond to one of four portions of one side of the other block.
12. Diagram the point correspondence across all the cuts in the twobody 0grid on p.34. Also give the relation between the indices of corresponding points on the cuts. Finally give the relation between the indices of points on a surrounding layer of points and points inside the field inside the cuts.
13. For one block of the system on p. 52, give the correspondence between indices of points on the surrounding layers and points inside adjacent blocks for the cuts.
14. Sketch the grid for a 2D composite system having two circular regions embedded in a grid which is generally rectangular. Let one of the circular regions have a cylindricalcoordinate type of grid and the other have a grid of the type with four "corners" on the circle as on p. 23.
15. Show that is is not possible to handle the point correspondence across the cuts in the embedded slab type system shown on p. 46 (a 2block system) by using an extra layer of points just inside the slab in the outer system. Also show that it is possible to represent the correspondence across the cuts using surrounding layers if a 4block composite system is used.