HPC MSU
    Cover
    Preface
  1. Introduction
  2. Boundary-Conforming Coordinate Systems
    1. Basic Concepts
    2. Generalization
      1. Boundary-value Problem -- Physical Region
      2. Boundary-value Problem -- Transformed Region
    3. Transformed Region Configurations
      1. Simply-connected Regions
      2. Multiply-connected Regions
      3. Embedded Regions
      4. Other Configurations
      5. Three-dimensional Regions
    4. Composite Grids
      1. Simply-connected Regions
      2. Multiply-connected Regions
      3. Embedded Regions
      4. Three-dimensional Regions
      5. Overlaid Grids
    5. Branch Cuts
      1. Point Correspondence
      2. Derivative Correspondence
    6. Implementation
    Exercises
  3. Transformation Relations
  4. Numerical Implementation
  5. Truncation Error
  6. Elliptic Generation Systems
  7. Parabolic and Hyperbolic Generation Systems
  8. Algebraic Generation Systems
  9. Orthogonal Systems
  10. Conformal Mapping
  11. Adaptive Grids
  12. Appendix A
    Appendix B
    Appendix C
    References

    Downloadable Version (PDF)
Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin

II. BOUNDARY-CONFORMING COORDINATE SYSTEMS

1. Basic Concepts

          To provide a familiar ground from which to view the general development to follow, consider first a two-dimensional cylindrical coordinate system covering the annular region between two concentric circles:

Here the curvilinear coordinates (r,) vary on the intervals [r1,r2] 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 = r1, and r = r2. 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 re-entrant 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 re-entrant 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 re-entrant boundaries will serve later for more general configurations.

          These simple concepts extend to more complicated two-dimensional 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 two-dimensional 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 one-dimensional 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 one-to-one. 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 x2 + y2 = and x2 + y2 = , subject to the boundary conditions

          Specified monotonic variation of over [0,1] on x2 + y2 = and on x2 + y2 = 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 x2( = 0) + y2,0) =
x(,1) and y(,1) specified on = 1 such that x2( = 1) + y2(,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 boundary-conforming 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 one-to-one, 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 two-dimensional 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. Boundary-value 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 re-entrant 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 re-entrant 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 re-entrant 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. Simply-connected 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 simply-connected region formed by four curves is logically treated by transforming to an empty rectangle:

Similarly, an L-shaped region could remain L-shaped in the transformed region:

Here, for instance, the cartesian coordinates of the desired points on the physical boundary segment 4-5 are specified as boundary conditions on the vertical line 4-5, 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 L-shaped 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 5-4-3 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 counter-example follows next:

In this case, the segment 1-2 on the physical boundary is a line of constant , while the segment 1-4 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 2-3-4-5 to a slit in the transformed region:

Here the physical cartesian coordinates are specified and are double-valued on the vertical slit, 2-9-5, in the transformed region. The cartesian coordinates of the desired points on the physical boundary 2-9 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 5-9 are used for generation to the right of the slit. Solution values in a numerical solution on such a coordinate system would also be double-valued on the slit, of course. This double-valuedness 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 double-valuedness, 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. Multiply-connected 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 simply-connected 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 1-2 and 4-3 form a branch cut, which becomes re-entrant 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 right-side boundary in the transformed region are the same as corresponding points at the same horizontal distance on the same horizontal line inside the left-side 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 re-entrant 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 O-type. Another possible configuration is as shown below, often called a C-type:

Opening the field at the cut we have, conceptually,

with 1-2-3-4 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 1-2, at a horizontal distance to the left of point 2, coincide with points at the same vertical distance above the segment 4-3, 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 C-type 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 simply-connected we have, using two branch cuts, a configuration related to the 0-type shown above for one internal boundary:

The conceptual opening here is as follows:

with segment 2-3-4-5-6-7 opening to the bottom. Here the pairs of segments (1-2,8-7) and (3-4,6-5) are the branch cuts, which form re-entrant 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 3-4, to the left of point 4, coincide with points above the bottom segment 6-5 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 two-dimensional 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 0-type 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 0-type system surrounding the interior boundary, we have the region inside the contour 12-13-6-9 opening as follows:

This then opens to the rectangular central portion of the transformed region shown above, with the inner boundary contour 8-7-8 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 12-13 and 9-6 as shown.

          Here points at a vertical distance below the segment 11-12 are coincident with points at the same vertical distance below the segment 10-9 on the same vertical line, and vice versa, with similar correspondence for the pair of segments 13-14 and 6-5. Points at a horizontal distance to the left of the segment B-12, at a vertical distance above point 8, coincide with points at the same horizontal distance to the right of the segment 8-9, at the same vertical distance below point 8. Similar correspondence holds for the pair 7-13 and 7-6. Boundary values are specified on the slit 8-7.

          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 re-entrant with the left half of the top boundary of the rectangle corresponding to the inner system, the upper side of the slit being re-entrant 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, 12-13, between the systems the lines are as follows:

while on the lower interface, 9-6 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 12-13 and 9-6, 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 9-6-13-12 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, 12-13, 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 six-sided 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, 7-8, in the outer system is re-entrant with the corresponding segment, 7-8, on a portion of the top of the inner system. The left side of the block, 6-7, and the bottom of the block, 6-5, are similarly re-entrant with single portions of the top of the inner system. Finally, the right side of the block, 5-12-8, is re-entrant with two portions, 5-12 and 12-8, 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 re-entrant 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 6-7 and 5-8, the line directions are as follows:

and on the top and bottom, segments 7-8 and 6-5, 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 C-type system embedded in another C-type 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 11-13. Finally, the outermost portion opens as follows:

which opens to a rectangle which is fitted to the first one along the segment 12-14.

          Again the embedded region inside the contour 14-12-11-13 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 16-12 coincide with points below segment 17-11 in this case. Points to the left of segment 15-12, above point 15, are coincident with points to the right of segment 15-11 below point 15. The slit here is formed of the segments 8-15 and 9-15. 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 3-4, to the left of point 4, coincide with points above segment 6-5, to the right of point 5. When calculations are made on or above the segment 12-14 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 13-11 on the larger block points above this segment coincide with points below the corresponding segment on the smaller block. Finally, points below the segment 7-8, to the left of point 8, on the smaller block are coincident with points above the segment 10-9, to the right of point 9.

          This configuration displays explicitly the correspondence of the embedded region inside the contour 14-12-11-13 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. Three-dimensional 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 three-dimensional 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 re-entrant 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.

          Three-dimensional grids may be constructed in some cases by simply connecting corresponding points on two-dimensional 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 two-dimensional systems. An example of such an analytical transformation of two-dimensional systems is the construciton of a three-dimensional grid for a curved pipe by rotating and translating (and scaling if the cross-sectional area of the pipe varies) two-dimensional grids generated for the pipe cross-section so as to place these transformed two-dimensional grids normal to the pipe axis at successive locations along the axis:

Another example is the rotation of a two-dimensional grid about an axis to produce an axi-symmetric 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 sub-regions, 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 sub-systems at the sub-region 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 sub-regions 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 sub-regions might be made to meet at the interface between with complete continuity:

---Interface

with some lesser degree of continuity, e.g., continuous line slope only:

---Interface

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 re-entrant 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 re-entrant 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 two-dimensional code discussed in Ref. [19] and in three recent three-dimensional codes, Ref. [13] Ref. [14], and [51].

A. Simply-connected Regions

          The first L-shaped simply-connected configuration on p. 21 can be interpreted as being composed of three blocks, with the sides of adjacent blocks forming pairs of re-entrant boundaries:

or two blocks with a portion of a side of one block re-entrant 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 re-entrant boundaries across adjacent sides.) In contrast, the L-shaped 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 re-entrant boundary:

