~ 

NASA Technical Memorandum 10221 3 

A Direct Procedure for 

Interpolation on a Structured 

Curvilinear Two-Dimensional 

Grid 

David W. Zingg, University of Toronto, Downsview, Ontario, Canada 

Maurice Yarrow, Sterling Software, Arnes Research Center, Moffett Field, California 

July 1989 

. 

National Aeronautics and 

Space Administration 

Ames Research Center 

Moffett Field, California 94035 
Summary 

A direct procedure is presented for locally bicubic interpolation on a structured, curvilinear, two-dimensional grid. The physical (Cartesian) space is transformed to a comput ational space in which the grid is uniform and rectangular by a generalized curvilinear 

coordinate transformation. Required partial derivative information is obtained by finite 

differences in the computational space. The partial derivatives in physical space are determined by repeated application of the chain rule for partial differentiation. A bilinear 

transformation is used to analytically transform the individual quadrilateral cells in physical space into unit squares. The interpolation is performed within each unit square using 

a piecewise bicubic spline. 

1. Introduction 

Problems in computational fluid dynamics (CFD) are often solved using a numerically generated grid which is structured but nonuniform, consisting of quadrilat'erals in 

two dimensions and hexahedra in three dimensions. The grid is structured in that each 

cell contains an implied set of coordinate directions. A uniform rectangular grid in computational space is obtained froin the physical space grid by a generalized curvilinear 

coordinate transformation (ref. 1). This facilitates the application of finite differences and 

finite volume methods. 

In many applications of CFD, such as multiblock algorithms using overlapping grids, 

surface definition, and graphical representation of flow fields, data which are known at grid 

points must be interpolated at intermediate points. In the overlapping grid technique described by Benek et al. (ref. 2), interpolation is accomplished using an iterative technique. 

The purpose of the present paper is to present an alternative approach which involves no 

iteration. 

The interpolation procedure presented uses a piecewise bicubic interpolant rather 

than the simpler bilinear interpolant. The bilinear interpolant has the often desirable 

property that it introduces no new extrema. However, the bicubic interpolant is inore 

accurate and produces a smooth interpolant. Furthermore, bicubic interpolation is capable 

of preserving derivative information generated by a finite-difference flow solver. Iniportant, 

considerations regarding monotonicity and conservative properties of the interpolant are 

not addressed in this paper. Examples of alternative techniques for bivariate interpolation 

on a rectangular grid are given by Carlson and Fritsch (ref. 3) and Schiess (ref. 4). 

The generalized curvilinear coordinate transformation is given in Sect,ion 2, and the 

iterative interpolation technique of Benek et al. is reviewed in Section 3. The direct, 

interpolation procedure is presented in Section 4. The Appendix contains the formulas for 

bicubic spline interpolation on a unit square. 

1 
I 

I 2. Generalized Curvilinear Coordinate Transformation 

The generalized curvilinear coordinate transformation maps a nonuniform, nonrectangular grid in physical (Cartesian) space to a uniform, rectangular computational grid. 

The Cartesian coordinates (x, y) are mapped to the curvilinear coordinates (t,~) by the 

following general transformation: 

where we consider two space dimensions only and ignore time dependence for simplicity. 

rectangular, and of unit length, i.e., 

, The mapping is chosen such that the resulting grid in computational space is uniform, 

Hence, finite-difference representations of tlt and a, are easily formulated. The mapping 

is one-to-one except at topologic,al singularities or cuts, where one physical point may be 

mapped to many computational points. 

The chain rule for partial differentiation gives, in matrix form, 

, The coordinate transformation is generally not known analytically and hence must be 

determined numerically. When the roles of the coordinate systems in equation (3) are 

reversed, t8he chain rule gives 

I 

Therefore, we must have 

