Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
VIII. ALGEBRAIC GENERATION SYSTEMS
As noted earlier, the problem of generating a curvilinear coordinate system can be formulated as a problem of generating values of the cartesian coordinates in the interior of the rectangular transformed region from specified values on the boundaries. This, of course, can be done directly by interpolation from the boundaries, and such coordinate generation procedures are referred to as algebraic generation systems. Thus (^{1},^{2},^{3}) is given as a specific function of the curvilinear coordinates. This function contains certain coefficients which are determined so that the function matches specified values of the cartesian coordinates, and perhaps derivatives also, on the boundary and perhaps elsewhere. Evaluation of this interpolation function at constant values of the curvilinear coordinates then defines the coordinate system. Algebraic grid generation is discussed in Ref. [31] and [8], as well as in the surveys, Ref. [l], [5] and [37], and in detail in Ref. [3236].
1. Unidirectional Interpolation
Unidirectional interpolation means the interpolation is in one curvilinear coordinate direction only. In this section the cartesian coordinate vector will be shown as a function of the coordinate involved in the interpolation, as the unidirectional interpolation is fundamentally between points. These points can, however, lie on boundary (and perhaps interior) curves or surfaces, and in this sense the unidirectional interpolation can be considered to be between these curves or surfaces. Therefore the singlevariable functional relationship () used in this section can be considered to represent dependence on all coordinates, the interpolation points _{i} being functions of the coordinates along the boundary curves or surfaces.
A. Lagrange interpolation
The simplest type of unidirectional interpolation is Lagrange interpolation, which is based on polynomials. In the linear form we have, with ,

(1) 
Here _{1}=(0) and _{2}=(I), so that () is defined in terms of the two boundary values, _{1} and _{2}. The grid points are located at the successive integer values of from 0 to I. One family of grid lines will be straight lines connecting corresponding boundary points with this linear interpolation.
The general form is

(2) 
with _{n} = (_{n}), and the functions _{n} being polynomials defined on the entire interval such that

(3) 
In the linear case given above we have, with N=2,

From Eq. (2) and (3), 

so that the interpolation function matches at the N points = _{1}, _{2}, ..., _{N}=I:
The specified interior points, _{n} for n=2,3,...N1, are not necessarily grid points, since the grid points are defined by evaluating the interpolation formula at successive integer values of , but are simply additional parameters that serve to control the distribution. It is possible to specify the locations of certain interior grids points, however, by taking the _{n} corresponding to the specified _{n} to be the value of at the grid point of interest.
The Lagrange interpolation polynomials, defined to satisfy by Eq. (3), are in general

(4) 
The quadratic forms thus are, with N=3, and _{2} = I/2.
for which () is defined in terms of the two boundary values, _{1} and _{3}, and one interior value, _{2}. It should be noted that the purpose of the inclusion of the interior points in grid generation is control of the grid point distribution, not to increase the accuracy of the interpolation as is normally the case. There is, in fact, no question of accuracy of the interpolation here, since the aim is just to generate a grid from the boundary values of the coordinates.
B. Hermite interpolation
Lagrange interpolation matches only function values. It is possible to match both function, , and firstderivative, ' = _{ }, values using Hermite interpolation defined by

(5) 
where the Hermite inerpolation polynomials are defined on and satisfy the conditions
These polynomials can be obtained from the Lagrange interpolation polynomials by

(6a) 

(6b) 
where the prime here indicates differentiation of the polynomial with respect to the argument, . With N=2 we have
and the function matches the two boundary values, _{1} and _{2}, and the first derivatives, ^{'}_{1} and ^{'}_{2}, at the two boundaries.
Extensions of polynomial interpolation to match higherorder derivatives is obviously possible, the degree of the polynomial increasing with each additional condition or point to be matched. The polynomials of high degree exhibit considerable oscillation, however, so such procedures are not of great importance to grid generation. The general form again includes matches at interior points, which can be used to control the coordinate line spacing, since the first derivative, ' = _{ }, is a measure of the grid point spacing here, with being unity between points by construction. As with Lagrange interpolation, these specified interior points may or may not be grid points.
It is also possible, of course, to omit points from either of the summations in Eq. (5), so that and its first derivatives are not both matched at all points (deficient Hermite interpolation). Thus, with N=2 and the n=1 term omitted from the second summation, the two boundary values would be matched, but the first dervative at only the =I boundary would be matched. Clearly, the Hermite interpolation form, Eq. (5), could be equivalently defined in terms of the Lagrangian interpolation form, Eq. (2), with 2N points, since both are polynomial representations. Obviously either approach can be used to control the grid point spacing in the field.
The capability of specifying _{ }, as well as , can be used to make the grid orthogonal at the boundary. From Eq. (III33) the unit normal to a ^{i}coordinate surface is given by