B. Multiply-connected 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 re-entrant boundaries:

          The multiply-connected 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 re-entrant boundaries:

or into four blocks, with entire sides as re-entrant boundaries in all cases:

          The double-body region on p. 34 opens to a single block as shown there, with portions of sides as re-entrant boundaries. A five block configuration would use only entire sides as re-entrant boundaries, however:

There is no real advantage, however, to the five-block 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 re-entrant with each other, while each of the remaining portions of these sides is re-entrant with half of one side of the smaller block. The left and right sides of the smaller block are re-entrant with each other. This configuration could also have been constructed with eight blocks:

with only entire sides being involved in re-entrant pairs as shown.

          With embedded systems the coordinate species often changes as the re-entrant 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. Three-dimensional Regions

          For general three-dimensional 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 sub-regions, each bounded by six curved surfaces, with each sub-region being transformed into a rectangular block. An individual grid is generated in each sub-region:

These sub-region grids are patched together to form the overall grads, as in the two-dimensional 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 sub-region interfaces by noting the correspondence of points exterior to one sub-region 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 sub-regions. 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 three-dimensional regions can be built up using sub-regions as follows: First, point distributions are specified on the edges of a curved surface forming one boundary of a sub-region:

and a two-dimensional coordinate system is generated on the surface:

When this has been done for all surfaces bounding the sub-region, the three-dimensional system within the sub-region 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 latitude-longitude 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 T-sections 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 sub-regions:

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 sub-region 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 re-entrant 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 re-entrant 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 re-entrant boundaries. There are essentially four types of pairs of re-entrant 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 multi-sheeted, 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 re-entrant 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 mirror-image 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 6-7, 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 6-7, 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 12-5 or 8-12, 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 5-6 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 7-8 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 re-entrant 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 six-sided 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,...,Ii over the grid points, where Ii 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 re-entrant 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 two-dimensional simply-connected 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, Ii, at the next corner, where Ii 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 I1 = I and I2 = J.

