Numerical Grid Generation
Foundations and Applications
##### By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin

APPENDIX C

CODE DEVELOPMENT AND COMPUTER EXERCISES

1. Code Development Exercises

1. Take the computational region to be a rectangle on which the curvilinear coordinates are defined to be equal to the indices of the field arrays, i.e., = I and = J, with x = X (I,J) and y = Y (I,J).

Make provision for reading in values of x and y on any segments of the boundary of the computational region. Generate x and y in the interior by interpolating linearly between the top and bottom boundaries. Plot the grid.

2. Modify the code to allow horizontal (in the computational plane) interpolation, as well, the choice being specified by input.

3. Now add the choice of interpolation from the four corners (tensor product interpolation).

4. Finally add the choice of transfinite interpolation.

5. Generalize the interpolation to cubic Hermite interpolation, with the grid being orthogonal at the boundary.

6. Generalize the interpolation to use the hyperbolic tangent distribution function, rather than being linear, with specified relative spacing on each end.

7. Modify the code to provide for reading in x and y on any segment of any horizontal or vertical line in the computational region. Also provide for the interpolation tion to be done on any rectangular segment of the computational region (including a segment that is only a line.)

8. Add another field array TYPE (I,J) which is a flag to identify each point as one for which the x,y values are (1) fixed, e.g., specified points on the physical boundary, (2) out of the computation, e.g., points inside a slab, or (3) to be generated. Provide for the designation (1) and (2) to be made by input for any rectangular segment of the computational region, the default being to the designation (3).

9. Modify the dimensions of the field arrays so that an extra layer of points surrounds the computational region. Also add two more field arrays, ILINK (I,J) and JLINK (I,J). Provide for any segments of any horizontal or vertical lines to be designated as image points in TYPE by input, i.e., points for which the values of x and y are set equal to those at some other point. Also provide for the indices of these other points to be put in ILINK and JLINK by input.

10. Add an elliptic generator, based on Laplace equation, to the code. Use the algebraic generator (the interpolation) to provide the initial guess for point SOR iteration.

11. Add control functions to the elliptic generator. Let the control function be evaluated on the boundaries and interpolated into the field by transfinite interpolation.

2. Computer Exercises

1. Generate an algebraic grid between two concentric circles. Use linear interpolation between the circles.

2. Generate an algebraic grid between two ellipses, both of which are centered at the origin but which may have different eccentricities, using interpolation between the ellipses. Compare grids generated using linear and Hermite interpolation, the latter being orthogonal at the boundaries.

3. Generate a C-type algebraic grid for an ellipse inside an outer boundary formed by a semicircle replacing one side of a rectangle:

Compare (1) vertical interpolation in the computational region boundary, (2) horizontal interpolation, (3) tensor product interpolation, and (4) transfinite interpolation, using linear interpolation in each case. Note that (2) and (3) are totally unreasonable.

4. Generate an algebraic grid for a circular simply-connected region by (1) unidirectional interpolation, (2) tensor product interpolation, and (3) transfinite interpolation. Note that here only (3)gives a reasonable grid. Compare linear and Hermite interpolation for (3).

5. Repeat Exercise 4 with a triangular boundary.

6. Using the boundary configuration of Exercise 3, but with a hyperbolic tangent point distribution on the right-hand boundary of the physical region with smaller spacing at the centerline than at the top and bottom. Compare algebraic grids generated using (1) linear interpolation between the inner and outer boundaries, (2) nonlinear interpolation, based on the hyperbolic tangent, between the inner and outer boundaries, (3) transfinite interpolation with linear blending functions, and (4) transfinite interpolation using the boundary point distribution (in terms of relative arc length) as the blending functions. Note that only (2) and (4) preserve the boundary point distribution in the field.

7. Generate an algebraic grid for a square inside a rectangle using linear interpolation between the inner and outer boundaries. Note the propagation of the boundary slope discontinuities into the field. Generate a grid from an elliptic generation system for the same boundary point distribution and note the difference.