Using Eq. (III10) this becomes 

The condition for orthogonally at the boundary then is that _{ i} be in the direction of the unit normal to the boundary:

(7) 
where is the spacing off the boundary to be specified. Since all the quantities with j and k subscripts can be evaluated from the points on the boundary, it remains only to specify the spacing, s^{i}, off the boundary and to use Eq. (7) for _{ i} on the boundary in the Hermite expressions.
C. Other forms of polynomial interpolation
As noted above, Hermite interpolation, which matches and _{ } at N points, can be equivalently constructed as an interpolant which matches at 2N points. Another form of expression of the polynomial interpolation uses the direct expression of the polynomial, so that

(8) 
Here we must have
where _{o} and _{I} are the boundaries. This form is not as straightforward as the Lagrange form for use in grid generation, since in the latter form certain grid point locations can be specified directly, while in the former the coefficients must be evaluated in terms of these specified points.
Still another form is that of Bezier, using Bernstein polynomials:
Thus the coefficients _{1} and _{N1} specify the slopes at the boundaries. An advantage of the Bezier form is that the coefficients define the vertices of an open polygon to which the curve is an approximation. Thus the general shape of the curve can be inferred by considering the coefficients to represent points in the field, with the lines from _{o} to _{1}, and from _{N1} to _{N}, defining the slopes at the two ends. The shape of the curve can then be designed by the placement of the vertices in the field as indicated below. Modifications of the curve can thus be made by adjusting the positions of these vertices.
Still another form can be defined using piecewise polynomials for the interpolation functions. Some degree of continuity must be lost in this case, of course. Continuity of the grid lines can be achieved using the piecewiselinear polynomials shown below (truncated versions apply at or near the end points):
while slope continuity can be gotten with the following piecewise polynomials:
Such piecewise polynomials allow a greater degree of local adjustment to be made, since the polynomial _{n} which multiplies the interpolation point _{n} vanishes except in the immediate vicinity of _{n}. By the conditions (3), any interpolation function _{n} must vanish at all the interpolation points except _{n}, but need not vanish between the points. Adding more interpolation points with global polynomials thus means increasing the degree of the polynomials, since the numbers of zeros must increase, and hence the polynomial becomes highly oscillitory.
D. Splines
The Lagrange and Hermite interpolation functions given above are completely continuous at all points. Complete continuity, however, may be attained at the price of oscillation. Both of these forms fit a single polynomial from one boundary to the other, matching specified values of the coordinates and perhaps the derivatives thereof (i.e., the point spacing). As more interior points are included, or as the first derivatives are included, the order of this global polynomial increases and thus oscillations become more likely. An alternative approach is to fit a loworder polynomial between each of the specified interior points, with continuity of as many derivatives as is possible enforced at the interior points. The interpolation function is then a piecewisecontinuous polynomial.
This type of interpolation function is called a spline and is formed as follows for the most common case of the cubic spline.
With a cubic polynomial fitted between points _{i} and _{i+1} we have a linear variation of the second derivative between these points and thus

(10) 
After two integrations and evaluation of the two constants ofintegration such that (_{i})=_{i} and (_{i+1})=_{i+1}, we have on _{i}_{i+1}

(11) 
Then, after differentiation and setting =_{i}, we have on _{i}_{i+1},

(12) 
Similar evaluation on the adjacent interval _{i1}_{i} gives

(13) 
Equating these two expressions in order to produce continuity of r' at the interior points, we have

(14) 
which is a tridiagonal equation for " at the interior points. It is necessary to set some conditions on " on the boundaries in order to solve this system, and the "natural" spline uses _{1}" = _{I}" = 0. This choice minimizes the total curvature, and thus the natural spline is the smoothest interpolant. This solution defines the _{i} in terms of the _{i}, so that substitution of these values for _{i}" into Eq. (11) then gives the spline in the general form of Eq. (2), except that the interpolation functions, _{n}, are, of course, different from the Lagrange interpolation polynomials. It should be recalled again that the interior points may or may not be grid points, the latter being defined by the interpolation formula evaluated at successive integer values of after the spline has been constructed over the entire field.
E. Tension splines
The spline tends to give a very smooth point distribution. Stronger localized curvature around the specified interior points can be obtained with the tension spline. Here Eq. (10) on _{i}_{i+1} is replaced by