The boundary specification is then completed by positioning I-2 points on the lower and upper boundary segments of the physical region as desired, and J-2 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 I-1, and those on the left and right segments are placed, respectively, in 1,j and I,j for j from 2 to J-1.

          This process of boundary specification can be most easily understood by viewing the rectangular boundary of the transformed region, with I equally-spaced points along two opposite sides and J equally-spaced 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 doubly-connected 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 re-entrant 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 J-1.

          The conceptual deformation of the rectangle for a C-type configuration is illustrated on p. 31. Here, with I1 the number of points on the segments 1-2 and 3-4 (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+I2-1.

          The first and last of these points must be coincident, i.e. I1,1 = I1 + I2-1,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 4-5 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 1-6 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 6-5 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(I1-1)+I2. Here the portions of the lower side of the rectangle, i.e., i from 2 to I1-1, and from I1+I2 to I-1 with j=1, are re-entrant 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 = I-i+1,1 for i=2 to I1-1 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 J2-J1+1 selected points on the segment 5-8 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 j-indices of the lower and upper sides, respectively, of the interior rectangle and I1 is the i-index of the left side of this rectangle. Similarly, J2-J1+1 points are positioned as desired on the segment 6-7 of the physical boundary and are placed in I2,j, where I2 is the i-index of the right side of the inner rectangle. Also I2-I1+1 points on the segments 5-6 and 8-7 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 J1-1 coordinate lines running from left to right below the inner boundary, and J-J2 lines running above the inner boundary. Similarly, there will be I1-1 lines running upward to the left of the interior boundary and I-I2 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 i-indices of the slit ends, 5 and 6, are I1 and I2, respectively, then the cartesian coordinates of I2-I1+1 points positioned as desired on the lower portion of the physical interior boundary, again proceeding from 5 to 6, are placed in a one-dimensional 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 one-dimensional 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 j-index 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., I1-1 lines on the left and I-I2 on the right. Similarly, there will be J1-1 lines below the interior boundary, and J-J1 above.

          Configurations, such as those illustrated on pp. 24-25, 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 i-indices I1 and I2, where I1-1 lines run vertically to the left of the slit and there are I2-I1+1 points on the slit. The vertical slit would be at i=I3 where there are I3-I2-1 vertical lines between this slit and the horizontal slit (and I-I2 lines to the right). Similarly, if the j-indices of the ends of the vertical slit or J1 and J2, there will be J1-1 horizontal lines below this slit and J-J2 lines above. With the j-index of the horizontal slit as J3, there will be J3-1 horizontal lines below this slit and J-J3 above. Provision will now have to be made in coding for two one-dimensional 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 2-3 and 6-7, 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, 3-4 and 5-6, 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 5-6, but proceeding from 6 to 5. Also values on the segments 2-1 and 7-8 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 3-block 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 8-15-9) are placed in i,1 of the first block, for i from the i-index at 8 to that at 9,proceeding clockwise from 8 to 9 on the physical boundary. (The difference between these two i-indices here is equal to the number of points on this interior boundary,less one.) Similary,the selected points on the left interior boundary (segment 4-5) are placed in i,1 of the second block for i from the i-index 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 i-indices 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 C-type lines occur in each block, and can be chosen independently. Here the segment 11-13 on the top of the first block interfaces with the corresponding segment on the top of the second block. The segment 12-14, 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 12-16, which forms the remainder of the bottom of the third block, interfaces with segment 11-17, which forms the remainder of the top of the second block. The segments 3-4 and 6-5 on the bottom of the second block interface with each other in the order indicated, as do also the segments 7-8 and 10-9 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 image-point 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 image-point 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 image-point array at the indices and block number of the other member of the pair. This correspondence is indicated on the block diagram on pp. 85-86 by the points enclosed in certain geometric symbols.

          Thus, for the 3-block configuration considered above, the indices (I1-i+1, 1) and block number 1, corresponding to a point on the segment 9-10 of the first block, would be placed in the image-point array at the point (i,1) on the segment 7-8 of this block, and vice versa. A similar pairing occurs for points on the segments 3-4 and 5-6 of the second block. The indices (I2-i+1, J2) and block number 2, corresponding to a point on the segment 11-13 of the second block would be placed in the image-point array at the point (i,Jl) of the first block on the segment 13-11 of that block, and vice versa. The indices (I3-I1+i, 1) and block number 3 (a point on the segment 12-14 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 12-14 of that block. Finally, the indices (i,1) and block number 3 (point on segment 16-12 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 17-11 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. 85-86 by geometric symbols.):

          As noted, all of this information would be input into the image-point 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 image-point 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 (I-1, 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. 35-37.

3. Sketch the grid for a O-grid, a C-grid, 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 C-type 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 cylindrical-coordinate 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 two-slit 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 two-block system for a circular pipe T-section. It is necessary to use a cross-sectional 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 six-block 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 hemi-spherical caps, these being coaxial with one inside the other. Sketch the grid and diagram the blocks for a three-block 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 two-body 0-grid 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 2-D composite system having two circular regions embedded in a grid which is generally rectangular. Let one of the circular regions have a cylindrical-coordinate 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 2-block 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 4-block composite system is used.