Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
VI. ELLIPTIC GENERATION SYSTEMS
As noted in Chapter II, the generation of a boundaryconforming coordinate system is accomplished by the determination of the values of the curvilinear coordinates in the interior of a physical region from specified values (and/or slopes of the coordinate lines intersecting the boundary) on the boundary of the region:
One coordinate will be constant on each segement of the physical boundary curve (surface in 3D), while the other varies monotonically along the segment (cf. Chapter II).
The equivalent problem in the transformed region is the determination of values of the physical (cartesian or other) coordinates in the interior of the transformed region from specified values and/or slopes on the boundary of this region, as discussed in Chapter II:
This is a more amenable problem for computation, since the boundary of the transformed region is comprised of horizontal and vertical segments, so that this region is composed of rectangular blocks which are contiguous, at least in the sense of being joined by reentrant boundaries (branch cuts), as described in Chapter II.
The generation of field values of a function from boundary values can be done in various ways, e.g., by interpolation between the boundaries, etc., as is discussed in Chapter VIII. The solution of such a boundaryvalue problem, however, is a classic problem of partial differential equations, so that it is logical to take the coordinates to be solutions of a system of partial differential equations. If the coordinate points (and/or slopes) are specified on the entire closed boundary of the physical region, the equations must be elliptic, while if the specification is on only a portion of the boundary the equations would be parabolic or hyperbolic. This latter case would occur, for instance, when an inner boundary of a physical region is specified, but a surrounding outer boundary is arbitrary. The present chapter, however, treats the general case of a completely specified boundary, which requires an elliptic partial differential system. Hyperbolic and parabolic generation systems are discussed in Chapter VII.
Some general discussion of elliptic generation systems has been given in Ref. [19], and numerous references to the application thereof appear in the surveys given by Ref. [1] and [5].
1. Generation Equations
The extremum principles, i.e., that extrema of solutions cannot occur within the field,that are exhibited by some elliptic systems can serve to guarantee a onetoone mapping between the physical and transformed regions (cf. Ref. [20] and [21]). Thus, since the variation of the curvilinear coordinate along a physical boundary segment must be monotonic, and is over the same range along facing boundary segments (cf. Chapter II), it clearly follows that extrema of the curvilinear coordinates cannot be allowed in the interior of the physical region, else overlapping of the coordinate system will occur. Note that it is the extremum principles of the partial differential system in the physical space, i.e., with the curvilinear coordinates as the dependent variables, that is relevant since it is the curvilinear coordinates, not the cartesian coordinates, that must be constant or monotonic on the boundaries. Thus it is the form of the partial differential equations in the physical space, i.e., containing derivatives with respect to the cartesian coordinates, that is important.
Another important property in regard to coordinate system generation is the inherent smoothness that prevails in the solutions of elliptic systems. Furthermore, boundary slope discontinuities are not propagated into the field. Finally, the smoothing tendencies of elliptic operators, and the extremum principles, allow grids to be generated for any configurations without overlap of grid lines. Some examples appear below:
There are thus a number of advantages to using a system of elliptic partial differential equations as a means of coordinate system generation. A disadvantage, of course, is that a system of partial differential equations must be solved to generate the coordinate system.
The historical progress of the form of elliptic systems used for grid generation has been traced in Ref. [1]. Consequently, references to all earlier work will not be made here. Numerous examples of the generation and application of coordinate systems generated from elliptic partial differential equations are covered in the above reference, as well as in Ref. [2].
A. Laplace system
The most simple elliptic partial differential system, and one that does exhibit an extremum principle and considerable smoothness is the Laplace system:

(1) 
This generation system guarantees a onetoone mapping for boundaryconforming curvilinear coordinate systems on general closed boundaries.
These equations can, in fact, be obtained from the Euler equations for the minimization of the integral

(2) 
as is discussed further in Chapter XI. Since the coordinate lines are located at equal increments of the curvilinear coordinate, the quantity ^{i} can be considered a measure of the grid point density along the coordinate line on which ^{i} varies, i.e., ^{i} must change rapidly in physical space where grid points are clustered. Minimization of this integral thus leads to the smoothest coordinate line distribution over the field.
With this generating system the coordinate lines will tend to be equally spaced in the absence of boundary curvature because of the strong smoothing effect of the Laplacian, but will become more closely spaced over convex boundaries, and less so over concave boundaries, as illustrated below. (In this and other illustrations and applications in two dimensions, ^{1} and ^{2} will be denoted and , respectively, while x and y will be used for x_{1} and x_{2}.)
In the left figure we have _{xx} > 0 because of the convex (to the interior) curvature of the lines of constant (lines). Therefore it follows that _{yy} < 0, and hence the spacing between the lines must increase with y. The lines thus will tend to be more closely spaced over such a convex boundary segment. For concave segments, illustrated in the right figure, we have _{xx} < 0, so that _{yy} must be positive, and hence the spacing of the lines must decrease outward from this concave boundary. Some examples of grids generated from the Laplace system are shown below. The inherent smoothness and the behavior near concave and convex boundaries are evident in these examples.
B. Poisson system
Control of the coordinate line distribution in the field can be exercised by generalizing the elliptic generating system to Poisson equations:

(3) 
in which the "control functions" P^{i} can be fashioned to control the spacing and orientation of the coordinate lines. The extremum principles may be weakened or lost completely with such a system, but the existence of an extremum principle is a sufficient, but not a necessary, condition for a onetoone mapping, so that some latitude can be taken in the form of the control functions.
Considering the equation ^{2} = Q and the figures above (P^{1} = P and P^{2} = Q in the illustrations here), since a negative value of the control function would tend to make _{yy} more negative, it follows that negative values of Q will tend to cause the coordinate line spacing in the cases shown above to increase more rapidly outward from the boundary. Generalizing, negative values of the control function Q will cause the lines to tend to move in the direction of decreasing , while negative values of P in ^{2} = P will cause lines to tend to move in the direction of decreasing . These effects are illustrated below for an line boundary:
With the boundary values fixed, the lines here cannot change the intersection with the boundary. The effect of the control function P in this case is to change the angle of intersection at the boundary, causing the lines to lean in the direction of decreasing .
These effects are illustrated in the following figures:
Here the lines are radial and the lines are circumferential. In the left illustration the control function Q is locally nonzero near a portion of the inner boundary as indicated, so the lines move closer to that portion of the boundary while in the right figure, P is locally nonzero, resulting in a change in intersection angle of the lines with that portion of the boundary. If the intersection angle, instead of the point location, on the boundary is specified, so that the points are free to move along the boundary, then the lines would move toward lines with lower values of :
In general, a negative value of the Laplacian of one of the curvilinear coordinates causes the lines on which that coordinate is constant to move in the direction in which that coordinate decreases. Positive values of the Laplacian naturally result in the opposite effect.
C. Effect of boundary point distribution
Because of the strong smoothing tendencies that are inherent in the Laplacian operator, in the absence of the control functions, i.e., with P_{i} = 0, the coordinate lines will tend to be generally equally spaced away from the boundaries regardless of the boundary point distribution. For example, the simple case of a coordinate system comprised of horizontal and vertical lines in a rectangular physical region, (cf. the right figure below) cannot be obtained as a solution of Eq. (3) with P=Q=0 unless the boundary points are equally spaced.
With _{yy} = _{xx} = 0, Eq. (3) reduces to
and thus P and Q cannot vanish if the point distribution is not uniform on the horizontal and vertical boundaries, respectively. With P=Q=0 the lines tend to be equallyspaced away from the boundary. These effects are illustrated further in the figures below. Here the control functions are zero in the left figure.
Although the spacing is not uniform on the semicircular outer boundary in this figure, the angular spacing is essentially uniform away from the boundary. By contrast, nonzero control functions in the right figure, evaluated from the boundary point distribution, cause the field spacing to follow that on the boundary. Thus, if the coordinate lines in the interior of the region are to have the same general spacing as the point distributions on the boundaries which these lines connect, it is necessary to evaluate the control functions to be compatible with the boundary point distribution. This evaluation of the control functions from the boundary point distribution is discussed more fully in Section 2 of this chapter.
D. General Poissontype systems
If a curvilinear coordinate system, ^{i} (i=1,2,3), which satisfies the Laplace system
is transformed to another coordinate system, (i = 1,2,3), then the new curvilinear coordinates, ^{i} satisfy the inhomogeneous elliptic system (cf. Ref. [19])