(15) 
where ^{2} is a constant to be specified. (The tension spline tends progressively toward a linear function for large values of , and toward a cubic spline for small values.) Integration and evaluation of constants then yields, on _{i}_{i+1},

(16) 
The requirement of continuity of first derivatives at the interior points then yields the tridiagonal equation

(17) 
where _{i} = _{i+1}  _{i} and _{i1} = _{i}  _{i1}. Some application of tension splines are given in Ref. [33].
F. BSplines
One further possibility is to use piecewise continuous functions which satisfy the cardinality conditions by vanishing identically outside some interval around _{n}, as discussed in Section C above. This type of function allows the interpolation to be modified locally without affecting the interpolation function elsewhere. The Bsplines are an example of this approach.
From Eq. (14) a cubic spline which matches the function at N points, with continuity of secondderivatives, requires N+2 items of data, i.e., the N values of _{n} (n=1,2,...N) and the values of " at each boundary. Therefore a cubic spline which has = ' = " = 0 at each boundary can be defined over five points if is specified at only a single interior point (since N+2=7 data items can be specified here). If such a spline over five points is joined to the line = 0 outside these five points, we have a function which is nonzero only over four intervals and yet which has continuous second derivatives everywhere. Such a function is called a Bspline, denoted N_{4N}(), where the endpoints of the nonzero interval are _{n4} and _{n}. Similarly, quadratic, linear and constant Bsplines are nonzero over three, two, and one intervals, respectively, and are denoted N_{qn}, where q = 3, 2, and 1. The endpoints of the interval of nonzero values for these splines are _{nq} and _{n}. The specification of a single value in this interval is usually replaced by the specification of the integral over the interval so that

(18) 
The practical importance of Bsplines is that any spline of order q (the cubic spline is of order 4) can be expressed as a sum of multiples of Bsplines. Thus the cubic spline can be written as

(19) 
Since the Bsplines are nonzero only over four intervals, the modification of one coefficient here only affects the function over four intervals, thus allowing more localized control of the resulting grid.
The Bsplines can be calculated from the recurrence relation

(20) 
Thus N_{4n}() requires the successive calculation of N_{1,n1}, N_{2,n1}, N_{2,n}, N_{3,n1}, N_{3,n}, and finally N_{4,n}. The constant Bspline, NN_{1,n1}, used to start this calculation, is given on the interval _{n2}_{n1} by N_{1,n1}=1 and vanishes elsewhere.
For the point _{n} we have, in view of the vanishing of the Bsplines outside four intervals,

(21) 
which is a tridiagonal relation (N+1 equations) for the coefficients _{n}, _{o}=_{o} and _{N}=_{N}. Thus, even though the modification of a single coefficient only affects four intervals, the modification of an interpolation point requires a redetermination of all the coefficients and thus affects the function over the entire range.
The coefficients, _{n}, in the Bspline representation may be interpreted as the vertices of an open polygon, to which the curve is an approximation, as for the Bezier form discussed above. The slopes at the ends are defined by the directions _{1}_{o} and _{N}_{N1}. The curve passes close to the midpoint of each side, with the exception of the first and last sides. The curve also passes through the points (_{k1} + 4_{k} + _{k+1})/6 for k = 2,3,...,N2. These points are onethird of the way along the straight line joining _{k} to the midpoint of the line joining _{k1} and _{k+1}. Since the Bsplines are nonzero only on four intervals, the alteration of one vertex only affects the curve in its immediate vicinity. An application of Bsplines in grid generation is given in Ref. [39].
G. Multisurface interpolation
The multisurface method, discussed in Ref. [32][36], is also a unidirectional interpolation procedure. This procedure is constructed from an interpolation of a specified vector field, followed by vector normalizations at each interpolation point in order to cause a desired telescopic collapse so that the boundaries are matched. The specified vector field is defined from piecewiselinear curves determined by the boundaries and successive intermediate control surfaces. Normals to such surfaces are special cases. Polynomial interpolants for the vector field yield all of the classical polynomial cases along with a rational method for avoiding disasters such as can occur with direct Hermite interpolation with excessively large or discontinuous derivatives. Here the immediate surfaces are not coordinate surfaces, but are used only to define the vector field. These vectors are taken to be tangents to the coordinate lines intersecting the surfaces, so that integration of this vector field produces the position vector field for the grid points.
A collection of subroutines which automatically perform the necessary parts of grid construction using this multisurface procedure has been written and is described in Ref. [34]. Some of the automation features of this collection are applicable to other grid construction procedures as well. These subroutines can rotate and move curves, prospect one curve from another, normalize and parameterize curves, cluster points on a curve, and perform other such utilitarian functions to aid in the setup of an overall configuration.
In the multisurface interpolation we have