=J[ ys 

-29 

All of the t,erins involving tlt and 

matrix form of the inetric relations. 

are evaluated as finite differences. Eyua,tion (6) is t81ie 
3. Iterative Interpolation Procedure 

The iterative scheme developed by Benek et al. (ref. 2) is applicable to threediiiiensional generalized curvilinear coordinates. Trilinear interpolation is used. Yarrow 

and Mehta (report to be published as a NASA Technical Memorandum) ut,ilize tricubic 

interpolants within tlie same iterative schenie. In this section, the iterative procedure is 

presented for a two-dimensional space. 

Since the grid is uniform and rectangular in computational space, the data can be 

readily interpolated as a function of tlie generalized coordinates by conventional niultivariate interpolation techniques. The interpolant for the function f can be written as 

The iterative procedure is used to deternine the generalized coordinates ((* , q*) which 

correspond to the physical point (z*,y*) at which tlie interpolation is to be performed. 

Since each point in computational space is mapped to only one physical point, we can form 

tlie following interpolants for the Cartesian coordinates: 

For a two-dimensional space, bivariate Newton-Raphson iteration is used to solve for 

((*, q*) such that 

The value of the function f at the point (z*,y*) is then determined from 

4. Direct Interpolation Procedure 

The direct interpolation procedure proceeds in three steps. First, required partial 

derivative information in physical space is determined using finite-difference approximations to the partial derivatives in coinputational space. The individual cells are then 

analytically transformed into unit squares by a bilinear transformation. Finally, a bicubic 

spline interpolant is determined within each unit square. 

For a function f known at the grid nodes, we require fz, fy , fzz, fzy , fyy . Using finite 

differences in computational space, we can determine ft , fv, ftt, f,tv, fvv. Applying the 

chain rule to equation (3) yields 

3 
Equations (3) and (11) may be assembled as 

where 

B= 

Similarly, 

8& = CdZY, 

where 

C= 

Therefore B = C-I. 

In block matrix form, B can be written as 

where 

Si lid arl y , 

4 
where 

Now B = C-' gives 

and 

B1 = L'I-' (the metric relations), 

B3 = C'3-l, 

B2 = -C'3-1C2C1-1. 

Substituting equations (16) and (17) into equation (18) gives 

Finally, use of the metric relations gives 

2 2 

XCC - 2Yll YCXCl, + Y1,YC2X,, - X,Y? YC€ + 2X,YllY€YCll - "CYt2Y?7?ll 13 

txx = -- [Yll 

72, = -- [-YtYll XCt + 2YC2Y, - Yt3X,, + XCY?l2YCC - 2XCYVYCYCll + XCYC Y,,] 

J3 

1 

J3 

1 

2 2 

t "y = -- J3 [-Y?12x,xCt + YdY?lrt + XllYC)XC, - YllYCXCX,, + Xl12Y1Y€€ 

%Y = -J3 [YtYllXllW - YdYllXC + XllYC)XC, + YC2XCX,11 - "€YvXC17YC€ 

tYY = -J3 [Y,XV2XC€ - 2YllXllXtXCll + YllXt2Xll, - Xl13YCt + 2XV2XtY€7l - X1,XC2Y?,] 

17YY = -- [-YCxa2x€€ + 2YCX,XtX€ll - YCXt2X,1, + ~C"Tl2Y€t - 2X,XC2YCll + XC3Yll?,] * 

- Xll(Y?lXC + XllYdYCll + ~,YC~€Ylllll 

+ XC(YllXt + XllYdY€ll - "€3Y€Ys,J 

1 

1 

1 

J3 

(20) 

Therefore, all of the elements of matrix B can be expressed in terms of the elements 

of matrix C. As a result, dZyf can be determined given atsf. The elements of C' and 

can be approximated using finite differences. An alternative procedure is to calculate fz 

and fy given f~ and f, from equation (3) and then to use finite differences to determine 

(.fz)t, (fy)t, (fZ),, (fv)ll. The required values of fzz, fxy, and fyy can then be calculated 

using equation (3) again. The first method involves more computer storage but it also 

involves a smaller stencil. 

Given fz, fy , fxz, fzy, fyy at the corners of a given cell, the problem remains to interpolate f on an arbitrary quadrilateral. We now consider the new two-dimensional space 

given by (p,q) shown in figure 1. The quadrilateral in physical space is related to a unit 

5 
square in (p, q) space by a bilinear mapping as follows. (Without loss of generality, we will 

assume the lower-left-hand corner to be at (O,O).) 

x = ap + bq + cpq 

Y = dP + eq + fpq 

Note that a different mapping is used for each quadrilateral, or cell, in the grid. This contrasts with the generalized curvilinear coordinate transformation in which a single smooth 

iiumerically generated mapping is applied to the entire grid. Consequently, analytical 

mappings can be utilized. 

Figure 1. - Two-dimensional space given by (p,q). 

The coefficients of the bilinear mapping are given by 

a = x4, b = 22, c = x3 - 24 - 22, 

d=y4, e=y2, f =y3-y4-y2* 

The required derivatives in (p, q) space are given (as in equation (12)) by 

ap = xpax + Ypdy 

aq = xqax + Yqdy 

a,, = xpqax + Ypqdy + xpxqaxx + (XpYq + x9Yp)dxy + YpYq d YY, 

where 

- and we have used xpp = xqg - ypp = yqq = 0. 

Froiii equation (21a), we have 

Given a location in physical space (20, yo), we must find the corresponding (PO, yo ). 

G 
Substituting this into equation (2lb) as follows: 

Yo = duo + (e + fpo)( 2o - UP0 ) b + cpo 

gives the following quadratic for po 

(25) 

po2(cd - uf) + po(-cyo + bd + .of - ae) + (-yob + exo) = 0. (26) 

The solution of equation (26) gives two values of po. The corresponding values of qo are 

determined from equation (24). One location (p0,qO) will not lie within the unit square 

and hence can be discarded. 

Using fz, fy , fzz , fry, fyy as calculated from equation (12)) we can determine fp, fq, 

and fpq from equat,ion (23) at each corner of the unit square in (p,q) space. We then 

use the standard formula for bicubic interpolation within a unit square (see Appendix). 

The bicubic is finally evaluated at (po,qo), as calculated from equations (26) and (24)) 

respectively. 

5. Conclusions 

A procedure has been presented for interpolation on a structured, curvilinear, twodimensional grid. In contrast to existing methods (ref. 2), the procedure avoids the need for 

iteration. A piecewise bicubic spline is used, leading to a smooth and accurate iiiterpolaiit 

which is capable of preserving derivative information generated by a flow solver based on 

finite differences. Potential applications include multiblock algorithms using overlapping 

grids in the solution of problems in computational fluid dynanlics, surface definition, and 

graphical representation of data. Extension to three dimensions is straightforward but 

nontrivial . 

7 
Appendix 

Locally Bicubic Spline Interpolation on a Unit Square 

In this appendix, we give a piecewise bicubic spline forinulation for interpolation 

within a unit square. De Boor (ref. 5) presents a technique for determining the required 

partial derivatives at the corners of each cell in a rectangular grid such that G2 continuity 

is obtained for the domain. In the formulation used here, the required partial derivatives 

at the corners of the unit square are given by centered difference formulas, leading to G1 

continuity, analogous to the univariate cubic Bessel spline given by de Boor (ref. 6). This 

formulation has the advantage of producing a local interpolant. 

The bicubic spline interpolant on a unit (p,q) square can be written as 

where 

bmn] = AKAt, 

1000 

A= [+ ;2 10 3 :1] 

-2 1 

and 

8 
Refer e nc es 

1. Prilliaiii, T. H.: Efficient, Solution Methods for the Navier-St,okes Equations. Lecture 

Notes for the Von Kariiian Institute for Fluid Dynaiilics Lecture Series: Numerical 

Techniques for Viscous Flow Computation in Turbomachinery Bladings, Brussels, 

Belgium, Jan. 1986. 

2. Benek, J. A.; Buning, P. G.; and Steger, J. L.: A 3-D Chimera Grid Embedding 

Technique. AIAA Paper 85-1523, June 1985. 

3. Carlson, R. E.; and Fritsch, F. N.: Monotone Piecewise Bicubic Interpolation. SIAM 

J. Numer. Anal., vol. 22, no. 2, April 1985. 

4. Schiess, J. R.: Two Algorithms for Rational Spline Interpolation of Surfaces. NASA 

TP 2536, 1986. 

5. De Boor, C.: Bicubic Spline Interpolation. J. Math. and Phys., vol. XLI, no. 3, 

Sept. 1962. 

6. De Boor, C.: A Practical Guide to Splines. Springer-Verlag, New York, 1978. 

c 

9 
Report Documentation Page 

1. Report No. 

NASA TM-102213 

2. Government Accession No. 

A Direct Procedure for Interpolation on a Structured Curvilinear 

Two-Dimensional Grid 

17. Key Words (Suggested by Author(s)) 

Interpolation 

Bicubic interpolation 

Generalized curvilinear coordinates 

7. Author@) 

David W. Zingg* and Maurice Yarrow** 

18. Distribution Statement 

Unclassified - Unlimited 

Subject Category - 64 

9. Performing Organization Name and Address 

Ames Research Center 

Moffett Field, CA 94035 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

3. Recipient's Catalog No. 

5. Report Date 

July 1989 

6. Performing Organization Code 

8. Performing Organization Report No. 

A-89210 

10. Work Unit No. 

505-60 

1 1. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 

15. Supplementaty Notes 

Point of Contact: Harvard Lomax, Ames Research Center, MS 202A-1, Moffett Field, CA 94035 

(4 15) 694-5 124 or FTS 464-5 124 

*University of Toronto, Downsview, Ontario, Canada M3H 5T6 

**Sterling Software, Ames Research Center, Moffett Field, CA 94035 

16. Abstract 

A direct procedure is presented for locally bicubic interpolation on a structured, curvilinear, twodimensional grid. The physical (Cartesian) space is transformed to a computational space in which the grid 

is uniform and rectangular by a generalized curvilinear coordinate transformation. Required partial 

derivative information is obtained by finite differences in the computational space. The partial derivatives 

in physical space are determined by repeated application of the chain rule for partial differentiation. A 

bilinear transformation is used to analytically transform the individual quadrilateral cells in physical space 

into unit squares. The interpolation is performed within each unit square using apiecewise bicubic spline. 

Structured grids 

19. Security Classif. (of this report) I 20. Security Classif. (of this page) I 21. NO. of Pages I 22. Price 

A02 I lo I Unclassified I Unclassified I 

ASA FORM 1626 OCT 86 

For sale by the National Technical Information Service, Springfield, Virginia 221 61