(4) 
where

(5) 
with the defined by the transformation from ^{i} to ^{i}:

(6) 
(It may be noted that if the subsequent transformation is onedimensional, i.e., if then only the three functions with i=1,2,3, are nonzero.)
These results show that a grid with lines concentrated by applying a subsequent transformation (often called a "stretching" transformation) to a grid generated as the solution of the Laplace system could have been generated directly as the solution of the Poisson system (4) with appropriate "control functions", derived from the subsequent concentrating transformation according to Eq. (6). Therefore, it is appropriate to adopt this Poisson system (4) as the generation system, but with the control functions specified directly rather than through a subsequent transformation.
Thus an appropriate generation system can be defined by Eqs. (4) and (5):

(7) 
with the control functions, considered to be specified. The basis of the generation system (7) is that it produces a coordinate system that corresponds to the subsequent application of a stretching transformation to a coordinate system generated for maximum smoothness. From Eq. (6), the three control functions (i = 1,2,3) correspond to onedimensional stretching in each coordinate direction and thus are the most important of the control functions. In applications, in fact, the other control functions have been taken to be zero, i.e., so that the generation system becomes

(8) 
It may be noted that, using Eq. (III37), Eq. (7) can be written as

(9) 
Actual computation is to be done in the rectangular transformed field, as discussed in Chapter II, where the curvilinear coordinates, ^{i}, are the independent variables, with the cartesian coordinates, x_{i}, as dependent variables. The transformation of Eq. (9) is obtained using Eq. (III71). Thus we have

(10) 
But ^{2}=0 and then using Eq. (7), we have

(11) 
This then is the quasilinear elliptic partial differential equation which is to be solved to generate the coordinate system. (In computation, the Jacobian squared, g, can be omitted from the evaluation of the metric coefficients, g^{ij} in this equation since it would cancel anyway, cf. Eq. III38.) As noted above, the more common form in actual use has been that with only three control functions, Eq. (8), which in the transformed region is

(12) 
Most of the following discussion therefore will center on the use of this last equation as the generation system. This form becomes particularly simple in one dimension, since then we have

(13) 
which can be integrated to give, with
The onedimensional control function corresponding to a distribution x() thus is given by

(14) 
In two dimensions, Eq. (11) reduces to the following form, using the twodimensional relations given in Section 8 of Chapter III (with ^{1} = and ^{2} = )

(15) 
where
with = x + y and
This corresponds to the following system in the physical space, from (7),

(19a) 

(19b) 
where g = (x_{ }y_{ }  x_{ }y_{ })^{2}.
The twodimensional form of the simpler generation system (12) with only two control functions is

(20) 
for which the system in the physical space is, from (8),

(21a) 

(21b) 
This generation system has been widely used, and a number of applications are noted in Ref. [1] and [5]. Several examples appear in Ref. [2].
Substitution of (3) in (10) gives the transformation of the original Poisson system (3) as

(22) 
This generation system has also been widely used, cf. Ref. [1] and [2], and the twodimensional form is

(23) 
corresponding in the physical space to

(24a) 

(24b) 
This system has also been widely used (cf. Ref. [1] and [5]), and its use predates that of Eq. (21). In general, however, the form of (12), corresponding to the system (8), is probably preferable over that of (22), which corresponds to (3), because of the simple form to which the former reduces in one dimension, and because the control functions in (8) are orders of magnitude smaller than those in (3) for similar effects.
E. Other systems
Other elliptic systems of the general form (4) have been considered, such as with P^{i}=gP_{i} where the P_{i} are the specified control functions, and with , where D is the control function. The latter form puts Eq. (4) in the form of a diffusion equation with the control function in the role of a variable diffusivity:

(25) 
This system also corresponds to the Euler equations for maximization of the smoothness, but now with the coefficient, D, serving as a weight function, i.e., multiplying the integrand in Eq. (2), so that the smoothness is emphasized where D is large. Both of these systems have actually been implemented only in two dimensions, although the formulations are general. Specific references to these and other related systems are given in Ref. [1] and [5].
Another elliptic system for the generation of an orthogonal grid has been constructed by combining the orthogonality conditions, ) with a specified distribution of the Jacobian over the field, . (This system is discussed further in Chapter IX.) Some twodimensional applications appear in Ref. [2], as noted in Ref. [5].
The secondorder systems allow the specification of either the point distribution on the boundary (Dirichlet problem):
or the coordinate line slope at the boundary (Neumann problem):
but not both. Thus it is not possible with such systems to generate grids which are orthogonal at the boundary with specified point distribution thereon. (This assumes that the control functions are specified. It is possible to adjust the control functions to achieve orthogonality at the boundary as is discussed in Section 2.)
A fourthorder elliptic system can be formulated by replacing the Laplacian operator, ^{2}, with the biharmonic operator, ^{4}. The analogous form to (4) then is

(26) 
which can be implemented as a system of two secondorder equations:

(27a) 

(27b) 
From (III71) and (22) above, the transformed system is

(28a) 

(28b) 
This generation system, being of higher order, allows more boundary conditions, so that the coordinate line intersection angles, as well as the point locations, can be specified on the boundary. It is therefore possible with this system to generate a coordinate system which is orthogonal at the boundary with the point distribution on the boundary specified, and for which the first coordinate surface off the boundary is at a specified distance from the boundary:
This allows segmented grids to be patched together with slope continuity as discussed in Chapter II.
In the above discussions, generation systems have been formulated based on linear differential operators in the physical space, e.g., the Laplacian with respect to the cartesian coordinates, resulting in quasilinear equations in the transformed space where the computation is actually performed. It is also possible to formulate the generation system using linear differential operators in the transformed space, e.g., the Laplacian with respect to the curvilinear coordinates:

(29) 
The use of some such generation systems is noted in Ref. [1], and such a biharmonic system is noted in Ref. [5]. Although this certainly produces simpler equations to be solved, since the computation is done in the transformed space, such systems transform to quasilinear equations in the physical space, and hence the extremum principles are lost in the physical space. This means that there is a possiblity of coordinate lines overlapping in general configurations. Therefore it is generally best to formulate the generation system using linear operators in the physical space.
As noted above, other variations of elliptic systems of the type discussed here are noted in Ref. [1] and [5]. Elliptic generation systems may also be produced from the Euler equations resulting from the application of variational principles to produce adpative grids, as is discussed in Chapter XI. Still another system, based on the successive generation of curved surfaces in the threedimensional region, is given in Section 3B of this chapter. Finally, quasiconformal mapping (Ref. [22] and [23]) is another example of an elliptic generation system.
2. Control Functions
For the elliptic generation system given by Eq. (12), the control functions that will produce a specified line distribution for a rectangular region, and for an annular region, are given as Eq. (14) and in Exercise 8, respectively. These functions could be used in other regions, of course, with the same general effect. In such extended use, the former would be more appropriate for simplyconnected regions, while the latter would be appropriate for multiplyconnected regions. Use of the rectangular function in a multiplyconnected region produces a stronger concentration than was intended because of the concentration over convex boundaries that is inherent in Poissontype generation systems (cf. Section 1A).
With generation systems of the Poisson type, negative values of the control function P^{i} > in Eq. (4), or P^{i} in Eq. (8) (since g^{ii} > 0), will cause the ^{i} coordinate lines to concentrate in the direction of decreasing ^{i} (cf. Section 1B). Several approaches to the determination of these control functions are discussed below.
A. Attraction to coordinate lines/points
This effect can be utilized to achieve attraction of coordinate lines to other coordinate lines and/or points by taking the form of the control functions to be, in 2D, (again with ^{1}=, ^{2}=, P^{1}=P, P^{2}=Q)