(22a) 
where 


(22b) 
and where the _{n} are specified points, with _{1}=(0) and _{N}=(I), on the boundary surfaces. (Recall the discussion at the beginning of this section, i.e., that the points can be considered to lie on curves or surfaces and thus the interpolation, while being fundamentally between points, can be considered to be between the surfaces on which those points lie.) Here the telescopic collapse for the series for =I matches the boundary at _{I}. The intermediate points here, _{2},_{3}....,_{N1}, are not grid points, but serve only to define the slopes _{ }, as given by Eq. (24) below. Eq. (22) is a polynomial if the functions _{n} are polynomials, but such is not required.
Since, by differentiation of Eq. (22),

(23) 
we have, for 0=_{1}<_{2}<....<_{n1}=I,

(24) 
if the functions _{n} satisfy the cardinality conditions

(25) 
The polynomials that satisfy these conditions are simply the Lagrange polynomials given by Eq. (4), here stated as

(26) 
Using Eq. (24), Eq. (23) can be written as

(27) 
This form thus is based on an interpolation of the firstderivatives _{ }, instead of , the interpolation expression for coming from an integration of the interpolation function for _{ }. Note, however, that this amounts to the specification of the slope _{ } at particular values of the curvilinear coordinate _{n}, and not at a specified position in space as is done in the Hermite form. It is clear from Eq. (24) that the intermediate points, _{n} for n=2,3,....,N2, serve to define the slopes _{ }(_{n}):
Because of the integration involved, the degree of the interpolation polynomial will be one greater than that of the functions _{n}.
Also with Eq. (24), the interpolation for , Eq. (22), can be written

(28) 
and thus is equivalent to a form of deficient Hermite interpolation. In implementation, however, it is the points _{n} that are specified, as in Eq. (22). Again it should be recalled that the point _{n}is not the grid point at _{n}, execept for the boundaries _{1}=0 and _{N}=I.
As has been noted, _{1} and _{n} are determined by the boundaries:

(29a) 
For N4, _{2} and _{N1} are determined by the intended values of _{ } at the boundaries through Eq. (24):

(29b) 
For N=3 only one of the above equations can be used, i.e., _{ } can be specified at either boundary but not at both. The use of the intermediate surfaces, instead of direct specification of the derivatives _{ } as in classical Hermite interpolation, provides a geometric interpretation that serves to help avoid the overlapping of grid lines that can occur if too large a value is given for _{ }.
Following Ref. [35], consider now for the _{n} the piecewiselinear functions diagramed below:
with the normalization
where _{n} = _{n+1}  _{n} so that each _{n} integrates to unity. These functions are given by, for n=2,3,....,N2,

(30) 
With these functions we have

(31) 
Note that G_{n}(I) = 1 here for all n and that G_{n}(_{1}) = 0,
These interpolation functions have the form
Now on the interval _{n}_{n+1} we have
G_{m}()=1 for m=1,2,....,n1 and G_{m}()=0 for m=n+2,n+3,....,N1. Therefore, on this interval
which, because of the telescopic collapse of the summation, reduces to, for _{n}_{n+1},

(32) 
Then 


(33a) 
and 


(33b) 
Also, from Eq. (32), for _{n}_{n+1}

(34) 
so that 


(35a) 
and 


(35b) 
We thus have on this interval
Thus on the interval _{n}_{n+1}, () is affected only by _{n}, _{n+1}, and _{n+2}. Conversely _{n} affects only the grid point locations on the interval _{n2}_{n+1}. Therefore local adjustments in the grid point locations can be made without affecting all of the points.
With the grid points located at unity increments of , so that I is the total number of points, we have, from Eq. (33), the grid points given by
The local control provided by these piecewise linear interpolants can be used to restrict undesirable mesh forms or to embed desirable ones within a global system with continuous first derivatives (cf. Ref. [34][35]). The second derivatives are, however, discontinuous. As examples, the propagation of boundary slope discontinuities can be arbitrarily restricted and general rectilinear Cartesian systems can be embedded to simplify problems over a large part of their domain.
In a further development (cf. Ref. [36]) the procedure is extended to use piecewise quadratic local interpolants, thus achieving continuity of second derivatives, with discontinuous third derivatives. The conceptual extension to higher order piecewise polynomial local interpolants, with consequent higher degree of continuity, is also discussed. Note that because of the integration of _{n}, the level of derivative continuity is always one greater than that of the piecewise polynomials.
H. Uniformity
It may be desirable for purposes of control of the grid point distribution to have a uniform distribution of the relative pro)ection of ()(0) along the straight line connecting the boundary points, i.e., (I)(0). This property has been called "uniformity" by Eiseman, cf. Ref. [32]  [36], and can be realized as follows: The unit vector along this straight line is
so that the relative pro)ection of ()(0) along this line is given by