8. Generate an algebraic grid for a square inside a circle using linear interpolation between the inner and outer boundaries. Show that it is possible to position the points on the circle such that the grid overlaps the corners of the square. Generate a grid from an elliptic generation system for the same boundary point distribution and note the difference.

3. Listing of Routine for Computer Exercises

```      SUBROUTINE INTERP
PARAMETER(NI=20,NJ=20,N=5)
COMMON/COORD/X(NI,NJ),Y(NI,NJ)
COMMON/CONST/CHOICE,IMAX,JMAX,NA,DS1,DS2
COMMON/ATTR/IAL(N),IAX(N),IAY(N),JAL(M),JAX(N),JAY(N)
COMMON/COEF/AI(N),BI(N),CI(N),DI(N),AJ(N),BJ(N)
COMMON/COEF/CJ(N),DJ(N)
DIMENSION P(NI,NJ),Q(NI,NJ),XX(0:NI,O:NJ),YY(0:MI,O:NJ)
DIMENSION X1(NI,NJ),X2(NI,NJ),Y1(NI,NJ),Y2(NI,NJ)
INTEGER CHOICE
C
C     BOUNDARY INTERPOLATION
C
C      X          X ARRAY OF XI-ETA COORDINATE
C      Y          Y ARRAY OF XI-ETA COORDINATE
C      IMAX       MAX. NUMBER OF GRID IN XI AXIS
C      JMAX       MAX. NUMBER OF GRID IN ETA AXIS
C      NA         MAX. NUMBER OF ATTRACTIONS
C      DSl        SPECIFIED LENGTH OF INITIAL INTERVAL
C      DS2        SPECIFIED LENGTH OF FINAL INTERVAL
C      ATTR       ARRAY OF ATTRACTION TO LINES/POINTS
C      COEF       ARRAY OF COEFFICIENT FOR ATTRACTION
C
C      CHOICE
C         1       VERTICAL INTERPOLATION
C         2       HORIZONTAL INTERPOLATION
C         3       TENSOR PRODUCT INTERPOLATION
C         4       TRANSFINITE INTERPOLATION
C         5       HERMITE CUBIC INTERPOLATION
C         6       HYPERBOLIC TANGENT INTERPOLATION
C         7       ELLIPTIC GRID GENERATION ( SOR ITERATION )
C         8       ATTRACTION TO COORDINATES
C
IF(CHOICE.EQ.1) GO TO 100
IF(CHOICE.EQ.2) GO TO 200
IF(CHOICE.EQ.3) GO TO 300
IF(CHOICE.EQ.H) GO TO 400
IF(CHOICE.EQ.5) GO TO 500
IF(CHOICE.EQ.6) GO TO 600
IF(CHOICE.EQ.7) GO TO 700
IF(CHOICE.EQ.S) GO TO 800
C
C      **** VERTICAL INTERPOLATION ****
C
100   DO 110 I 1,IMAX
DO 110 J 1,JMAX
RJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)
RJ2=FLOAT(J-1)/FLOAT(JMAX-1)
C      *** ( EQ. 8-1 )
X(I,J)=RJ1*X(I,l)+RJ2*X(I,JMAX)
110   Y(I,J) RJl*Y(I.1)+RJ2*Y(I,JMAX)
RETURN
C
C      **** HORIZONTAL INTERPOLATION ****
C
200   DO 210 I=l,JMAX
DO 210 J=1,IMAX
RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)
RI2=FLOAT(I-1)/FLOAT(IMAX-1)
C      *** ( EQ. 8-1 )
X(I,J)=RI1*X(1,J)+RI2*X(IMAX,J)
210   Y(I,J)=RI1*Y(1,J)+RI2*Y(IMAX,J)
RETURN
C
C      **** TENSOR PRODUCT INTERPOLATION ****
C
300   DO 310 I=1,IMAX
DO 310 J=l,JMAX
RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)
RI2=FLOAT(I-1)/FLOAT(IMAX-1)
RJl=FLOAT(JMAX-J)/FLOAT(JMAX-1)
RJ2=FLOAT(J-1)/FLOAT(JMAX-1)
C      *** ( EQ. 8-69 )
X(I,J)=RI1*RJ1*X(1,1)+RI1*RJ2*X(1,JMAX)
*+RI2*RJ1*X(IMAX,1)+RI2*RJ2*X(IMAX,JMAX)
Y(I,J)=RI1*RJ1*Y(1,1)+RI1*RJ2*Y(1,JMAX)
*+RI2*RJ1*Y(IMAX,1)+RI2*RJ2*Y(IMAX,JMAX)
310   CONTINUE
RETURN
C
C      **** TRANSFINITE INTERPOLATION ****
C
400   DO 410 I=l,IMAX
DO 410 J=1,JMAX
RI1=FLOAT(I-1)/FLOAT(IMAX-1)
RI2=FLOAT(IMAX-I)/FLOAT(IMAX-1)
X1(I,J)=RI1*X(IMAX,J)+RI2*X(1,J)
410   Yl(I,J)=RI1*Y(IMAX,J)+RI2*Y(1,J)
DO 420 I=1,IMAX
DO 420 J=1,JMAX
RJ1=FLOAT(J-1)/FLOAT(JMAX-1)
RJ2=FLOAT(JMAX-J)/FLOAT(JMAX-1)
X2(I,J)=RJ1*(X(I,JMAX)-X1(I,JMAX))+RJ2*(X(I,1)-X1(I,1))
420   Y2(I,J)=RJ1*(Y(I,JMAX)-Y1(I,JMAX))+RJ2*(Y(I,l)-Y1(I,1))
C      *** ( EQ. 8-73 )
DO 430 I=1,IMAX
DO 430 J=1,JMAX
X(I,J)=X1(I,J)+X2(I,J)
430   Y(I,J)=Y1(I,J)+Y2(I,J)
IF(CHOICE.NE.4) GO TO 740
RETURN
C
C      ***  HERMITE CUBIC INTERPOLATION (ORTHOGONAL BOUNDARY) ***
C
500   DO 510 I=1,IMAX
DO 510 J=1,JMAX
XX(I,J)=X(I,J)
510   YY(I,J)=Y(I,J)
DO 520 J=1,JMAX
XX(O,J)=XX(IMAX-1,J)
YY(O,J)=YY(INAX-1,J)
XX(IMAX+1,J)=XX(2,J)
520   YY(IMAX+1,J)=YY(2,J)
DO 530 I=1,IMAX
DO 530 J=1,JMAX
RJJ=FLOAT(J-1)/FLOAT(JMAX-1)
C      *** ( EQ. 8-6 a and b, n=2 )
PHI1=(1.+2.*RJJ)*(1.-RJJ)*(1.-RJJ)
PHI2=(3.-2.*RJJ)*RJJ*RJJ
PSI1=(1.-RJJ)*(1.-RJJ)*RJJ
PSI2=(RJJ-1.)*RJJ*RJJ
C
C      ** CAL. NORMAL DERIV. **
C
XXI1=.5*(XX(I+1,1)-XX(I-1,1))
XXI2=.5*(XX(I+1,JMAX)-XX(I-1,JMAX))
YXI1=.5*(YY(I+1,1)-YY(I-1,1))
YXI2=.5*(YY(I+1,JMAX)-YY(I-1,JMAX))
UNIT1=SQRT(XXI1*XXI1+YXI1*YXI1)
UNIT2=SQRT(XXI2*XXI2+YXI2*YXI2)
C      *** ( EQ. 3-108 )
XNl=-YXI1/UNIT1*DSl
XN2=-YXI2/UNIT2*DS2
YN1=XXI1/UNIT1*DSl
YN2=XXI2/UNIT2*DS2
C      *** ( EQ. 8-5 )
XX(I,J)=PHI1*XX(I,l)+PHI2*XX(I,JMAX)+PSI1*XN1+PSI2*XN2
530   YY(I,J)=PHI1*YY(I,1)+PHI2*YY(I,JMAX)+PSI1*YN1+PSI2*YN2
DO 540 I=1,IMAX
DO 540 J=1,JMAX
X(I,J)=XX(I,J)
540   Y(I,J)=YY(I,J)
RETURN
C
C      **** HYPERBOLIC TANGENT SPACING INTERPOLATION ****
C
600   TOL=1.0E-10
C      *** ( EQ. 8-49, 50 and 51
A=SQRT(DS2/DS1)
B=1./(FLOAT(JMAX-1)*SQRT(DS1*DS2))
C      *** INITIAL GUESS BY SERIES EXPANSION
DELTA=SQRT(6.*(B-1.))
DO 610 IT=1,20
RESID=SINH(DELTA)/(DELTA*B)-1.
IF(ABS(RESID)LT.TOL) GO TO 630
610   CALL AITKEN(DELTA,RESID,DELTO,RO,RSO)
PRINT 620, RESID,DELTA,IT-1
620   FORMAT(//, 5X, 'DELTA IS NOT CONVERGE ?', 5X, 2E15.5,
*5X, I3, //)
GO TO 660
630   CONTINUE
C      *** ( EQ. 8-52, 53 and 54 )
DO 650 I=1,IMAX
DO 650 J=2,JMAX-1
RATIO FLOAT(J-1)/FLOAT(JMAX-1)
U=.5*(1.+TANH(DELTA"(RATIO-.5))/TANH(.5*DELTA))
S=U/(A+(1.-A)*U)
X(I,J)=X(I,1)+(X(I,JMAX)-X(I,1))*S
650   Y(I,J)=X(I,1)+(Y(I,JMAX)-Y(I,1))*S
660   RETURN
C
C      **** ELLIPTIC GRID GENERATION ( SOR ITERATION ) ****
C      *** CAL. P AND Q ON THE BOUNDARY **
C
700   DO 710 I=1,IMAX
DO 710 J=1,JMAX
XXI=.5*(X(I+1,J)-X(I-1,J))
XXIXI=X(I+1,J)-2.*X(I,J)+X(I-1,J)
XETA=.5*(X(I,J+1)-X(I,J-1))
XETA2=X(I,J+1)-2.*X(I,J)+X(I,J-1)
YXI=.5*(X(I+1,J)-Y(I-1,J))
XXIXI=Y(I+1,J)-2.*X(I,J)+Y(I-1,J)
YETA=.5*(Y(I,J+1)-Y(I,J-1))
YETA2-Y(I,J+1)-2.*Y(I,J)+Y(I,J-1)
IF(ABS(XETA2).LT.10E-3) XETA2=0.
IF(ABS(YETA2).LT.10E-3) YETA2=0.
RXI2=XXI*XXI+YXI*YXI
RETA2=XETA*XETA+YETA*YETA
C      *** ( EQ. 8-70 )
P(I,J)=(XXI*XXIXI+XXI*YXIXI)/RXI2
Q(I,J)=(XETA*XETA2+YETA*YETA2)/RETA2
710   CONTINUE
C
C      **    INTERPOLATE P AND Q BETWEEN BOUNDARY **
C            P : VERTICAL, Q : HORIZONTAL
C
DO 720 I=1,IMAX
DO 720 J=l,JMAX
RJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)
RJ2=FLOAT(J-1)/FLOAT(JMAX-1)
RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)
RI2=FLOAT(I-1)/FLOAT(IMAX-1)
P(I,J)=RJl*P(I,l)+RJ2*P(I,JMAX)
Q(I,J)=RI1*Q(1,J)+RI2*Q(IMAX,J)
720   CONTINUE
C
C      **    INITIAL GUESS WITH TRANSFINITE INTERPOLATION **
C
GO TO 400
740   CONTINUE
C
C      *** ITERATION ( SOR ) ***
C
ITMAX=200
TOL=10.E-5
W=1.8
DO 760 IT=1,ITMAX
ERRX=0.
ERRY=0.
DO 750 J=2,JMAX-1
DO 750 I=2,IMAX-1
XXI=.5*(X(I+1,J)-X(I-1,J))
YXI=.5*(Y(I+1,J)-Y(I-1,J))
XXIXI=X(I+1,J)+X(I-1,J)
YXIXI=Y(I+1,J)+Y(I-1,J)
XETA=.5*(X(I,J+1)-X(I,J-1))
YETA .5*(Y(I,J+1)-Y(I,J-1))
XXIETA=.25*(X(I+1,J+1)-X(I+1,J-1)-X(I-1,J+1)+X(I-1,J-l))
YXIETA=.25*(Y(I+1,J+1)-Y(I+1,J-1)-Y(I-1,J+1)+Y(I-1,J-1))
XETA2=X(I,J+1)+X(I,J-1)
YETA2=Y(I,J+1)+Y(I,J-1)
C      *** ( EQ. 6-18 and 6-20 )
G11=XXI*XXI+YXI*YXI
G22=XETA*XETA+YETA*YETA
G12=XXI*XETA+YXI*YETA
XTEMP=.5*(G22*(P(I,J)*XXI+XXIXI)+G11*(Q(I,J)*XETA+XETA2)
*-2.*G12*XXIETA)/(G11+G22)
YTENP=.5*(G22*(P(I,J)*YXI+YXIXI)+G11*(Q(I,J)*YETA+YETA2)
*-2.*Gf2*YXIETA)/(Gll+G22)
XTENP=W*XTEMP+(1.-W)*X(I,J)
YTEMP=W*YTEMP+(1.-W)*Y(I,J)
ERRX=AMAXO(ERRX,ABS(XTEMP-X(I,J)))
ERRY=AMAXO(ERRY,ABS(YTEMP-Y(I,J)))
X(I,J)=XTEMP
Y(I,J)=YTEMP
750   CONTINUE
IF(ERRX.LT.TOL.AND.ERRY.LT.TOL) GO TO 78O
760   CONTINUE
PRINT 770,ERRX,ERRY,IT-1
770   FORMAT(//, 5X, 'X AND Y ARE NOT CONVERGE ?', 2E15.5,
*5X, I5, //)
780   CONTINUE
IF(CHOICE.EQ.8) GO TO 830
RETURN
C
C      **** ATTRACTION TO COORDINATE LINE/POINT ****
C
800   DO 810 I=1,IMAX
DO 810 J=1,JMAX
P(I,J)=0.
810   Q(I,J)=0.
DO 820 NS=1,NA
DO 820 I=1,IMAX
DO 820 J=1,JMAX
XL=FLOAT(I-IAL(NS))
XI=FLOAT(I-IAX(NS))
XJ=FLOAT(J-IAY(NS))
YL=FLOAT(J-JAL(NS))
YI=FLOAT(I-JAX(NS))
YJ=FLOAT(J-JAY(NS))
C      *** ( EQ. 6-30 )
P(I,J)=P(I,J)-AI(NS)*(XL/ABS(XL))*EXP(-CI(NS)*ABS(XL))
*-BI(NS)*(XI/ABS(XI))*EXP(-DI(NS)*SQRT(XI*XI+XJ*XJ))
Q(I,J)=Q(I,J)-AJ(NS)*(YL/ABS(YL))*EXP(-CJ(NS)*ABS(YL))
*-BJ(NS)*(YJ/ABS(YJ))*EXP(-DJ(NS)*SQHT(YI*YI+YJ*YJ))
820   CONTINUE
GO TO 400
830   CONTINUE
RETURN
END
```

4. Examples for Computer Exercises