(30) 
and an analogous form for Q(,) with and interchanged. (Here the subscripts identify particular lines and are not to be confused with the superscripts used to refer to the curvilinear coordinates in general.) In this form, the control functions are functions only of the curvilinear coordinates.
In the P function, the effect of the amplitude a_{i} is to attract lines toward the _{i}line:
while the effect of the amplitude b_{i} is to attract lines toward the single point (_{i},_{i}):
Note that this attraction to a point is actually attraction of lines to a point on another line, and as such acts normal to the line through the point. There is no attraction of lines to this point via the P function. In each case the effect of the attraction decays with distance in  space from the attraction site according to the decay factors, c_{i} and d_{i}. This decay depends only on the distance from the _{i}line in the first term, so that entire lines are attracted to the entire _{i} line. In the second term, however, the decay depends on both the and distances from the attraction point (_{i},_{i}), so that the effect is limited to portions of the lines. With the inclusion of the sign changing function, the attraction occurs on both sides of the line, or the (_{i},_{i}) point, as the case may be. Without this function, attraction occurs only on the side toward increasing , with repulsion occuring on the other side. A negative amplitude simply reverses all of these effects, i.e., attraction becomes repulsion, and vice versa. The effect of the Q function on lines follows analogously.
In the case of a boundary that is an line, positive amplitudes in the Q function will cause lines off the boundary to move closer to the boundary, assuming that increases off the boundary. The effect of the P function will be to alter the angle at which the lines intersect the boundary, if the points on the boundary are fixed, with the lines tending to lean in the direction of decreasing . These effects have been noted in figures above, and further examples are given below:
The first two figures here show the result of attraction to the two circled points, in comparison with the case with no control function. The last figure illustrates strong attraction to the coordinate line coincident with the inner boundary and the branch cut in this Ctype system. If the boundary is such that decreases off the boundary, then the amplitudes in the Q function must be negative to achieve attraction to the boundary. In any case, the amplitudes a_{i} cause the effects to occur all along the boundary (as in the last figure above), while the effects of the amplitudes b_{i} occur only near selected points on the boundary (second figure above).
In configurations involving branch cuts, the attraction lines and/or points in this type of evaluation of the control function strictly should be considered to exist on all sheets. In the 0type configuration shown on p. 29, where the two sides of the cut are on opposite sides of the transformed region, the control function P for attraction to the _{i}line must be constructed as follows: In the figure below,
when the attraction line is the _{i}=2 line, the =I1 line experiences a counterclockwise attraction to this line at a distance of (I1)2. However, the _{i}=2 attraction line also appears as a I+(_{i}1)=I+(21)=I+1 attraction line on the next sheet as the cut is crossed. Therefore, the =I1 line also experiences a clockwise attraction to this I+1 line at a distance of (I+1)(I1)=2, and this attraction is, of course, stronger than the first mentioned. In fact, since the attraction line is repeated on all sheets there strictly must be a summation over all sheets in Eq. (30), i.e., a summation over k, with _{i} replaced by _{i}+k where is the jump in at the cut (=I1 in the above figure). Thus _{i} in Eq. (30) would be replaced by the _{i}+k, and the rightside would be summed from k= to +=. However, because of the exponential decay, the terms decrease rapidly as k increases, so that only the term with the smallest distance in the k summation really needs to be included, i.e., only the term giving clockwise attraction at a distance of 2 from the attraction line for the =I1 line in the above figure. Since there is no jump in across the cut in this configuration, the evaluation of Q is affected by this out only through the replacement of _{i} as above in the term for the point attraction, with summation over k of only this part of the right side. Again only the term with the smallest distance need actually be included.
For the Ctype configuration on p. 30, with the two sides of the out on the same side of the transformed region, is reflected in the cut, and the construction of the control function Q is as follows. With reference to the figure below,
the attraction line, _{i}=2, is located on both sides of the cut in this configuration. Now the =3 line above the cut experiences a downward attraction toward the _{i}=2 attraction line at a distance of 32=1. Strictly speaking, this line above the cut should also experience a downward attraction toward the portion of the _{i}=2 attraction line below the cut as it appears on the next sheet (and, in fact, on all other sheets), i.e., at _{cut}(_{i}_{cut}), where _{cut} is the value of on the cut (_{cut}=1 here). This attraction line on the next sheet is at a distance [_{cut}(_{i}_{cut})] from the line of interest, i.e., at 3[1(21)]=3 from the =3 line above the cut. This attraction line on the next sheet is therefore farther away and hence its effect can perhaps be neglected. However, for lines between the attraction line and the cut, the effect of the attraction line on the next sheet should be considered. In any case it is necessary to take into account the attraction lines appearing on the next sheet, those on all other sheets being too far away to be of consequence. Here the evaluation of the control function P is affected by the cut only through the point attraction part, with _{i} replaced as above.
The third type of cut, illustrated on p. 40, for which the two sides of the cut face across a void of the transformed region, is treated by replacing with _{i} in both the control functions, where 1 is the number of lines in the void. There is no additional summation in this case.
The case on p.52, where the coordinate species changes sign at the cut, requires individual attention at each cut. For example, the contribution to the control functions in region A at a point (,) from an attraction site (_{i},_{i}) in region B would be evaluated using distances of (_{max})+(_{max}_{i}) and (_{i}) in place of _{i} and _{i}, respectively.
B. Attraction to lines/points in space
If the attraction line and/or points are in the field, rather than on a boundary, then the above attraction is not to a fixed line or point in space, since the attraction line or points are themselves determined by the solution of the generation system and hence are free to move. It is, of course, also possible to take the control functions to be funtions of x and y instead of and , and thus achieve attraction to fixed lines and/or points in the physical field. This case becomes somewhat more complicated, since it must be ensured that coordinate lines are not attracted parallel to themselves.
With the attraction discussed in the previous section, lines are attracted to other lines, and lines are attracted to other lines. It is unreasonable, of course, to attempt to attract lines to lines, since that would have the effect of collapsing the coordinate system. When, however, the attraction is to be to certain fixed lines in the physical region, defined by curves y=f(x), care must be exercised to avoid attempting to attract coordinate lines to specified curves that cut the coordinate lines at large angles. Thus, in the figure below,
it is unreasonable to attract glines to the curve y f(x),
while it is natural to attract the qlines to this curve.
However, in the general situation, the specified line y=f(x) will not necessarily be aligned with either a or line along its entire length. Since it is unreasonable to attract a line tangentially to itself, some provision is necessary to decrease the attraction to zero as the angle between the coordinate line and the given line y=f(x) approaches 90°. This can be accomplished by multiplying the attraction function by the cosine of the angle between the coordinate line and the line y=f(x). It is also necessary to change the sign on the attraction function on either side of the line y=f(x). This can be done by multiplying by the sine of the angle between the line y=f(x) and the vector to the point on the coordinate line.
These two purposes can be accomplished as follows. Let a general point on the line be located by the vector (x,y), and let the attraction line y=f(x) be specified by the collection of points S(x_{i},y_{i}) i=1,2,...N. Let the unit tangent to the attraction line be (x_{i},y_{i}), and the unit tangent to a line be ^{ ()}. Then, with the unit vector normal to the twodimensional plane, and with reference to the following figure,
the control functions, P(x,y) and Q(x,y), may logically be taken as