(36) 
where 


(37) 
Uniformity is then achieved by choosing the interpolation parameters such that S() is linear. This does not completely determine all the interpolation parameters, however, so that some remain to be specified as desired. Uniformity is trivially assumed for linear interpolation of course.
For Lagrange interpolation we have, from Eq. (2),

(38) 
so that uniformity is achieved by selecting the _{n}, for n=2,3,....,N1, to cause all of the terms in S that are quadratic or higher to vanish. For Hermite interpolation, (5),

(39) 
For the multisurface interpolation defined by Eq. (22) we have

(40) 
For S() to be linear we must have


or 


(41) 
But, using Eq. (25), we then have 

(42) 
as the uniformity condition on the _{n}'s (cf. Ref. [35]). Both the _{n} (for n=2,3,....,N1) and the _{n} (for n=2,3,....,N2) are free to be chosen in order to satisfy the uniformity conditions (42). Thus a oneparameter family of cubic forms (N=4) results, a twoparametric family of quartic forms, etc. Substitution of Eq. (42) back into (41) yields a restriction on the choice of the functions _{n} since these must satisfy the relation

(43) 
Uniformity is particularly useful when the distribution function, such as those discussed in the next section, is used to redistribute the points on the grid lines set up by the interpolation (cf. Ref. [34][36]). Thus the interpolation is first applied with I=1 and with the uniformity conditions enforced. The final grid points then are placed according to the distribution function on the grid lines set up by the interpolation variable in place of the arc length, s, in the distribution function s().
I. Functions other than polynomials
The interpolation functions in the general forms given by Eq. (2), (5), and (22) do not have to be polynomials, and, in fact, if the variation in spacing over the field is large, other functions are better suited for grid generation. With N=2, Eq. (2) can be written in the form

(44) 
where can be any function such that (0)=0 and (1)=1. Here we have taken _{1}=1 and _{2}=. The linear polynomial case is obtained here with () = . The function in this form may contain parameters which can be determined so as to match the slope at the boundary, or to match interior points and slopes.
The interpolation function, , in this form is often referred to as a "stretching" function, and the most widely used function has been the exponential:

(45) 
where is a parameter that can be determined to match the slope at a boundary. Thus, since, from Eq. (44)

(46) 
we can determine from the equation

(47) 
with (_{ })_{1} specified.
As noted in Chapter V, the truncation error is strongly affected by the point distribution, and studies of distribution functions have been made in that regard. The exponential, while reasonable, is not the best choice when the variation of spacing is large, and polynomials are not suitable in this case. The better choices are the hyperbolic tangent and the hyperbolic sine. The hyperbolic sine gives a more uniform distribution in the immediate vicinity of the minimum spacing, and thus has less error in this region, but the hyperbolic tangent has the better overall distribution (cf. Section 3 of Chapter 5). These functions are implemented as follows (following Ref. [18]), with the spacing specified at either or both ends, or a point in the interior, of a point distribution on a curve.
Let arc length, s, vary from 0 to 1 as varies from 0 to I: s(0)=0, s(I)=1. Then let the spacing be specified at =0 and =I:

(48) 
The hyperbolic tangent distribution is then constructed as follows.
First,

(49) 

(50) 
Then the following nonlinear equation is solved for :

(51) 
The arc length distribution then is given by

(52) 
where 


(53) 
If this is applied to a straight line on which varies from _{0} to _{I} we have for the point locations:
The points are then located by taking integer values of :
Clearly the arc length distribution, s(), here is the function of Eq. (44).
With the spacing hs specified at only =0, the construction proceeds as follows. First B is calculated from

(55) 
and Eq. (51) is solved for . The arc length distribution then is given by

(56) 
With the spacing specified only at =I the procedure is the same, except that Eq. (56) is replaced by

(57) 
If the spacing s is specified at only an interior point s=, B is again calculated from Eq. (55), and then is determined as the solution of

(58) 
The value of at which s = is obtained by solving the nonlinear equation

(59) 
The arc length distribution then is given by

(60) 
This last distribution is based on the hyperbolic sine. From this a distribution baaed on the hyperbolic sine with the spacing specified at one end can be derived. Here B is evaluated from Eq. (55), and then is determined as the solution of

(61) 
The arc length distribution then is given by

(62) 
if the spacing is specified at =0. With the specification at =I, the distribution is

(63) 
It is also possible to construct a distribution based on the hyperbolic sine with specified spacing on each end. Here A and B are again calculated from Eq. (49) and (50), but is determined from

(64) 
The distribution is then given by Eq. (52), but with

(65) 
Finally, a procedure for incorporating the effect of curvature into the distribution function is given in Ref. [38], where the arc length distribution is given in the inverse form by

(66) 
where

(67) 
s_{t} is the total arc length, K() is the curvature, and A() is any distribution function (in the inverse form) that would be used without consideration of curvature.
2. MultiDirectional Interpolation
A. Transfinite interpolation
In two directions we may write a linear Lagrange interpolation function individually in each curvilinear direction:

(68a) 
and

(68b) 
This interpolation is now called "transfinite" since it matches the function on the entire boundary defined by =0 and =I in the first equation, or by =0 and =J in the second, i.e., at a nondenumerable number of points, cf. Ref. [40] and [41]).
The tensor product form

(69) 
where _{nm} = (_{n}, _{m}) matches the function at the four corners:
It does not, however, match the function on all the boundary.
The sum of Eq. (68a) and (68b),

(70) 
when evaluated on the =0 boundary gives

(71) 
This does not match the function on the =0 boundary because of the second term on the right, which is an interpolation between the ends of this boundary:
Similar effects occur on all the other boundaries, and the discrepancy on the =I boundary is
The discrepancies on both of these boundaries can be removed by subtracting from (,) a function formed by interpolating the discrepancies between the two boundaries:

(72) 
But this is simply the tensor product form given by Eq. (69), which matches the function at the four corners.
The function  then matches the function on all four sides of the boundary, so that we have the transfinite interpolation form,

(73) 
which matches the function on the entire boundary. By contrast, the tensor product form

(74) 
matches the function only at the four corners on the boundary. This generalizes to the interpolation from a set of N+M intersecting curves for which the univariate interpolation is given by

(75a) 
and

(75b) 
where now the "blending" functions, _{n} and _{m}, are any functions which satisfy the cardinality conditions

(76) 
The general form of the transfinite interpolation then is

(77) 
while the tensorproduct form is

(78) 
Eq. (77) can be written in the form
But here the first term is the result at each point in the field of the unidirectional interpolation in the direction, and the bracket is the difference between the specified values on the =_{n} lines and the result of the unidirectional interpolation on those lines. The twodirectional transfinite interpolation can thus be implemented in two unidirectional interpolation steps by first performing the unidirectional interpolation in one direction, say , over the entire field, calling the result _{1}(,):

(79a) 
then interopolating the discrepancy on the =_{n} lines over the entire field in the other direction, here, calling the result _{2}(,):

(79b) 
and then adding _{1} and _{2}.

(79c) 
The transfinite interpolation form given by Eq. (77) is the algebraically best approximation, while the tensor product from of Eq. (78) is the algebraically worst (cf. Ref. [40]). The difference between these two forms should be fully understood. The transfinite interpolation form, Eq. (77), interpolates to the entirety of a set of intersecting arbitrary curves, while the tensor product form, Eq. (78), interpolates only to the intersections of these curves. The interpolation function defined by Eq. (77) with N=M=2, using the Lagrange interpolation polynomials as the blending functions, is termed the transfinite bilinear interpolant. With N=M=3, this form is the transfinite biquainterpolant. Other immediate candidates for the blending functions are the Hermite interpolation polynomials and the splines, since these all can be expressed in the form of Eq. (75). The splineblended form gives the smoothest grid with continuous second derivatives.
B. Projectors
Now let P_{ }() be a onedimensional interpolation function in the direction which matches on the N lines, =_{n} (n=1,2,...N),:
(Note that the subscript here does not denote differentiation.) Similarly, let P_{ }() match on the M lines, =_{m} (m=1,2,...M). These interpolations are performed by projectors, P_{ } and P_{ }, which are assumed to be idempotent linear n operators. Protectors are discussed in more detail in Ref. [40]. Some discussion is also given in Ref. [37]. The product projector, P_{ }[P_{}()], then matches the function P_{ } (), instead of , on the N lines, =_{n}:
Then, since P_{ }() matches on the M lines, =_{m}, it follows that the product projector will match at the NxM points (_{n},_{m}):
Clearly the same conclusion is reached for the product projector P_{ }[P_{ }()], so that the protectors P_{ } and P_{ } commute.
The sum projector, P_{ } ()+P_{ }() matches +P_{ }() on the N lines =_{n}, and matches +P_{ }() on the M lines =_{m}. It should be clear then that the projector, P_{ } ()+P_{ }()P_{ }[P_{ }()] will match on the N lines =_{n}, since P_{ }[P_{ }()] matches P_{ }() on these lines. Similarly, the projector P_{ }()+P_{ }()P_{ }[P_{ }()] matches on the M lines =_{m}. Therefore, since P_{ } + P_{ } = P_{ }P_{ }, the Boolean sum projector, P_{ }P_{ } = P_{ }+P_{ }P_{ }P_{ }, will match on the entirety of the N+M lines =_{n} and =_{m} which includes, of course, the entire boundary of the region.
In summary, the individual protectors, P_{ } and P_{ }, interpolate undirectionally between two opposing boundaries:
The product progeotor, P_{ }P_{ }, interpolates in two directions from the four corners:
The Boolean sum projector, P_{ }P_{ }, interpolates from the entire boundary:
In three dimensions, the individual protectors, P_{ }, P_{ }, and P_{ }, interpolate undirectionally between two opposing faces of the sixsided region:
(matching on each of the two faces in each case). The double product projector, P_{ }P_{ }, interpolates in two directions from the four edges along which and are constant:
(matching on each of these edges). The Boolean sum projector