(31) 
The equation for Q simply has replaced by in the above. These functions depend on x and y through both and ^{ ()} or ^{ ()}, and thus must be recalculated at each point as the iterative solution proceeds. This form of coordinate control will therefore be more expensive to implement than that based on attraction to other coordinate lines.
There is no real distinction between "line" and "point" attraction with this type of attraction. "Line" attraction here is simply attraction to a group of points that form a line, y=f(x). If line attraction is specified then the tangent to the line y=f(x) is computed from the adjacent points on the line. If point attraction is specified, then the "tangent" must be input for each point. The unit tangents to the coordinate lines are computed from Eq. (III3):
The presence of branch cuts introduces no complication with this type of attraction since the distances involved are in terms of the cartesian coordinates, rather than the curvilinear coordinates. This form of attraction makes the control functions dependent on both the curvilinear and cartesian coordinates, and thus attraction to space lines and/or points involves more complicated equations in the transformed region than does attraction to other coordinate lines and/or points, since for the former, coefficients of the first derivatives are functions of the dependent variables. Attraction to lines and/or points in space has not been widely used, and the use of Eq. (31) has not been fully tested.
C. Evaluation along a coordinate line
As has been noted above, if it is desired that the spacing of the coordinate lines in the field generally follow that of the points on the boundary, the control functions must be evaluated so as to correspond to this boundary point distribution. This can be accomplished as follows. (The developments in this and the next two sections are generalizations of that given in Ref. [12], and other works cited therein and in Ref. [5].)
The projection of Eq. (12) along a coordinate line on which ^{l} varies is found by forming the dot product of this equation with the base vector , which is tangent to the line.
Thus we have

(32) 
Now assume for the moment that the two coordinate lines crossing the coordinate line of interest do so orthogonally. Then on this line we have
and
which leads to an explicit equation for P_{l} on the coordinate line of interest:

(33) 
If it is further assumed for the moment that the two coordinate lines crossing the coordinate line of interest are also orthogonal to each other, i.e., complete orthogonality on the line of interest,
we have on this line g^{ij} = _{ij}g_{ii} and g_{ij} = _{ij}gii Also, from Eq. (III38),
since . But also by Eq. (III16) g = g_{ll}g_{mm}g_{nn} so that g^{ll}g_{ll}=1. Then Eq. (33) becomes

(34) 
which can also be written, using Eq. (III3), as

(35) 
By Eq. (III7) the derivative of arc length along the coordinate line on which ^{l} varies is

(36) 
Then

(37) 
so that the logarithmic derivative of arc length along this coordinate line is given by

(38) 
which is exactly the i=1 term in the summation in Eq. (34).
The unit tangent to a coordinate line on which ^{m} varies is

(39) 
and the derivative of this unit tangent with respect to arc length is a vector that is normal to this line, the magnitude of which is the curvature, K, of the line. The unit vector in this normal direction is the principal normal, , to the line.
Thus, using Eq. (36),

(40) 
Then

(41) 
so that the curvature is

(42) 
The component of K^{m}^{m} along the coordinate line on which ^{l} varies
is given by

(43) 
since for . Then the two terms of the summation in Eq. (34) for which can be written as

(44) 
Thus Eq. (34) can be written

(45) 
where (l,m,n) are cyclic, and using Eq. (III3), we have

(46) 
and

(47) 
with an analogous equation for (K^{n}^{n})^{(1)}. The arc length in the expressions (45) for the control function P_{l} along the coordinate line on which ^{l} varies can be determined entirely from the grid point distribution on the line using Eq. (46). The other two terms in P_{l}, however, involve derivatives off this line and therefore must either be determined by specifications of the components of the curvature, K, of the crossing lines along the line of interest, or by interpolation between values evaluated on coordinate surfaces intersecting the ends of this line.
If it is assumed that the curvatures of these crossing lines vanish on the coordinate line of interest, then the last two terms in Eq. (45) are zero, and the control function becomes simply

(48) 
and then can be evaluated entirely from the specified point distribution on the coordinate line of interest.
The neglect of the curvature terms, however, is illadvised since the elliptic system already has a strong tendency to concentrate lines over a convex boundary, as has been discussed earlier in this chapter. Therefore neglect of the curvature terms will result in control functions which will produce a stronger concentration than intended over convex boundaries (and weaker over concave). When interpolation from the end points is used to determine the curvature term, the entire term (K) should be interpolated, since individual interpolation of the vectors _{l} and (_{m})_{ m} can give an inappropriate value for the dot product.
It should be noted that the assumptions of orthogonality, and perhaps vanishing curvature, that were made in the course of the development of these expressions for the control functions on a coordinate line are not actually enforced on the resulting coordinate system, but merely served to allow some reasonable relations for these control functions corresponding to a specified point distribution on a coordinate line to be developed. This should not be considered a source of error since the control functions are arbitrary in the generation system (12).
D. Evaluation on a coordinate surface
In a similar fashion, expressions for the control functions on a coordinate surface on which ^{l} is constant can be obtained from the projections of Eq. (12) along the two coordinate lines lying on the surface, i.e., the lines on which ^{m} and ^{n} vary, (l,m,n) being cyclic.
These projections are given by Eq. (32) with l replaced by m and n, respectively. If it is assumed for the moment that the coordinate line crossing the coordinate surface of interest is orthogonal to the surface
then , so that P_{l} is removed from both of these two equations to yield the equation
and an analogous equation with m and n interchanged. Solution of these two equations for P_{m} and P_{n} then yields

(49) 
with an analogous equation for P_{n} with m and n interchanged. Since g_{lm} = g_{ln} = 0, we have by Eq. (III38), g^{lm} = g^{ln} = 0. Therefore only the five terms, ll, mm, nn, mn, nm, are nonzero in the summation. Also from Eq. (III38) we have
since here . An analogous equation for g^{nn} is obtained by interchanging m and n.
Then Eq. (49) can be rewritten as

(50) 
and an analogous equation for P_{n} with m and n interchanged. But, again using Eq. (III38), we have
Therefore

(51) 
and the analogous equation with m and n interchanged. This can also be written as

(52) 
and the analogous equation for P_{n}.
All of the terms, except the first, in the above equations can be evaluated completely from the point distribution on the coordinate surface of interest.
From Eq. (47) the first term in (52) can be written

(53) 
where (K^{l}^{l})^{m} and (K^{l}^{l})^{n} are the components of the curvature K for the coordinate line crossing the coordinate surface of interest along the two coordinate lines on the surface.
These quantities must be either specified on the surface or interpolated from values evaluated on its intersections with the other coordinate surfaces. If it is further assumed that the curvature of the crossing line vanishes at the surface, then this first term in Eq. (52) vanishes also.
As was noted for the control functions on a line, the curvature terms should not be neglected, however, else the concentration will be stronger than intended over convex boundaries and weaker over concave. Also, it is the entire term K which should be interpolated, not the individual vectors involved, else the dot product can have inappropriate values.
E. Evaluation from boundary point distribution
Using the relations developed in the previous two sections for the control functions on a coordinate line and on a coordinate surface, an interpolation procedure can be formulated for evaluation of the control functions in the entire field. If the point distribution is specified on all the boundary surfaces of a threedimensional field, the control functions can be evaluated on these boundaries using the relations in Section D, and then the control functions in the entire threedimensional field can be interpolated from these values on the bounding surfaces using transfinite interpolation (discussed in connection with algebraic grid generation in Chapter VIII.)
To be definite, consider a general threedimensional region bounded by six curved sides:
with curvilinear coordinates as shown, which transforms to a rectangular block. From Eq. (52) the two control functions, P_{j} and P_{k}, can be evaluated from the specified boundarypoint distribution on the two faces on which ^{i} is constant, i.e., the left and right faces in the figure. Similar evaluations yield two control functions on each face, with the result that the control function P_{k} will be known on the four faces on which ^{k} varies, i.e., the front, back, left, and right faces in the figure. Thus, in general, interpolation for the control function P_{k} in the interior of the region is done from the boundary values on the four faces on which ^{k} varies:

(54) 
Here I^{i} is the maximum value of ^{i}, etc., i.e., , i=1,2,3. In an analogous manner all three control functions can be determined in the interior of the region.
It may be desirable in some cases to generate a twodimensional coordinate system on a curved surface, as discussed in Section 3, rather than specifing the point distribution on the surface. The two control functions needed on the surface for this purpose can be determined by interpolation from values evaluated on the four edges of the surface:
Eq. (45) allows the control function P_{i} to be evaluated on the edges on which ^{i} varies, i.e., the top and bottom edges in the figure. This control function on the surface can then be evaluated by interpolation between these two edges:

(55) 
Both of the necessary control functions on the surface can thus be determined from the specified boundary point distributions on the edges of the surface.
F. Iterative determination
As noted above, a secondorder elliptic generation system allows either the point locations on the boundary or the coordinate line slope at the boundary to be specified, but not both. It is possible, however, to iteratively adjust the control functions in the generation system of the Poisson type discussed above until not only a specified line slope but also the spacing of the first coordinate surface off the boundary is achieved, with the point locations on the boundary specified.
In three dimensions the specification of the coordinate line slope at the boundary requires the specification of two quantities, e.g., the direction cosines of the line with two tangents to the boundary.
The specification of the spacing of the first coordinate surface off the boundary requires one more quantity,
and therefore the three control functions in the system (12) are exactly sufficient to allow these three specified quantities to be achieved, while the one boundary condition allowed by the secondorder system provides for the point locations on the boundary to be specified.
The capability for achieving a specified coordinate line slope at the boundary makes it possible to generate a grid which is orthogonal at the boundary, with a specified point distribution on the boundary, and also a specified spacing of the first coordinate surface off the boundary. This feature is important in the patching together of segmented grids, with slope continuity, as discussed in Chapter II, for embedded systems.
An iterative procedure can be constructed for the determination of the control functions in two dimensions as follows (cf. Ref. [25]): Consider the generation system given by Eq. (20). On a boundary segment that is a line of constant we have _{ } and _{ } known from the specified boundary point distribution
also , the spacing off this boundary, is specified
as is the condition of orthogonality at the boundary, i.e., ,
But specification of , together with the condition provides two equations for the determination of x_{ } and y_{ } in terms of the already known values of the x_{ } and y_{ }. Therefore _{ } is known on the boundary.
Because of the orthogonality at the boundary, Eq. (20) (Eq. (23) is used instead in Ref. [25]) reduces to the following equation on the boundary:
Dotting _{ } and _{ } into this equation, and again using the condition of orthogonality, yields the following two equations for the control functions on the boundary:

(56a) 

(56b) 
All of the quantities in these equations are known on the boundary except _{ }. (On a boundary that is a line of constant , the same equations for the control functions result, but now with _{ } the unknown quantity.)
The iterative solution thus proceeds as follows:
(1). Assume values for the control function on the boundary.
(2). Solve Eq. (20) to generate the grid in the field.
(3). Evaluate _{ } on line boundaries, and _{ } on line boundaries, from the result of Step (2), using onesided difference representations. Then evaluate the control functions on the boundary from Eq. (56).
Evaluate the control functions in the field by interpolation from the boundary values.
Steps (2) and (3) are then repeated until convergence.
This type of iterative solution has been implemented in the GRAPE code of Ref. [24]  [26], some results of which are shown below:
These grids are orthogonal at the boundary, and the spacing of the first coordinate surface (line in 2D) off the boundary is specified at each boundary point, the locations of which are specified.
An iterative solution procedure for the determination of the three control functions for the general threedimensional case can be constructed as follows. Eq.(52) gives the two control functions, P_{m} and P_{n}, for a coordinate surface on which ^{l} is constant (1,m,n cyclic) for the case where the coordinate line crossing the surface is normal to the surface. Taking the projection of the generation equation (12) on the coordinate line along which ^{l} varies, we have on this same surface,
since g_{lm} = g_{ln} = 0 on the surface. Using the relations for the metric components obtained for this situation in Section D, this equation reduces to

(57) 
Since the coordinate line intersecting the surface is to be normal to the surface, we may write

(58) 
since
using the identity (III9). Eq. (57) can then be written
With the spacing along the coordinate line intersecting the surface specified at the surface, we have known on the surface. Since all the quantities subscripted m or n in Eq. (52) and (59) can be evaluated completely from the specified point distribution on the surface, we then have all quantities in these equations for the three control functions on the surface known except for (g_{ll})_{ l} and (_{l})_{ l}. These two quantities are not independent, and using Eq. (58), we have

(60) 
Recall also that (_{l}) = _{ l}_{ l}.
Therefore, with the control functions in the field determined from the values on the boundary by interpolation, as discussed in the preceeding section, Eq. (52) and (59) can be applied to determine the new boundary values of the control functions in terms of the new values of (_{l})_{ l} in an iterative solution. Upon convergence, the coordinate system then will have the coordinate lines intersecting the boundary normally at fixed locations and with the specified spacing on these lines off the boundary.
A similar iterative determination of the two control functions for use in generating a coordinate system on a surface can be set up using only Eq. (52), and the analogous equation for P_{n}, with the first term either omitted, amounting to the assumption of vanishing curvature of the crossing line at the surface or with this term considered as specified on the surface, either directly or by interpolation from the edges of the surface (the edges assumed to be on coordinate lines) to provide two equations for the two control functions, P_{m} and P_{n}, on these edges. Here the two dimensional surface coordinate system is to be orthogonal on the bounding edges of the surface with the spacing off the edges, and the point distribution thereon, specified on these edges.
Since the coordinate system is to be orthogonal on the edges, g_{mn} = 0 there so that the last term in the bracket in Eq. (52) vanishes. Eq. (52) then reduces to the following expression

(61) 
and the analogous equation for P_{n} is

(62) 
These equations can also be written

(63) 
and