(80) 
interpolates in two directions from the four faces on which either or is constant:
(matching on all of these faces).
The Boolean sum projector

(81) 
interpolates in three directions, matching on the four edges on which and are constant and also on the two faces on which is constant:
The Boolean sum projector

(82) 
interpolates in three directions, with matched on all twelve edges:
The triple product projector, , interpolates from the eight corners:
Finally the Boolean sum projector

(83) 
matches on the entire boundary.
Much cancellation occurs in the algebraic manipulation of the projectors involved in developing the above relations, since P_{ }P_{ } = P_{ }, etc. Thus, for example,
This is to be expected since interpolation by P_{ } matches the function on all of the two sides on which is constant, while P_{ }P_{ } matches the function on the four edges on which and are constant. But these edges are contained on the two sides cited, so that nothing is changed by adding P_{ }P_{ } to P_{ } in the Boolean sense. The projector formed as the Boolean sum of all three of the individual projectors is algebraically maximum, while the triple product projector is algebraically minimal.
The importance of the projectors is that the structure given above allows multidirectional interpolation to be constructed systematically from unidirectional forms. With onedimensional interpolation of the form of Eq. (75) we have

(84a) 

(84b) 
so that

(85) 
which is just the tensor product form given previously in Eq. (78), so that the twodirectional transfinite interpolation corresponding to the projector P_{ }P_{ } is just that given by Eq. (77). As noted above, spline interpolation also falls directly into this form, so that the multidirectional transfinite interpolation based on splines requires only the determination of the splines separately in the individual directions.
Although Hermite interpolation can be defined in terms of additional points, and thus be put in this same form also, the use of projectors allows a more direct statement as follows. For the projectors we have, following Eq. (5),

(86a) 
and

(86b) 
Now

(87) 
Then the twodirectional transfinite interpolation can be constructed by substitution of Eq. (86) and (87) into the projector P_{ }P_{ }. Here the tensor product form, P_{ }P_{ }, interpolates from the values of the function, its two first derivatives, and the crossderivative at the four corners of the boundary. The transfinite interpolation form, P_{ }P_{ }, however, interpolates from the value of the function and its normal derivative on the entire boundary.
The triple product corresponding to Eq. (84) is simply

(88) 
Recall that with the unidirectional form given by Eq. (44), we have in these relations L=M=N=2 and

The above evaluations of the product projectors serve to illustrate the evaluation of such products for general projectors, i.e., that the effect of the product protectors is simply an interpolation in onedirection of an interpolant in another direction. This allows the multidirectional transfinite interpolation to be constructed from the Boolean sums of the protectors given above using any appropriate unidirectional interpolation forms as the basis projectors. It should also be noted that the unidirectional interpolation does not have to be of the same form in all the directions. Thus Lagrange interpolation could be used in one of the directions while Hermite is used in the other direction of a twodirectional construction. As noted above, the blending functions do not have to be polynomials. In fact, all of the unidirectional interpolation that was discussed earlier in this chapter can be applied in the context of multidirectional interpolation based on the protectors. This freedom to combine different types of univariate interpolation gives considerable flexibility to transfinite interpolation based on the projector structure, and allows attention to be focused on developing appropriate unidirectional interpolations, the multidirectional format then following automatically.
The protectors allow the transfinite interpolation to be easily set up as a sequence of unidirectional interpolations, in the manner discussed above. Thus in the twodirectional case, Eq. (80) can be written as

(89) 
where I indicates the identity operation. But here the first term is clearly the unidirectional interpolation in the direction, while the parenthesis (IP_{ }) is the discrepancy on the lines on which is specified that results from this interpolation. The second term then is the unidirectional interpolation of this descrepancy interpolated in the direction.
The twodirectional interpolation thus can be implemented by: (1) interpolating in the direction, (2) calculating the discrepancy between and this result on the lines that are to be used in the interpolation, (3) interpolating this discrepancy in the direction, and (4) adding the result of this interpolation to that of the interpolation. Symbolically these steps can be stated as the following:

(90) 
Obviously, the order of the unidirectional interpolation is immaterial.
Similarly, Eq. (83) can be written as

(91) 
The threedirectional interpolation thus can be implemented by (1) interpolating in the direction, (2) calculating the discrepancy between and this result on the surfaces and surfaces that are to be used in the interpolation in those directions, (3) interpolating this discrepancy by twodirectional interpolation, and (4) adding this result to that of the interpolation. These operations can be stated as

(92) 
Exercises
1. Show that with N=2 and _{n} constant, the multisurface interpolation is equivalent to the linear Lagrange interpolation. Mote that the Lagrange polynomials here satisfy Eq. (25) with . For other choices of _{2}, other quadratic polynomials result from Eq. (25), so that there exists a oneparameter family of cubic forms of the multisurface interpolation. Similarly, a twoparameter family of quadratic forms exists, etc.
2. Show that the quadratic form of the multisurface interpolation is given by
3. Show that the quadratic forms of the multisurface and Bezier interpolations are equivalent.
4. Show that with N=4 and _{n} the quadradic Lagrange interpolation polynomials given on p. 282, the interpolation functions for the multisurface interpolation are given by
the multisurface interpolation is equivalent to the cubic Hermite interpolation.
5. Show that S() is given by Eq. (38) for Lagrange interpolation. (Hint: If all the _{n} in Eq. (2) are the same, the interpolation must reproduce this value, hence the Lagrange interpolation polynomials satisfy
6. Show that for quadratic Lagrange interpolation, uniformity requires that _{2} be selected such that
Note that this does not completely determine _{2}.
7. Show that uniformity is achieved with cubic Hermite interpolation with with orthogonality at the boundaries if the spacings at the boundaries are given by
where is the unit normal to the boundary. This completely specifies the interpolation in this case. However, as noted in Exercise 1, _{2} is a free parameter.
8. Show that for multisurface interpolation with N=4 and uniformity is achieved with
9. Show also that with orthogonality at the boundary the result of Exercise 8 completely determines all of the interpolation parameters, i.e., that
where is the unit normal to the boundary. Hint: Use Eq. (7) and (29). For general _{2} the 1/6 is replaced by 1/2  1/[6(_{2}/I)] and 1/2  1/[6(1_{2}/I)] in the above expressions involving _{2} and _{3}, respectively. Some effects of the choice of _{2} are shown in Ref. [32].
10. Show that local uniformity on the interval _{n}_{n+1} for the multisurface interpolation based on piecewiselinear functions requires that
where , with
11. Consider a rectangular physical region with equallyspaced points on the bottom and top, but with unequal spacing on the left and right sides (but with the same point distribution on both of these sides). Show that horizontal interpolation will reflect the unequal spacing of the horizontal grid lines in the field, but that vertical interpolation will not. Show also that the unequal spacing is reflected with transfinite interpolation.
12. Show that transfinite interpolation based on linear blending functions will reflect the unequal boundary point spacing in the field for the rectangular physical region of Exercise 11, but will not for a Cgrid. From the consideration of transfinite interpolation as a sequence of unidirectional interpolations, explain why this is so.
13. Show that with cubic Lagrange interpolation the locations of the two intermediate surfaces, _{2} and _{3}, are related to the slopes at both ends. Note the contrast between this and the multisurface interpolation where each of the intermediate surfaces depends on only the slope at one end.
14. Give the cubic form of Lagrange interpolation.
15. Show that in two dimensions transfinite interpolation is equivalent to a generation system based on the fourth  order partial differential equation
(This is also equivalent to the quadralaterial isoparemetric elements often used to construct finite element meshes.)