(64) 
If the point distribution is specified along an edge on which ^{m} varies, then g_{mm}, _{m}, and (_{m}_{ m} can all be calculated on this edge. The specification of the spacing from this edge to the first coordinate line off the edge determines g_{nn} on this edge. Also, because of the orthogonality on the edge, we have

(65) 
where is the unit normal to the surface. Note that will vary along the edge if the surface is curved. Since the surface normal, , will be known, all quantities in Eq. (63) and (64) are known except (g_{nn})_{ n} and (_{n}_{ n}. These two quantities are not independent and, in fact,

(66) 
On edges along which ^{n} varies, Eq. (65) and (66) are replaced by

(67) 
and

(68) 
and it is (g_{mm})_{ m} and (_{m})_{ m} that are not known.
The iterative solution then proceeds as described above, with the new control functions being determined from Eq. (63) and (64), together with Eq. (65) and (66) on edges along which ^{m} varies, or with Eq. (67) and (68) on edges along which ^{n} varies.
3. Surface Grid Generation Systems
The grid generation systems discussed in the preceeding sections of this chapter have been for the generation of curvilinear coordinate systems in general threedimensional regions. Twodimensional forms of these systems serve to generate curvilinear coordinate systems in general twodimensional regions in a plane. It is also of interest, however, to generate twodimensional curvilinear coordinate systems on general curved surfaces.
Here the surface is specified, and the problem is to generate a twodimensional grid on that surface, the third curvilinear coordinate being constant on the surface. The configurations of the transformed region will be the same as described in Chapter II for twodimensional systems in general, i.e., composed of contiguous rectangular blocks in a plane, with point locations and/or coordinate line slopes specified on the boundaries. These boundaries now correspond to bounding curves on the curved surface of the physical region. The problem is thus essentially the same as that discussed above for twodimensional plane regions, except that the curvature of the surface must now enter the partial differential equations which comprise the grid generation system.
As for general regions, algebraic generation systems based on interpolation can be constructed, and such systems are discussed in Chapter VIII. The problem can also be considered as an elliptic boundaryvalue problem on the surface with the same general features discussed above being exhibited by the elliptic generation system.
A. Surface grid generation
An elliptic generation system for surface grids can be devised from the formulae of Gauss and Beltrami, of. Ref. [27]. Some related, but less general, developments are noted in Ref. [9] and [5]. The starting point is the set formed by the formulae of Gauss for a surface, which for a surface, ^{ } = constant, ( = 1,2, or 3) are given by Eq. (34) of Appendix A:

(69) 
where the variation of the indicies , and is over the two coordinate indices different from . (Greek coordinate indices are used here to set apart the coordinates generated on a surface from those generated in a threedimensional region in general).
The unit normal ^{ }, the coefficients b_{ } and the surface Christoffels () have all been defined in Eq. (15), (20), and (33) of Appendix A, respectively. The indices , each assume the two values different from . For each , with (,,) taken in cyclic order, we have
with assuming the two values, ,.
A surface grid generation system that is analogous in form to that based on Poissontype equations in a plane given earlier can be constructed by multiplying Eq. (70a,b,c), respectively, by G_{ }g^{}, 2G_{ }g^{}, G_{ }g^{} and adding. This given, after some algebra,

(71) 
where
The quantities and are the local principal curvatures of the surface _{ }=constant. It must be noted here the R^{ } as defined in (74) is based on the intrinsic values of b_{ }. That is, the b_{ } are solely determined by the data and coordinates as available in the surface. If, however, it is desired to use Eq. (71) for generating a series of surfaces in a threedimensional space, as in the following section, from the data of a given surface, then it is desirable to have an extrinsic form for b_{ }. To obtain the extrinsic form, we use Eq. (29) of Appendix A, i.e.,

(76) 
Equating the right hand sides of Eq. (69) and (76), taking the dot product with ^{ } on both sides, and noting that , we get

(77a) 
where

(77b) 
Thus

(78) 
The operator _{2} is called the Beltrami secondorder differential operator, and in general is defined as

(79) 
Thus

(80a) 

(80b) 
The generation system is now formed by taking, in analogy with the system (7),

(81a) 

(81b) 
where and each assume the two values , in the summation. Here the are the symmetric control functions. Thus the equations for the generation of surface grids are (with = x + y + z)

(82) 
where
The lefthand side of Eq. (82) here corresponds exactly to that of Eq. (15) for the plane. However, here we have in place of (18) the relations
The effect of the surface curvature enters through the inhomogenous term, in particular through R^{ ()} which is, in fact, equal to twice the product of and the mean curvature of the surface. Here, as for the plane, the control functions, , are considered to be specified. This system corresponds to the following system in the physical space, from (81),

(87a) 

(87b) 
Thus the Beltrami operator on the general surface replaces the Laplacian operator in the plane. If the surface is a plane, the Beltrami operator reduces to the Laplacian.
If only the two control functions and are included, the surface grid generation system reduces to the more practical system

(88) 
corresponding to the plane system given by Eq. (20). In the physical plane this system is

(89a) 

(89b) 
Clearly, we could also replace the system (89) with the simplier system

(90a) 

(90b) 
in analogy with the sytem (24) in the plane, to obtain the surface grid generation system

(91) 
which is analogous to the plane system (23).
Equation (71) is the basic equation for the generation of curvilinear coordinates in a given surface. From (74) the function R^{()} depends on the principal curvatures and . The sum is twice the mean curvature of the surface, and its value is invariant to the coordinates introduced in the given surface. If the equation of the surface in the form x_{3} = f(x_{1},x_{2}) is available, then from elementary differential geometry

(92) 
where
For arbitrary surfaces it is always possible to use a nummerical method, e.g., the least square method, to fit anequation in the form x_{3}=f(x_{1},x_{2}) or F(x_{1},x_{2},x_{3})=0 and to obtain the needed partial derivatives to find + as a function of x_{1},x_{2},x_{3}. (Surface grids have been obtained for simply and double connected regions in a surface using the above method.)
It may be desired in some applications to generate a new coordinate system based on an already existing coordinate system in a given surface. In the formulation of this problem Eq. (71) can have the form of R^{()} given in (74), (75) or (78). Let the surface on which the new grid is to be generated be specified parametrically by

(93) 
(For example, the parameters (u,v) might be latitude and longitude on a spherical surface.) If the specified cartesian coordinates on the surface form a finite set of discrete points, a smooth interpolation scheme is needed to recover the differentiable functions in (93). To attain the desired smoothness in the parametric representation (93), it is generally preferable to divide the given surface into a suitable number of patches such that each patch is representable by a bicubio spline with suitable blending functions. Having once established the smooth parametric functions (93), it is now possible to introduce any other desired coordinate system, say (^{ }, ^{ }) on the surface. For example, a surface coordinate system ^{ }, ^{ } of the configuration
might be generated on a surface defined by the parametric coordinates (u,v) in a latitudelongitude configuration:
Alternatively, a surface may be defined in terms of crosssections, in which case one of the parametric coordinates (u,v) runs around the section and the other connects the sections:
The fundamental equations for the generation of (^{ },^{ }) on the surface ^{ }=constant can be obtained from Eq. (75) in the form

(94) 
Here we have taken

(95a) 

(95b) 
Now using the chain rule of differentiation, we can write _{ }, _{ } _{ } etc., in terms of _{u}, _{v}, _{uu}, etc. Thus,
for example,
with the quantities as defined below. Substituting these derivatives in (94), and also in the expressions for ^{ } and ^{ }, we get

(99) 
where
To isolate the differential equations for u and v as dependent variables from (99), we take the dot product of Eq. (99) with _{u}, and then with _{v}, and use the conditions
Writing
the required equations are

(104a) 

(104b) 
where

(105a) 

(105b) 

(105c) 

(106a) 

(106b) 
Note that the metric quantities with an overbar relate to the surface definition in terms of the parametric coordinates and therefore can be calculated directly from the surface specification, Eq. (92).
Clearly we could redefine the control functions so that (104) is replaced by the following system, which is analogous in form to the plane system (20):

(107a) 

(107b) 
B. Threedimensional grids
As mentioned earlier, the system of Eq. (71), or Eq. (82), is also capable of generating threedimensional grids. This capability in the set of equations is incorporated through ^{ ()} as defined in (78).
The strategy of the method is to generate a series of surfaces on each of which two curvilinear coordinates vary while the third remains fixed. The variation along the third coordinate is specified as a surface derivative condition, which in turn depends on the given boundary data.
A study of Eq. (82)  (84) shows immediately that for the solution of Eqs. (82) we need to specify the values of and _{ } on certain curves ^{ }=constant. To fix ideas, let us consider the problem of coordinate generation between two given surfaces and as shown below:
The coordinates on these surfaces are ^{ } and ^{ }. To start solving these equations we need the values of and _{ } on the surfaces and . These values of are the input conditions for the solution of Eqs. (82), and are either prescribed analytically or numerically. On the other hand, the values of _{ } at B and are available easily, based on the values of , simply by numerical differentiation. The values of r_{ } in the field for each surface to be generated are then obtained by interpolation between the available values of (_{ })_{B} and (_{ })_{}. A simple formulae which has been used with success is

(108) 
where
4. Implementation
The setup of the transformed region configuration is done as described in Chapter II. This includes the placing of the cartesian coordinates of the selected points on the boundary of the physical region into _{ijk} for each block and the setting of the interface correspondence between points on the surrounding layer for each block and points inside the same, or another, block via input to an imagepoint array as described in Section 6 of Chapter II.
A. Difference equations
Implementation of an elliptic generation system then is accomplished by devising an algorithm for the numerical solution of the partial differential equations comprising the generation system. Recall that the use of the surrounding layer for each block, as described in Section 6 of Chapter II, allows the same difference representations that are used in the interior to be used on the interfaces. The usual approach is to replace all derivatives in the partial differential equations by secondorder central difference expressions, as given in Chapter IV, and then to solve the resulting system of algebraic difference equations by iteration. As noted above, most generation systems of interest are quasilinear, so that the difference equations are nonlinear.
A number of different algorithms have been used for the soltuion of these equations, including point and line SOR, ADI, and multigrid iteration (cf. Ref. [1] and [5]). For general configurations, point SOR is certainly the most convenient to code and has been found to be rapid and dependable, using overrelaxation, for a wide variety of configurations. The optimum acceleration parameters and the convergence rate decrease as the control functions increase in magnitude. Some consideration has been given to the calculation of a field of locallyoptimum acceleration parameters (cf. Ref. [1]), but the predicted values generally tend to be too high, and the desired increases in convergence rate were not obtained.
Since the system is nonlinear, convergence depends on the initial guess in iterative solutions. The algebraic grid generation procedures discussed in Chapter VIII can serve to generate this initial guess, and transfinite interpolation generally produces a more reliable initial guess than does unidirectional interpolation because of the reduced skewness in the former. In fact with strong line concentration, convergence may not be possible from an initial guess constructed from unidirectional interpolation, while rapid convergence occurs from an initial guess formed with transfinite interpolation. With the slab and slit configurations, the interpolation must be unidirectional between the closest facing boundary segments as illustrated below:
In a block structure, however, the slab/slit configuration can be avoided so that transfinite interpolation can be used.
Since the coordinate lines tend to concentrate near a convex boundary, very sharp convex corners may cause problems with the convergence of iterative solutions of the generation equations. These equations are nonlinear, and therefore convergence of an iterative procedure requires that the initial guess be within some neighborhood of the solution. With control functions designed to cause attraction to the boundary, it is possible for the coordinate lines to overlap a very sharp convex corner during the course of the iteration, even though a solution with no overlap exists:
This problem may be handled by first converging the solution with the coordinate lines artificially locked off the corner. Thus, if newly calculated values of the cartesian coordinates at a point during the iteration would cause this point to move farther from its present location than the distance to the adjacent point on the curvilinear coordinate line running to the corner, then these new values are replaced by the average of the coordinates of the old point and the adjacent point. After convergence, this lock is removed and final convergence to the solution is obtained. Note that this problem does not arise when the curvilinear coordinate line emanating from the corner is the same as that on the boundary, as in the Ctype configuration on p. 30 since then the lines do not wrap around the corner.
With very large cell aspect ratio, e.g., for g_{11}>>g_{22}, the generation equation is dominated by the term containing the second derivative along the curvilinear coordinate line on which the shorter are length lies. This causes the cartesian coordinates to tend strongly toward averages of adjacent points on this line during the course of the iteration. Therefore, when strong control functions are used to attract coordinate lines to the boundary in a Ctype configuration,
the points on the cut are very slow to move from the initial guess during the iteration. Convergence in such a case is very slow, and it is expedient to artificially fix the points on the cut as if it were a boundary. This will cause the coordinate lines crossing the cut to have discontinuous slopes at the cut, but since the spacing along these crossing lines is very small, the error thus incurred in difference solutions on the coordinate system is small.
B. Control functions
Several types of control functions have been discussed in Section 2 which serve to control the coordinate line spacing and orientation in the field. Most of these functions are set before the solution algorithm begins, either directly through input or by calculation from the boundary point distributions that have been input.
For the attraction to other coordinate lines/points, described in Section 2A, it is necessary to input the indices of the lines/points, i.e., the _{i} and _{i} of Eq. (30), to which attraction is to be made. In the case of attraction to lines, the line is identified by the single index which is constant thereon, while a point requires the specification of two indices (in 2D, with analogous generalization to 3D). The attraction amplitude and decay factor in Eq. (30) must also be input for each line/point. The control functions are then calculated at each point in the field (,) by performing the summations in Eq. (30), those summations being over all the attraction lines/points that have been input. As noted in Section 2A, these summations must also extend over some lines/points on other sheets across branch cuts in some cases.
This type of control function was used in the original TOMCAT code (cf. Ref. [1]), but is not really suitable as a primary means of control function definition because it only provides controlnot control to achieve a specified spacing distribution, since the appropriate values of the various parameters involved can only be determined by experimentation. This form does, however, still serve as a useful addition to other types of control, in that it allows particular ad hoc concentrations or adjustments of line spacing and orientation to be made. This can be particularly useful near the special points discussed in Chapter II where the grid line configuration departs locally from the usual simple coordinate line intersections.
The attraction to lines/points in space, implemented through Eq. (31), requires input similar to that just described, except that here the location of the attraction lines must be defined in the physical region by inputing a set of points along the line sufficient for its definition in discrete form. For attraction to a point a unit vector must also be input with each point. Again, attraction amplitude and decay factors must be input.
More important is the evaluation of the control functions from the boundary point distribution that has been input, as described in Section 2E. With the point distribution specified on a boundary line, the control functions on this line can be evaluated from Eq. (45)(47). Here the derviatives in Eq.(46) are best calculated from Eq. (36) and (37), using secondorder, central difference expressions along the line:
(Recall that .) The curvature terms given by Eq. (47), if included, must either be input at each point on the line, or, as is more likely, must be interpolated from values on the ends of the line. In this latter case, the ^{m} and ^{n} derivatives are off the line and are evaluated from the point distribution on the other coordinate lines intersecting the line of interest at its ends, using firstorder onesided difference expressions along these intersecting lines:
Onedimensional linear interpolation in ^{l} then serves to define the curvature term quantities at each point on the line of interest. Recall that it is the entire curvature term, rather than the individual vectors involved, that should be interpolated.
This evaluation determines the P_{l} control function on a boundary line on which ^{l} varies. Such an evaluation can be made on each edge of a surface, corresponding to one face of a block in three dimensions (cf. Section 6 of Chapter II). If it is desired to generate a twodimensional grid on this surface, control functions on the surface can be evaluated by interpolation from the function values on the edges, using linear interpolation between the two edges on which ^{i} is constant to evaluate P_{j}, and between the two edges on which ^{j} is constant to evaluate P_{i} (cf. the figure on p. 227). With the control functions thus defined on the surface, a twodimensional grid on the surface can now be generated using a surface grid generation system described as in Section 3. If the surface is a portion of the physical boundary, then a parametric definition of the surface will need to be input, so that the system defined by Eq. (107) can be applied. If, however, the surface is simply an interface between blocks, then its position is arbitrary and either a plane twodimensional generation system, such as Eq. (20),can be used, or surface curvature values could be input at each point on the surface and the surface system Eq. (82) used. The former is the more likely choice.
With the grid points on all the block faces defined, either by surface generation systems or by direct input, two control functions on each face can be evaluated from the surface point distribution using Eq. (52). Here the m and n derivatives are along coordinate lines on the surface and thus can be represented by secondorder central differences between points on the surface:
(Recall that (_{m})_{ m} can be expanded to _{ mm} for evaluation.) The lderivatives are off the surface and must either be specified by input at each point on the surface, or, as is more likely, must be interpolated from values evaluated along the coordinate lines intersecting the surface at its edges using firstorder, onesided difference expressions. The interpolation would here properly be twodimensional transfinite interpolation discussed in Chapter VIII.
This then serves to determine the two control functions P_{j} and P_{k}, on a surface on which ^{i} is constant (cf. the figure on p. 226), so that each control function will be defined on four faces of the block. Transfinite interpolation among these four faces then determines this control function in the interior of the block (cf. p. 227).
Another possibility is to evaluate the radius of curvature, , of the surface and to replace the curvature terms in Eq. (45) with (cf. Exercise 9). Here the radius of curvature should be interpolated unidirectionally between facing surfaces, and the same twodirectional transfinite interpolation used for the first term of the control function should be used for the spacing .
Still another approach is to solve the three generation system equations for the three control functions at each point using an algebraic grid, but with the offdiagonal metric elements set to zero. This will produce a grid which will have a greater degree of smoothness and orthogonality than the algebraic grid and yet has the same general spacing distribution. Here the result of the Computer Exercise 6 in Appendix C must be considered since the algebraic grid influences the spacing distribution.
In generation systems that iteratively adjust the control functions during the course of the solution of the difference equations (Section 2F) to achieve a specified spacing and angle of intersection, e.g., orthogonality, at the boundary, this spacing and intersection angle are input for each boundary point and it is, of course, not necessary to calculate the control functions beforehand. Several references to discussion of such systems are given in Ref. [5]. The GRAPE code is based on this approach, cf. the users manual Ref. [24].
C. Surface generation systems
A boundary surface in the physical region will typically be input by giving the cartesian coordinates of points on a series of crosssections, or other set of space curves:
These input points may then be splined to provide a functional definition of these curves. These curves are then parameterized in terms of normalized arc length thereon, i.e., so that this normalized parameter varies over the same range on each curve.
This normalized arc length then provides one parametric coordinate on the surface. The other coordinate is defined by connecting points at the same value of the first coordinate on the successive curves, again using a spline fit:
This second coordinate is then also expressed in terms of normalized arc length
(On a sphere these two parametric surface coordinates could correspond to longitude and latitude, the latter arising from the crosssections and the former from the connecting thereof.)
There are other techniques of surface definition and parameterization, cf. especially works on computeraided design, but the above decription is representative. The end result of this stage in any case is (u,v), i.e., the cartesian coordinates on the surface in terms of two surface parametric coordinates.
The two parametric coordinates (u,v) used to define the surface can also be adopted as the curvilinear coordinates defining the surface grid. However it is more likely that these coordinates were selected for convenience of input definition of the surface than for the definition of an appropriate grid thereon. This is particularly true when two such intersection surfaces, e.g., a wingbody, are input, each with its own set of parametric coordinates. Therefore, the surface grid generation system defined by Eq. (107) or (104) is used to generate a new surface coordinate system (^{ },^{ }) by generating values of the parametric coordinates (u,v) as functions of the curvilinear coordinates (^{ },^{ }), analogous to the plane generation systems which generate values of the cartesian coordinates as functions of the curvilinear coordinates. In fact, as noted above, the surface generation system degenerates to the plane system when the surface curvature vanishes.
With (u,v) now available, as described above, the metric elements with overbars can be calculated from the definitions in Eq. (103), using secondorder central differences for all derivatives as in the plane case. The quantities _{2}u and _{2}v are then calculated in the same manner from Eq. (106). Also the control fucntions are evaluated from the same relations given above for the plane case. All derivatives in the system (107) or (104) are represented by secondorder central difference epxressions, and the resulting nonlinear difference equations are solved as in the plane case.
Exercises
1. Demonstrate the validity of Eq. (4)  (6).
2. For plane polar coordinates (r,) defined as
show that the curvilinear coordinates
are solutions of the Laplace equations ^{2}=0, ^{2}=0.
3. Show that the onedimensional control function in Eq. (13) that is equivalent to the use of a subsequent exponential stretching transformation by the function given by Eq. (VIII26) is P = /I. Hint: x/L = (/I)
4. Show that the onedimensional control function in Eq. (13) that corresponds to a hyperbolic tangent stretching transformation by the function given by Eq. (VIII32) is
where u is given by Eq. (VIII33) and
5. Show that the onedimensional control function for the generation system given by Eq. (23) corresponding to a distribution x() is
Note that this control function will be considerably larger than that for Eq. (12) because of the higher inverse power of x
6. Show that a solution of Eq. (20), with P = P() and Q = Q(), for a rectangular region with and , and , is given by
where
7. Show that a solution of Eq. (20), with P = P() and Q(), for an annular region between two concentric circles of radius r_{1} and r_{2} is given, with and , by
where
with F() and G() given in the preceding exercise. Show also that for P = p_{ }/p_{ } and Q = q_{ }/q_{ }, R() and () become
8. From the result of the preceding exercise show that the control function Q() required to produce a specified radial distribution r() is given by
9. Show that the first term in the control function Q() given in the preceding exercise arises from the first term in Eq. (45), and that the second term arises from curvature term in (45).
10. Consider the generating system (23) for plane curvilinear coordinates. Let the control functions P and Q be defined as follows:
where k > 0 is a constant. Let it be desired to solve Eq. (23) for the generation of coordinates in the region of a circular annulus with = _{i}(r = 1) as the inner circle and = _{o}(r = R) as the outer circle. Considering the clockwise traverse in the direction as positive, set
where
in Eq. (23), and show that
where
11. Show that the control function P in Exercise 4 has the following values at the boundaries:
Note that an iterative procedure could be set up in the manner of Section 2F in which and A are determined from P(0) and P(I) and then these are used in the P() of Exercise 4 to define the control function in the field, rather than interpolating from the boundary values.
12. Show that the Beltrami operator reduces to the Laplacian for a plane surface.
13. Verify Eq. (71).
14. Verify Eq. (77a).
15. Consider a sphere of unit radius in which it is desired to introduce a coordinate system (,) in such a way that (i) is orthogonal, and (ii) the resulting metric coefficients g_{33} and g_{11} are equal. (Such systems are known to be isothermic.)
(a) Verify by inspection that for isothermic coordinates Eq. (80) is identically satisfied.
(b) To obtain the isothermic coordinates on a sphere set
and show that
(c) Show that the relation between the standard longitude and latitude surface coordinates and where 0 < 2 and 0 < < , is
16. Using Eq. (15), (20) and (21) of Appendix A, show that the sum of the principal curvatures, + , of a prolate ellipsoid defined as
is
17. Verify the correspondence between Eq. (90) and (91).
18. Verify Eq. (92).
19. Let (,) be the surface coordinates in the surface on which = constant. Then as shown in Appendix A, Eq. (21),
Let a new coordinate system (,) be introduced in the same surface such that = (,), and = (,) are admissible transformation functions.
(a) Use the chain rule of differentiation to show that the components of the normal to the surface are coordinate invariants, i.e.,
(b) Also show that on coordinate transformation
20. Let it be desired to obtain the 3D curvilinear coordinates in the region bounded by a prolate ellipsoid as an inner boundary ( = _{i}) and a sphere as an outer boundary ( = _{o}). The (x,y,z) for both the inner and outer bodies are given below in which and _{i} are the parameters of the ellipsoid:
(a) First write Eq. (91) as three equations in x, y and z for the generation of those surfaces on which = constant. Also set P = Q = 0 and transform the three equations mentioned above from to , where = _{i} + ().
(b) Assume the solution is
and compute all the needed derivatives to find g_{11}, g_{12}, g_{22} while keeping fixed. Also using Eq. (15), (20) and (21) of Appendix A obtain the expressions for the components of i^{(3)} and R^{(3)}.
(c) Use all the quantities obtained in (b) in the equations written in (a), and show that
where
21. Let a surface in the xyzspace by given as
Show the following:
(a) The components of the unit normal vector to the surface area are
(b) The element of area dA on the surface is
(c) The element of length ds of a surface curve is given by
(d) The sum of the principal curvatures is given by
22. (a) Show that the unit tangent vector to the curve of intersection of two surfaces F(x,y,z) = 0, G(x,y,z) = 0 is
where
and m,n,k are in the cyclic per mutations of 1,2,3.
(b) Using the formula for the normal vector to F constant, .i,e.,
find the Cartesian components of .
23. Verify Eq. (104).