NASA Technical Memorandum 103766

//_'-_ /

GRID3D-v2: An Updated Version of

the GRID2D/3D Computer Program for

Generating Grid Systems in Complex-Shaped

Three-Dimensional Spatial Domains

E. Steinthorsson and T.I-P. Shih

Carnegie Mellon University

Pittsburgh, Pennsylvania

and

R.J. Roelke

Lewis Research Center

Cleveland, Ohio

June 1991

N_A

1.0

2.0

3.0

4.0

CONTENTS

INTRODUCTION

THEORY AND METHODOLOGY

2.1 Boundary Curves and Distribution of Grid Points

2.2 Connecting Curves Based on Tension Spline Interpolation

2.3 Specification of Derivatives at Boundaries

USING GRID2D/3D-v2

3.1 Generation of a 3-D Grid Using GRID3D-v2

3.2 Example: A Spatial Domain With Irregular Boundaries

SUMMARY

APPENDIXES

A Support Programs

A. 1 Description of Support Programs

A.2 Listing of PRSURF

A.3 Listing of 3DPREP

A.4 Listing of EDGE

A.5 Listing of EDGPREP

A.6 Listing of GRIDTST

B Listing of GRID3D-v2

REFERENCES

Page

1

3

3

7

10

12

12

13

15

15

15

2O

28

42

5O

54

5B

93

GRID3D-v2: An Updated Version of the GRID3D

Computer Program for Generating Grid Systems in

Complex-Shaped Three-Dimensional Spatial Domains

E. Steinthorsson and T. I-P. Shih

Department of Mechanical Engineering

Carnegie Mellon University

Pittsburgh, Pennsylvania 15213

R. J. Roelke

National Aeronautics and Space Administration

Lewis Research Center

Cleveland, Ohio 44135

1.0 INTRODUCTION

GRID2D/3D (refs. 1 and 2) is a powerful grid-generation package, capable of generating grid

systems for complicated geometries in both two and three dimensions. This package, which

employs algebraic grid-generation techniques, is computationally efficient and easy to use.

Nonetheless, when the geometry is unusually complex (e.g., see fig. 1-1), the partitioning of the

geometry into zones or blocks that are suitable for GRID2D/3D becomes very cumbersome. In order

to make GRID2D/3D more versatile and readily applicable to geometries like the one shown in figure

1-1, a number of modifications were made to the package. These modifications are as follows:

(1) specification of boundary curves has been made more flexible;

(2) control over grid-point distribution has been increased;

(3) new interpolating functions based on tension splines have been added; and

(4) control over orthogonality at boundary surfaces has been increased.
GRID2D/3D is madeupof two programs,GRID2DandGRID3D. GRID2D generatesgrid

systemsfor two-dimensional(2-D) spatialdomains,andGRID3Dgeneratesgrid systemsfor threedimensional(3-D) domains.Theaforementionedmodificationsweremadeto GRID3Donly. In this

report, the original version of GRID3D asreportedin references1 and 2 will bereferred to as

GRID3D-v1. Thenewmodifiedversionwill bereferredto asGRID3D-v2.

In theremainderof thisreport,thetheoryandmethodbehindthemodificationsaredescribed,

andtheuseof GRID3D-v2is explainedandillustratedby anexample.
2.0 THEORY AND METHODOLOGY

In this section, the theory and methodology behind the modifications that were incorporated

into GRID3D-v2 are described. First, specification of boundary curves is discussed, and a new

approach for controlling distribution of grid points is explained. Next, the functional forms of new

tension-spline based connecting curves are derived, and the properties of these functions examined.

Finally, control of orthogonality of grid lines at boundaries is discussed.

Note that throughout this section, the reader is assumed to be familiar with the theory behind

GRID3D-v 1 which is described in reference 1.

2.1 Boundary Curves and Distribution of Grid Points

When generating 3-D grid systems with GRID3D-vl or GRID3D-v2, we assume that the

geometry of the spatial domain for which a grid is to be generated is completely described by the

edge curves of the boundary surfaces (as used here, the edge curves are the four plane or twisted

curves that define the boundary of a surface). The algorithms used in GRID3D-vl and GRID3D-v2

were designed to construct grid systems from the edge curves; however, these algorithms differ in

the way they specify edge curves and control grid-point distributions. This difference is described in

this section.

In GRID3D-vl, the following approach is used:

(1) Each edge curve is given by specifying a set of node points that lie on the curve. From

these node points, a parametric description of the curve is constructed by using tensionspline interpolation.

(2) The distribution of grid points within the domain (including edge curves) is controlled by

specifying three one-dimensional stretching functions -- one for each family of grid lines

(i.e,, one for _-grid lines, one for rl-grid lines, and one for _-grid lines).

Although this approach provides flexibility in generating grid systems within complexshaped spatial domains, it is inadequate in some cases. For example, it does not allow for edge

curves to have derivative discontinuities such as those at cusps. Also, it does not provide adequate

control over distribution of grid points in regions where geometry changes appreciably. As a

specific example, the distribution of grid points on all constant-_ surfaces must be the same and

cannot, even with partitioning, be made to vary from one section of the grid system to the next.
To overcome the shortcomingsof GRID3D-vl, the following approach was used in

GRID3D-v2:

(1) Each edge curve is defined by specifying the location of the grid points on the curve.

This can be done by either one of the following two procedures: (a) Specify a set of node

points that lie on the edge curve for interpolation with a tension spline, and specify a

stretching function that controls where on the edge curve the grid points lie. (b) Specify

the grid points directly (in this case, a stretching function is not specified by the user,

rather it is calculated based on arc length as will be shown herein).

(2) The grid point distribution along grid lines of a given family is that obtained by the

bilinear interpolation of the stretching functions used for the edge curves belonging to that

family.

The specification of edge curves as described under (l(a)) is self explanatory, but the

calculation of stretching functions as described under (1 (b)) and (2) requires further explanation.

For illustration, consider the _-family of grid lines. Suppose all edge curves belonging to this

family of grid lines are specified by using procedure (l(a)), and suppose the stretching functions

_00(;), _10(;), _01(;) and _ll(;) describe the distribution of the grid points on the edge curves

located at (_ = 0,rl = 0), (_ = 1,1"!= 0), (_ = 0,rl = 1) and (_ = 1,1"1= 1), respectively. In GRID3Dv2, the stretching function for a _-grid line at any _-lq location is given by the bilinear interpolation

of the stretching functions at the edge curves; namely

_n(;) = 00(;) (l-k) + ;10(;) (1-rl) + 01(_) (l-f) + _11(_) 11 (9"°1)

Now, instead of all edge curves being defined by (l(a)), suppose that the edge curve at _ = 0

and 1"1-0 is defined by using procedure (l(b)); that is, by specifying the grid point coordinates

directly. In order to use equation (2.1) for this case, a stretching function must be calculated for the

edge curve. Since the stretching function is needed only at the grid point-locations along the edge

curve, it can be calculated by using approximate arc length as follows:

and

where

and

_oo(_k)= 0 k=l (2.2,_)

_00(_k) = d---k-k k=2,3,4 ..... KL (2.2b)

dKL

k

dk : E [(Xn'Xn-1)2 +(Y n'yn-1)2 + (Zn-Zn-1)2] 1/2

n=2

(2.2c)
tk = (k-l) At At -- I/(KL-I) (2.2d)

In equation (2.2), k = 1,2,3 ..... KL denotes the grid points on the edge curve; KL is the total

number of grid points on the curve; and Xn, Yn and Zn are the x-, y- and z-coordinates of the n-th grid

point on the curve.

The same approach as that just described for the t-family of grid lines is employed to

determine the stretching functions for the _- and rl-families of grid lines. This approach gives a

smooth distribution of grid points throughout the domain.

Finally, note that all stretching functions available in GRID3D-vl for controlling the

distribution of grid points in the entire domain are available in GRID3D-v2 for controlling the

distribution of grid points along edge curves defined by using procedure (1 (a)) (i.e., by specifying a

set of node points that lie on the edge curve and interpolating with a tension spline). In GRID3D-v2,

a stretching function that allows asymmetric clustering of grid points along the edge curve was

added. The new stretching function, which was developed by Vinokur (ref. 3), is given here.

Let t e [0,1] be normalized distance or any monotonic parameter along a curve, and let _

[0,1] be the computational coordinate along which grid points are equally spaced. Two user

controlled parameters, so and Sl, and two secondary parameters, A and B, are defined as

d_(t : I) (2.3) d_(t = 0) and sl= so,Sl > 0

so =- dt dt

A = _]- and B = "/sosl (2.4)

In terms of these parameters, the functional form of the stretching function can be written as

t(_) =

A + (l-A) u(¢) (2.5)

where the function u(_) depends on the value of the parameter B, as shown in the following

equations.

IfB > 1.001, then

tanh[Ay (¢-_)] (2.6a) u(¢) = 1 +

2 tanh[Ay/2]

5
whereAy is obtainedfrom therelation

B - sinh(Ay)

Ay

If B < 0.999, then

tan[_ (_- _)]

2 tan[Ard2]

(2.7a)

where Ax is obtained from the relation

B = sin(Ax) (2.7b)

Ax

Finally, if 0.999 < B < 1.001, then

u(_)=_[1 +2(B-1)(_-1⁄2)(1-_)]

The amount of clustering produced by the stretching function is controlled by the parameters

so and Sl which are defined by equation (2.3). If so and sl are greater than one, then grid points are

clustered near the boundaries where t = 0 and t = 1. The greater so and sl, the greater is the

clustering of grid points near the t = 0 and t = 1 boundaries, respectively. If so and Sl are less than

one, then the grid spacing is larger near the boundaries than in the interior; the smaller so and Sl, the

greater the grid spacing near the boundaries.

Frequently, when using the stretching function given by equations (2.3) to (2.8), we must

either solve equation (2.6(b)) for Ay, or solve equation (2.7(b)) for Ax. Vinokur (ref. 3) developed

approximate analytical relations for both of these inversion problems as follows:

For equation (2.6(b)), which is used when B > 1.001, the approximate inverse when 1.001

< B < 2.7829681 is

Ay ='Vr6-_-(1 - 0.1513 + 0.057321429_ 2- 0.024907295133

+ 0.0077424461_ 4- 0.0010794123_ 5)

(2.0-,)

where

13= B- 1 (2.0b)
where

and

When B > 2.7829681,

Ay = v + (1 + I/v) ln(2v) - 0.02041793 + 0.24902722w

+ 1.9496443w 2 - 2.6294547w 3 + 8.56795911w 4

(2. lOa,)

v = In (B) (2.10b)

w = (l/B) - 0.028527431 (2.10c)

For equation (2.7(b)), which is used when B < 0.999, the approximate inverse when 0 < B

< 0.26938972 is

Ax = rt [1 - B + B 2- (1 + _2/6) B 3 + 6.794732B 4

- 13.205501B 5 + 11.726095B6)]

(2.11)

where

When 0.26938972 < B < 0.999,

Ax = _ (1 + 0.1513 + 0.057321429132 + 0.048774238133

- 0.053337753134 + 0.075845134135)

(2.12e,)

13= 1- B (2.12b)

2.2 Connecting Curves Based on Tension Spline Interpolation

Experience has shown that Hermite interpolation (cubic polynomials) as used in GRID3D-v 1

sometimes results in connecting curves with too much curvature. In GRID3D-v2, the Hermite

interpolation is replaced by tension-spline interpolation. The most attractive feature of tension-spline

interpolation is that as the tension parameter is increased from zero to infinity, the interpolation

function varies from being a cubic polynomial to being a linear polynomial. Thus, tension-spline

interpolation offers increased control over the shape of the grid lines in the grid system. In this

section, a derivation of the tension-spline interpolation function is given for the two- and four-
boundarymethods.First,theinterpolationfor anarbitraryvariableis derived. Then, the application

to algebraic grid generation is illustrated.

The tension-spline interpolation function is derived as follows: suppose the variable X is a

function of the parameter s on an interval [0,1], but only X(0), X(1), X'(0) and X'(1) (X' denotes

dX/ds) are known. A tension-spline interpolation of X(s) on the interval [0,1] is sought. A tensionspline interpolation of X(s) is traditionally written in terms of X(0), X(1), X"(0) and X"(1), where

X" = d2X/ds 2, as follows (see, e.g., ref. 4):

X(s) =

X"(0) sinh[t_( l-s)]

_-(X(O) _ X"(O)/(l_s)

0 2 sinh[G] G 2 I (2.15)

+

X"(1) sinh[as] + (X(1)X"(1)_) s

a 2 sinh[c] c 2

where t_ is the tension parameter. By differentiating equation (2.13) and evaluating the resulting

equation at the end points s = 0 and s = 1, we obtain

( )X l, /

X'(0) = X"(0) cosh[c] X(0) X"(0) + t- X(1) - (2.14_)

a sinh[_] _2 _ sinh[_] _2

and

X'(1) = - (s sinh[o] (x o,x o, +x" l,cosh,o, c2 ! c sinh[c_] _ (X(1) - (I 2

The above two simultaneous equations can be solved to give expressions for X"(0) and X"(1) in terms

of X(0), X(1), X'(0) and X'(1). Substituting the resulting expressions into equation (2.13) gives

X(s) = X(0)hl(S) + X(1)h2(s) + X'(0)h3(s) + X'(1)h4(s)

/sinh[c(1-s)] - sinh[os] 1

hl(s) = Cl(1-s) + c2s + c2/ si---_[_] J

(9.. 15)

(2.16a)

/sinh[o(1-s)] - sinh[cs] ) h2(s) = ClS + c2 (l-s) - c2 / s_nh[_] (2.16b)
h3(s) = c3 ((l-s)-sinh[a(1-s)] )+c4 (s h,os ) ssinh[c] sinh[a]

(2.16e)

h4(s)= -C4((1-S)-sinh[cr(1-s)] sinh[o] )-c3 (s -sinh[o] sinh[os] ) (_.xea)

where

cl = 1 - c2 (2.17n)

sinh[o] (2.17b)

c2 = 2 sinh[o] - o cosh[o] - o

= " _ sinh[o] (2.17e)

= 13 sinh[o] (2.17d)

c4 {132_(_2)

O_= 0 cosh[o] - sinh[o] (2.17e)

13= sinh[_] - o (2.17f)

Equations (2.15) to (2.17) can be used to interpolate any function on an interval [0,1], when

the function's values and its first derivatives are known at the end points of the interval. The

application of these equations to algebraic grid generation is straightforward. Consider, for example,

the two-boundary technique (ref. 1). Suppose we want to generate a grid system between two

constant-rl boundary surfaces. The formulation in this case can be written as follows:

r (_,rl,_) = r (_,rl = 0,_) hl(rl) + r (_,q = 1,_) h2(rl)

(3r(_,,r I = 1,_)

+ _r(_,rl = 0,4) h3(q) + h4(rl)

0n 0n

(2.18)

where
x(_,_,4)

r (_,rl,4) =/y(_,rl,4)

Lz(_,n,4)

(_.ll)

The above equation has the same form as if Hermite interpolation were used, except for the definition

of the functions hl, h2, h3, and h4, which in the case of Hermite interpolation are cubic polynomials

(ref. 1). Note, however, that the functions used in tension-spline interpolation (eqs. (2.16) and

(2.17)) have the property that as o--o0, hi, h2, h3, and h4 approach the cubic polynomials used in

Hermite interpolation (see ref. 1). Also, as o--o-o, hl(s)--@(1-s), h2(s)---_s, h3(s)---_0, and h4(s)--o0,

giving rise to linear connecting functions (Lagrange interpolation; see ref. 1).

2.3 Specification of Derivatives at Boundaries

When the two- and four-boundary methods, as described in reference 1, are used to generate

grid systems, the first-order derivatives involved (e.g., _r(_,rl = 0,4)/3'q and _r(_,ri = 1,4)/i9rl in

equation (2.18)) need to be specified. In GRID3D-vl and GRID3D-v2, these derivatives are

specified such that grid lines intersect boundary surfaces orthogonally. In this section, the methods

used in GRID3D-vl and GRID3D-v2 to calculate these derivatives will be explained. The twoboundary technique given by equation (2.18) will be used to illustrate the concepts.

In GRID3D-vl, the derivatives Or({,rl = 0,4)//)rl and _gr(_,rl = 1,4)/_)rl (eq. (2.18)) are

chosen as follows:

ar(_,lq =0,4) = K_1(_,4) t_l ar(_'rl = 1'4) = K_2(_,_) t_2 (2.20)

0q _n

where

br(_,,T1 = O,r,) br(_,rl = 0,4)

tql = - x (2.21a)

¢3r(_,q = 1,_) _r(_,rl = 1,_)

tq2 = - × (2.21b)

In equation (2.20), the terms Kql(_,_) and Krl2(_,_) -- that is,the K-factors -- are specified by the

user and are intended to control the magnitude of the derivatives. However, the magnitude of the

10
derivatives will also depend on the magnitude of the vectors %1 and tn2, which in turn depend both

on the geometry of the boundary surfaces at r I = 0 and rI = 1, respectively, and on the grid spacing

on the surfaces. Thus, full control cannot be exerted over the magnitudes of the derivative terms

0r(__,rl=0,t_)/011 t,.. and 0r(_,q=l,_)/0q, and for complex-shaped geometries this may pose a problem.

In order to overcome the aforementioned problems, in GRID3D-v2 the derivative terms

0r(__,q = O,_)/o_q and ar(_,q = 1,_)hgq were defined as follows:

ar( ,n = (2.22) _gr(__,H = 0,_) _ Knl(_,_) - = Kr12(7_,_) erl2

0_1 __ e_l O_l -

where en t and en2 are unit vectors that are normal to the boundary surfaces at r I = 0 and rI = 1,

respectively: that is

tnl t_2 (2.23)

This approach allows total control over the magnitude of the derivatives Or(_,q = O,_)/Oq and

Or(_,q = 1,_)/Oq through the K-factors alone.

Finalls'. note that in GRID3D-vl, the K-factors were taken to be constants on each boundary

surface, even though the authors realized that they could be allowed to vary. In GRID3D-v2, the Kfactors are allowed to vary from point to point and can be controlled by the user.

In the next section, we show how GRID3D-v2 is used to generate grid systems. Several

auxiliary programs were written to assist grid generation with GRID3D-v2. A description of these

programs is given in appendix A. A complete listing of the GRID3D-v2 computer program is given

in appendix B.

11
3.0 USING GRID3D-v2

In this section, the use of GRID3D-v2 is described. The section consists of two parts: first,

an explanation of generating 3-D grid system with GRID3D-v2; and second, an example of such a grid

system generated by using GRID3D-v2.

3.1 Generating a Three-Dimensional Grid System with GRID3D-v2

When using GRID3D-v2, the user must answer the following questions:

(1) Should the two boundary technique or the four boundary technique be

used?

(2) Should any edge curve be defined by specifying the grid points

directly?

(3) How many grid points are desired in each of the _, 1"1,and _ directions?

(4) Is any clustering of grid points needed?

(5) What K-factors should be used (see section 2.3 and pp. 18 to 20 in ref.

1) and are constant K-factors sufficient?

Once these questions have been answered, an input file for GRID3D-v2 must be constructed. The

form of the input file is shown in figure 3-1. The various parameters in the input file are explained in

table 3.1. The surface and edge curve numbering scheme embedded in GRID3D-v2 is shown in

figure 3-2. Note that some of the edge curves are identical; that is, curves 3 and 9 are identical and

so are curves 4 and 13, curves 7 and 10, and curves 8 and 14 (for further explanation of the edge

curve numbering scheme, see p. 5 of ref. 2)

In the input file are values for K-factors at boundary surfaces -- a single value is assigned for

all grid points on each surface. The user can modify the K-factor for any individual grid point or

groups of grid points on each surface by adding FORTRAN statements into subroutine KFACTOR

(see listing in Appendix B). Note that in addition to being used to generate grid points in the interior

of the spatial domain, K-factors are also used to generate grid points on boundary surfaces

themselves. For this latter case, the relevant K-factors are those at grid points that lie on the four

edge curves bounding the surface. These K-factors can also be modified in subroutine KFACTOR.

Once an input file has been prepared (and an executable file created for GRID3D-v2 if

subroutine KFACTOR was modified), GRID3D-v2 can be executed. Note, GRID3D-v2 reads the

input file from unit 7 and writes output to unit 8.

12
Grid generationis aniterativeprocess;thatis, anacceptablegrid systemis generatedby trial

anderror after generatinga seriesof unacceptablegrid systems.Someobservationsandrulesof

thumbthatmightbeusefulwhengeneratinggrid systemswith GRID3Dareasfollows:

(1) If it is necessaryto partition aspatialdomain,thenselectpartitioning surfaces

that intersectboundarysurfacesasorthogonallyaspossible. This approach

minimizesskewnessbothatboundariesandin theinteriorof thedomain.

(2) A commonboundarybetweentwopartitionsshouldbespecifiedin exactly the

sameway in theinput files in orderto guaranteeacontinuousgrid acrossthe

commonboundary.

(3) For improved flexibility in modifying a grid system,the boundarycurves

shouldbedefinedby usingnodepoints for tension-splineinterpolationand a

stretchingfunction(seeSection2.1),ratherthanspecifyinggrid pointsdirectly.

(4) When generating a grid system for a complicated geometry, first find stretching

functions and K-factors for edge curves that give the desired distribution of grid

points on the boundary surfaces. Afterwards, try to optimize grid-point

distribution in the interior by modifying K-factors on boundary surfaces and/or

the amount of tension in the connecting curves. If this two-step process does

not yield an acceptable grid system, then try to modify grid-point distribution on

boundary surfaces and/or re-partition the domain.

(5) Start with low tension and low K-factors. Slowly increase K-factors to

improve orthogonality at boundaries and eliminate overlapping grid lines.

Increasing the amount of tension in the connecting curves also can eliminate

overlapping grid lines by straightening grid lines.

(6) Increased tension tends to straighten out grid lines whereas increased K-factors

tend to increase the curvature of grid lines.

(7) If K-factors are too high, then grid lines can overlap.

(8) K-factors affect the grid spacing near the boundary surfaces. Increasing the Kfactor at a boundary increases the grid spacing adjacent to that boundary. This

effect can be beneficial in some instances but detrimental in others.

(9) K-factors should vary smoothly from grid line to grid line. An exception to

this general rule is when a grid line intersects a boundary surface at a cusp in the

surface.

3.2 Example: A Spatial Domain With Irregular Boundaries

Figure 1-1 shows a cooling passage in a radial turbine blade. Figure 3-3 shows how this

cooling passage was partitioned into blocks or zones for the purpose of grid generation. The

partitioning that is shown was deemed necessary in order to get an acceptable grid system.

13
Figure3-4showsthegrid systemfor partitionnumber 18. Theinputfile for thispartitionis

givenin table 3-2. Theentiregrid systemgeneratedby usingGRID3D-v2 is shownin figures3-5

and3-6. Theplottingin figures3-4 to 3-6wereobtainedby using3DSURF-- theplottingpackage

suppliedwithGRID2D/3D.

4.0 SUMMARY

A new version of the grid generation program GRID3D, which is a part of the grid

generation package GRID2D/3D, has been developed. The new program is referred to as GRID3Dv2. This report describes GRID3D-v2 and how to use it. The capability of the program was

demonstrated by generating a grid system for a very complicated geometry, namely a cooling

passage inside a radial turbine blade.

14
Appendix A-- Support Programs

A.I Description of Support Programs

To support the task of grid generation with GRID3D-v2, several auxiliary programs have

been developed. Two of these programs, namely PRSURF and 3DSURF, were developed to allow

the user to view grid systems generated with the package. Three more programs, namely 3DPREP,

EDGE, and EDGPREP, were developed to aid in the preparation of input files for GRID3D-v2.

Finally, one program, GRIDTST, was written to test the grid system for problems such as

overlapping grid lines. Each of these programs and their use will be described in the following

pages.

3DSURF

The program 3DSURF was developed alongside the original version of GRID2D/3D to allow

the user to plot grid systems on the computer screen. It was designed for IBM PC, XT and AT

computer systems and compatibles. Two-dimensional grid output files from GRID2D/3D are

already in the proper format for use with 3DSURF. Three-dimensional grid output files must, on the

other hand, be processed using PRSURF (described next) to create input files for 3DSURF. A

general description of 3DSURF and its use is given in Section 3.1 in reference 2.

PRSUrRF

The PRSURF program was written to process output files from GRID3D and create input

files for 3DSURF. This program allows the user to select surfaces or parts of surfaces from the grid

system (i.e., constant-_, constant-q, and constant-_ surfaces) and to store them in a format

compatible with 3DSURF. The program is interactive and self-explanatory.

PRSURF uses three files. The file containing the input to PRSURF (i.e., the output file

from GRID3D) is read from unit 7. The output file from PRSURF (which becomes the input file

for 3DSURF) is written to unit 9. Last, a file used for temporary storage of data is accessed as unit

8.

15
Finally, we mentionthat since 3DSURFis limited to handlingsurfacesthat have40 grid

pointsper sideor less,PRSURFautomaticallybreaksanygrid surfaceinto sectionsthatare40grid

pointsby 40grid pointsor smaller.

3DPREP

The user can apply 3DPREP to create input files for GRID3D-v2. The program prompts the

user for all control parameters such as stretching functions and k-factors. It reads from files the

coordinates of points defining the edge curves. These files must have the following format:

NP (number of points)

xl, Yl, Zl

x2, Y2, z2

x3, Y3, z3

XNp, YNP, zNP

The program can read data points from these files in both forward and reverse order, and the user

can let the program read data from several files (the whole file or only a part of the file) to put

together a single edge curve. 3DPREP is useful, primarily, when many input files for GRID3D-v2

need to be created from the same set of data.

3DPREP is designed to run on the IBM PC, XT, AT, and compatibles, but it can easily be

modified to run on other computer systems. The only modification that should be needed in such a

case is the insertion of open statements into the program so that the user can interactively specify

which files the program must access.

EDGE

The EDGE program was written to aid the user in generating grid points along an edge curve

that cannot be represented by a single spline curve (e.g., edge curves possessing derivative

discontinuities such as cusps) and which must, therefore, be defined in the input files for GRID3Dv2 by giving the grid point coordinates directly (see Section 2.1). The program was designed to

16
generategrid points on anedgecurvethat is composedof severalsections,whereeachsectionis

definedby a setof nodalpointsthatareinterpolatedby asplinecurve. Thenumberof gridpointson

eachsectionandtheirdistributionwithinthesectioniscontrolledindependently.

EDGE, which is designed to run on the IBM PC, XT, AT, and compatibles, reads input data

from UNIT 1 and writes the output to UNIT 20. The input file must have the following format:

NS (number of sections)

Data for1 section

Data for section 2

I Data forNSsection

where the data for an arbitrary section number i are the fo110wing:

IPi (number of grid points on section i )

_i (tension parameter for the spline interpolation)

NNi (number of node points given on section i )

xl, Yl, Zl

x2, Y2, z2

x3, Y3, z3

XNNi, YNNi, ZNNi

Su-etchTypei

Betali, Bet_i

The meaning of the p:u'ameters StretchType, Betal, and Beta2 is explained in table 3-1. Note that in

EDGE it is assumed that the curve being generated is continuous; that is, the last grid point on one

17
sectionmustbethesameasthefirst grid point onthenextsection.Thustheoutputfrom EDGEhas

thefollowing format:

IL (numberof gridpointsontheedgecurve)

Xl,1, Yl,1, Zl,1

Xl,2, Yl,2, Zl,2

X1,3, yl,3, Zl,3

Xl,IPl-1, Yl,IPI-1, ZI,IPI-1

X2,1, y2,1, Z2,1

X2,2, Y2,2, Z2,2

X2,3, Y2,3, Z2.3

X2,1P2-1, Y2,1P2-1, Z2,IP2-1

XNS,1, YNS,1, ZNS,I

XNS.2, YNS,2, ZNS,2

XNS.3, YNS,3, ZNS,3

XNS,IPNs, YNS,IPNS, ZNS,IPNS

where xi.j, Yi.j, and Zi, j

number of grid points is

are the coordinates of grid point number j on section number i. The total

NS

IL=I +_{IPi-1}

i=l

18
Intentionally, this is the sameformat asfor the input files for 3DPREPso theoutput files from

EDGEcanbereadby 3DPREP.

GRIDTST

The GRIDTST program is used to check grid systems of 3-D spatial domains for defects,

such as over lapping grid lines, that result in a negative Jacobian, where the Jacobian is defined as

follows (for further explanation of the Jacobian, see ref. 1):

J = x_(y_z; - y;z,q) - xrl(y_z; - y;z_) + x;(y_z_ - yrlz_)

GRIDTST evaluates the Jacobian at every grid point, estimating the derivatives (i.e., x_, y_,

z_,, etc.) by using central difference approximations for grid points that do not lie on the boundary

surfaces of the domain and by using second order accurate one-sided difference formulas where

necessary on the boundary surfaces. If any negative Jacobians are found, GRIDTST prints a

message on the computer screen, and the Jacobians and the grid point locations are written into a file.

The input into GRIDTST is the grid system generated by GRID3D-v2. The input is read

from UNIT 1 whereas the output (if any) is written into UNIT 2. GRIDTST is written to run on the

IBM PC, AT, XT, and compatibles, but it can be used on any computer system.

19
A.2 Listing of PRSURF

PROGRAM PRSURF

C

C This program writes out user picked surfaces from a 3-D grid I

C

PARAMETER (IM=II,JM=51,KM=I51)

INTEGER i, j, k, IL, JL, KL

C

5

5001

C

C

C

C

C

C

C

REAL X (IM, JM, KM), Y (IM, JM, KM), Z (IM, JM, KM)

IERR=0

WRITE (*, 5001)

READ (*, *) IPICK

FORMAT('2',///, ' Please enter:',//,

$' 1 - if you want the boundary surfaces of the grid',/,

$' to be saved',/,

$' 2 - if you want to select surfaces to be saved',///)

IF (IPICK.EQ. i)THEN

CALL PRGRID (X, Y, Z, IM, JM, KM)

ELSE IF(IPICK.EQ.2)THEN

CALL PRSRFS (X, Y, Z, IM, JM, KM)

ELSE IF(IERR.EQ.0)THEN

WRITE (*, *) '

WRITE (*, *) '

WRITE (*, *) '

READ (*, *)

IERR= 1

GOTO 5

ELSE

STOP

ENDIF

You will get one more chance to make a'

selection - Please enter any character'

and then press RETURN'

STOP

END

SUBROUTINE PRGRID(XPnt,YPnt,ZPnt,IL2,JL2,KL2)

C This subroutine reads in grid point coordinates and writes out the

C coordinates of the grid points which lie along specified planes.

INTEGER i, j, k, !L, JL, KL, ILl, JLI, KLI

REAL

$

$

XPnt (IL2, JL2, KL2),

YPnt (IL2,JL2,KL2),

ZPnt (IL2, JL2, KL2)

2O
C Read in the grid size.

READ (7, *) IL

READ(7, *) JL

READ (7,*) KL

C Read in the grid point locations.

DO 7 i=l, IL

DO 6 j=l, JL

DO 5 k=I,KL

READ (7,*) XPnt(i, j,k),XPnt(i, j,k),ZPnt (i, j,k)

CONTINUE

CONTINUE

CONTINUE

C Calculate the number of sections the grid must be split into

C for plotting purposes. (Plotting routines can handle only a grid

C with a maximum dimension of 40. Here it is assumed that only the

C zeta-coordinate direction can involve more grid points than that)

KS= (KL-I)/39

JS= (JL-I) /39

C Print out the number of surfaces.

IF((KS*39+I).LT.KL)THEN

NSK=KS+I

ELSE

NSK=KS

ENDIF

IF((JS*39+I).LT.JL)THEN

NSJ=JS+I

ELSE

NSJ=JS

ENDIF

NOSURF=2*NSJ+2*NSK+3*(NSK*NSJ)

WRITE(9,*) NOSURF

C Print out the grid points.

DO 21 m=I,NSJ

j0=39" (m-l) +I

jI=MIN (j0+39, JL)

JLm=jl- j0+l

WR!TE (9, *) !L

WRITE (9, *) JLm

i0

20

DO 20 i=l, IL

DO I0 j=j0, jl

7_....r rn•E(9 t 25) XPnt(i, j,KL),YPnt(i, j,KL),ZPnt(i, j,KL)

CONTINUE

CONT INUE

21
21

25

27

28

29

CONTINUE

FORMAT(IX, FI0.6,3x, FI0.6, 3X,FI0.6)

DO 29 n=I,NSK

k0=39" (n-l) +I

kI_MIN (k0+39,KL)

KLn=kl-k0+l

WRITE(9, *) IL

WRITE(9, *) KLN

DO 28 i=l,IL

DO 27 k=k0,kl

WRITE(9,25)

CONTINUE

CONTINUE

CONTINUE

XPnt(i, JL, k), YPnt(i, JL, k), ZPnt (i, JL, k)

42

43

44

DO 44 m=I,NSJ

j0=39*(m-l)+l

jI=MIN(j0+39,JL)

JLm=jl-j0+l

DO 44 n=I,NSK

k0=39*(n-l)+l

kI=MIN(k0+39,KL)

KLn=kl-k0+l

WRITE(9, *) JLm

WRITE(9,*) KLn

DO 43 j=j0, jl

DO 42 k=k0,kl

WRITE(9, 25)

CONTINUE

CONTINUE

CONTINUE

XPnt(IL, j, k), YPnt(IL, j, k), ZPnt(IL, j, k)

DO I!0 m=I,NSJ

j0=39 * (m-l) +I

jI=MIN (j0+39, JL)

JLm=jl-j 0+I

DO ii0 n=I,NSK

k0=39" (n-l) +I

kI=MIN (k0+39,KL)

KLn=kl-k0+l

WRITE(9, *) JLm

WRITE(9,*) KLn

22
108

109

II0

DO 109 j=j0,jl

DO 108 k=k0, kl

WRITE(9,25)

CONTINUE

CONTINUE

CONTINUE

XPnt(l, j,k),YPnt(l, j,k),ZPnt (I, j,k)

125

130

131

132

133

134

DO 131 m=I,NSJ

j0=39 * (m-l) +I

j I=MIN (j0+39, JL)

JLm=j i- j 0+I

WRITE (9,*) IL

WRITE (9, *) JLm

DO 130 i=l,IL

DO 125 j=j0, jl

WRITE(9,25)

CONTINUE

CONTINUE

CONTINUE

DO 134 n=I,NSK

k0=39" (n-l) +I

kI=MIN (k0+39, KL)

KLn=kl-k0+l

WRITE (9, *) IL

WRITE (9, *) KLn

DO 133 i=l,IL

DO 132 k=k0,kl

WRITE (9,25)

CONT I NU E

CONTINUE

CONTINUE

XPnt (i, j, i), YPnt (i, j, I), ZPnt (i, j, i)

XPnt (i, I, k), YPnt (i, i, k), ZPnt (i, i, k)

DO 144 m=I,NSJ

j0=39" (m-l) +I

jI=MIN (j0+39, JL)

JLm= j i- j0+l

DO 144 n=I,NSK

k0=39*(n-l)+l

kI=MIN(k0+39,KL)

KLn=kl-k0+l

WRITE (9, _) JLm

WRITE(9, *) KLn

!H=(IL+I)/2

23
142

143

144

DO 143 j=j0, jl

DO 142 k=k0,kl

WRITE(9,25) XPnt(IH, j,k),YPnt(IH, j,k),ZPnt(IH, j,k)

CONTINUE

CONTINUE

CONTINUE

C

C

STOP

END

SUBROUTINE PRSRFS (X, Y, Z, IM, JM, KM)

C This subroutine reads in grid point coordinates and writes out the

C coordinates of the gri d points which lie along planes specified

C by the user.

INTEGER i, j, k, IL, JL, KL

REAL X (IM, JM, KM), Y (IM, JM, KM), Z (IM, JM, KM)

C

WRITE (*, 5001)

C

C Read in the grid size.

READ (7, *) IL

READ(7,*) JL

READ (7, *) KL

C Read in the grid point locations.

DO 3 i=l, IL

DO 2 j=I,JL

DO 1 k=I,KL

READ (7, *) X(i, j, k) ,Y (i, j, k), Z (i, j, k)

CONTINUE

CONTINUE

CONTINUE

NOSURF=0

WRITE (*, 5002) IL, JL, KL

READ (*, *) IPICK

IF (IPICK.EQ. i) THEN

WRITE (*, 5003)

READ (*, *) IPICK2

IF (IPICK2 .EQ. I) THEN

WRITE (*, *)

WRITE(*,*) 'Please enter the value of i'

READ(*, *) I

JFIRST=I

JLAST=JL

KFIRST=I

KLAST=KL

ELSE

WRITE (*, *)

WRITE(* l *)'Please enter the value of i'

READ (*, *) I

24
30

4O

5O

60

WRITE(* *)

WRITE(*,*)' Please enter the lower and upper limit'

WRITE(*,*)' for the j coordinate (JFIRST,JLAST)'

READ (*, *) JFIRST, JLAST

WRITE (* *)

WRITE(*,*)' Please enter the lower and upper limit'

WRITE(* *)' for the k coordinate (KFIRST, KLAST)'

READ (*, * )KF IRST, KLAST

ENDIF

NSJ= (JLAST-JFIRST)/39

IF (JFIRST+NSJ*39. LT. JLAST) NSJ=NSJ+I

NSK= (KLAST-KFIRST)/39

IF (KFIRST+NSK*39. LT. KLAST) NSK=NSK+I

NOS URF=NOSURF+NS J*NSK

DO 60 m=I,NSJ

j0--39" (m-l) +JFIRST

jI=MIN (j0+39, JLAST)

JLm=jl- j 0+I

DO 50 n=I,NSK

k0=39 * (n-l) +KFIRST

kI=MIN (k0+39, KLAST)

KLn=kl-k0+l

WRITE (8, *) JLm

WRITE (8, *) KLn

DO 40 j=j0, jl

DO 30 k=k0, kl

WRITE (8, 25) X(I, j, k), Y (I, j,k), Z (I, j, k)

CONTINUE

CONTINUE

CONT I NU E

CONTINUE

IERR=0

GOTO 5

ELSEIF (IPICK. EQ. 2) THEN

WRITE (*, 5003)

READ(* *) IPICK2

IF (IPICK2 .EQ. i) THEN

WRITE (*, 2)

WRITE(*, 2) 'Please enter the value of j'

READ (*, *) J

IFIRST=I

ILAST=IL

KFIRST=I

KLAST=KL

ELSE

WRITE(* *)

WRITE(*,*) 'Please enter the value of j'

READ (*, *) J

WRITE (2, .)

WRITE(*,*)' Please enter the lower and upper limit'

Z_RITE(.,.) , for the i coordinate (IFIRST, ILAST) '

READ (2, 2) IFIRST, ILAST

_R: _E , )

WRITE(*,*) ' Please enter the lower and upper limit'

Z_.ITE(*,*) ' for the k coordinate (KFIRST,KLAST) '

READ (*, _)KFIRST,KLAST

25
70

8O

90

I00

C

NSI-- (ILAST-IFIRST)/39

IF (IFIRST+NSI*39. LT. ILAST)NSI=NSI+I

NSK= (KLAST-KFIRST)/39

IF (KF IRST+NSK* 39. LT. KLAST )NSKmNSK+I

NO SURF=NOSURF+NS I*NSK

DO 100 m=I,NSI

i0=39" (m-l) +IFIRST

iI=MIN (j0+39, ILAST)

ILm=il-i0+l

DO 90 n--I,NSK

k0=39 * (n-I) +KFIRST

kI=MIN (k0+39, KLAST)

KLn--kl-k0+l

WRITE (8,*) ILm

WRITE (8, *) KLn

DO 80 i=i0,il

DO 70 k=k0,kl

WRITE (8, 25)

CONTINUE

CONTINUE

CONTINUE

CONTINUE

IERR=0

GOTO 5

X (i, J, k), Y (i, J, k) ,Z (i, J, k)

ELSEIF (IPICK.EQ. 3) THEN

WRITE (*, 5003)

READ (*, *) IPICK2

IF (IPICK2 .EQ. i) THEN

WRITE (*, *)

WRITE(*,*) 'Please

READ (*, *) K

IFIRST=I

ILAST=IL

JFIRST=I

JLAST=JL

ELSE

WRITE (*, *)

WRITE(*,*) 'Please

READ (*, *) K

WRITE (*, *)

enter the value of k'

enter the value of k'

WRITE(*,*)' Please enter the lower and upper limit'

Z_ITE(*,*) ' for the i coordinate (IFIRST, ILAST)'

READ (*, *) IFIRST, ILAST

WRITE (*, *)

WRITE(*,*)' Please enter the lower and upper limit'

_ITE(*,*)' for the j coordinate (JFIRST, JLAST)'

READ(*,*)JFIRST, JLAST

ENDIF

NSI=(ILAST-IFIRST)/39

IF(IFIRST+NSI*39.LT.ILAST)NSI=NSI+I

NSJ=(JLAST-JFIRST)/39

IF(JFIRST+NSJ*39.LT.JLAST)NSJ=NSJ+I

NOSURF=NOSURF+NSI*NSJ

DO 140 m:I,NSI

i0=39*(m-I)+IFIRST

26
ii0

120

130

140

999

i000

C

C

25

5001

5002

50O3

iI=MIN(j0+39,ILAST)

ILm=il-i0+l

DO 130 n=I,NSJ

j0=39*(n-I)+JFIRST

jI=MIN(j0+39,JLAST)

JLn=jl-j0+l

WRITE(8,*) ILm

WRITEI8,*) JLn

DO 120 i=i0,il

DO ii0 j=j0, jl

WRITE(8,25) X(i,j,K),Y(i,j,K),Z(i,j,K)

CONTINUE

CONTINUE

CONTINUE

CONTINUE

IERR=0

GOTO5

ELSEIF(IPICK.EQ.0 .OR. IERR.EQ.I)THEN

REWIND(8)

WRITE(9, *) NOSURF

DO i000 N=I,NOSURF

READ(8, *) II

READ(8, *) I2

WRITE(9, *) II

WRITE(9, *) I2

DO 999 I=l,II*I2

READ(8, *) Xl, X2,X3

WRITE(9, 25)Xl, X2,X3

CONTINUE

CONTINUE

ELSEIF(IERR.EQ.0)THEN

IERR=I

WRITE(*,*)' INVALID SELECTION'

WRITE(*,*)' You will get one more chance to make a'

WRITE(*,*)' selection - Please enter any character'

WRITE(*,*)' and then press RETURN'

READ(*,*)

GOTO5

ENDIF

RETUP_

FORMAT(IX, FI0.6,3x, FI0.6, 3X,FI0.6}

FORMAT(' ',///,

$' Reading in the grid. Please wait',///)

FORMAT(' ',///,' Please enter:',//,

$' 1 - if you want to save a constant-i surface',/,

$' 2 - if you want to save a constant-j surface',/,

$' 3 - if you want to save a constant-k surface',/,

$' 0 - if you want to QUIT',/,

$' (Recall: IL=',I3, ', JL=', I3, ' KL=',I3, ') ',//)

FORW-AT( ' ',///, ' Please enter: ',//,

$' i - if you want the whole surface saved',/,

$' 2 - if you want to specify a part of the',/,

$' surface to be saved',//)

END

27
A.3 Listing of 3DPREP

c

c

c

c

c

c

c

PROGRAM PREPARE

DIMENSION TENSION (16) ,BETA1 (16) ,BETA2 (16) ,

/ Xl (i00) , Y1 (100) , Z1 (100) ,X(16,100) , Y (16, I00) , Z (16, I00)

INTEGER STRTYPE(16),N2B(4),TYPE(16),NODES(16)

REAL KXII, KXI2, KETAI, KETA2, KZETAI, KZETA2

WRITE (*, * )

WRITE (*, * )

WRITE (*, *)

WRITE (*, * )

WRITE (*, *) '

WRITE(*, *) '

WRITE (*, *) '

WRITE (*, *)

WRITE (*, * )

This program prepares input files for GRID3D by'

reading the necessary information from the screen'

and from files.'

WRITE(*,*)'What technique is to be used.'

WRITE(*,*)'Enter 2 for the two-boundary technique'

WRITE(*,*)' or 4 for the four-boundary technique'

READ(*,*)ITECH

WRITE(*,*) 'Enter IL, JL and KL'

READ (*, *) IL, JL, KL

WRITE(*,*) 'Enter SigmaXi, SigmaEta, and SigmaZeta'

READ (*, *) SigXi, SigEt, SigZt

WRITE(*,*) 'Enter kXIl and kXI2'

READ (*, *) KXII, KXI2

WRITE(*,*) 'Enter kETAI and kETA2'

READ (*, * )KETAI, KETA2

WRITE(*,*)'Enter kZETA! and kZETA2'

READ (*, *) KZETAI, KZETA2

WRITE(*,*) 'The

WRITE (20, *) ITECH, '

WRITE (20, *) IL, '

WRITE (20, *) JL, '

WRITE(20,*)KL, '

WRITE (20, *) SigXi, '

WRITE (20, *

WRITE (20, *

WRITE (20, *

WRITE (20, *

WRITE (20, *

WRITE (20, *

WRITE (20, *

_RITE (20, *

output file will be UNIT 20'

Technique'

IL'

JL'

KL '

SigEt,'

SigZt,'

KXII, '

KXI2, '

KETAI, '

KETA2, '

KZETAI,'

KZETA2,'

SigmaXi'

SigmaEta'

SigmaZeta'

kXIl'

kXI2'

kETAI'

kETA2'

kZETAI'

kZETA2'

DO 200 NSRF=I,2

NEI=I+4* (NSP_C-I)

28
C

C

C

C

5

C

I0

C

C

C

C

2O

C

C

C

C

C

NE2=NEI+I

NE3=NEI+2

NE4=NEI+3

Get data for edge NEI:

/

!

Get

WRITE(*,2002)NEI

WRITE(*,2001)NEI

READ(*,*)ITYPE

IF(ITYPE.NE.I .AND. ITYPE.NE.2)GOTO 5

IF(ITYPE.EQ.I)THEN

CALL GETGRP(XI,YI,ZI,NOP)

IF(NOP.NE.KL)WRITE(*,*)

'WARNING --- NUMBER OF GRID POINTS

' EDGE',NEI

WRITE(20,3001)ITYPE,NEI

DO 10 K=I,KL

WRITE(20,3004)XI(K),YI(K),ZI(K),K

CONTINUE

XlI=XI(1)

YII=YI(1)

ZlI=ZI(1)

XlL=XI(KL)

YIL=YI(KL)

ZIL=ZI(KL)

INCONSISTENT -',

ELSEIF (ITYPE.EQ. 2) THEN

CALL GETNODES(XI, YI,ZI,NOP)

WRITE(*,2006)NEI

READ(*,*)TENSN

CALL GETSTR(NEI,ISTRI,BETAII,BETA21)

WRITE (20,3001) ITYPE, NEI

WRITE (20,3002) TENSN

WRITE (20,3003)NOP

DO 20 I=I,NOP

WRITE (20,3004)XI(I),YI(I),ZI(I),I

CONTINUE

WRITE (20,3005)ISTRI

IF (ISTRI .NE. 4)WRITE (20, 3006) BETA11

IF (ISTRI .EQ. 4) WRITE (20,3007) BETA11, BETA21

Xll=Xl (I)

YII=YI (I)

ZlI=ZI (I)

XIL=XI (NOP)

YIL=YI (NOP)

ZlL=ZI(NOP)

ENDIF

data for edge NE2:

WRITE(*,2002)NE2

2g
25

C

30

C

C

C

C

4O

C

C

C

C

C

45

C

WRITE (*, 2001)NE2

READ (*, * )ITYPE

IF(ITYPE.NE.I .AND. ITYPE.NE.2)GOTO 25

IF (ITYPE.EQ. I) THEN

CALL GETGRP(XI,YI,ZI,NOP)

IF (NOP .NE. KL) WRITE (*, *)

'WARNING --- NUMBER OF GRID POINTS

' EDGE',NE2

WRITE (20,3001) ITYPE, NE2

DO 30 K=I,KL

WRITE (20, 3004) Xl (K) ,Y1 (K) ,Zl (K) ,K

CONT INU E

X21=Xl (i)

Y21=YI (i)

Z21=ZI (I)

X2L=Xl (KL)

Y2L=YI (KL)

Z2L=Zl (KL)

INCONSISTENT -',

ELSEIF(ITYPE.EQ.2)THEN

CALL GETNODES(XI,YI,ZI,NOP)

WRITE (*, 2006)NE2

READ (*, * )TENSN

CALL GETSTR(NEI,ISTR2,BETAI2,BETA22)

WRITE (20,3001) ITYPE, NE2

WRITE (20,3002) TENSN

WRITE (20,3003) NOP

DO 40 I=I,NOP

WRITE (20,3004)XI(I),YI(I),ZI(I),I

CONTINUE

WRITE (20,3005) ISTR2

IF (ISTR2 .NE. 4) WRITE (20,3006) BETA12

IF (ISTR2 .EQ. 4) WRITE (20,3007) BETA12, BETA22

X21=Xl (i)

Y21=YI (i)

Z21=ZI (I)

X2L=Xl (NOP)

Y2L=YI (NOP)

Z2L=Zl (NOP)

ENDIF

Get data for edge NE3:

WRITE (*, 2002) NE3

WRITE (*, 2003) NE3

READ (*, *) ITYPE

IF (ITYPE.NE. 1 .AND. ITYPE.NE.2 .AND. ITYPE.NE. 3) GOTO 45

IF(ITYPE.EQ.I)THEN

CALL GETGRP(XI,YI,ZI,NOP)

IF(NOP.NE.IL)WRITE(*,*)

'WA_ING NUMBER OF GRID POINTS INCONSISTENT -',

30
5O

C

C

C

6O

C

C

C

' EDGE',NE3

WRITE (20,3001) ITYPE, NE3

DO 50 I=I,IL

WRITE (20,3004)XI(I),YI(I),ZI(I),I

CONTINUE

ELSEIF (ITYPE .EQ. 2) THEN

CALL GETNODES(XI, YI,ZI,NOP)

WRITE(*, 2006)NE3

READ (*, * )TENSN

CALL GETSTR (NE3, ISTR3, BETA13, BETA23)

WRITE (20, 300 i) ITYPE, NE3

WRITE (20,3002) TENSN

WRITE (20, 3003) NOP

DO 60 I=I,NOP

WRITE (20,3004)XI(I),YI(I),ZI(I),I

CONTINUE

WRITE (20,3005) ISTR3

IF (ISTR3 .NE. 4)WRITE (20, 3006) BETA13

IF (ISTR3. EQ. 4) WRITE (20,3007) BETA13, BETA23

ELSEIF (ITYPE.EQ. 3) THEN

ITYPE=2

NOP=2

Ii=l

I2=2

WRITE (*, 2006)NE3

READ (*, * )TENSN

CALL GETSTR (NE3, ISTR3, BETA13, BETA23)

WRITE (20, 3001) ITYPE, NE3

WRITE (20,3002) TENSN

WRITE (20, 3003)NOP

WRITE (20,3004) XII, YII, ZII, II

WRITE (20,3004) X21, Y21, Z21, I2

Xl (1)=Xll

X1 (2)=X21

Y1 (1 )=YI 1

Y1 (2)=Y21

Z1 (1)=ZII

Z1 (2)=Z21

WRITE (20,3005) ISTR3

IF (ISTR3.NE. 4)WRITE (20, 3006)BETA13

IF (ISTR3.EQ. 4) WRITE (20,3007) BETA13, BETA23

ENDIF

TYPE (8+NSRF)=ITYPE

NODES (8+NSRF )=NOP

TENS ION (8+NSRF )=TENSN3

STRTYPE (8+NSRF) =ISTR3

BETA1 (8+_{SRF) =BETA13

BETA2 (8+NSRF)=BETA23

IF (ITYPE. EQ. i) IMAX=IL

IF (ITYPE. EQ. 2) IMAX=NOP

DO 70 I=l, IF_AX

31
70

C

C

C

C

75

8O

C

C

C

C

9O

C

Get

X(8+NSRF,I)=Xl (I)

Y(8+NSRF,I) =YI (I)

Z(8+NSRF,I)=ZI (I)

CONTINUE

data for edge NE4:

WRITE(*, 2002)NE4

WRITE(*, 2003)NE4

READ(*, *) ITYPE

IF(ITYPE.NE.I .AND. ITYPE.NE.2 .AND. ITYPE.NE.3) GOTO 75

IF (ITYPE.EQ. I) THEN

CALL GETGRP(XI,YI,ZI,NOP)

IF (NOP.NE. IL) WRITE (*, *)

'WARNING --- NUMBER OF GRID POINTS

' EDGE',NE4

WRITE (20, 3001) ITYPE, NE4

DO 80 I=I,IL

WRITE (20,3004) X1 (K), Y1 (K), Z1 (K), K

CONTINUE

INCONSISTENT -'

ELSEIF (ITYPE.EQ. 2) THEN

CALL GETNODES(×I,Y!,ZI,NOP)

WRITE (*, 2006)NE4

READ (*, * )TENSN

CALL GETSTR(NE4,ISTR4,BETAI4,BETA24)

WRITE (20,3001) ITYPE, NE4

WRITE (20,3002) TENSN

WRITE (20, 3003)NOP

DO 90 I=I,NOP

WRITE (20,3004)XI(I),YI(I),ZI(I),I

CONTINUE

WRITE (20,3005)ISTR4

IF (ISTR4 .NE. 4) WRITE (20,3006) BETA14

IF (ISTR4 .EQ. 4) WRITE (20,3007) BETA14, BETA24

ELSEIF (ITYPE.EQ. 3) THEN

ITYPE=2

NOP = 2

II=l

I2=2

WRITE (*, 2006)NE4

READ (*, *) TENSN

CALL GETSTR(NE4,ISTR4,BETAI4,BETA24)

WRITE(20,3001) ITYPE,NE4

WRITE (20,3002) TENSN

WRITE (20,3003) NOP

WRITE (20,3004)X!L, YIL, ZIL, II

WRITE (20,3004)X2L, Y2L, Z2L, I2

X1 (1)=XIL

X1 (2) =X2L

Y1 (i)=YiL

Y1 (2)=Y2L

32
C

C

I00

C

C

200

C

C

201

C

C

C

C

C

C

210

zl (1)=ZIL

Zl (2) =Z2L

WRITE (20, 3005) ISTR4

IF (ISTR4 .NE. 4)WRITE (20,3006) BETA14

IF (ISTR4 .EQ. 4) WRITE (20,3007) BETA14, BETA24

ENDIF

TYPE (12+NSRF)=ITYPE

NODES (12+NSRF) =NOP

TENSION (12+NSRF) =TENSN4

STRTYPE (12+NSRF) =ISTR4

BETAI (12 +NSRF) =BETA14

BETA2 (12 +NSRF) =BETA24

IF (ITYPE.EQ. I) IMAX=IL

IF (ITYPE .EQ. 2) IMAX=NOP

DO I00 I=I,IMAX

X (16+NSRF, I)=Xl (I)

Y (16+NSRF, I)=YI (I)

Z (16+NSRF, I)=ZI (I)

CONTINUE

CONTINUE

IF (ITECH.EQ. 2) THEN

WRITE(*,' (////////)')

WRITE(*,*) ' Specify

WRITE(*,*) ' ii, 12,

WRITE(*, ' (////) '}

N2B(1)=II

N2B(2}=I2

N2B(3)=I5

N2B(4)=I6

DO 201 IE=I,4

CALL GETSTR(N2B(IE),ISTR, BETAIIE,BETA2IE)

WRITE (20,3010) ISTR,N2B(IE)

IF (ISTR.NE. 4) WRITE (20, 3006) BETAIIE

IF (ISTR.EQ. 4) WRITE (20,3007) BETAIIE, BETA21E

CONTINUE

stretching parameters for

15 and 16'

ELSEIF (ITECH.EQ. 4) THEN

DO 500 NSRF=3,4

NEI=I+4*(NSRF-I)

NE2=NEI+I

NE3=NEI+2

NE4=NEI+3

Write data for edge NEI:

WRITE (2,0,3001) TYPE (NEI), NEI

IF (TYPE(NE!) .EQ.I)THEN

DO 210 I=I,IL

WRITE (20,3004)X(NEI, I),Y(NEI,I),Z(NEI,I),I

CONTINUE

edges'

33
C

C

220

C

C

C

C

C

230

C

C

240

C

C

C

C

C

3O5

C

ELSEIF (TYPE (NEI) .EQ, 2) THEN

WRITE (20, 3002)TENSION (NEI)

WRITE (20, 3003) NODES (NEI)

DO 220 I--I,NODES(NEI)

WRITE (20, 3004) X (NEI, I), Y (NEI, I), Z (NEI, I), I

CONTINUE

WRITE (20, 3005)STRTYPE (NEI)

IF (STRTYPE (NEI) .NE. 4) WRITE (20, 3006) BETA1 (NEI)

IF (STRTYPE (NEI) .EQ. 4)

WRITE (20, 3007)BETA1 (NEI), BETA2 (NEI)

ENDIF

Write data for edge NE2:

WRITE (20,3001) TYPE (NE2) ,NE2

IF (TYPE (NE2) .EQ. I) THEN

DO 230 I=I,IL

WRITE (20, 3004) X (NE2, I), Y (NE2, I), Z (NE2, I), I

CONT INUE

ELSEIF (TYPE (NE2) .EQ. 2) THEN

WRITE (20, 3002) TENSION (NE2)

WRITE (20, 3003) NODES (NE2)

DO 240 I=i,NODES(NE2)

WRITE (20,3004)X (NE2, I), Y (NE2, I), Z (NE2, I), I

CONT INU E

WRITE (20,3005) STRTYPE (NE2)

IF(STRTYPE(NEI) .NE.4)WRITE(20,3006)BETAI (NE2)

IF (STRTYPE (NEI) .EQ. 4)

WRITE (20,3007) BETA1 (NE2), BETA2 (NE2)

ENDIF

Get data for edge NE3:

WRITE (*, 2002) NE3

WRITE(*,*) Enter

WRITE *, *) 1

WRITE * I *) 2

WRITE *, *) 3

WRITE *f *)

WRITE *, *) 4

WRITE * ¢ *)

WRITE * *t

WRITE * *

Z<R!TE (* r *

WRITE (* *

_,qR!TE (* f *

WRITE( _,* 'Enter your

READ (*, * )ITYPE

the following:'

if grid points are specified'

if nodes for splining are specified'

if you want to let the end points of'

edges already entered define the curve'

if you want to use the end points of'

edges already defined but add (by '

typing in directly) some points in '

between'

choice'

Z4
C

310

C

C

C

C

320

C

C

C

C

C

IF(ITYPE.NE.I .AND. ITYPE.NE.2 .AND. ITYPE.NE.3

.AND. ITYPE.NE.4)GOTO 305

IF (ITYPE.EQ. i) THEN

CALL GETGRP(XI,YI, ZI,NOP)

IF (NOP. NE. JL) WRITE (*, *)

'WARNING --- NUMBER OF GRID POINTS INCONSISTENT

' EDGE',NE3

WRITE (20,3001) ITYPE, NE3

DO 310 J=I,JL

WRITE (20, 3004) Xl (J) ,Y1 (J) ,Z1 (J) ,J

CONTINUE

ELSEIF(ITYPE.EQ.2)THEN

CALL GETNODES(XI,YI,ZI,NOP)

WRITE (*, 2006)NE3

READ (*, *) TENSN

CALL GETSTR(NE3,ISTR3,BETAI3,BETA23)

WRITE (20, 3001) ITYPE, NE3

WRITE (20, 3002) TENSN

WRITE (20,3003) NOP

DO 320 I=I,NOP

WRITE (20, 3004)XI(I),YI(I),ZI(I),I

CONTINUL

WRITE (20,3005) ISTR3

IF (ISTR3 .NE. 4)WRITE (20, 3006) BETA13

IF (ISTR3.EQ. 4)WRITE (20, 3007)BETA13, BETA23

ELSEIF(!TYPE.EQ.3)THEN

ITYPE=2

NOP=2

II=I

!2=2

WRITE (*, 2006)NE3

READ (*, *) TENSN

CALL GETSTR(NE3, ISTR3,BETAI3,BETA23)

WRITE (20,3001) ITYPE, NE3

WRITE (20,3002)TENSN

WRITE (20,3003) NOP

9_RITE (20,3004) X (NEI, I), Y (NEI, I), Z (NEI, I), Ii

9_qRITE (20,3004) X (NE2, i), Y (NE2, I), Z (NE2, i), I2

WRITE (20,3005) ISTR3

IF (ISTR3.NE. 4)WRITE (20,3006)BETA13

IF (ISTR3. EQ. 4) WRITE (20,3007) BETA13, BETA23

ELSEIF(!TYPE.EQ.4)THEN

WRITE(_,*) ' How many nodes do

READ (*, *) NO._

DO 330 I=I,NOP

you want to add?'

D!

35
330

C

340

C

C

C

C

C

4O5

C

C

C

C

410

C

CALL GETNEWND(XI,YI,ZI,I)

CONTINUE

ITYPE=2

NOP=NOP+2

II=l

WRITE(*, 2006)NE3

READ(*, *)TENSN

CALL GETSTR(NE3,ISTR3,BETA13,BETA23)

WRITE(20,3001) ITYPE,NE3

WRITE(20, 3002)TENSN

WRITE(20, 3003)NOP

WRITE(20,3004) X(NEI, I), Y(NEI, i), Z(NEI, i), Ii

DO340 I=2,NOP-I

WRITE(20, 3004)Xl (I-l), Y1(I-l), Z1(I-l), I

CONTINUE

WRITE(20, 3004)X(NE2,I), Y(NE2,I), Z(NE2,i) ,NOP

WRITE(20,3005) ISTR3

IF (ISTR3.NE.4)WRITE(20, 3006)BETA13

IF (ISTR3.EQ.4)WRITE(20, 3007)BETA13,BETA23

ENDIF

Get data for edge NE4:

WRITE* 2002)NE4

WRITE*,*) Enter

WRITE*, *) 1

WRITE(*, *) 2

WRITE(*,*) 3

WRITE(*, *)

WRITE(*,*) 4

WRITE(*, *)

WRITE(*, *)

WRITE(*, *)

WRITE(*, *)

WRITE(*, *)

WRITE(*, *)

Z_RITE(*, *) 'Enter your choice'

READ(*, *) ITYPE

the following:'

if grid points are specified'

if nodes for splining are specified'

if you want to let the end points of'

edges already entered define the curve'

if you want to use the end points of'

edges already defined but add (by '

typing in directly) some points in '

between'

IF(ITYPE.NE.I .AND. ITYPE.NE.2 .AND. ITYPE.NE.3

.AND. ITYPE.NE.4)GOTO405

IF (ITYPE.EQ.I) THEN

CALL GETGRP(XI,YI,ZI,NOP)

IF (NOP.NE.JL)Z_RITE(*, *)

'WARNING--- NUMBEROF GRIDPOINTSINCONSISTENT

' EDGE',NE4

WRITE(20, 3001)ITYPE,NE4

DO 410 J=!,JL

WRITE(20,3004)XI(J),YI(J),ZI(J),J

COh_TINUE

_ !

36
C

C

C

420

C

C

C

C

C

C

430

C

C

ELSEIF(ITYPE.EQ.2)THEN

CALL GETNODES(XI,YI,ZI,NOP)

WRITE(*, 2006)NE4

READ(*, *)TENSN

CALL GETSTR(NE4,ISTR4,BETA14,BETA24)

WRITE(20, 3001)ITYPE,NE4

WRITE(20, 3002)TENSN

WRITE(20, 3003)NOP

DO 420 I=I,NOP

WRITE(20, 3004)X1(I), Y1(I), Z1(I), I

CONTINUE

WRITE(20, 3005)ISTR4

IF (ISTR4.NE.4)WRITE(20, 3006)BETA14

IF (ISTR4.EQ.4)WRITE(20, 3007)BETA14,BETA24

ELSEIF(ITYPE.EQ.3)THEN

ITYPE=2

NOP=2

Ii=I

12=2

WRITE(*, 2006)NE4

READ(*, *) TENSN

CALL GETSTR(NE4,ISTR4,BETAI4,BETA24)

WRITE(20, 3001) ITYPE,NE4

WRITE(20,3002) TENSN

WRITE(20, 3003)NOP

IMAX=IL

IF (TYPE(NEI) .EQ.2) IMAX=NODES(NEI)

WRITE(20, 3004)X (NEI, IMAX), Y(NEI, IMAX), Z(NEI, IMAX), Ii

IMAX=IL

IF (TYPE(NE2).EQ.2) IMAX=NODES(NE2)

WRITE(20,3004) X(NE2,IMAX),Y(NE2,IMAX), Z(NE2,IMAX), I2

WRITE(20,3005)ISTR4

IF (ISTR4.NE.4) WRITE (20, 3006) BETA14

IF (ISTR4 .EQ. 4)WRITE (20,3007) BETA14, BETA24

ELSEIF (ITYPE.EQ. 4) THEN

WRITE(*,*) ' How many nodes do

READ (*, * )NOP

DO 430 I=I,NOP

CALL GETNEWND(XI,YI,ZI, I)

CONTINUE

ITYPE=2

NOP=NOP+2

II=l

you want to add?'

WRITE (*, 2006)_E4

READ (*, *) TE_4SN

CALL GETSTR(NE4,ISTR4,BETAI4,BETA24)

WRITE (20,3051) ITYPE, NE4

37
44O

C

5OO

C

2001

2002

2003

2006

C

3001

3002

3003

3004

3005

3006

3007

3010

C

C

C

C

WRITE (20,3002) TENSN

WRITE (20, 3003) NOP

IMAX=IL

IF (TYPE (NEI) .EQ. 2) IMAX=NODES (NEI)

WRITE (20,3004) X (NEI, IMAX) ,Y (NEI, IMAX) , Z (NEI, IMAX) ,Ii

DO 440 I=2,NOP-I

WRITE (20, 3004) X1 (I-l) ,Y1 (I-l) ,Zl (I-l) ,I

CONTINUE

IMAX=IL

IF (TYPE (NE2) .EQ. 2) IMAX=NODES (NE2)

WRITE (20,3004)X (NE2, IMAX), Y (NE2, IMAX), Z (NE2, IMAX), I2

WRITE (20, 3005) ISTR4

IF (ISTR4 .NE. 4) WRITE (20,3006) BETA14

IF (ISTR4 .EQ. 4)WRITE (20, 3007)BETA14, BETA24

ENDIF

CONTINUE

ENDIF

FORMAT(/////,' Enter the TYPE for edge ',I2,

$///,' Enter:',//,

$' 1 if grid points along the edge are to be ',

$'specified',/,' 2 if nodes for splining are to be',

$' specified',/)

FORMAT(////////////,' SPECIFYING EDGE ',I2,///)

FORMAT(/////,' Enter the TYPE for edge ',I2,

$///,' Enter:',//,

$' 1 if grid points along the edge are to be ',

$'specified',/,' 2 if nodes for splining are to be',

$' specified',/,' 3 if end nodes of edges already ',

$'entered ',/,' define the curve completely',/)

FORMAT(///////,' Enter the TENSION parameter for curve ',I2)

FORMAT ' ',I3,' Type - EDGE NO: ',I2,

$ .... ,)

FORMAT ' ',F6.2, ' Tension parameter')

FORMAT ' ',I3, ' Number of nodes')

FORY_kT (3X, 3 (F8.5, 3X) , ' --', I2)

FORMAT(' ',I3, ' StretchType')

FOP_AT (' ', F8.4, ' Stretching parameter BETA' )

FOR_nkT(' ',2(FS.4,3X),' Stretching parameters BETA1 and',

$ ' BETA2' )

FORMAT(' ',I3, ' StretchType ',I2, ' -')

STOP

END

C

SUBROUTINE GETNODES(XI,YI,ZI,NOP)

C

C This subroutine reads in and arranges nodal points to define

C an edge.

C

DIME_SIO_ XI(i00),YI(100),ZI(100)

C

!NNUMO = iO

38
C

C

C

C

C

I0

C

20

C

3O

C

4O

C

i00

C

200

201

202

C

C

CNOP=0

WRITE(*,*) 'How many sections is the edge composed of?'

READ (*, * )NOSECT

DO i00 ISECT=I,NOS ECT

INNUM = INNUM 0+ ISECT

WRITE (*, 200) ISECT, INNUM

READ (INNUM, *)NP

WRITE (*, 201) NP

READ (*, *)NI

WRITE (*, 202)

READ (*, *) N2

NI2=IABS (N2-NI) +I

NOP=NOP+NI 2

IF(N2.GE.NI)THEN

DO I0 J=I,NI-I

READ(INNUM,*)

CONTINUE

I=NOP-NI 2

DO 20 J=NI,N2

I=I+l

READ (INNUM, *) Xl (I), Y1 (I), Zl (I)

CONTINUE

ELSE

DO 30 J=I,N2-1

READ(INNUM,*)

CONTINUE

I=NOP+I

DO 40 J--N2,NI

I=I-i

READ (INNUM, *) X1 (I), Y1 (I), Zl (I)

CONTINUE

ENDIF

CLOSE(INNUM)

CONTINUE

FO_/AT(' Section ',I2,' will be read in from UNIT',I3)

FORb_T(' There are ',I2,' points in the file.',/

/ ,' Enter the number of the point that is to be the',/

/ ,' the first on the current section.')

FO_:AT( ' Enter the number of the point that is to be the',/

/ , ' the last on the current section.')

RETURN

END

C

SUBROUTINE GETGRP(XI,YI,ZI,NOP)

39
C

C This subroutine reads in grid point coordinates for an

C edge and stores in either forward or reversed order.

C

DIMENSIONXI(100),YI(100),ZI(100)

C

C

C

C

C

C

i0

C

2O

C

6O

C

70

C

5OO

C

200

201

INNUM0=I0

NOP=0

WRITE(*,*) 'HOWmany sections is the edge composed of?'

READ(*, *)NOSECT

DO500 ISECT=I,NOSECT

INNUM=INNUM0+I SECT

WRITE(*, 200)ISECT,INNUM

READ(INNUM,*) NP

WRITE(*, 201)NP

READ(*, *)NI

WRITE (*, 202)

READ(*, *)N2

NI2=IABS(N2-NI)+I

NOP=NOP+NI2

IF(N2.GE.NI)THEN

DO I0 J=I,NI-I

READ(INNUM,*)

CONTINUE

I=NOP-NI2

DO 20 J=NI,N2

l=I+l

READ(INNUM,*)XI(I),YI(I),ZI(I)

CONTINUE

ELSE

DO 60 J=I,N2-1

READ(INNUM,*)

CONTINUE

I=NOP+I

DO 70 J=N2,NI

I=I-I

READ (INNUM, *)Xl (I), Y1 (I), Z1 (I)

CONTINUE

ENDIF

CLOSE(INNUM)

CONTINUE

FORb_T(' Section ',I2,' will be read in from UNIT',I3)

FOR<_T( ' There are ',I2, ' grid points in the file.',/

/ ,' Enter the number of the grid point that is to be the',/

/ ,' the first on the current section.')

4O
202

C

C

C=- -

FORMAT(' Enter the number of the grid point that is to be the',/

/ ,' the last on the current section.')

RETURN

END

C

C

SUBROUTINE GETNEWND (Xl, YI, ZI, I)

C

C This subroutine reads in from the screen new points that are

C to be included on the edge.

C

C

C

2001

C

C

CDIMENSION XI(100),YI(100),ZI(10O)

WRITE (*, 2001) I

READ(*, *)Xl (I) ,YI (I), Zl (I)

FORMAT(' Please enter the x, y, and z coordinates for inter-'

/ ,/, ' mediate point number ', I2)

RETURN

END

C

SUBROUTINE GETSTR (NEI ,ISTR, BETA1, BETA2 )

C

C This subroutine reads information regarding stretching function

C along edge NEI.

C

WRITE (*, 2011) NEI

WRITE (*, 2020)

I00 READ (*, *) ISTR

IF(ISTR.LT.4 .AND. ISTR.GT.0)THEN

WRITE (*, 2031) NEI

READ (*, * )BETAI

ELSEIF (ISTR.EQ. 4) THEN

WRITE (*, 2036) NEI

READ (*, * )BETA1, BETA2

ELSEIF (ISTR. EQ. 0) THEN

BETA1= 1.1

ELSE

WRITE(*,*)' Please enter a number from 0 to 4'

GOTO i00

ENDIF

C

2011

2020

2031

2036

C

C

FORMAT(/, ' Enter:

/ 1 for

/ 2 for

/ 3 for

/

/ 4 for

/

FORMAT(' Enter the STRETCH TYPE for edge ',I2)

0 for no stretching',/,

concentration near lower boundary',/,

concentration near upper boundary',/,

concentration near both boundaries',/,

(one parameter stretching function)',/,

concentration near both boundaries',/,

(two-parameter stretching function)')

FO_,_T(' Enter the stretching parameter (BETA) for edge ',

/ I2)

FO_,_?(' Enuer :he stretchinc parameters (BETA1 and BETA2) ',

/ 'fcr ed:e ',i2)

41
A.4 Listing of EDGE

_.

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

CPROGRAM EDGE

This program generates grid points along an edge which is given I

by a set of discrete nodes. The edge can be made up of up to I

several sections. A parametric tension spline is fit through I

each section. Each section has its own control parameters I

such as number of grid points, tension, and stretching functions. I

NOTE: Subroutines have been adopted from GRID3D (and modified

slightly) to generate the grid points on the curves.

PARAMETERS:

MxBPts - Maximum number of nodes per section

MxGSiz - Maximum number of grid points per section

MxSect - Maximum number of sections

MxBCvs - Should not be modified

I

I

I

I

I

I

I

I

I

I

I

I

C

C

C

C

C

PARAMETER (MxBPts=31, MxGSiz=101, MxSect=5, MxBCvs =I)

DIMENSION x(MxSect,MxBPts), y(MxSect,MxBPts),

$ z (MxSect, MxBPts), zx (MxSect, MxBPts),

$ zy(MxSect,MxBPts), zz(MxSect,MxBPts),

$ s(MxSect,MxBPts), Tensn(MxSect)

DIMENSION Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts)

DIMENSION XB(MxGSiz,MxSect), YB(MxGSiz,MxSect),

$ ZB(MxGSiz,MxSect), StrB(MxGSiz,MxBCvs)

INTEGER NDPts(MxSect), ILS(MxSect), StrTp

WRITE(*,*) ' Running PROGRAM EDGE'

WRITE (*, *)

WRITE (*, *)

WRITE(*,*)' The input data will be read in from UNIT I'

WRITE(*,*)' The output will be written into UNIT 20'

WRITE (*, *)

WRITE (*, *)

READ (I, *)NOSECT

IL=0

DO 200 iS=I,NOSECT

CALL RdSctIn (x, y, z, NDPts, is, Tensn, l,MxBPts,MxSect,

StrTp, Betal,Beta2, ILS(is))

CALL PTSpln (x, y, z, s, zx, zy, zz,Diag,OfDiag,Right,NDPts,

Tensn (is) ,is, MxBPt s, MxSect)

CALL CalcStr2(l,ILS(is),StrTp, Betal,Beta2,

StrB, MxBCvs,MxGSiz)

42
200

300

301

C

C

3001

CALL EdgGPts(is,I,ILS(is),XB,YB, ZB, StrB, x,Y,Z,S,

zx, zy, zz,NDPts,Tensn(is),

MxBCvs,MxBPts,MxGSiz,MxSect)

IL=IL+ILS(is)

CONTINUE

IL=IL-(NOSECT-I)

WRITE(20,*)IL,' Number of grid points'

IP=0

DO 301 is=I,NOSECT

DO 300 i=l,ILS(is)-I

IP=IP+I

WRITE(20,3001)XB(i,is),YB(i,is),ZB(i,is),IP

CONTINUE

CONTINUE

is=NOSECT

i=ILS(NOSECT)

WRITE(20,3001)XB(i, is),YB(i,is),ZB(i, is),IL

FORMAT ( '

STOP

END

',3X,3(F9.6,3X),' ',I3)

C

C..... C

C

SUBROUTINE RdSctIn (x,y,z,NDPts,CrvNum, Tensn, InNum, MxBPts,MxSect,

$ StrTp, Betal, Beta2, IL)

C

C This SUBROUTINE reads in the information concerning discrete points on

C the boundaries. This information is used for generating spline-fitted

C boundary approximation curves.

C

C

C

C

i0

C

C

INTEGER CrvNum, i, NDPts(MxSect), InNum,

REAL x(MxSect,MxBPts), y(MxSect,MxBPts),

$ z(M×Sect,MxBPts), Tensn(MxSect)

READ (InNum, * )

READ (InNum, *)

READ (InNum, *)

IL

Tensn (CrvNum)

NDPts (CrvNum)

DO I0 i=l,NDPts(CrvNum)

READ(InNum,*) x(CrvNum, i), y(CrvNum, i),

CONT INUE

READ (InNum, * )StrTp

IF (StrTp.NE. 4) THEN

READ (InNum, *) Betal

ELSE

READ (InNum, *) Betal, Beta2

ENDIF

RETU_

END

StrTp, IL

z(CrvNum, i)

C

C--

C

SUBROUTINE Ca!cS (x,y,z,s,NDPts,CrvNum, MxBPts,MxSect}

43
C This SUBROUTINEcalculates the spline parameter, s, as

C arc length.

C

INTEGER NDPts(MxSect), CrvNum, i

REAL x(MxSect,MxBPts), y(MxSect,MxBPts),

$ z(MxSect,MxBPts), s(MxSect,MxBPts)

s (CrvNum, i )=0.0

I0

C

RETURN

END

C

C ......

DO I0 i=2,NDPts(CrvNum)

s (CrvNum, i)--s (CrvNum, i-l)

$ +SQRT ( (x (CrvNum, i) -x (CrvNum, i-l) )**2

$ + (y (CrvNum, i) -y (CrvNum, i-l) )**2

$ + (z (CrvNum, i) -z (CrvNum, i-l) )**2)

CONT INUE

an approximate

C

SUBROUTINE SplMat (Diag, OfDiag,Right,w,s,NDPts,T, CrvNum,

tension

$ MxBPts,MxSect)

C

C This SUBROUTINE forms the parametric

C particular boundary curve data set.

C

C

C

spline matrix

INTEGER i, NDPts(MxSect), CrvNum

REAL Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts),

$ w(MxSect,MxBPts), s(MxSect,MxBPts), T, h, hm

Diag (I) =!. 0

Of Diag (I) =0.0

Right (i) =0.0

for a

C

$

i0

c

C

DO i0 i=2,NDPts(CrvNum)-I

h=s (CrvNum, i+l) -s (CrvNum, i)

hm=s (CrvNum, i)-s (CrvNum, i-l)

Diag (i) = (T*COSH (T*hm)/SINH (T*hm) -I/hm+T*COSH (T'h)/SINH (T'h)

-I/h)/T**2

ofmiag(i)=(i/h-T/SINH(T*h))/T**2

Right(i) = (w(CrvNum, i+l)-w(CrvNum, i))/h

- (w (CrvNum, i) -w (CrvNum, i-l) )/hm

CONTINUE

Diag (NDPts (CrvNum)) =i. 0

OfDiag (NDPt s (CrvNum) -i) =0.0

Right (NDPts (CrvNum)) =0.0

RETU_{

END

SUBROUTINE SplSlv

$

This SUBROUTINE solves

(Diag,OfDiag,Right,Derv2,NDPts,CrvNum,

MxBPts,MxSect)

the diagonally dominant parametric tension

44
C spline matrix for a given data set using the Gauss-Seidel iteration.

C Convergence is assumed after 20 iterations.

C

INTEGER i, j, NDPts(MxSect), CrvNum

C

20

3O

C

RETURN

END

C

C .......

REAL Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts),

$ Derv2(MxSect,MxBPts)

C

C Initialize the second derivative matrix to all zeroes.

C

DO I0 i=l,NDPts(CrvNum)

Derv2(CrvNum, i)=0.0

i0 CONTINUE

C

C Calculate the second derivative values using 20 iterations of

C the Gauss-Seidel method.

C

DO 30 j=l,20

DO 20 i=2,NDPts(CrvNum)-i

Derv2 (CrvNum, i)= (Right (i) -OfDiag (i) *Derv2 (CrvNum, i+l)

$ -OfDiag (i-l) *Defy2 (CrvNum, i-l) )

$ /Diag (i)

CONTINUE

CONTINUE

FUNCTION SplVal (s,w,Derv2,sval,T,n,CrvNum, MxBPts,MxSect)

C

C This real function finds the w-value (x-value or y-value) corresponding

C to a specified s-value using the parametric tension spline curve

C generated for a particular boundary curve data set.

C

INTEGER n, CrvNum

C

C

C

C

C

C

REAL s(MxSect,MxBPts), w(MxSect,MxBPts), Derv2(MxSect,MxBPts),

$ sval, T, h, Interim, Templ, Temp2

Templ=sval-s (CrvNum, n)

h=s (CrvNum, n+l) -s (CrvNum, n)

Temp2=s (CrvNum, n+l )-sval

Interim=Derv2 (CrvNum, n) /T**2*SINH (T*Temp2)/SINH (T'h)

$ + (w (CrvNum, n) -Derv2 (CrvNum, n)/T**2) *Temp2/h

SplVal=Interim+Derv2 (CrvNum, n+l )/T* *2 *S INH (T* Templ )

$ /SINH (T'h) + (w (CrvNum, n+l)

$ -Derv2 (CrvNum, n+l )/T**2} *Templ/h

RETUF2_

END

u._.,: PTSp!n(x,y,z,s,XDerv2 YDerv2,ZDerv2,Diag,OfDiag,

$ Right,NDPts,Tensn,CrvNum, MxBPts,MxSect)

C

C This SUBROUTINE forms the main routine for the parametric tension

45
C spline process.

C

C

C

C

C

INTEGER NDPts(MxSect), CrvNum

REAL

$

$

$

$

Diag (MxBPts), OfDiag (MxBPts), Right (MxBPts),

XDerv2 (MxSect,MxBPts), YDerv2 (MxSect,MxBPts),

ZDerv2(MxSect,MxBPts), Tensn,

x (MxSect, MxBPt s), y (MxSect, MxBPt s),

z(MxSect,MxBPts), s(MxSect,MxBPts)

CALL

CALL

SCALL

CALL

$

CALL

CALL

$

CALL

CalcS(x,y,z,s,NDPts,CrvNum, MxBPts,MxSect)

SplMat(Diag, OfDiag, Right,x,s,NDPts,Tensn,CrvNum,

MxBPts,MxSect)

SplSlv(Diag, OfDiag, Right,XDerv2,NDPts,CrvNum, MxBPts,MxSect)

SplMat(Diag,OfDiag, Right,y,s,NDPts,Tensn,CrvNum,

MxBPts,MxSect)

SplSlv(Diag, OfDiag,Right,YDerv2,NDPts,CrvNum, MxBPts,MxSect)

SplMat(Diag, OfDiag, Right,z,s,NDPts,Tensn,CrvNum,

MxPts,MxSect)

SplSlv(Diag,OfDiag,Right, ZDerv2,NDPts,CrvNum, MxBPts,MxSect)

RETURN

END

C

i0

C

C

SUBROUTINE SplInt(n,s,SValue,NDPts,CurCrv,MxBPts,MxSect)

This SUBROUTINE finds the proper interval in which a point on a specified

boundary lies. The interval indicates which initial data points the

point in question lies between and thus which spline coefficients to

use.

INTEGER i, n, CurCrv, NDPts(MxSect)

REAL Temp, SValue,

n=l

i=NDPt s (CurCrv)

s(MxSect,MxBPts)

IF ((n. EQ.I) .AND. (i.GT.I)) THEN

I=I-I

Temp=SValue-s (CurCrv, i )

IF (Temp.GT. 0.0)

n=i

ENDIF

THEN

GOTO I0

ENDIF

RETURN

END

C

SUBROUTI_,_E FAINew(AiNew,Alpha,B, Str)

46
C This SUBROUTINE computes the new Alpha value after stretching as

C AiNew. Alpha is a dummy variable representing either Xi, Eta or Zeta.

C

INTEGER Str

C

C

C

C

C

C

C

----.

REAL Alpha, Templ, Temp2, B2, AiNew, B

AiNew=Alpha

Templ=(B+l)/(B-l)

IF (Str.EQ.l) THEN

Temp2=Templ**(l-Alpha)

AiNew=((B+l)-(B-l)*Temp2)/(Temp2+l)*l

ENDIF

IF (Str.EQ.2) THEN

B2=0

Temp2=Templ** ((Alpha-B2) / (I-B2))

AiNew= ((B+2*B2) *Temp2-B+2*B2) / ((2"B2+I) * (l+Temp2))

ENDIF

IF (Str.EQ.3) THEN

B2=0.5

Temp2=Templ**((Alpha-B2)/(l-B2))

AINew=((B+2*B2)*Temp2-B+2*B2)/((2*B2+I)*(I+Temp2))

ENDIF

RETURN

END

SUBROUTINE CalcStr2(EdgNum, NGPts, StrTp,Betal,Beta2,

$ StrB,MxBCvs,MxGSiz)

C

C This subroutine calculates the distribution function base on the

C stretching parameters 'StrTp' and 'Beta'

C

INTEGER NGPts, StrTp, EdgNum, i

REAL StrB(MxGSiz,MxBCvs), Betal, Beta2, A, B, DZ

I0

2o

C

StrB (i, EdgNum) =0.

IF (StrTp. LE. 3) THEN

DO i0 i=I,NGPts-I

Alpha = (i-l.) / (NGPts-I.)

CALL FAiNew(AiNew,Alpha,Betal,StrTp)

StrB (i, EdgNum) =AiNew

CONTINUE

ELSEIF (StrTp. EQ. 4) THEN

CALL Str4Prm(Betal,Beta2,A,B, DZ)

DO 20 i=2,NGPts-I

Alpha = (i-l.) / (NGPts-I.)

CALL Str4 (AINew, AIpha,A,B, DZ)

StrB (i, EdgNum) =AiNew

CONTINUE

ENDIF

StrB (NGPt s, EdgNum) =i.

47
C

Co

C

C

C

C

C

C

C

C

RETURN

END

SUBROUTINE

$

$

===============================================

EdgGPts (CrvNum, EdgNum, NGPts, XB, YB, ZB, StrB,

x,y,z,s,zx, zy, zz,NDPts, Tensn,

MxBCvs,MxBPts,MxGSiz,MxSect)

This subroutine calculates the grid point location along an edge

based on a spline curve fitted through specified nodal points and a

given distribution function.

INTEGER CrvNum, EdgNum, NGPts, NDPts(MxSect), i, n

REAL

$

$

$

XB(MxGSiz,MxSect), YB(MxGSiz,MxSect), ZB(MxGSiz,MxSect),

StrB(MxGSiz,MxBCvs), x(MxSect,MxBPts), y(MxSect,MxBPts),

z(MxSect,MxBPts), zx(MxSect,MxBPts), zy(MxSect,MxBPts),

zz(MxSect,MxBPts), s(MxSect,MxBPts), Tensn

SRa=S (CrvNum, NDPts (CrvNum))

I0

C

RETURN

END

C

-- ......................

DO i0 i=I,NGPts

SB=SRa*StrB (i, EdgNum)

CALL Splint (n,s, SB,NDPts,CrvNum, MxBPts,MxSect)

XB (i, CrvNum) =SplVal (s, x, zx, SB, Tensn, n, CrvNum, MxBPts, MxSect)

YB (i, CrvNum) =SpIVal (s, y, zy, SB, Tensn, n, CrvNum, MxBPts, MxSect)

ZB (i, CrvNum) =SplVal (s, z, zz, SB, Tensn, n, CrvNum, MxBPt s, MxSect)

CONTINUE

C

C

SUBROUTINE Str4Prm(S0,SI,A,B, DZ)

REAL SO, Sl, A, B, DZ, Y, PI

C

C This subroutine calculates the parameters A, B, and DZ for the twoC sided Vinokur stretching function.

C

C

C

PI=ACOS (-i.)

A=SQRT(S0/SI)

B=SQRT(S0*SI)

IF (B.GT. 1 .001) THEN

IF (B.LE. 2 .7829681) THEN

Y=B-I

DZ=SQRT (6. *Y) * (i .-0.15*Y+0. 057321429* (Y**2}

$ -0. 024907295* (Y**3) +0. 0077424461" (Y**4)

$ -0. 0010794123* (Y**5))

ELSEIF (B. GT. 2. 7829681) THEN

V=LOG (B)

W=Io/B - 0.028527431

DZ=V* (i. +!./V) *LOG (2. *V) -0. 02041793+0. 24902722"W

S tl.9496443*(W**2)-2.6294547*(W**3)+8.56795911*(W**4)

ENDIY

ELSE IF (B.LT. 0 .999) THEN

48
C

C

C

C--

IF (B. LE. 0. 26938972) THEN

DZ=PI* (i .-B+B**2- (i.+ (PI**2)/6. )* (B**3) +6. 794732* (B**4)

$ -13. 205501* (B**5) +II. 726095* (B**6))

ELSE

Y=B-I

DZ=SQRT (6.*Y) * (I .+0.15*Y+0. 057321429* (Y**2)

$ +0. 048774238* (Y*'3)-0.053337753" (Y**4)

$ +0. 075845134* (Y**5))

ENDIF

ENDIF

RETURN

END

SUBROUTINE Str4(AiNew, Alpha,A,B,DZ)

REAL AiNew, Alpha, A, B, DZ, U, T

This subroutine calculates the value of the two-sided Vinokur

stretching function based on the value of the parameters A, B,

and DZ, and on the value of the "computational" coordinate Alpha.

IF (B.GT. 1. 001) THEN

U=0.5+TANH (DZ* (Alpha-0.5)) / (2.*TANH (DZ/2.) )

ELSEIF (B.LT. 0. 999) THEN

U=0.5+TAN (DZ* (Alpha-0.5)) / (2. *TAN (DZ/2.) )

ELSE

U=Alpha * (I .+2. * (B-l) * (Alpha-0.5) * (1-Alpha))

ENDIF

T=U/ (A+ (I.-A) *U)

AINew=T

RETURN

END

4g
A.5 Listing of EDGPREP

C

C

C

C

PROGRAM EDGPREP

PARAMETER (MxNode=100)

INTEGER ISTR, ILS,NOP

REAL TENSION,BETAI,BETA2,

/ Xl(MxNode),Yl(MxNode),Zl(MxNode)

C This program prepares input files for EDGE by reading the

C necessary information from the screen and from files.

C

WRITE(*,*)

WRITE(*,*)

WRITE(*,*)

WRITE(*,*)'This program prepares an input file for the program'

WRITE(*,*)'EDGE. EDGE generates the grid point coordinates on'

WRITE(*,*)'a curve given by a set of nodes through which'

WRITE(*,*)'a spline curve can be fitted. The spline curve can'

WRITE(*,*)'be made up from several sections.'

WRITE(*,*)

WRITE(*,*) 'Enter the number of sections'

WRITE (*, *)

READ (*, * )NOSECT

C

C

C

C

C

2O

WRITE (*, *)

WRITE (*, *)

WRITE(*,*) 'The data will be written into UNIT 20'

WRITE(20,*)NOSECT, ' Number of sections'

DO I00 is=I,NOSECT

WRITE (*, *)

WRITE (*, *)

WRITE (*, *)

WRITE (*, *)

WRITE(*,*) 'Now enter data for section number',is

WRITE (*, * )

WRITE(*,*) 'Enter number of grid points to be on the section'

WRITE (*, *)

READ(*, *) ILS

CALL GETNODES(XI,YI,ZI,NOP,MxNode)

WRITE (*, 2001) is

READ (*, *)TENSION

CALL GETSTR(is, ISTR, BETAI,BETA2)

WRITE (20, 3001) ILS, is

}_ITE (20,3002) TENSION

WRITE (20,3003) NOP

DO 20 I=I,L_OP

WRITE (20,3004)XI (I),YI(I),ZI(I),I

CONTINUE

WRITE (20, 3005) ISTR

5O
C

I00

C

C

C

2001

3001

3002

3003

3004

3005

3006

3007

C

C

IF (ISTR.NE. 4) WRITE (20, 3006) BETA1

IF (ISTR.EQ. 4) WRITE (20,3007) BETA1, BETA2

CONT INUE

FORMAT(//,' Enter the TENSION parameter for section ',I2)

FORMAT(' ',I3,' No of gridpts: Section ',I2,

$, ......... ,)

FORMAT(' ',F6.2, ' Tension parameter')

FORMAT(' ',I3, ' Number of nodes')

FORMAT(3X,3(F8.5,3X),' ---',I2)

FORMAT (' ', I3, ' StretchType' )

FORMAT(' ',F8.4,' Stretching parameter BETA')

FORMAT(' ',2(F8.4,3X),' Stretching parameters BETA1 and',

$ ' BETA2' )

STOP

END

C

C- =........ #

C

SUBROUTINE GETNODES(XI,YI, ZI,NOP,MxNode)

C

C This subroutine reads in and arranges nodal points to define

C two parallel curve sections (one on each blade surface).

C

DIMENSION Xl(MxNode),Yl(MxNode),Zl(MxNode)

C

C

C

INNUM0=I0

NOP=0

PI=ACOS(-I.)

WRITE(*,*)'The data for the section can be read in several'

WRITE(*,*)'parts, where each part is a single node point or a'

WRITE(*,*)'series of node points. Note, the node points can'

WRITE(*,*)'be read in both forward and reverse order from'

WRITE(*,*)'the input files.'

WRITE(*,*)

WRITE(*,*)'How many parts is the section composed of?'

READ(*,*)NOSECT

DO I00 ISECT=I,NOSECT

INNUM=INNUM0+ISECT

Z_RITE (*, 200) ISECT, INNUM

READ (INNUM, *)NP

WRITE (*, 201) NP

READ (*, *)NI

h_R!TE (*, 202)

READ (*, *)N2

N!2=IARSt_o _,.-:__-_]+I

NOP=NOP+Ni2

IF (_$2.GE .NI) THE_

51
I0

C

20

C

C

30

C

40

C

I00

C

2OO

201

202

C

DO I0 J=I,NI-I

READ (INNUM, * )

CONTINUE

I=NOP-NI 2

DO 20 J=NI,N2

I=I+l

READ (INNUM, *)Xl (I), Y1 (I), Zl (I)

CONTINUE

ELSE

DO 30 J=I,N2-1

READ (INNUM, * )

CONTINUE

I=NOP + 1

DO 40 J=N2,NI

I=I-!

READ (INNUM, *) Xl (I), Y1 (I), Z1 (I)

CONTINUE

ENDIF

CLOSE (INNUM)

CONTINUE

FORMAT(' Part ',I2, ' will be read in from UNIT',I3)

FORMAT(' There are ',I2,' node points in the file.',/

/ ,' Enter the number of the node point that is to be the',/

.the' first on the current part. ')

/FORMAT(' Enter the number of the node point that is to be the',/

/ , ' the last on the current part.')

RETURN

END

C

C ................

C

SUBROUTINE GETSTR(NEI,ISTR, BETAI,BETA2)

C

C This subroutine reads information regarding the stretching function

C along curve section NEI.

C

WRITE (*, 2011) NEI

WRITE (*, 2020)

I00 READ (*, *) ISTR

IF(ISTR.LT.4 .AND. ISTR.GT.0)THEN

WRITE (*, 2031) NEI

READ (*, * )BETA1

ELSEIF (ISTR.EQ. 4) THEh _

WRITE (*, 2036) NEI

READ (*, *) BETA1, BETA2

ELSEIF (ISTR.EQ. 0) THEN

BETA1=1.1

ELSE

WRITE(*,*)' Please enter a number from 0 to 4'

GOTO I00

E_,_DIF

52
C

2011

2020

2031

2036

C

FORMAT('Enter the STRETCHTYPEfor section',I2)

FORMAT(/,' Enter: 0 for no stretching',/,

/ ' 1 for concentration near lower boundary',/,

/ ' 2 for concentration near upper boundary',/,

/ ' 3 for concentration near both boundaries',/,

/ ' 4 for two-parameter stretching function')

FORMAT(' Enter the stretching parameter (BETA) for section ',

/ I2)

FORMAT(' Enter the stretching parameters (BETA1 and BETA2) ',

/ 'for section ',I2)

RETURN

END

53
A.6 Listing of GRIDTST

PROGRAM GRIDTST

PARAMETER (IM=II, JM=51,KM=I51)

C

:::::::::::::::::::::

C This program is used to test whether

C system are positive.

C

C

C

C

C

C

C

5

C

C

C

C

C

C

I0

2O

3O

all Jacobians for a grid

DIMENSION X (IM, JM, KM) ,Y (IM, JM, KM) , Z (IM, JM, KM)

DIMENSION DXDXI(IM),DXDET(IM),DXDZET(IM),DYDXI(IM),DYDET(IM)

/ , DYDZET (IM) ,DZDXI (IM) ,DZDET (IM) , DZDZET (IM)

Read in the grid point coordinates.

READ(I,*) IL

READ (1, * )JL

READ (I, *) KL

DO 5 I=l, IL

DO 5 J=l, JL

DO 5 K=I,KL

READ (i, *) X (I, J, K), Y (I, J,K), Z (I, J, K)

CONTINUE

DXI=I ./FLOAT (IL)

DET=I ./FLOAT (JL)

DZET=I ./FLOAT (KL)

Calculate the metric coefficents at the regular grid points.

INEG=0

DO 30 K=I,KL

DO 20 J=I,JL

CALL DDXIJK(DXI, IL, JL, KL,X,Y, Z,J,K

/ , DXDXI, DYDXI, DZDXI, IM, JM, KM)

CALL DDETJK(DET, IL, JL, KL,X,Y,Z,J,K

/ ,DXDET, DYDET, DZDET, IM, JM, KH)

CALL DDZETJK(DZET, IL, JL, KL,X,Y,Z,J,K

/ ,DXDZET, DYDZET, DZDZET, IM, JM, KM)

DO I0 I=I,IL

RJACB =

(DXDXI (I) * (DYDET (I) *DZDZET (I) -DZDET (I) *DYDZET (I))

+DYDXI (I) * (DZDET (I) *DXDZET (I) -DXDET (I) *DZDZET (I))

+DZDXI (I) _ (DXDET (I) *DYDZET (I) -DYDET (I) *DXDZET (I)) )

IF (RJACB. LE. 0) ZIRITE (2,50) I, J, K, RJACB

IF (RJACB. LE. 0) INEG=INEG+I

CONTINUE

CONTI_[UE

CONTINUE

IF(INEG.GT.0)WRiTE(*,*)' NEGATIVE JACOBIANS FOUND'

54
C

C

5O

C

C =

C

C

C

C

C ....

RETURN

FORMAT ('

END

', '(i, j,k)=(',2(i3, ', '),i3, ') ', J=' f el6 8)

SUBROUTINE DDXIJK (DXI, IL, JL, KL, X, Y, Z, J, K, DXDXI

/ ,DYDXI,DZDXI, IM, JM, KM)

This subroutine calculates the derivatives of X, Y, and Z I

(i.e., the coordinates of the Cartesian coordinate system) I

with respect to the coordinate XI of the transformed I

coordinates. I

C

C

C

C

C

i0

C

C

C

C

C

C

DIMENSION X(IM, JM, KM),Y(IM, JM,KM),Z(IM, JM, KM)

/ ,DXDXI (IM), DYDXI (IM) ,DZDXI (IM)

DXI2=2. *DXI

DXDXI (1 )= (-X (3, J, K) +4. *X (2, J, K) -3. *X (i, J, K) )/DXI 2

DYDXI(1)=(-Y(3,J,K)+4.*Y(2,J,K)-3.*Y(I,J,K))/DXI2

DZDXI (I) = (-Z (3, J, K) +4 .*Z (2, J, K) -3. *Z (i, J, K) )/DXI2

DO I0 I=2,1L-I

DXDXI (I) = (X (I+l, J, K) -X (I-l, J, K) )/DXI2

DYDXI (I) = (Y (I+l, J, K) -Y (I-l, J, K) )/DXI2

DZDXI(I)=(Z(I+I,J,K)-Z(I-1,J,K))/DXI2

CONTINUE

DXDXI(IL)=(3.*X(IL, J,K)-4.*X(IL-I,J,K)

/ + X(IL-2,J,K))/DXI2

DYDXI (IL) = (3. *Y (IL, J, K) -4. *Y <IL-I, J, K)

/ + Y(IL-2,J,K) )/DXI2

DZDXI(IL)=(3.*Z(IL, J, K) -4 .*Z (IL-I, J, K)

/ + Z(IL-2,J,K) )/DXI2

RETURN

END

SUBROUTINE DDETJK(DET, IL, JL, KL, X, Y, Z, J, K, DXDET

/ ,DYDET,DZDET, IM, JM,KM)

C

C

C

C

C _-

This subroutine calculates the derivatives of X, Y, and Z

(i.e., the coordinates of the Cartesian coordinate system)

with respect to the coordinate ETA of the transformed

coordinates.

D=kENS:O., X(!M, JM,}E<),Y(IM, JM, KM) Z(IM, JM, KM)

/ ,DXDET (IM), DYDET (IM) ,DZDET (IM)

55
C

C

I0

C

20

C

30

C

C

C_

I

I

I

DET2=2. *DET

IF (J.EQ.I) THEN

DO I0 I=I,IL

DXDET (I) = (-X (I, J+2,K) +4 .*X (I, J+I,K)

-3.*X (I, J, K) )/DET2

DYDET (I)-- (-Y (I, J+2, K) +4 .*Y (I, J+I,K)

-3.*Y (I, J, K) )/DET2

DZDET (I) = (-Z (I, J+2, K) +4. *Z (I, J+l, K)

-3.*Z (I, J,K) )/DET2

CONTINUE

ELSE IF (J.EQ.JL) THEN

DO 20 I=I,IL

DXDET (I) = (3. *X (I, J, K)

DYDET (I) = (3. *Y (I, J, K)

DZDET (I) = (3. *Z (I, J, K)

CONTINUE

- 4.*X(I, J-l, K)

+X(I, J-2,K) )/DET2

- 4.*Y(I, J-l, K)

+Y(I, J-2, K) )/DET2

- 4.*Z (I, J-I,K)

+Z (I, J-2, K) )/DET2

ELSE

DO 30 I=I,IL

DXDET (I) = (X (I, J+l, K)-X (I, J-l, K) )/DET2

DYDET (I) = (Y (I, J+l, K)-Y (I, J-l, K) )/DET2

DZDET (I) = (Z (I, J+l, K)-Z (I, J-l, K) )/DET2

CONTINUE

ENDIF

RETURN

END

SUBROUTINE DDZETJK(DZET, IL, JL, KL, X,Y,Z,J,K, DXDZET

/ ,DYDZET,DZDZET,IM, JM, KM)

C

C

C

C

CThis subroutine calculates the derivatives of X, Y, and Z

(i.e., the coordinates of the Cartesian coordinate system)

with respect to the coordinate ZETA of the transformed

coordinates.

DIMENSION X (IM, JM, KM), Y (IM, JM, KM), Z (IM, JM, KM)

/ ,DXDZET (IM), DYDZET (IM), DZDZET (IM)

DZET2=2 .*DZET

IF (K.EQ. I) THEN

DO I0 I=I,IL

DXDZET (I) = (-X (I, J, K+2) +4. *X (I, J, K+I)

-3. *X (I, J, K) )/DZET2

DYDZET (1) = (-Y (I, J, K+2) +4. *Y (I, J, K+I)

-3 .*Y (I, J, K) )/DZET2

56
I0

C

2O

C

30

C

DZDZET (I)= (-Z (I, J, K+2) +4 .*Z (I, J, K+I)

-3 .*Z (I, J,K) )/DZET2

CONTINUE

ELSE IF (K.EQ.KL) THEN

DO 20 I=I,IL

DXDZET (I) = (3.*X (I, J, K)

DYDZET (I) = (3.*Y (I, J,K)

DZDZET (1) = (3.*Z (I, J, K)

CONTINUE

- 4.*X (I, J,K-I)

+X (I, J, K-2) )/DZET2

- 4.*Y (I, J,K-I)

+Y (I, J, K-2) )/DZET2

- 4.*Z (I, J,K-I)

+Z (I, J, K-2) )/DZET2

ELSE

DO 30 I=I,IL

DXDZET (I) = (X (I, J, K+I) -X (I, J, K-I) )/DZET2

DYDZET (I) = (Y (I, J, K+I) -Y (I, J, K-l) )/DZET2

DZDZET (I) = (Z (I, J, K+I) -Z (I, J, K-I) )/DZET2

CONTINUE

ENDIF

RETURN

END

57
Appendix B -- LISTING OF PROGRAM GRID3D-v2

PROGRAM GRID3D

C

C

PARAMETER (MxSrfs =4, MxBCvs=16, MxBPts=21, MxGSiz=31)

C

C This SUBROUTINE generates a three-dimensional grid system using the

C "two-boundary" or "four-boundary" algebraic grid generation techniques.

C Boundary surface edge curves are formed from sets of nodal points by

C using parametric tension splines. Boundary surfaces are formed by

C using the .,bi-directional 3-D Hermite interpolation" technique.

C

INTEGER CrvNum, SrfNum, NSurfs, InNum, OutNum,

$ StrXi, StrEt, StrZt, StrAA, StrBB, II, JJ, KK, AL, BL,

$ i, j, k, NDPts(4), AAL(MxSrfs), BBL(MxSrfs),

$ NGPts(MxBCvs), Type, StrTp, EdgNum, ZoneNo

C

C

REAL EtStep, XiStep, ZtStep, AAStep, BBStep,

$ SigmaXi, SigmaEt, SigmaZt,

$ KXil, KXi2, KEtal, KEta2, KZetal, KZeta2,

$ BetaXi, BetaEt, BetaZt, BetaAA, BetaBB,

$ hl (MxGSiz), h2(MxGSiz), h3(MxGSiz), h4 (MxGSiz),

$ h5(MxGSiz), h6(MxGSiz), h7(MxGSiz), h8(MxGSiz),

$ kS(MxSrfs,MxGSiz,MxGSiz), kl(MxSrfs,MxGSiz),

$ k2(MxSrfs,MxGSiz), k3(MxSrfs,MxGSiz), k4(MxSrfs,MxGSiz),

$ SigmaAA(MxSrfs), SigmaBB(MxSrfs),

$ XB (MxGSiz, 4), YB (MxGSiz, 4), ZB (MxGSiz, 4),

$ Xl(MxGSiz), X2(MxGSiz), X3(MxGSiz), X4(MxGSiz),

$ YI(MxGSiz) , Y2(MxGSiz), Y3(MxGSiz), Y4 (MxGSiz),

$ ZI(MxGSiz), Z2(MxGSiz), Z3(MxGSiz), Z4 (MxGSiz),

$ StrB (MxGSiz,MxBCvs) ,

$ EtStll (MxGSiz), EtStl2 (MxGSiz), EtStl5 (MxGSiz),

$ EtStl6 (MxGSiz), ZtStl (MxGSiz), ZtSt2 (MxGSiz),

$ ZtSt5(MxGSiz),ZtSt6(MxGSiz)

REAL

$

$

$

$

$

$

REAL

$

$

$

$

$

$

$

REAL

$

$

$

$

$

PXSIPE (MxGSiz, MxGSiz) ,

PYSIPE (MxGSiz, MxGSiz),

PZSIPE (MxGSiz, MxGSiz) ,

PXS3Zt (MxGSiz, MxGSiz) ,

PYS3Zt (MxGSiz,MxGSiz) ,

PZS3Zt(MxGSiz,MxGSiz),

Tensn (4) ,

PXS2PE(MxGSiz,MxGSiz),

PYS2PE(MxGSiz,MxGSiz),

PZS2PE(MxGSiz,MxGSiz),

PXS4Zt(MxGSiz,MxGSiz),

PYS4Zt(MxGSiz,MxGSiz),

PZS4Zt(MxGSiz,MxGSiz)

Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts),

XDerv2(4,MxBPts), YDerv2(4,MxBPts),

ZDerv2(4,MxBPts),

x(4,MxBPts), y(4,MxBPts),

z(4,MxBPts), s(4,MxBPts),

zx(4,MxBPts), zy(4,MxBPts),

zz(4,MxBPts)

PXlPBB (MxGSiz

PYIPBB (MxGSiz

PZIPBB (MxGSiz

PXIPAA (M>:GSiz

PYIPAA (MxGSiz

PZ!PAA (MxGSiz

PX2PBB(MxGSiz),

PY2PBB(MxGSiz),

PZ2PBB(MxGSiz),

PX2PAA(MxGSiz),

PY2PAA(MxGSiz),

PZ2PAA(MxGSiz),

S8
PX3PBB(MxGSiz), PX4PBB(MxGSiz),

PY3PBB(MxGSiz), PY4PBB(MxGSiz),

PZ3PBB(MxGSiz), PZ4PBB(MxGSiz),

PX3PAA(MxGSiz), PX4PAA(MxGSiz),

PY3PAA (MxGSiz), PY4PAA(MxGSiz),

PZ3PAA(MxGSiz), PZ4PAA(MxGSiz),

XS (MxSrfs, MxGsiz, MxGsiz),

YS (MxSrfs, MxGsiz,MxGsiz),

ZS (MxSr fs, MxGsiz, MxGsi z ),

XPnt (MxGSiz, MxGSiz,MxGSiz),

YPnt (MxGSiz, MxGSiz, MxGSiz ),

ZPnt (MxGSiz, MxGSiz,MxGSiz)

EQUIVALENCE (XB(I,I),XI(1)), (YB(I,I),YI(1)), (ZB(I,I),ZI(1)),

$ (XB(I,2),X2(1)), (YB(I,2),Y2 (I)), (ZB(I,2),Z2(1)),

$ (XB(I,3),X3(1)), (YB(I,3),Y3(1)), (ZB(I,3),Z3(1)),

$ (XB(I,4),X4 (I)), (YB(I,4),Y4 (i)), (ZB(I,4),Z4 (i))

EQUIVALENCE (StrB(l,l),ZtStl(1)), (StrB(l,2),ZtSt2(1)),

$ (StrB(l,5),ZtSt5(1)), (StrB(l,6),ZtSt6(1)),

$ (StrB(l,ll),EtStll(1)), (StrB(l,12),EtSt12(1)),

$ (StrB(l,15),EtStl5(1)), (StrB(l,16),EtStl6(1))

C

C Specify input and output device unit numbers for Region i. This is

C convenient for running the program on a PC. For a mainframe, you will

C need to use the FORTRAN OPEN and CLOSE statements or alter the input

C to use a namelist.

C

I nNum= 7

OutNum=8

C

C

C

C

C Read in the grid control information.

C

CALL RdGrIn(II,JJ, KK, NSurfs,SigmaXi,SigmaEt,SigmaZt,

$ kXil, kXi2, kEtal, kEta2, kZetal, kZeta2, InNum)

C

C

C Set

C

C

various parameters for the grid generation routines.

CALL KFctrs(l,kS,kl,k2,k3,k4,kXil,kXi2,kEtal,kEta2,

$ kZetal, kZeta2, MxSrfs, MxGSiz)

AAL (I

AAL (2

AAL (3

AAL (4

BBL (1

BBL (2

BBL (3

BBL (4

,=KK

=KK

=II

=II

=If

=II

=JJ

=gJ

S igma._.h (1 )=Si_maZt

Si97rLaAA (2) =Si_maZt

S igmaL.A (3 )=Sig_.,aXi

SigmaAA (4 ) =SigmaXi

Sigw_aBB (1 )=Si_maXi

59
C

C

C

C

C

C

C

C

C

C

C

30

C

C

C

C

SigmaBB (2 )=SigmaXi

SigmaBB (3) =SigmaEt

SigmaBB (4 )=SigmaEt

NGPts (I)=KK

NGPts

NGPts

NGPts

NGPts

NGPts

NGPts

NGPts

NGPts

NGPt s

NGPt s

NGPts

NGPt s

NGPts

NGPt s

NGPts

(2 )=KK

(3) =II

(4)=II

(5) =KK

(6 )=KK

(7)--II

(8)=II

(9) =II

(i0 =II

(ii =JJ

(12 =JJ

(13 =II

(14 =II

(15 =JJ

(16 =JJ

Calculate the grid point spacings in the transformed domain.

XiStep=l .0/(II-I)

EtStep=l.0/(JJ-l)

ZtStep=l. 0/ (KK-I)

Calculate the boundary surface grid point locations for each surface.

DO 40 SrfNum=l,NSurfs

Read in the edge curve nodal points and form the boundary surface

edge curves for surface SrfNum by sp!ining.

DO 30 CrvNum=l,4

EdgNum=(SrfNum-l)*4

READ(InNum,*)Type

IF(Type.EQ.I)THEN

CALL

CALL

ELSE

CALL

CALL

CALL

CALL

$

$

ENDIF

CONTINUE

+ CrvNum

RdGrPIn (NGPts (EdgNum) ,XB, YB, ZB, CrvNum, MxGSiz, InNum)

CalStl (NGPt s (EdgNum) ,XB, YB, ZB, CrvNum, EdgNum,

StrB, MxBCvs, MxGSiz)

RdCv!n(x,y,z,NDPts,CrvNum, Tensn, InNum, MxBPts,

StrTp, Betal,Beta2)

PTSpln(x,y,z,s, zx, zy, zz,Diag, OfDiag, Right,NDPts,

Tensn(CrvNum),CrvNum, MxBPts)

CalSt2(EdgNum, NGPts(EdgNum),StrTp, Betal,Beta2,

StrB,MxBCvs,MxGSiz)

EdgGPts(CrvNum, EdgNum, NGPts(EdgNum),XB,YB, ZB, StrB,

x,y,z,s, zx, zy, zz,NDPts,Tensn(CrvNum),

MxBCvs,MxBPts,MxGSiz)

Calculate the boundary

CALL

surface edge derivative values for surface SrfNum.

EdgDer(PX!P_.,PX2P_.,PYIP_A,PY2PAA, PZlPAA, PZ2PAA,

6O
C

C

C

C

C

4O

C

5O

60

C

C

C

C

C

C

PX3PBB, PX4PBB, PY3PBB, PY4PBB, PZ3PBB, PZ4PBB,

Xl, X2, X3, X4, YI, Y2, Y3, Y4, ZI, Z2, Z3, Z4,

AAL (Sr fNum), BBL (SrfNum), MxGSiz)

Calculate the boundary surface grid point locations for surface

CALL TwoBnd(XS, YS, ZS, SrfNum, AAL(SrfNum),BBL(SrfNum),

SigmaBB(SrfNum),kl,k2,StrB,

hl,h2,h3,h4,Xl,X2,X3,X4,

YI,Y2,Y3, Y4,ZI,Z2, Z3,Z4,PXlPBB, PX2PBB,

PYIPBB, PY2PBB, PZlPBB, PZ2PBB, PXlPAA, PX2PAA,

PYIPAA, PY2PAA, PZIPAA, PZ2PAA, PX3PBB,PX4PBB,

PY3PBB, PY4PBB,PZ3PBB, PZ4PBB,

MxBCvs,MxGSiz,MxSrfs)

$

$

$

$

$

$

$

$

$

CALL ForBnd (XS, YS, ZS, SrfNum, AAL (SrfNum), BBL (SrfNum),

SigmaAA (SrfNum) ,SigmaBB (SrfNum) ,

k3, k4, StrB,

hl, h2, h3, h4, h5, h6, h7, h8,

Xl, X2, X3, X4, YI, Y2, Y3, Y4, Zl, Z2, Z3, Z4,

PXIPBB, PX2PBB, PYIPBB, PY2PBB, PZIPBB, PZ2PBB,

PXlPAA, PX2PAA, PYIPAA, PY2PAA, PZIPAA,PZ2PAA,

PX3PBB, PX4PBB, PY3PBB,PY4PBB, PZ3PBB,PZ4PBB,

PX3PAA, PX4PAA, PY3PAA, PY4PAA, PZ3PAA, PZ4PAA,

MxBCvs,MxGSiz,MxSrfs)

CONTINUE

IF(NSurfs.EQ.2)THEN

DO 60 SrfNum=3,4

DO 50 CrvNum=3,4

EdgNum=(SrfNum-i)*4 + CrvNum

READ(InNum,*)StrTp

IF(StrTp. NE.4)THEN

READ (InNum, *) Betal

ELSE

READ (InNum, *) Betal, Beta2

ENDIF

CALL CalSt2(EdgNum, NGPts(EdgNum),StrTp, Betal,Beta2,

StrB,MxBCvs,MxGSiz)

CO_TINUE

CONTINUE

ENDIF

Calculate the interior grid point locations.

CALL

$

$

$

$

TwoSrf(XPnt,YPnt, ZPnt,II,JJ, KK, SigmaEt,kS,

EtStll,EtStl2,EtStl5,EtStl6,

XiStep,EtStep, ZtStep, XS, YS,ZS,hl,h2,h3,h4,

PXS!PE, PXS2PE,PYSIPE,PYS2PE,PZSIPE,

PZS2PE,MxGSiz,MxSrfs)

IF (NSurfs.EQ.4) THEN

CALL ForSrf(XPnt,YPnt,ZPnt,II,JJ, KK,

$ SigmaEt,Sig_r, aZt,kS,

$ EtStll,EtStl2,EtStl5,EtStl6,

SrfNum.

6]
C

C

$

$

$

$

$

$

ENDIF

ZtStl, ZtSt2, ZtSt5, ZtSt6,

XS,YS, ZS,XiStep, EtStep, ZtStep,

hl,h2,h3,h4,h5,h6,h7,h8,

PXSIPE,PXS2PE, PYSIPE,PYS2PE,PZSIPE, PZS2PE,

PXS3Zt,PXS4Zt,PYS3Zt,PYS4Zt,PZS3Zt,PZS4Zt,

MxGSiz,MxSrfs)

CALL PrGrid (XPnt, YPnt, ZPnt, II, JJ, KK, OutNum, MxGSiz )

C

C

C

C

C

C

C

C

C

SUBROUTINE

$

$

$

$

TwoSrf(XPnt,YPnt, ZPnt, II,JJ, KK, SigmaEt,kS,

EtStll,EtStl2,EtStl5,EtStl6,

XiStep,EtStep, ZtStep, XS, YS, ZS,

hl,h2,h3,h4,PXSIPE, PXS2PE,PYSIPE,PYS2PE, PZSIPE,

PZS2PE,MxGSiz,MxSrfs)

This SUBROUTINE calculates the grid point locations between two specified

surfaces using the "two-boundary technique".

INTEGER i, j, k, StrXi, StrEt, StrZt, II, JJ, KK

REAL Xi, Eta, Zeta, XiNew, EtaNew, ZtaNew, LLI, LL2,

$ PXSIXi, PXS2Xi, PYSIXi, PYS2Xi, PZSIXi, PZS2Xi,

$ PXSIZt, PXS2Zt, PYSIZt, PYS2Zt, PZSIZt, PZS2Zt,

$ Beta×i, BetaEt, BetaZt, XiStep, EtStep, ZtStep,

$ EtStll(MxGSiz),EtStl2(MxGSiz),

$ EtStl5(MxGSiz),EtStl6(MxGSiz),

$ kS(MxSrfs,MxGSiz,MxGSiz),

$ hl(MxGSiz), h2(MxGSiz), h3(MxGSiz), h4(MxGSiz),

$ PXSIPE(MxGSiz,MxGsiz), PXS2PE(MxGSiz,MxGsiz),

$ PYSIPE(MxGSiz,MxGsiz), PYS2PE(MxGSiz,MxGsiz),

$ PZSIPE(MxGSiz,MxGsiz), PZS2PE(MxGSiz,MxGsiz),

$ XS(MxSrfs,MxGsiz,MxGsiz),

$ YS(MxSrfs,MxGsiz,MxGsiz),

$ ZS(MxSrfs,MxGsiz,MxGsiz),

$ XPnt(MxGSiz,MxGsiz,MxGSiz),

$ YPnt(MxGSiz,MxGsiz,MxGSiz),

$ ZPnt(MxGSiz,MxGsiz,MxGSiz)

Calculate the derivative values along the constant Xi/Zeta

boundaries.

PXSIXi=

PXS2Xi =

PYSIXi=

PYS2Xi= YS

PZSIXi=ZS

PZS2Xi = ZS

PXSIZt=(XS

PXS2Zt=(XS

PYSIZt=(YS

PYS2Zt=(YS

XS I 1,2)-XS

XS 2 1,2)-XS

YS I I, 2 -YS

2 i, 2 -YS

I I, 2 -ZS

2 1,2 -ZS

1,2, ! -XS

2,2, 1 -XS

i, 2, 1 -YS

2,2, 1 -YS

(i,I,

(2,1,

(I,I,

(2,1,

(I,I

(2,1

(I,I

(2,1

(I,i

(2,1

l))/XiStep

l))/XiStep

l))/XiStep

l))/XiStep

l))/XiStep

l))/XiStep

l))/ZtStep

l))/ZtStep

l))/ZtStep

l})/ZtStep

62
C

C

C

C

PZSlZt = (ZS (i, 2, l)-ZS (i, i, I) )/ZtStep

PZS2Zt= (ZS (2,2, I) -ZS (2, I, I) )/ZtStep

LLI = ((PYSIXi*PZSIZt-PZSIXi*PYSIZt) **2

/ +(PXSIXi*PZSIZt-PZSIXi*PXSIZt)**2

/ +(PXSIXi*PYSIZt-PYSIXi*PXSIZt)**2)**0.5

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

/ +(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

/ +(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIPE (I, I)=-kS (I, I, I) *

PXS2PE (I, l)=-kS (2, I, I) *

PYSIPE (i, I) = kS (i, i,I)*

PYS2PE (i, I) = kS (2, I,I)*

PZSIPE(I,I)=-kS(I,I,I)*

PZS2PE (I, I) =-kS (2, i, I) *

(PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

(PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LLI

(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

(PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

PXSIXi = (XS (I, i, II) -XS (I,

PXS2Xi= (XS (2, I, II) -XS (2,

PYSIXi = (YS (I, i, II)-YS (i,

PYS2Xi =(YS(2,I,II)-YS(2,

PZSIXi =(ZS(I,I,II)-ZS(I,

PZS2Xi =(ZS(2,I,II)-ZS(2,

PXSiZt = (XS (i, 2, II)-XS (i,

PXS2Zt =(XS (2, 2,II)-XS(2,

PYSIZt =(YS(I,2,II)-YS(I,

PYS2Zt = (YS (2, 2, II)-YS (2,

PZSIZt = (ZS (i, 2, II)-ZS (i,

PZS2Zt= (ZS (2, 2, II)-ZS (2,

i, II-l) )/XiStep

i, II-l) )/XiStep

I, II-l) )/XiStep

I, II-l) )/XiStep

I, II-l) )/XiStep

i, II-l) )/XiStep

l, II

l, II

l,II

l,II

l,II

i, II

)/ZtStep

)/ZtStep

)/ZtStep

)/ZtStep

)/ZtStep

)/ZtStep

LLI=((PYSIXi*PZSIZt-PZSIXi*PYSIZt)**2

/ +(PXSIXi*PZS!Zt-PZSIXi*PXSIZt)**2

/ +(PXSIXi*PYSIZt-PYSIXi*PXSIZt)**2)**0.5

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

/ +(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

/ +(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIPE(II,I)=-kS(I,II, I)* (PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

PXS2PE (II, I)=-kS (2, II, i) * (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

PYSIPE(II,I)= kS(I, II, I)*(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LL!

PYS2PE(II,I)= kS (2,II,I)*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

PZSIPE (II, I) =-kS (I, II, i) * (PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

PZS2PE(II,I)=-kS(2,II, I)* (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

DO 55 i=2,II-i

PXSIXi = XS 1

PXS2Xi= XS 2

PYSIXi = YS 1

PYS2Xi= YS 2

PZSIXi = ZS 1

PZS2Xi = ZS 2

PXSIZt = XS (I

PXS2Zt = XS(2

PYSIZt = YS(I

PYS2Zt= YS(2

PZSIZt = ZS (i

PZS2Zt= ZS(2

I, i+l -XS

I, i+l -XS

i, i+l -YS

I, i+l -YS

i, i+l -ZS

I, i+l -ZS

2, i) -XS (I

2, i) -XS (2

2,i)-YS(I

2, i) -YS (2

2, i) -ZS (I

2,i)-ZS(2

(i, I, i-i

(2, l,i-I

(I, l,i-I

(2, I, i-I

(i, I, i-i

(2, I, i-I

l,i)

l,i)

l,i)

l,i)

l,i)

l,i)

/2/XiStep

/2/XiStep

/2/XiStep

/2/XiStep

/2/XiStep

/2/XiStep

/ZtStep

/ZtStep

/ZtStep

/ZtStep

/ZtStep

/ZtStep

LLi=((PYSIXi*PZSIZt-PZSIXi*PYSlZt)**2

t(PXSIXi*PZSiZt-PZSIXi*PXSIZt)**2

+(PXS!Xi*PYS!Zt-PYSIXi*PXSIZt)**2)**0.5

63
55

C

C

C

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

+(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

+(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIPE(i, I) =-kS (i, i, i)* (PYSIXi*PZSIZt-PZSlXi*PYSIZt)/LLI

PXS2PE (i, i) =-kS (2, i, i) * (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

PYSIPE (i, I)= kS (I, i, I) * (PXSIXi*PZSIZt-PZSIXi*PXSIZt) /LLI

PYS2PE(i,l)= kS(2,i,I)*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

PZSIPE(i,I)=-kS(I,i,I)*(PXSIXi*PYSlZt-PYSIXi*PXSIZt)/LLI

PZS2PE(i,I)=-kS(2,i,I)*(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

CONTINUE

DO 70 k=2, KK-I

PXSIXi = (XS (I,

PXS2Xi = (XS (2,

PYSIXi= (YS (i,

PYS2Xi = (YS (2,

PZSlXi =(zS(I,

PZS2Xi= (ZS (2,

PXSIZt= (XS (i,

PXS2Zt = (XS (2,

PYSIZt = (YS (I,

PYS2Zt = (YS (2,

PZSIZt = (ZS (i,

k, 2) -XS (I, k, i) )/XiStep

k, 2) -XS (2, k, I) )/XiStep

k, 2) -YS (I, k, I) )/XiStep

k, 2)-YS (2, k, I) )/XiStep

k, 2) -ZS (i, k, i) )/XiStep

k, 2) -ZS (2, k, I) )/XiStep

k+l, i) -XS (I, k-l, i) )/2/ZtStep

k+l, i) -XS (2, k-l, I) )/2/ZtStep

k+l, I)-YS (I, k-l, I) )/2/ZtStep

k+l, i) -YS (2, k-l, i) )/2/ZtStep

k+l, I)-ZS (I, k-l, i) )/2/ZtStep

PZS2Zt=(ZS(2,k+I,I)-ZS(2,k-I,I))/2/ZtStep

LLI=((PYSIXi*PZSIZt-PZSIXi*PYSIZt)**2

+(PXSIXi*PZSIZt-PZSIXi*PXSIZt)**2

+(PXSIXi*PYSIZt-PYSIXi*PXSIZt)**2)**0.5

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

+(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

+(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIPE (i, k)

PXS2PE (l,k)

PYSIPE (I, k)

PYS2PE (I, k)

PZSlPE (l,k)

PZS2PE (i, k)

=-kS(I,I,k)*(PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

=-kS(2,I,k)*(PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

= kS(I,I,k)*(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LLI

= kS(2,I,k)*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

=-kS(I,I,k)*(PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

=-kS(2,!,k)*(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

PXSlXi = (XS (I,

PXS2Xi = (XS (2,

PYSIXi =(Ys(I,

PYS2Xi= (YS (2,

PZSIXi = (ZS (I,

PZS2Xi = (ZS (2,

PXSIZt = (XS (I,

PXS2Zt = (XS (2,

PYSIZt = (YS (I,

PYS2Zt= (YS (2,

PZSIZt = (ZS (I,

k, II

k,II

k, II

k, II

k, II

k, ll

k+l, II

k+l, II

k+l, II

k+l, II

k+l, II

-XS (i, k, II-l) )/XiStep

-XS (2, k, II-l) )/XiStep

-YS (I, k, II-l) )/XiStep

-YS (2, k, II-l) )/XiStep

-ZS (i, k, II-l) )/XiStep

-ZS (2, k, II-l) )/XiStep

-XS(l,k-l,II) /2/ZtStep

-XS(2,k-l,II) /2/ZtStep

-YS(l,k-l,II) /2/ZtStep

-YS(2,k-l,II) /2/ZtStep

-ZS (I, k-l, II) /2/ZtStep

PZS2Zt =(ZS(2,k+I,II -ZS(2,k-l,II)

LLI = ( (PYSIXi*PZSIZt-PZSIXi*PYSIZt

+ (PXS IXi*PZSIZt-PZS IXi*PXSIZt

+ (FXS IXi*PYSIZr-PYSIXi*PXS 1Zt

LL2 = ((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt

+ (PXS2Xi*PZS2 Zt-PZS2Xi*PXS2Zt

+ (PXS2Xi* PYS2 Zt-PYS2Xi*PXS2Zt

/2/ZtStep

**2

**2

**2)**0.5

**2

**2

**2)**0.5

64
C

60

70

C

PXSIPE (II, k) =-kS (I, II, k) * (PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

PXS2PE (II, k)=-kS (2, II, k) * (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

PYSIPE(II,k) = kS(I,II,k)*(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LLI

PYS2PE(II,k)-- kS (2,II,k)*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

PZSIPE (II, k) =-kS (i, II, k) * (PXSIXi*PYSIZt-PYSIXi*PXSlZt)/LLI

PZS2PE (II, k) =-kS (2, II, k) * (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

DO 60 i--2, II-i

PXSlXi= (XS (i, k, i+l) -XS

PXS2Xi = (XS

PYSIXi = (YS

PYS2Xi = (YS

PZSIXi = (ZS

PZS2Xi = (ZS

PXSI Zt-- (XS

PXS2Zt= (XS

PYSIZt= (YS

PYS2Zt= (YS

PZSIZt = (ZS

(2, k, i+l) -XS

(I, k, i+l) -YS

(2, k, i+l) -YS

(i, k, i+l) -ZS

(2, k, i+l) -ZS

(i, k+l, i) -XS

(2, k+l, i) -XS

(I, k+l, i)-YS

(2, k+l, i) -YS

(I, k+l, i) -ZS

(I, k, i-l) )/2/XiStep

(2, k, i-l) )/2/XiStep

(I, k, i-l) )/2/XiStep

(2, k, i-l) )/2/XiStep

(i, k, i-l) )/2/XiStep

(2, k, i-l) )/2/XiStep

(i, k-l, i) )/2/ZtStep

(2, k-l, i) )/2/ZtStep

(I, k-l, i) )/2/ZtStep

(2, k-l, i))/2/ZtStep

(l,k-l,i))/2/ZtStep

PZS2Zt = (ZS (2, k+l, i) -ZS (2, k-l, i) )/2/ZtStep

LLI = ((PYSIXi*PZSIZt-PZSIXi*PYSIZt) **2

+ (PXSIXi*PZSIZt-PZSIXi*PXSIZt) **2

+ (PXSIXi*PYSIZt-PYSIXi*PXSIZt) **2) **0.5

LL2 = ((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt) **2

+ (PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt) **2

+ (PXS2Xi*PYS2 Zt-PYS2Xi*PXS2Zt) **2) **0 .5

PXSIPE (i, k) =-kS

PXS2PE (i, k) =-kS

PYSIPE (i, k) = kS

PYS2PE (i,k) = kS

PZSIPE (i, k) =-kS

PZS2PE (i, k) =-kS

CONT I NUE

CONTINUE

i, i, k) * (PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

2,i,k)*(PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

I,i,k)*(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LLI

2, i, k) * (PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

I, i, k) * (PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

2, i, k) * (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt) /LL2

PXSIXi= (XS (I, KK, 2) -XS

PXS2Xi = (XS (2, KK, 2) -XS

PYSlXi = (YS (I, KK, 2)-YS

PYS2Xi =(Ys(2 KK, 2)-YS

PZSIXi =(ZS(I KK,2)-Z$

PZS2Xi:(ZS(2 KK, 2)-ZS

PXSlZt =(XS(I KK, I)-XS

PXS2Zt =(XS(2 KK, I)-XS

PYSIZt= (YS (i KK, I)-YS

PYS2Zt =(YS(2 KK, I)-YS

PZSlZt=(ZS(I KK, I)-ZS

PZS2Zt =(ZS (2 KK, I)-ZS

I,KK, i) )/XiStep

2,KK, I) )/XiStep

l,KK, i))/XiStep

2,KK, I) )/XiStep

I, KK, I) )/XiStep

2, KK, I) )/XiStep

I, KK-I, I) )/ZtStep

(2,KK-l,l))/ZtStep

(I, KK-I, I) )/ZtStep

(2, KK-I, i) )/ZtStep

(l,KK-l,l))/ZtStep

(2,KK-I,I))/ZtStep

LLI=((PYSIXi PZSIZt-PZSIXi*PYSIZt)**2

/ +(PXSIXi*PZSIZt-PZSlXi*PXSlZt)**2

/ +(PXSIXi*PYSIZt-PYSIXi*PXSIZt)**2)**0.5

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

/ +(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

/ +(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIFE(I,KK)=-kS(I, I,KK)*(PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

PXS2FE(!,KK)=-kS(2, I,KK)* (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

PYS!FE(I,KK): kS (i, I,KK) _ (PXSIXi*PZSIZt-PZSIXi*PXSIZt) /LLI

68
C

PYS2PE(i, KK)= kS(2, I, KK)

PZSIPE(I, KK)=-kS (I, I, KK)

PZS2PE(i, KK)---kS(2, I, KK)

*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

*(PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

*(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

PXSIXi= (XS(i, KK,II) -XS(i, KK,II-l) )/XiStep

PXS2Xi= (XS(2, KK,II) -XS(2, KK,II-l) )/XiStep

PYSIXi--(YS(i, KK,II) -YS(i, KK,II-l) )/XiStep

PYS2Xi--(YS(2, KK,II) -YS(2, KK,II-l) )/XiStep

PZSIXi= (ZS(i, KK,II) -ZS(i, KK,II-l) )/XiStep

PZS2Xi= (ZS(2, KK,II) -ZS(2, KK,II-l) )/XiStep

PXSIZt--(XS(I, KK,II) -XS(I, KK-I, II) )/ZtStep

PXS2Zt--(XS(2, KK,II) -XS(2, KK-I, II) )/ZtStep

PYSIZt= (YS(i, KK,II) -YS(i, KK,I, II) )/ZtStep

PYS2Zt=(YS(2, KK,II) -YS(2, KK-I, II) )/ZtStep

PZSIZt= (ZS(i, KK,II) -ZS(i, KK-I, II) )/ZtStep

PZS2Zt=(ZS(2,KK, II)-ZS(2,KK-I,II))/ZtStep

LLI=((PYSIXi*PZSIZt-PZSIXi*PYSIZt)**2

/ +(PXSIXi*PZSIZt-PZSIXi*PXSIZt)**2

/ +(PXSIXi*PYSIZt-PYSIXi*PXSIZt)**2)**0.5

LL2=((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)**2

/ +(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)**2

/ +(PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)**2)**0.5

PXSIPE (II, KK)---kS (i, II, KK) * (PYSIXi*PZSIZt-PZSIXi*PYSIZt)/LLI

PXS2PE (II, KK) =-kS (2, II, KK) * (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt)/LL2

PYSIPE(II,KK)= kS(I,II,KK)*(PXSIXi*PZSIZt-PZSIXi*PXSIZt)/LLI

PYS2PE(II,KK)= kS (2,II,KK)*(PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt)/LL2

PZSiPE (II, KK) =-kS (I, II, KK) * (PXSIXi*PYSIZt-PYSIXi*PXSIZt)/LLI

PZS2PE (II, KK) =-kS (2, II, KK) * (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt)/LL2

/

!

DO 75 i=2,II-I

PXSIXi = (XS

PXS2Xi = (XS

PYSIXi= (YS

PYS2Xi= (YS

PZSIXi= (ZS

PZS2Xi = (ZS

PXSIZt = (XS

PXS2Zt = (XS

PYSIZt = (YS

PYS2Zt = (YS

PZSIZt = (ZS

i, KK, i+l) -XS (I, KK, i-l) )/2/XiStep

2, KK, i+l) -XS (2, KK, i-l) ) /2/XiStep

i, KK, i+l) -YS (i, KK, i-l) )/2/XiStep

2, KK, i+l) -YS (2, KK, i-l) ) /2/XiStep

i, KK, i+l) -ZS (I, KK, i-l) )/2/XiStep

2, KK, i+l) -ZS (2, KK, i-l) )/2/XiStep

i, KK, i) -XS (I, KK-I, i) )/ZtStep

2, KK, i) -XS (2, KK-I, i) )/ZtStep

I, KK, i) -YS (I, KK-I, i) )/ZtStep

2, KK, i) -YS (2, KK-I, i) )/ZtStep

I, KK, i) -ZS (i, KK-I, i) )/ZtStep

PZS2Zt=(ZS 2,KK, i)-ZS(2,KK-I,i))/ZtStep

LLI = ((PYSIXi*PZSIZt-PZSIXi*PYSIZt) **2

+ (PXS IXi*PZSIZt-PZSIXi*PXSIZt) **2

+ (PXS IXi*PYSIZt-PYSIXi*PXSIZt } **2) **0.5

LL2= ((PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt) **2

+ (PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt) **2

+ (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt) **2) **0.5

75

PXSIPE (i

PXS2PE (i

PYSIPE (i

PYS2PE (i

PZSIFE(i

PZS2PE(i

CONT INUE

KK) =-kS

KK) =-kS

KK)= kS

W_K)= kS

KK) =-kS

KK) =-kS

i, i, KK) * (PYSIXi*PZSIZt-PZSIXi*PYSIZt

2, i, KK) * (PYS2Xi*PZS2Zt-PZS2Xi*PYS2Zt

I, i, KK) * (PXSIXi*PZSIZt-PZSIXi*PXSlZt

2, i, KK) * (PXS2Xi*PZS2Zt-PZS2Xi*PXS2Zt

I, i, KK) * (PXSlXi*PYSIZt-PYSIXi*PXSIZt

2, i, KK) * (PXS2Xi*PYS2Zt-PYS2Xi*PXS2Zt

/LLI

/LL2

/LLI

/LL2

/ELl

/LL2

66
C

C

8O

90

I00

C

C

C _-

Calculate the interior grid point locations.

DO I00 k=I,KK

Zeta = (k-l.) / (KK-I.)

DO 90 i=l,II

Xi=(i-l.)/(II-l.)

DO 80 j=l,JJ

EtaNew= (EtStll (j) * (l.-Xi) +EtStl2 (j) *Xi) * (l.-Zeta)

$ + (EtStl5 (j) * (l.-Xi) +EtStl6 (j) *Xi) *Zeta

CALL FindHs(hl(j),h2(j),h3(j),h4(j),EtaNew, SigmaEt)

XPnt (i, j, k)=hl (j)

$ *XS (i, k, i) +h2 (j) *XS (2,k, i)

$ +h3 (j) *PXSIPE (i, k)

$ +h4 (j) *PXS2PE (i, k)

YPnt(i, j,k)=hl(j)

$ *YS (I, k, i) +h2 (j) *YS (2, k, i)

$ +h3 (j) *PYSIPE (i, k)

$ +h4(j)*PYS2PE(i,k)

ZPnt (i, j,k)=h! (j)

$ *ZS (i, k, i) +h2 (j) *ZS (2, k, i)

$ +h3 (j) *PZSIPE (i, k)

$ +h4 (j) *PZS2PE (i, k)

CONTINUE

CONT INU E

CONTINUE

RETURN

END

C

C

C

C

C

C

SUBROUTINE

$

$

$

$

$

$

$

ForSrf (XPnt, YPnt, ZPnt, II, JJ, KK, SigmaEt, SigmaZt, kS,

EtSt !I, EtStl2, EtStl5, EtStl6,

ZtStl, ZtSt2, ZtSt5, ZtSt6,

XS, YS, ZS, XiStep, EtStep, ZtStep,

hl, h2, h3, h4, h5, h6, h7, h8,

PXSlPE, PXS2PE, PYSIPE, PYS2PE, PZSlPE, PZS2PE,

PXS3Zt, PXS4Zt, PYS3Zt, PYS4 Zt, PZS3Zt, PZS4Zt,

MxGSiz,MxSrfs)

This SUBROUTINE adjusts the grid so that the other two surfaces of

region are mapped correctly using the "four-boundary technique".

INTEGER i, j, k, StrXi, StrEt, StrZt, II, JJ, KK

REAL Xi, Eta, Zeta, XiNew, EtaNew, ZtaNew, LL3, LL4,

$ hl(MxGSiz), h2(MxGSiz), h3(MxGSiz), h4(MxGSiz),

$ h5(MxGSiz), h6(MxGSiz), h7(MxGSiz), hS(MxGSiz),

$ PXS3Xi, PXS4Xi, PYS3Xi, PYS4Xi, PZS3Xi, PZS4Xi,

$ PXS3PE, PXS4PE, PYS3PE, PYS4PE, PZS3PE, PZS4PE,

$ P2X00, P2X01, P2XI0, P2Xll, P2Y00, P2Y01, P2YI0, P2YII,

$ P2Z00, P2Z01, P2ZI0, P2ZII

REAL BetaXi, BetaEt, BetaZt, XiStep, EtStep, ZtStep,

$ EtStll(MxGSiz),EtStl2(MxGSiz),

$ EtSt!5(MxGSiz),EtStl6(MxGSiz),

$ Z_S_i(MxGSiz),ZtSt2(MxGSiz),

$ ZtS:5(MxGSiz),ZtSt6(MxGSiz),

the

67
C

C

C

C

C

C

C

kS(MxSrfs,MxGSiz,MxGSiz),

PXSIPE (MxGSiz .MxGsiz),

PXS3Zt (MxGSiz .MxGsiz),

PYSIPE (MxGSiz .MxGsiz),

PYS3Zt (MxGSiz .MxGsiz),

PZSIPE (MxGSiz .MxGsiz),

PZS3Zt (MxGSiz .MxGsiz),

PXS2PE(MxGSiz,MxGsiz),

PXS4Zt(MxGSiz,MxGsiz),

PYS2PE(MxGSiz,MxGsiz),

PYS4Zt(MxGSiz,MxGsiz),

PZS2PE(MxGSiz,MxGsiz),

PZS4Zt(MxGSiz,MxGsiz),

XS (MxSr f s, MxGSiz, MxGsiz ),

YS (MxSrfs, MxGSiz, MxGsiz ),

ZS (MxSrfs, MxGSiz, MxGsiz),

XPnt (MxGSiz, MxGsiz, MxGSiz ) ,

YPnt (MxGSiz, MxGsiz, MxGSiz ) ,

ZPnt (MxGSiz, MxGsiz, MxGSiz )

Calculate the derivative values along the constant Xi/Eta

boundaries.

PXS3Xi= (XS

PXS4Xi = (XS

PYS3Xi= (YS

PYS4Xi= (YS (4

PZS3Xi= (ZS (3

PZS4Xi = (ZS (4

PXS3PE= (XS (3

PXS4PE = (XS (4

PYS3PE = (YS (3

PYS4PE = (YS (4

PZS3PE= (ZS (3

PZS4PE= (ZS (4

LL3 = ((PYS3Xi

(3,2, !) -XS (3, i, i) )/XiStep

(4,2,1)-XS(4,1,1))/XiStep

(3, 2, I) -YS (3, i, i) )/XiStep

2, I) -YS (4, i, i) )/XiStep

2, I) -ZS (3, i, i) )/XiStep

2,1)-ZS(4,1,1))/XiStep

i, 2) -XS (3, i, i) )/EtStep

1,2)-XS(4,1,1))/EtStep

1, 2) -YS (3,1,1) )/EtStep

1,2)-YS(4,1, l))/EtStep

I, 2) -ZS (3, i, I) )/EtStep

1, 2) -ZS (4,1,1) )/EtStep

*P ZS 3PE-PZS3Xi*PYS3PE) **2

/ +(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

/ +(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

/ +(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

/ +(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (I, I)=-kS (3

PXS4Zt(l,l)=-kS(4

PYS3Zt (I,I)= kS (3

PYS4Zt (I,I)= kS (4

PZS3Zt (l,l)=-kS(3

PZS4Zt (l,l)=-kS(4

I)*(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

I)*(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

I)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

I)*(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

I)*(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

I)*(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

PXS3Xi =(xs(3,II,l

PXS4Xi =(XS (4, II,l

PYS3Xi = (YS (3, II, 1

PYS4Xi =(YS(4,II,I

PZS3Xi= (ZS (3

PZS4Xi =(zs(4

PXS3PE= (XS (3

PXS4PE =(xS (4

PYS3PE = (YS (3

PYS4FE = (YS (4

PZS3PE =(ZS(3

PZS4PE =(ZS(4

LL3 = ((PYS3Xi

/ + (PXS3Xi

II,l

II,l

11,2

11,2

11,2

II,2

II,2

II,2

-XS (3, II-l, I)

-XS (4, II-l, i)

-YS (3, II-l, i)

-YS (4, II-l, I)

-ZS(3, II-l,l)

-ZS(4,II-I i)

-XS(3,II, 1

-XS (4, II, 1

-YS(3,II,I

-YS (4, II, 1

-ZS(3,II,I

-ZS (4, II, 1

/XiStep

/XiStep

/XiStep

/XiStep

/XiStep

/XiStep

/EtStep

/EtStep

/EtStep

/EtStep

/EtStep

/EtStep

*PZS3PE-PZS3Xi*PYS3PE)**2

_PZS3PE-PZS3Xi*PXS3PE)**2

68
C

C

45

C

/ +(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

/ +(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

/ +(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (II, I)=-kS (3, II, I) *

PXS4Zt (II, i) =-kS (4, II, i) *

PYS3Zt(II,I)= kS (3,II, l)*

PYS4Zt (II, i) = kS (4, II, i)*

PZS3Zt (II, i) =-kS (3, II, i) *

PZS4Zt (II, l)---kS (4, II, i) *

(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

DO 45 i=2,II-i

PXS3Xi = (XS (3, i+l, I)-XS (3, i-l, I) )/2/XiStep

PXS4Xi = (XS (4, i+l, I)-XS (4, i-l, I) )/2/XiStep

PYS3Xi = (YS (3, i+l, I)-YS (3, i-l, i) )/2/XiStep

PYS4Xi = (YS (4, i+l, I)-YS (4, i-l, I) )/2/XiStep

PZS3Xi =(ZS(3,i+l,l)-ZS(3,i-l,l))/2/xiStep

PZS4Xi =(ZS(4,i+l,l)-ZS(4,i-l,l))/2/xiStep

PXS3PE = (XS (3, i, 2) -XS (3, i, I) )/EtStep

PXS4PE =(xS(4,i,2)-xS(4,i,I))/EtSteP

PYS3PE = (YS (3, i, 2) -YS (3, i, i) )/EtStep

PYS4PE= (YS (4, i, 2) -YS (4, i, i) )/EtStep

PZS3PE = (ZS (3, i, 2) -ZS (3, i, I) )/EtStep

PZS4PE= (ZS (4, i, 2)-ZS (4, i, I))/EtStep

LL3 = ((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE) **2

+ (PXS3Xi*PZS3PE-PZS3Xi*PXS3PE) **2

+ (PXS3Xi*PYS3PE-PYS3Xi*PXS3PE) **2) **0 .5

LL4 = ((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE) **2

+ (PXS 4Xi*PZS4PE-PZS4Xi*PXS4PE) **2

+ (PXS4Xi*PYS4PE-PYS4Xi*PXS4PE) *'2) **0.5

PXS3Zt (i, I) =-kS (3, i, I) * (PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

PXS4Zt (i, I) =-kS (4, i, I) * (PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

PYS3Zt (i, I)= kS (3,i,I)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

PYS4Zt (i,l)= kS (4,i,I)*(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

PZS3Zt (i, i) =-kS (3, i, I) * (PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

PZS4Zt (i,I)=-kS(4,i,I)*(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

CONTINUE

DO 60 j=2, JJ-I

PXS3Xi = (XS

PXS4Xi = (XS

PYS3Xi = (YS

PYS4Xi = (YS

PZS3Xi = (ZS

PZS4Xi = (ZS

PXS3PE = (XS

PXS4PE = (XS

PYS3PE = (YS

PYS4PE= (YS

PZS3PE= (ZS

3,2

4,2

3,2

4,2

3,2

4,2

3,1

4,1

3,1

4,1

3,1

-XS (3, I, j) )/XiStep

-XS(4,1, j))/XiStep

-YS(3,1, j))/XiStep

-YS(4,1, j))/XiStep

-ZS(3,1, j))/XiStep

-ZS(4,1, j))/XiStep

j+l)-XS(3,1, j-1))/2/EtStep

j+I)-XS (4, !, j-1))/2/EtStep

j+I)-YS (3,1, j-l))/2/EtStep

j+I)-YS(4,1, j-1))/2/EtStep

j+I)-ZS (3,1, j-1))/2/EtStep

PZS4PE=(ZS 4,1 j+l)-ZS(4,l,j-l))/2/EtStep

LL3:((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)**2

.(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

+(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4:((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

+(PXS4Xi-PZS4PE-PZS4Xi*PXS4PE)**2

6g
C

C

C

C

+(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (i, j)=-kS (3, i, j) *

PXS4Zt(I, j)=-kS (4,1, j)*

PYS3Zt (i, j) = kS (3,1, j)*

PYS4Zt (i, j) = kS (4,1, j)*

PZS3Zt (i, j)---kS (3, I, j) *

PZS4Zt (I, j)=-kS (4,1, j)*

(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4 °

(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

PXS3Xi=

PXS4Xi =

PYS3Xi =

PYS4Xi =

PZS3Xi=

PZS4Xi =

PXS3PE=

PXS4PE =

PYS3PE=

PYS4PE =

PZS3PE=

(XS (3, II, j)-XS (3

(XS (4, II, j)-XS (4

(YS (3, II, j)-YS (3

(YS(4,II, j)-YS(4

(ZS(3, II, j)-ZS(3

,II-l,j))/XiStep

,II-l,j))/XiStep

,II-l,j))/XiStep

,II-l,j))/XiStep

,II-l,j))/XiStep

(ZS (4, II, j)-ZS (4, II-l, j) )/XiStep

(XS (3, II, j+l) -XS (3, II, j-l) )/2/EtStep

(XS (4, II, j+l) -XS (4, II, j-l) )/2/EtStep

(YS (3, II, j+l) -YS (3, II, j-l) )/2/EtStep

(YS(4,II, j+I)-YS(4,II, j-I))/2/EtStep

(ZS(3, II, j+l)-ZS(3, II, j-l))/2/EtStep

PZS4PE=(ZS(4,II, j+I)-ZS(4,II, j-I))/2/EtStep

LL3=((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)**2

+(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

+(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

+(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

+(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt

PXS4Zt

PYS3Zt

PYS4Zt

PZS3Zt

PZS4Zt

II

II

II

II

II

II

j)=-kS (3, II,

j)=-kS(4,II,

j)= kS (3, II,

j)= kS(4, II,

j)=-kS(3, II,

j)=-kS(4, II,

j)*(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

j)*(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

j)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

j)*(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

j)*(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

j)*(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

DO 50 i=2, II-i

PXS3Xi = (XS (3, i+l, j) -XS

PXS4Xi =(xS(4,i+l, j)-XS

PYS3Xi =(YS(3,i+I, j)-YS

PYS4Xi =(YS(4,i+I, j)-YS

PZS3Xi =(ZS(3,i+I, j)-ZS

PZS4Xi =(ZS(4,i+I, j)-ZS

PXS3PE =(xS(3,i, j+I)-XS

PXS4PE=(XS(4,i, j+I)-XS

PYS3PE =(YS(3,i, j+I)-YS

PYS4PE =(YS(4,i, j+I)-YS

PZS3PE =(ZS(3,i, j+I)-ZS

3, i-l, j) )/2/XiStep

4,i-l, j))/2/XiStep

3,i-l, j))/2/XiStep

4,i-l, j))/2/XiStep

3,i-l, j))/2/XiStep

4, i-l, j) )/2/XiStep

3, i, j-l) )/2/EtStep

4,i, j-1))/2/EtStep

3, i, j-l) )/2/EtStep

4,i, j-1))/2/EtStep

3,i, j-1))/2/EtStep

PZS4PE=(ZS(4,i,j+I)-ZS 4,i,j-l))/2/EtStep

LL3=((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)**2

+(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

+(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

+(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

+(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (i, j)=-kS(3, i, j)*(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

PXS4Zt (i, j)=-kS (4, i, j)*(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

PYS3Zt (i, j)= kS (3, i, j)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL5

7O
5O

60

C

C

C

C

C

PYS4Zt (i, j)= kS(4,i, j)*(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

PZS3Zt (i, j)=-kS (3, i, j) * (PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

PZS4Zt (i, j)=-kS (4, i, j) * (PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

CONTINUE

CONTINUE

PXS3Xi = (XS (3,

PXS4Xi = (XS (4,

PYS3Xi = (YS (3,

PYS4Xi= (YS (4,

PZS3Xi= (ZS (3,

PZS4Xi = (ZS (4,

PXS3PE = (XS (3,

PXS4PE = (XS (4,

PYS3PE = (YS (3,

PYS4PE = (YS (4,

PZS3PE= (ZS (3,

2, JJ)-XS (3, I, JJ) )/XiStep

2, JJ)-XS (4, I, JJ) )/XiStep

2, JJ)-YS (3, i, JJ) )/XiStep

2, JJ)-YS (4, i, JJ) )/XiStep

2, JJ)-ZS (3, I, JJ) )/XiStep

2, JJ)-ZS (4, l,JJ) )/XiStep

I, JJ)-XS (3, l,JJ-l) )/EtStep

i, JJ) -XS (4, i, JJ-l) )/EtStep

i, JJ) -YS (3, i, JJ-l) )/EtStep

i, JJ) -YS (4, i, JJ-l) )/EtStep

I, JJ) -ZS (3, I, JJ-l) )/EtStep

PZS4PE=(ZS(4,I,JJ)-ZS(4, I,JJ-I))/EtStep

LL3=((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)**2

/ +(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

/ +(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

/ +(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

/ +(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (I, JJ) =-kS (3

PXS4Zt (l, JJ)=-kS (4

PYS3Zt(I,JJ)= kS (3

PYS4Zt(I,JJ)= kS (4

PZS3Zt(I,JJ)=-kS(3

PZS4Zt(I,JJ)=-kS(4

l,JJ) *

l,JJ) *

l,JJ) *

I, JJ) *

l,JJ)*

l,JJ) *

(PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

(PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

PXS3Xi = XS(3

PXS4Xi = XS(4

PYS3Xi= YS(3

PYS4Xi = YS(4

PZS3Xi = ZS(3

PZS4Xi = ZS(4

PXS3PE= (XS (3

PXS4PE = (XS (4

PYS3PE= (YS (3

PYS4PE = (YS (4

PZS3PE= (ZS (3

PZS4PE = (ZS (4

LL3 = ((PYS3Xi

/ + (PXS3Xi

/ + (PXS3Xi

LL4 = ((PYS4Xi

/ + (PXS4Xi

/ + (PXS4Xi

,II,JJ)-XS(3

II, JJ) -XS (4

II, JJ) -YS (3

II, JJ) -YS (4

II, JJ) -ZS (3

II, JJ) -ZS (4

II JJ)-XS(3

II JJ)-XS(4

II JJ)-YS(3

II JJ}-YS(4

II-l, JJ) )/XiStep

II-l, JJ) )/XiStep

II-l, JJ) )/XiStep

II-l, JJ) )/XiStep

II-l, JJ) )/XiStep

II-l, JJ) )/XiStep

II,JJ-l))/EtStep

II,JJ-l))/EtStep

II, JJ-l) )/EtStep

II, JJ-l) )/EtStep

II JJ)-ZS(3,II,JJ-l))/EtStep

II JJ)-ZS(4,II,JJ-l))/EtStep

*PZS3PE-PZS3Xi*PYS3PE)**2

*PZS3PE-PZS3Xi*PXS3PE)**2

*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

*PZS4PE-PZS4Xi*PYS4PE)**2

*PZS4PE-PZS4Xi*PXS4PE)**2

*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt (II

PXS4 Zt (I I

PYS3Zt (ii

PYS4 Zt (I i

PZS3Zt (!I

PZS4Zt (II

J J}

JJ)

JJ)

JJ)

J J}

JJ)

=-kS

=-kS

= kS

: kS

=-kS

=-kS

3, II, JJ) * (PYS3Xi*PZS3PE-PZS3Xi*PYS3PE

4, II, JJ) * (PYS4Xi*PZS4PE-PZS4Xi*PYS4PE

3, II,JJ)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE

4, II, JJ)* (PXS4Xi*PZS4PE-PZS4Xi*PXS4PE

3,!I,JJ)*(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE

4, II,JJ)* (PXS4Xi*PYS4PE-PYS4Xi*PXS4PE

/LL3

/LL4

/LL3

/LL4

/LL3

/LL4

71
C

C

C

C

65

DO 65 i=2, II-I

PXS3Xi=(XS(3,

PXS4Xi=(XS(4,

PYS3Xi=(YS(3,

PYS4Xi=(YS(4,

PZS3Xi=(ZS(3,

PZS4Xi=(ZS(4,

PXS3PE= (XS(3,

PXS4PE=(XS(4,

PYS3PE= (YS(3,

PYS4PE= (YS(4,

i+l, JJ)-XS (3, i-l, JJ) )/2/XiStep

i+l, JJ) -XS(4, i-l, JJ) )/2/XiStep

i+l, JJ) -YS(3, i-l, JJ) )/2/XiStep

i+l, JJ) -YS(4, i-l, JJ) )/2/XiStep

i+l, JJ) -ZS(3, i-l, JJ) )/2/XiStep

i+l, JJ) -ZS(4, i-l, JJ) )/2/XiStep

i,JJ)-XS (3, i, JJ-l) )/EtStep

i, JJ)-XS (4, i, JJ-l) )/EtStep

i,JJ)-YS (3, i, JJ-l) )/EtStep

i, JJ)-YS (4, i, JJ-l) )/EtStep

PZS3PE=(ZS(3,i,JJ)-ZS(3, i,JJ-I))/EtStep

LL3=((PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)**2

+(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)**2

+(PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)**2)**0.5

LL4=((PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)**2

+(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)**2

+(PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)**2)**0.5

PXS3Zt(i, JJ)=-kS (3, i, JJ) * (PYS3Xi*PZS3PE-PZS3Xi*PYS3PE)/LL3

PXS4Zt(i, JJ)=-kS (4, i, JJ) * (PYS4Xi*PZS4PE-PZS4Xi*PYS4PE)/LL4

PYS3Zt(i,JJ)= kS(3, i, JJ)*(PXS3Xi*PZS3PE-PZS3Xi*PXS3PE)/LL3

PYS4Zt(i,JJ)= kS(4,i, JJ)*(PXS4Xi*PZS4PE-PZS4Xi*PXS4PE)/LL4

PZS3Zt(i, JJ)=-kS (3, i, JJ) * (PXS3Xi*PYS3PE-PYS3Xi*PXS3PE)/LL3

PZS4Zt(i, JJ) =-kS (4, i, JJ) * (PXS4Xi*PYS4PE-PYS4Xi*PXS4PE)/LL4

CONTINUE

P2X00=0 0

P2XI0=0 0

P2X0!=0 0

P2XI!=0 @

P2Y00=0 0

P2YI0=0 0

P2Y01=0 0

P2YII=0 0

P2Z00=0 0

P2ZI0=0 0

P2Z01=0 0

P2ZII=0 @

Calculate the grid point locations everywhere.

DO 90 k=l, KK

Zeta = (k-I.) / (KK-I.)

DO 80 i=l,II

Xi =(i-l.)/(II-l.)

DO 70 j=l, JJ

Eta =(j-l.)/ (OJ-I

EtaNew =(EtStll (j

+(EtStl5 (j

ZetaNew = (ZtSt i (k

+ (ZtSt 5 (k

CALL FindHs(hl (j

CALL FindHs(h5(k

XPnt (i, j,k)=XPnt

+(XS(3,i, j)-hl(j)*XS(l,l,i)

-h2(j)*XS(2, l,i)

)

*(l.-Xi)+EtStl2(j)*Xi)*(l.-Zeta)

*(l.-Xi)+EtStl6(j)*Xi)*Zeta

* (i .-Xi) +ZtSt2 (k) *Xi) * (i .-Eta)

* (i .-Xi) +ZtSt6 (k) *Xi) *Eta

,h2(j),h3(j),h4(j),EtaNew,SigmaEt)

,h6 (k), h7 (k), h8 (k), ZetaNew, SigmaZt)

i, j,k)

72
-h3 (j) *PXSIPE (i, i)

-h4 (j) *PXS2PE (i, i) )*h5 (k)

+ (XS (4, i, j) -hl (j) *XS (I, KK, i)

-h2 (j) *XS (2, KK, i)

-h3 (j) *PXSIPE (i, KK)

-h4 (j) *PXS2PE (i, KK) )*h6 (k)

+ (PXS3Zt (i, j) - (hl (j) *PXS3Zt (i, I)

+h2 (j) *PXS3Zt (i, JJ)

+h3 (j) *P2X00+h4 (j) *P2X01) )*h7 (k)

+ (PXS4Zt (i, j) - (hl (j) *PXS4Zt (i, i)

+h2 (j) *PXS4Zt (i, JJ)

+h3 (j) *P2Xl0+h4 (j)*P2XlI))*h8 (k)

YPnt (i, j, k)=YPnt (i, j, k)

+(YS(3,i, j)-hl(j)*YS(l,l,i)

-h2 (j) *YS (2, i, i)

-h3 (j) *PYSIPE (i, I)

-h4 (j) *PYS2PE (i, I) )*h5 (k)

+(YS(4,i, j)-hI(j)*YS(I,KK, i)

-h2 (j) *YS (2, KK, i)

-h3 (j) *PYSIPE (i, KK)

-h4 (j)*PYS2PE(i,KK))*h6(k)

+ (PYS3Zt (i, j) - (hl (j) *PYS3Zt (i, I)

' +h2 (j) *PYS3Zt (i, JJ)

+h3(j)*P2Y00+h4(j)*P2Y01))*h7(k)

+(PYS4Zt(i, j)-(hI(j)*PYS4Zt(i,I)

+h2 (j) *PYS4Zt (i, JJ)

+h3(j)*P2Yl0+h4(j)*P2Yll))*hS(k)

ZPnt (i, j,k)=ZPnt(i, j,k)

+(ZS(3,i, j)-hl(j)*ZS(l,l,i)

-h2 (j) *ZS (2, i, i)

-h3 (j) *PZSIPE (i, i)

-h4 (j) *PZS2PE (i, I) )*h5 (k)

+(ZS(4,i, j)-hI(j)*ZS(I,KK, i)

-h2 (j) *ZS (2,KK, i)

-h3 (j) *PZSIPE (i, KK)

-h4 (j)*PZS2PE(i,KK))*h6(k)

+(PZS3Zt(i, j)-(hl(j)*PZS3Zt(i,l)

+h2 (j) *PZS3Zt (i, JJ)

+h3 (j) *P2Z00+h4 (j) *P2Z01) )*h7 (k)

+(PZS4Zt (i, j)-(hl(j)*PZS4Zt (i, i)

+h2(j)*PZS4Zt(i, JO)

+h3(j)*P2Zl0+h4(j)*P2Zll))*h8(k)

70 CONTINUE

80 CONTINUE

90 CONTINUE

C

C

C

RETURN

END

=C

C

SUBROUTINE PrGrid (XPnt,YPnt, ZPnt,II,JJ, KK, OutNum, MxGSiz)

C

C This SUBROUTINE prints (to output) the grid point x, y, and z

C

INTEGER i, j, k, If, JJ, KK, OutNum

C

REAL XFnt(MxGSiz,MxGsiz,MxGSiz),

coordinates.

73
C

C

i0

20

30

C

35

C

C

C-

$

$

YPnt (MxGSiz, MxGsiz, MxGSiz ),

ZPnt (MxGSi z, MxGsiz, MxGSi z)

WRITE (OutNum, *) II

WRITE (OutNum, *) JJ

WRITE (OutNum, *) KK

DO 30 i=l,II

DO 20 j=l,JJ

DO i0 k=I,KK

WRITE(OutNum, 35)

CONTINUE

CONTINUE

CONTINUE

XPnt (i, j, k), YPnt (i, j, k), ZPnt (i, j, k)

FORMAT (IX, Fi0.6, 3x, FI0.6, 3X, FI0.6)

RETURN

END

C

C

C

C

C

C

C

C

C

SUBROUTINE RdGrIn(II,JJ, KK, NSurfs, SigmaXi,SigmaEt,SigmaZt,

$ kXil,kXi2,kEtal,kEta2,kZetal,kZeta2,InNum)

This SUBROUTINE reads in the desired grid information for grid control.

INTEGER StrXi, StrEt, StrZt, InNum, II, JJ, KK

REAL kXil, kXi2, kEta], kEta2, kZetal, kZeta2,

$ BetaXi, BetaEt, BetaZt, SigmaXi, SigmaEt, SigmaZt

READ(InNum,*) NSurfs

READ (InNum, *) II

READ (InNum, *) JJ

READ (InNum, *) KK

READ

READ

READ

READ

READ

READ,

READ

READ,

READ

InNum, *) SigmaXi

InNum,*} SigmaEt

InNum,*) SigmaZt

InN_m, *) kXil

InN_m, *) kXi2

InNum,*) kEtal

InNum,*) kEta2

InNum,*) kZetal

InNum,*) kZeta2

RE TUF_

END

C

C

C

SUBROUYI}_E RdCvIn

$

This SUBROUT!_E reads

(x, y, z, NDPts, CrvNum, Tensn, InNum, MxBPts,

StrTp, Betal, Beta2)

in the information concerning discrete points on

74
C

C

C

C

C

C

I0

C

C

C

the boundaries. This information

boundary approximation curves.

INTEGER CrvNum, i, NDPts(4),

is used for generating

InNum, StrTp

REAL x (4, MxBPts), y (4,MxBPts),

$ z(4,MxBPts), Tensn(4)

READ(InNum,*) Tensn(CrvNum)

READ(InNum,*) NDPts(CrvNum)

DO i0 i=l,NDPts(CrvNum)

READ(InNum,*) x (CrvNum,i) ,

CONTINUE

y (CrvNum,i) , z (CrvNum,i )

READ(InNum,*) StrTp

IF (StrTp.NE. 4)THEN

READ(InNum,*)Betal

ELSE

READ(InNum,*) Betal, Beta2

ENDIF

RETURN

END

spline-fitted

C

C

C

C

C

C

C

C

I0

C

C

C

SUBROUTINECalcS (x,y,z,s,NDPts,CrvNum, MxBPts)

This SUBROUTINEcalculates the spline parameter, s,

arc length.

INTEGER NDPts(4 , CrvNum, i

as an approximate

REAL x (4,MxBPts

$ z (4,MxBPts

S(CrvNum,i) =0.0

, y (4,MxBPts) ,

, s (4,MxBPts)

DO i0 i=2,NDPts(CrvNum)

s (errS'Sum,i)=s (CrvNmm,i-l)

$ +SQRT( (x(CrvNum,i) -x (CrvNum,i-1))**2

$ + (y (CrvNum,i) -y (CrvNum,i-l) )**2

$ + (z (CrvNum,i) -z (CrvNum,i-l) )**2)

CONTINUE

RETU_$

END

SUBROUTINESplMat (Diag,OfDiag,Right,w,s,NDPts,T,CrvNum, MxBPts)

C

C This SUBROUT!_[Eforms the parametric tension spline matrix for a

C particular boundary curve data set.

C

INTEGER _,_ _,_ts(4), CrvNum

75
C

C

i0

C

C

C=- -

REAL Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts),

$ w(4,MxBPts), s(4,MxBPts), T, h, hm

Diag (I) =i. 0

Of Diag (I) =0.0

Right (I) =0.0

DO I0 i=2,NDPts(CrvNum)-I

h=s (CrvNum, i+l) -s (CrvNum, i)

hm=s (CrvNum, i) -s (CrvNum, i-l)

Diag (i) = (T*COSH (T*hm)/SINH (T*hm) -I/hm+T*COSH (T'h)/SINH (T'h)

$ -l/h)/T**2

OfDiag (i) = (I/h-T/SINH (T'h))/T**2

Right (i )= (w (CrvNum, i+ 1 )-w (CrvNum, i ) )/h

$ - (w (CrvNum, i) -w (CrvNum, i-l) )/hm

CONTINUE

Diag (NDPt s (CrvNum)) =I. 0

OfDiag (NDPt s (CrvNum) -I) =0.0

Right {NDPts {CrvNum) ) =0.0

RETURN

END

C

SUBROUTINE SplSlv (Diag,OfDiag,Right,Derv2,NDPts,CrvNum, MxBPts)

C

C This SUBROUTINE solves the diagonally dominant parametric tension

C spline matrix for a given data set using the Gauss-Seidel iteration.

C Convergence is assumed after 20 iterations.

C

INTEGER i, j, NDPts(4), CrvNum

C

REAL Diag(MxBPts), OfDiag(MxBPts), Right(MxBPts),

$ Defy2 (4,MxBPts)

C

C Initialize the second derivative matrix to all zeroes.

C

DO I0 i=l,NDPts(CrvNum)

Defy2 (CrvNum, i) =0 .0

I0 CONTINUE

C

C Calculate the second derivative

C the Gauss-Seidel method.

C

DO 30 j=i,20

DO 20 i=2,NDPts(CrvNum)-I

20

3O

C

C

Cvalues using 20 iterations of

Derv2(CrvNum, i)=(Right(i)-OfDiag(i)*Derv2(CrvNum, i+l)

$ -OfDiag(i-l)*Derv2(CrvNum, i-l))

$ /Diag(i)

CONTI_IUE

CONTINUE

RETURn[

END

76
C

C

C

C

C

C

C

C

C

FUNCTION SplVal (s,w,Derv2,sval,T,n, CrvNum,MxBPts)

This real function finds the w-value (x-value or y-value) corresponding

to a specified s-value using the parametric tension spline curve

generated for a particular boundary curve data set.

INTEGER n, CrvNum

REAL s(4,MxBPts), w(4,MxBPts), Derv2(4,MxBPts),

$ sval, T, h, Interim, Templ, Temp2

Templ=sval-s (CrvNum, n)

h=s (CrvNum, n+l)-s (CrvNum, n)

Temp2=s (CrvNum, n+l )-sval

Interim=Derv2 (CrvNum, n)/T**2*SINH (T*Temp2)/SINH (T'h)

$ + (w (CrvNum, n) -Derv2 (CrvNum, n)/T**2) *Temp2/h

SplVal=Interim+Derv2 (CrvNum, n+l )/T**2*SINH (T*Templ)

$ /SINH (T'h} + (w (CrvNum, n+l)

$ -Defy2 (CrvNum, n+l) /T**2) *Templ/h

RETURN

END

C

C

C

C

C

C

C

C

C

C--

SUBROUTINE PTSpln(x,y,z,s,XDerv2,YDerv2,ZDerv2,Diag, Of Diag,

$ Right, NDPt s, Tensn, CrvNum, MxBPt S )

This SUBROUTINE forms the main routine for the parametric tension

spline process.

INTEGER NDPts(4), CrvNum

REAL Diag(MxBPts), Of Diag(MxBPts), Right (MxBPts),

$ XDerv2(4,MxBPts), YDerv2(4,MxBPts),

$ ZDerv2(4,MxBPts), Tensn,

$ x (4,MxBPts) , y (4,MxBPts) ,

$ z(4,MxBPts), s (4,MxBPts)

CALL

CALL SplMat(Diag

CALL SplSlv(Diag

CALL SplMat(Diag

CALL SplSlv(Diag

CALL SplMat(Diag

CALL Sp!Slv(Diag

CalcS (x, y, z, s,NDPts, CrvNum, MxBPts)

Of Diag, Right, x, s, NDPt s, Tensn, CrvNum, MxBPt s)

OfDiag, Right, XDerv2, NDPt s, CrvNum, MxBPt s)

Of Diag, Right, y, s, NDPts, Tensn, CrvNum, MxBPts)

OfDiag, Right, YDerv2, NDPt s, CrvNum, MxBPt s )

Of Diag, Right, z, s, NDPts, Tensn, CrvNum, MxBPts)

Of Diag, Right, ZDerv2, NDPts, CrvNum, MxBPt s)

C

C

C This

SU_io_.I.,- Dr _T_c? _J? FindHs(hi r h2,h3,h4,n t s)

__-_._I_[E computes the h factors used in Hermite interpolation.

77
C

REAL hl, h2, h3, h4, n, s

/ al, a2, a3, a4, a, b, bbaa, sh, ch, shsn, shsnl

IF(s.NE.0)THEN

sh=sinh(s)

ch=cosh(s)

a2=sh/{2.*sh-s*ch-s)

al=l-a2

a=s*ch-sh

b=sh-s

bbaa=b*b-a*a

a3=-a*sh/bbaa

a4=b*sh/bbaa

shsn=sinh (s'n)/sh

shsnl=sinh (s* (i .-n) )/sh

hl=a2* (shsnl-shsn) +al* (I .-n) +a2*n

h2=a2* (shsn-shsnl) +a2* (I .-n) +al*n

h3=a3 * ((l.-n)-shsnl)+a4*(n-shsn)

h4=a4 * (shsnl-(l.-n))+a3*(shsn-n)

ELSE

h!= 2*n**3-3*n**2+l

h2=-2*n**3+3*n**2

h3 = n**3-2*n**2+n

h4 = n**3-n**2

ENDIF

RETURN

END

C

_-- ..........................

C

SUBROUTINE SplInt(n,s,SValue,NDPts,CurCrv,MxBPts)

C

C This SUBROUTINE finds the proper interval in which a point on a specified

C boundary lies. The interval indicates which initial data points the

C point in question lies between and thus which spline coefficients to

C use.

C

C

C

C

i0

C

C

C

INTEGER i, n, CurCrv, NDPts(4)

REAL Temp, SValue, s(4,MxBPts]

n=l

i=NDPts(CurCrv)

IF ((n.EQ.I).AND. (i.GT.I)) THEN

I=I-I

Temp=SValue-s(CurCrv, i)

IF (Temp. GT.0.0) THEN

n=i

ENDIF

GOTO !0

ENDIF

78
C

C=

C

RETURN

END

SUBROUTINE FAINew (AiNew, Alpha, B, Str)

C

C This SUBROUTINE computes the new Alpha value after stretching as

C AiNew. Alpha is a dummy variable representing either Xi, Eta or Zeta.

C

INTEGER Str

C

C

C

C

C

CREAL Alpha, Templ, Temp2, B2, AINew, B

AiNew=Alpha

Templ=(B+l)/(B-l)

IF (Str.EQ.I) THEN

Temp2=Templ**(l-Alpha)

AiNew=((B+l)-(B-l)*Temp2)/(Temp2+l)*l

ENDIF

IF (Str.EQ.2) THEN

B2=0

Temp2=Templ**((Alpha-B2)/(l-B2))

AINew=((B+2*B2)*Temp2-B+2*B2)/((2*B2+I)*(I+Temp2))

ENDIF

IF (Str.EQ.3) THEN

B2=0.5

Temp2=Templ ** ((Alpha-B2) / (I-B2))

AiNew = ((B+2*B2) *Temp2-B+2*B2) / ((2"B2+I) * (l+Temp2))

END IF

RETURN

END

C

SUBROUTINE EdgPts(XI,X2,X3,X4,YI,Y2,Y3,Y4,ZI,Z2,Z3, Z4,AL, BL,

$ AAStep, BBStep, x, y, z, s, zx, zy, zz, NDPts, Tensn,

$ StrAA, StrBB, BetaAA, BetaBB, MxBPts, MxGSiz )

C

C This SUBROUTINE calculates the grid point locations along the surface

C edges.

C

INTEGER ACt, BCt, nl, n2, n3, n4,

$ AL, BL, StrAA, StrBB, NDPts(4)

C

C

REAL

$

$

$

S

S

$

S

z_., BB, AANew, BBNew, SI, $2, $3, $4, BBStep, AAStep,

SIAAR, S2_.AR, S3BBR, S4BBR,

X! (MxGSiz), X2(MxGSiz), X3(MxGSiz), X4 (MxGSiz),

Y1 (MxGSiz), Y2(MxGSiz), Y3(MxGSiz), Y4 (MxGSiz),

Z1 (MxGCiz), Z2(MxGSiz), Z3(MxGSiz), Z4 (MxGSiz),

x(4,MxBPts), y (4,MxBPts) , z (4,MxBPts) ,

s(4,MxBPts), zx(4,MxBPts), zy(4,MxBPts),

zz (4,MxBPts), Tensn (4), BetaAA, BetaBB

79
C

C

C

C

C

C

i0

20

C

SIAAR=s (i, NDPts (I))

S2AAR=s (2, NDPt s (2) )

S3BBR--s (3, NDPts (3))

S4BBR=s (4, NDPts (4))

Calculate the grid point locations along boundaries 1 and 2.

AA-0.0

C

C'--

DO I0 ACt=I,AL

CALL FAiNew (AANew, AA, BetaAA, StrAA)

S I=AANew*S IAAR

S2=AANew*S2AAR

CALL SplInt (nl, s, Sl,NDPts,l,MxBPts)

CALL SplInt (n2, s, S2,NDPts, 2,MxBPts)

Xl (ACt) =SplVal (s, x, zx, Si,Tensn(1),nl, l MxBPts)

X2 (ACt) =SplVal (s, x, zx, $2, Tensn (2), n2, 2 .MxBPts)

Y1 (ACt)--SplVal (s, y, zy, SI, Tensn (I), nl, 1 .MxBPts)

Y2 (ACt) =SplVal (s, y, zy, $2, Tensn (2), n2, 2 .MxBPts)

Zl (ACt) =SplVal (s, z, zz, SI, Tensn (I), nl, 1 .MxBPts)

Z2(ACt)=SplVal(s,z,zz,S2,Tensn(2),n2,2.MxBPts)

AA=AA+AASt ep

CONTINUE

Calculate the grid point locations along boundaries 3 and 4.

BB=0.0

DO 20 BCt=I,BL

CALL FAINew(BBNew,BB,BetaBB, StrBB)

S3=BBNew*S3BBR

S 4=BBNe w* S 4BBR

CALL Splint (n3,s,S3,NDPts, 3,MxBPts)

CALL Splint (n4 s,S4,NDPts, 4,MxBPts)

X3(BCt)=SplVal s,x,zx,S3,Tensn(3),n3,3,MxBPts)

X4(BCt)=SpIVal s,x,zx,S4,Tensn(4),n4,4,MxBPts)

Y3(BCt)=SpIVal s,y, zy, S3,Tensn(3),n3,3,MxBPts)

Y4 (BCt)=SpIVal s,y, zy,S4,Tensn(4),n4,4,MxBPts)

Z3(BCt)=SpIVal s,z,zz,S3,Tensn(3),n3,3,MxBPts)

Z4(BCt)=SplVal s, z,zz,S4,Tensn(4),n4,4,MxBPts)

BB=BB+BBStep

CONTINUE

RETURN

END

SUBROUTINE EdgDer

$

$

$

PXIPAA, PX2PAA, PYIPAA, PY2PAA, PZIPAA, PZ2PAA,

PX3PBB, PX4PBB, PY3PBB, PY4PBB, PZ3PBB, PZ4PBB,

XI, X2, X3, X4, YI, Y2, Y3, Y4, Zl, Z2, Z3, Z4,

AL, BL, MxGSiz)

INTEGER ACt, BCt, AL, BL

REAL AAStep, BBStep,

$ PY3FBB(MxGSiz),

$ PZ4PBB(MxGSiz),

PX3PBB (MxGSiz),

PY4PBB (MxGSiz) ,

PXlPAA (MxGSiz) ,

PX4PBB(MxGSiz),

PZ3PBB(MxGSiz),

PX2PAA(MxGSiz),

8O
C

C

C

C

C

C

C

C

i0

C

C

C

C

C

$ pYIPAA(MxGSiz), PY2PAA(MxGSiz), PZlPAA(MxGSiz),

$ PZ2PAA(MxGSiz), XI(MxGSiz), X2(MxGSiz), X3(MxGSiz),

$ X4(MxGSiz), YI(MxGSiz), Y2(MxGSiz), Y3(MxGSiz),

$ Y4(MxGSiz), Zl(MxGSiz), Z2(MxGSiz), Z3(MxGSiz),

$ Z4 (MxGSiz)

Calculate step size in the AA and BB directions.

AAStep=l. / (AL-I.)

BBStep=l. / (BI-I.)

Calculate the derivative values along the constant AA boundaries.

PXIPAA (i

PX2PAA (I

PYIPAA (i

PY2PAA (1

PZ IPAA (1

PZ2PAA (i

= (Xl (2)-Xl (i))/AAStep

= (X2 (2)-X2 (I))/AAStep

= (YI (2)-YI (I))/AAStep

= (Y2 (2)-Y2 (I))/AAStep

=(ZI(2)-ZI(1))/AAStep

= (Z2 (2)-Z2 (i))/AAStep

PXlPAA (AL)

PX2PAA (AL)

PYIPAA (AL)

PY2PAA (AL)

P Z IP_-A (AL)

PZ2PAA (AL)

=(XI(AL) -XI(AL-I

= (x2 (AL) -X2 (AL-I

=(YI(AL) -YI(AL-I

=(Y2(AL) -Y2(AL-I

=(ZI(AL) -ZI(AL-I

=(Z2(AL) -Z2(AL-!

/AAStep

/AAStep

/AAStep

/AAStep

/AAStep

/AAStep

DO

PXIPAA (ACt

PX2PAA (ACt

PYIPAA (ACt

PY2PAA (ACt

PZIPAA (ACt

PZ2PAA (ACt

CONT I NUE

I0 ACt=2,AL-I

(Xl (ACt+i

(X2 (ACt+l

(YI (ACt+l

(Y2 (ACt+l

(Zl (ACt+I

(Z2 (ACt+I

-Xl (ACt-l))/2/AAStep

-X2 (ACt-l)) /2/AAStep

-YI (ACt-l)) /2/AAStep

-Y2 (ACt-l))/2/AAStep

-ZI (ACt-l))/2/AAStep

-Z2 (ACt-l)) /2/AAStep

Calculate the derivative values along the constant BB boundaries.

PX3PBB (i) = (X3 (2

PX4PBB(1) = (X4 (2

PY3PBB (i) = (Y3(2

PY4PBB (I) = (Y4(2

PZ3PBB (I) = (Z3(2

PZ4PBB (i) = (Z4(2

-X3 (I))/BBStep

-X4 (i))/BBStep

-Y3 (i))/BBStep

-Y4 (1))/BBStep

-Z3 (I))/BBStep

-Z4 (i)) /BBStep

PX3PBB (BL) = (X3

PX4PBB (BL) = (X4

PY3PBB (BL)= (Y3

PY4PBB (BL) = (Y4

PZ3FBB (BL)= (Z3

PZ4PBB (BL)= (Z4

(BL)

(BL)

(BL)

(BL)

(BL)

(BL)

-X3 (BL-I

-X4 (BL-I

-Y3 (BL-I

-Y4 (BL-I

-Z3 (BL-I

-Z4 (BL-I

/BBStep

/BBStep

/BBStep

/BBStep

/BBStep

/BBStep

DO 20 BCt=2,BL-I

FX3FBB (Btt) =

PX4PBB (Bet) =

PY3PBB (Bet) =

FY4PBE (BCt] =

PZ3._BB (BCt) =

(X3 (Bet+l) -X3 (BCt-I))/2/BBStep

(X4 (BCt+I)-X4(BCt-I))/2/BBStep

(Y3(BCt+I)-Y3(BCt-I))/2/BBStep

(Y4 (BCt+I)-Y4 (BCt-I))/2/BBStep

(Z3 (BCt+I)-Z3 (BCt-I)) /2/BBStep

81
20

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

PZ4PBB(BCt)= (Z4(BCt+I)-Z4 (BCt-l))/2/BBStep

CONTINUE

RETURN

END

SUBROUTINE

$

$

$

$

$

$

TwoBnd (XS, YS, ZS, SrfNum, AL, BL, SigmaBB, kl, k2,

StrB, hl, h2, h3, h4, Xl, X2, X3, X4,

YI, Y2, Y3, Y4, Zl, Z2, Z3, Z4, PXlPBB, PX2PBB,

PYIPBB, PY2PBB, PZIPBB, PZ2PBB, PXlPAA, PX2PAA,

PYIPAA, PY2PAA, PZIPAA, PZ2PAA, PX3PBB, PX4PBB,

PY3PBB, PY4PBB, PZ3PBB, PZ4PBB,

MxBCvs, MxGSiz, MxSr fs )

This SUBROUTINE calculates the interior grid point locations between

two specified boundaries (I and 2) by using transfinite Hermite

interpolation.

INTEGER ACt, BCt, AL, BL, StrAA, StrBB, SrfNum, Edgl, Edg2,

$ Edg3, Edg4

REAL

$

$

$

$

$

$

$

$

REAL

$

$

$

$

$

$

$

$

$

$

$

$

AA, BB, AANew, BBNew, LLI, LL2,

Boxli, Boxlj, Boxlk, Box2i, Box2j, Box2k,

Templi, Templj, Templk, Temp2i, Temp2j, Temp2k,

kl(MxSrfs,MxGSiz), k2(MxSrfs,MxGSiz),

BetaAA, BetaBB, BBStep, AAStep,

hl(MxGSiz), h2(MxGSiz), h3(MxGSiz), h4 (MxGSiz),

Xl(MxGSiz), X2(MxGSiz), X3(MxGSiz), X4(MxGSiz),

YI(MxGSiz), Y2(MxGSiz), Y3(MxGSiz), Y4(MxGSiz),

Zl(MxGSiz), Z2(MxGSiz), Z3(MxGSiz), Z4(MxGSiz)

PXIPBB MxGSiz , PX2PBB(MxGSiz),

PYIPBB MxGSiz

PZIPBB MxGSiz

PXIPAA MxGSiz

PYIPAA MxGSiz

PZIPAA MxGSiz

PX3PBB MxGSiz

PY2PBB(MxGSiz),

PZ2PBB(MxGSiz),

PX2PAA(MxGSiz),

PY2PAA(MxGSiz),

PZ2PAA(MxGSiz),

PX4PBB(MxGSiz),

PY3PBB MxGSiz PY4PBB(MxGSiz),

PZ3PBB MxGSiz , PZ4PBB(MxGSiz),

XS(MxSrfs,MxGSiz,MxGSiz),

YS(MxSrfs,MxGSiz,MxGSiz),

ZS(MxSrfs,MxGSiz,MxGSiz),

StrB(MxGSiz,MxBCvs)

Calculate the step size in the AA and BB directions.

AAStep=I./(AL-I.)

BBStep=i./(BL-I.)

Calculate the edge numbers for the surface 'SrfNum'

Edg!=(SrfNum-!)*4 + 1

Edg2=(Srf_$um-l)*4 + 2

Edg3=(Srfhum-l)*4 + 3

Edg4=(Srf_um-l)*4 + 4

82
C

C

C

C

C

C

C

Calculate the derivative values for grid line orthogonality.

AA=0.0

DO 20 ACt=I,AL

Boxli-- AA*(PYIPAA(AL)*PZ4PBB(1)

$ -PZIPAA (AL) *PY4PBB (1 ) )

$ +(I-AA)*(PYIPAA(1) *PZ3PBB(1)

$ -PZIPAA(1) *PY3PBB(1))

Boxlj= AA*(PXlPAA(AL)*PZ4PBB(1)

$ -PZIPAA (AL) *PX4PBB (i))

$ +(I-AA)*(PXIPAA(1) *PZ3PBB(1)

$ -PZIPAA(1) *PX3PBB(1))

Boxlk = AA*(PXlPAA(AL)*PY4PBB(1)

$ -PYIPAA (AL) *PX4PBB (I))

$ +(I-AA)*(PXlPAA(1) *PY3PBB(1)

$ -PYIPAA(1) *PX3PBB(1))

Box2i = AA*(PY2PAA(AL)*PZ4PBB(BL)

-PZ2PAA (AL) *PY4PBB (BL))

+(I-AA)*(PY2PAA(i) *PZ3PBB(BL)

-PZ2PAA(1) *PY3PBB(BL))

Box2j= AA*(PX2PAA(AL)*PZ4PBB(BL)

-PZ2PAA (AL) *PX4PBB (BL))

+(I-AA)*(PX2PAA(1) *PZ3PBB(BL)

-PZ2PAA(1) *PX3PBB(BL))

Box2k= AA*(PX2PAA(AL)*PY4PBB(BL)

-PY2PAA (AL) *PX4PBB (BL))

+(I-AA)*( PX2PAA_i) *PY3PBB(BL)

-PY2PAAI

Templi=PYIP_(ACt)*Boxl

Templj=PXIPAA(ACt)*Boxl

Templk=PXIPAA(ACt)*Boxl

Temp2i=PY2PAA(ACt)*Box2

Temp2j=PX2P_A(ACt)*Box2

i) *PX3PBB(BL))

k+PZIPAA(ACt)*Boxlj

k-PZIPAA(ACt)*Boxli

j+PYIPAA(ACt)*Boxli

k+PZ2PAA(ACt)*Box2j

k-PZ2PAA(ACt)*Box2i

Temp2k=PX2PAA(ACt)*Box2j+PY2PAA(ACt)*Box2i

LLl=(Templi**2+Templj**2+Templk**2)**0.5

LL2=(Temp2i**2+Temp2j**2+Temp2k**2)**0.5

PXIPBB (ACt

PX2PBB (ACt

PYIPBB (ACt

PY2PBB (ACt

PZlPBB (ACt

PZ2PBB (ACt

=-kl(SrfNum, ACt)*Templi/LLl

=-k2(SrfNum, ACt)*Temp2i/LL2

= kl (SrfNum, ACt) *TempIj/LLI

= k2 (SrfNum, ACt) *Temp2j/LL2

= kl (SrfNum, ACt) *Templk/LLl

= k2(SrfNum, ACt)*Temp2k/LL2

AA=AA+AAStep

20 CONTINUE

Calculate the interior grid point locations.

DO 40

DO

ACt=I,AL

30 BCt=I, BL

_n=(ACt-I.)/ (AL-!.)

BBNew=StrB (BCt, Edg3) * (I. -AA) +StrB (BCt, Edg4) *AA

CALL FindHs(hl(BCt),h2(BCt),h3(BCt),h4 (BCt),BBNew, SigmaBB)

XS (SrfNum, AC_,BCt] = hi (BCt)*Xl (ACt)+h2 (BCt)*X2(ACt)

+h3 (BCt) *PXlPBB (ACt) +h4 (BCt) *PX2PBB (ACt)

83
30

40

C

YS(SrfNum, ACt,BCt)= hI(BCt)*YI(ACt)+h2(BCt)*Y2(ACt)

$ +h3 (BCt) *PYIPBB (ACt) +h4 (BCt) *PY2PBB (ACt)

ZS(SrfNum, ACt,BCt) = hI(BCt)*ZI(ACt)+h2(BCt)*Z2(ACt)

$ +h3 (BCt) *PZIPBB (ACt)+h4 (BCt) *PZ2PBB (ACt)

CONTINUE

CONTINUE

RETURN

END

C

C

SUBROUTINE ForBnd (XS, YS, ZS, SrfNum, AL, BL, SigmaAA, SigmaBB, k3, k4,

$ StrB, hl, h2, h3, h4, h5, h6, h7, h8,

$ Xl, X2, X3, X4, YI, Y2, Y3, Y4, ZI, Z2, Z3, Z4,

$ PXIPBB, PX2PBB, PYIPBB, PY2PBB, PZIPBB, PZ2PBB,

$ PXlPAA, PX2PAA, PYIPAA, PY2PAA, PZlPAA, PZ2PAA,

$ PX3PBB, PX4PBB, PY3PBB, PY4PBB, P Z3PBB, PZ4PBB,

$ PX3PAA, PX4PAA, PY3PAA, PY4PAA, PZ3PAA, PZ4PAA,

$ MxBCvs, MxGSiz, MxSr fs )

This SUBROUTINE adjusts the grid so that the other two boundaries

surface are mapped correctly using transfinite Hermite

C

C

C

C

C

C

C

C

(3 and 4) of the

interpolation.

INTEGER ACt, BCt, AL, BL, StrAA, StrBB, i, j, SrfNum,

$ Edgl, Edg2, Edg3, Edg4

REAL

$

$

$

$

$

$

$

$

$

$

$

REAL

$

$

$

$

$

$

$

$

$

$

$

$

$

$

$

Calculate

AA, BB, AANew, BBNew, LL3, LL4,

Box3i, Box3j, Box3k, Box4i, Box4j, Box4k,

Temp3i, Temp3j, Temp3k, Temp4i, Temp4j, Temp4k,

P2Y00, P2Y01, P2YI0, P2YII, P2X00, P2X01, P2XI0,

P2Z00, P2Z01, P2ZI0, P2ZII,

k3(MxSrfs,MxGSiz), k4(MxSrfs,MxGSiz),

BetaAA, BetaBB, BBStep,

hl(MxGSiz), h2(MxGSiz),

hS(MxGSiz), h6(MxGSiz),

Xl(MxGSiz), X2(MxGSiz),

YI(MxGSiz), Y2(MxGSiz),

ZI(MxGSiz), Z2(MxGSiz),

PX!PBB MxGSiz),

PYIPBB MxGSiz),

PZIPBB MxGSiz),

PXIPAA MxGSiz),

PYIPAA MxGSiz),

PZIPAA MxGSiz),

PX3PBB MxGSiz),

PY3PBB MxGSiz),

PZ3PBB MxGSiz),

PX3PAA MxGSiz),

PY3PAA MxGSiz),

PZ3PAA MxGSiz),

AAStep,

h3 (MxGSiz),

h7 (MxGSiz),

X3 (MxGSiz) ,

Y3 (MxGSiz),

Z3 (MxGSiz),

PX2PBB(MxGSiz),

py2PBB(MxGSiz),

pZ2PBB(MxGSiz)

PX2PAA(MxGSiz)

py2PAA(MxGSiz)

PZ2PAA(MxGSiz)

PX4PBB(MxGSiz)

py4PBB(MxGSiz}

PZ4PBB(MxGSiz),

PX4PAA(MxGSiz),

py4PAA(MxGSiz),

PZ4PAA(MxGSiz),

XS (MxSr fs, MxGSiz, MxGSiz) ,

YS(MxSrfs,MxGSiz,MxGSiz),

ZS (MxSrfs, MxGSiz, MxGSiz) ,

StrB (MxGSiz, MxBCvs)

h4 (MxGSiz) ,

h8 (MxGSiz) ,

X4 (MxGSiz) ,

Y4 (MxGSiz) ,

Z4 (MxGSiz)

P2XII,

the step size for directions AA and BB.

84
C

C

C

C

C

C

C

C

C

AAStep=l. / (AL-I.)

BBStep=l. / (BL-I.)

Calculate edge numbers for the

Edgl=(SrfNum-l)*4 + 1

Edg2=(SrfNum-l)*4 + 2

Edg3=(SrfNum-l)*4 + 3

Edg4=(SrfNum-l)*4 + 4

surface 'SrfNum'

Calculate the derivative values for grid line orthogonality.

BB=0.0

DO 20 BCt=I,BL

Box3i = BB*

+(I-BB) *

Box3j = BB*

+(I-BB)*

Box3k = BB*

+(I-BB)*

Box4i = BB*

+(I-BB)*

Box4j= BB _

+(I-BB)*

Box4k = BB*

+(I-BB)*

Temp3i=PZ3PBB(BCt

Temp3j=PZ3PBB(BCt

Temp3k:PY3PBB(BCt

Temp4i=PZ4PBB(BCt

Temp4j=PZ4PBB(BCt

Temp4k=PY4PBB(BCt *Box4

LL3=(Temp3i**2+Temp3j **

LL4=(Temp4i**2+Temp4j**

PY2PAA (1 )*PZ3PBB (BL)

-PZ2PAA (I) *PY3PBB (BL))

PYIPAA (I) *PZ3PBB (I)

-PZIPAA (I) *PY3PBB (I))

PX2PAA (i) *PZ3PBB (BL)

-PZ2PAA(I *PX3PBB(BL))

PXIPAA (i) *PZ3PBB (I)

-PZIPAA(I *PX3PBB(1))

PX2PAA (i) *PY3PBB (BL)

-PY2PAA(I *PX3PBB(BL))

PXlPAA (I) *PY3PBB (i)

-PYIPAA(I *PX3PBB(1))

PY2PAA (I) *PZ4PBB (BL)

-PZ2PAA(I *PY4PBB(BL))

PYIPAA (i) *PZ4PBB (I)

-PZIPAA(I *PY4PBB(1))

PX2PAA (I) *PZ4PBB (BL)

-PZ2PAA (I *PX4PBB (BL))

PXIPAA (i) *PZ4PBB (I)

-PZIPAA (i) *PX4PBB (i))

PX2PAA (I) *PY4PBB (BL)

-PY2PAA (I) *PX4PBB (BL))

PXIPAA (I) *PY4PBB (i)

-PYIPAA (i) *PX4PBB (I))

*Box3 j+PY3PBB (BCt) *Box3k

*Box3i-PX3PBB (BCt) *Box3k

*Box3i+PX3PBB (BCt) *Box3 j

*Box4 j+PY4PBB (BCt) *Box4k

*Box4i-PX4PBB (BCt) *Box4k

i+PX4PBB (BCt) *Box4 j

2+Temp3k**2) **0.5

2+Temp4k**2) **0.5

PX3PAA(BCt)= k3(SrfNum

PX4PAA(BCt) = k4 (SrfNum

PY3PAA(BCt) = k3(SrfNum

PY4PAA(BCt)= k4 (SrfNum

PZ3P;_A (ECt) =-k3 (SrfNum

PZ4PAA(BCt)=-k4 (SrfNum

BCt

BCt

BCt

BCt

BCt

BCt

*Temp3i/LL3

*Temp4i/LL4

*Temp3j/LL3

*Temp4j/LL4

*Temp3k/LL3

*Temp4k/LL4

BB=BB+BBSte_ ,

85
20

C

C Set

C

CONTINUE

the cross-derivative terms equal to zero.

P2X00=0 0

P2XI0=0 0

P2X01=0 0

P2XlI=0 0

P2Y00=0 0

P2YI0=0 0

P2Y01=0 0

P2YII=0.0

P2Z00=0.0

P2ZI0=0.0

P2Z01=0.0

P2ZlI=0.0

C

C Calculate the

C

DO 40 i=l

grid point locations everywhere.

,AL

$

$

$

$

$

$

$

$

$

$

$

$

$

$

DO 30 j=I,BL

AA= (i-l.) / (AL-I.)

BB =(j-I.)/(BL-I.)

AANew=StrB (i, Edgl) * (i .-BB) +StrB (i, Edg2) *BB

BBNew=StrB (j, Edg3) * (i .-AA) +StrB (j, Edg4) *AA

CALL FindHs(hl(j),h2(j),h3(j),h4(j),BBNew, SigmaBB)

CALL FindHs(h5(i),h6(i),h7(i),h8(i),AANew, SigmaAA)

XS (SrfNum, i, j)=XS (SrfNum, i, j)

+(X3 (j) -hl (j) *Xl (I)

-h2(j)*X2(1)

-h3 (j) *PXlPBB (i)

-h4 (j)*PX2PBB(1))*h5(i)

+ (X4 (j)-hl (j) *XI (AL)

-h2 (j) *X2 (AL)

-h3 (j) *PXlPBB (AL)

-h4 (j) *PX2PBB (AL)) *h6 (i)

+(PX3PAA(j)-(hI(j)*PX3PAA(1)

+h2 (j) *PX3PAA(BL)

+h3 (j) *P2X00+h4 (j) *P2X01) )*h7 (i)

+ (PX4PAA(j)- ( hl (j) *PX4PAA (i)

+h2 (j) *PX4PAA (BL)

+h3(j)*P2Xl0+h4 (j)*P2Xll))*h8(i)

YS (SrfNum, i, j)=YS (SrfNum, i, j)

+(Y3(j)-hI(j)*YI(1)

-h2(j)*Y2(1)

-h3 (j) *PYIPBB (i)

-h4 (j)*PY2PBB(1))*h5(i)

+(Y4 (j)-hI(j)*YI(AL)

-h2 (j) *Y2 (AL)

-h3 (j) *PYIPBB (AL)

-h4(j)*PY2PBB(AL))*h6(i)

+(PY3PAA(j)-(hI(j)*PY3PAA(1)

+h2 (j) *PY3PAA (BL)

+h3(j)*P2Y00+h4(j)*P2Y01))*hT(i}

+ (PY4PAA(j)- ( hl (j) *PY4PAA(1)

+h2 (j) *PY4PAA (BL)

+h3(j)*P2Yl0+h4(j)*P2Yll))*h8 (i)

ZS (Sr f_,_u."_, i, j)=ZS (SrfNum, i, j)

86
30 CONTINUE

40 CONTINUE

C

C

C

RETURN

END

+(Z3 (j) -hl (j) *ZI (I)

-h2 (j) *Z2 (I)

-h3 (j) *PZIPBB (I)

-h4 (j)*PZ2PBB(1))*h5(i)

+ (Z4 (j) -hl (j) *ZI (AL)

-h2 (j) *Z2 (AL)

-h3 (j) *PZIPBB (AL)

-h4 (j) *PZ2PBB (AL)) *h6 (i)

+(PZ3PAA(j)-(hI(j)*PZ3PAA(1)

+h2 (j) *PZ3PAA (BL)

+h3(j)*P2Z00+h4(j)*P2Z01))*h7(i)

+ (PZ4PAA(j)- ( hl (j) *PZ4PAA (I)

+h2 (j) *PZ4PAA (BL)

+h3(j)*P2Zl0+h4(j)*P2Zll))*h8(i)

SUBROUTINE XiEtFl (XPnt, YPnt,ZPnt,II,JJ, KK, JI,J2,NKCsps,KCusp,

$ NoIts,MxGSiz)

C

C This SUBROUTINE smooths 3D, constant Zeta grid planes which have been

C disturbed by a constant Zeta cusp in a Xi-Zeta boundary surface. The

C process produces smoother grid lines in the Zeta direction.

C

INTEGER kc, i, j, k, i, II, JJ, KK, NICsps, NoIts,

$ KCusp(MxGSiz), Jl, J2

C

C

C

REAL XPnt (MxGSiz,MxGSiz,MxGSiz), YPnt (MxGSiz,MxGSiz,MxGSiz),

$ ZPnt (MxGSiz,MxGSiz,MxGSiz)

REAL XPnt(II,JJ, KK),YPnt (II,JJ, KK),ZPnt(II,JJ,KK)

DO

$

$

$

$

$

$

30 k=l, NKCsps

kc=KCusp (k)

DO 20 i=2,II-!

DO I0 j=Jl+l,J2-1

DO 5 l=l,No!ts

XPnt (i, j,kc)=0.5* (XPnt(i,

+XPnt (i,

YPnt (i, j,kc)=0.5*(YPnt(i,

+YPnt (i,

Zmnt (i, j,kc)=0.5*(ZPnt(i,

+ZPnt (i,

XPnt (i, j kc-I)=0.25 * (XPnt

+XPnt

YPnt (i, j kc-I)=0.25 * (YPnt

+YPnt

ZPnt (i, j kc-l)=0.25*(ZPnt

+ZPnt

j, kc+l) -2*XPnt (i, j, kc)

j, kc-l) )+ XPnt (i, j, kc)

j,kc+l)-2*YPnt (i, j,kc)

j,kc-l))+ YPnt (i, j,kc)

j,kc+l)-2*ZPnt (i, j,kc)

j,kc-l))+ ZPnt (i, j,kc)

(i, j,kc)-2* XPnt(i,j,kc-l)

(i, j,kc-2))+XPnt (i, j,kc-l)

(i, j,kc)-2* YPnt(i, j,kc-l)

(i, j,kc-2))+YPnt(i, j,kc-l)

(i, j,kc)-2* ZPnt (i, j, kc-l)

(i, j,kc-2))+ZPnt(i, j,kc-l)

XPnt (i, j

YPnt (i, j

ZPnt (i, j

kc+I)=0.25" (XPnt

+XPnt

kc+l =0.25* (YPnt

+YPnt

kc+l =0.25*(ZPnt

(i, j, kc+2)-2*XPnt (i, j, kc+l)

(i, j,kc))+ XPnt (i, j, kc+l)

(i, j,kc+2)-2*YPnt (i, j,kc+l)

(i, j,kc))+ YPnt(i, j,kc+l)

(i, j,kc+2)-2*ZPnt(i, j,kc+l)

87
C

C_- -

C

$

5 CONTINUE

i0 CONTINUE

20 CONTINUE

30 CONTINUE

+ZPnt(i, j,kc))+

RETURN

END

ZPnt (i, j, kc+l)

SUBROUTINE RdGrPIn (NGPt s, XB, YB, ZB, CrvNum, MxGSiz, InNum)

CrvNum

YB(MxGSiz,4), ZB(MxGSiz, 4)

the coordinates of the grid

INTEGER NGPts,

C

REAL XB(MxGSiz,4),

C

C This subroutine reads in

C on one edge.

C

DO I0 i=l,NGPts

READ(InNum,*)XB(i,CrvNum),YB(i,CrvNum),ZB(i,CrvNum)

i0 CONTINUE

C

C

points

RETURN

END

SUBROUTINE CalStI(NGPts,XB, YB, ZB, CrvNum, EdgNum,

$ StrB,MxBCvs,MxGSiz)

C

C This subroutine calculates the distribution function

C to the given distribution of grid points along the edge

C

INTEGER NGPts, CrvNum, EdgNum, i

C

REAL XB(MxGSiz,4), YB(MxGSiz,4), ZB(MxGSiz,4),

$ StrB (MxGSiz, MxBCvs)

C

C

I0

C

20

C

C

Cthat corresponds

'EdgNum'

StrB (I, EdgNum) =0.

DO i0 i=2,NGPts

StrB(i,EdgNum)=StrB(i-l,EdgNum) +

$ SQRT ((XB (i, CrvNum) -XB (i-l, CrvNum) )*'2

$ (YB (i, CrvNum) -YB (i-l, CrvNum) )**2

$ (ZB (i, CrvNum) -ZB (i-l, CrvNum) )*'2)

CONTINUE

SMax=StrB(NGPts,EdgNum)

DO 20 i=2,NGPts

StrB(i,EdgNum)=StrB(i,EdgNum)/SMax

CONTINUE

RETU_[

END

88
SUBROUTINECalSt2(EdgNum,NGPts,StrTp, Betal,Beta2,

$ StrB,MxBCvs,MxGSiz)

C

C This subroutine calculates the distribution function based

C stretching parameters 'StrTp' and 'Beta'

C

INTEGERNGPts, StrTp, EdgNum, i

C

C

i0

2o

C

C

__--.

REAL StrB(MxGSiz,MxBCvs), Betal, Beta2, A, B, DZ

_trB(l,EdgNum)=0.

IF(StrTp.LE.3)THEN

DO i0 i=I,NGPts-I

Alpha=(i-l.)/(NGPts-l.)

CALL FAiNew(AiNew,Alpha,Betal,StrTp)

StrB(i,EdgNum)=AiNew

CONTINUE

ELSEIF(StrTp. EQ.4)THEN

CALL Str4Prm(Betal,Beta2,A,B,DZ)

DO 20 i=2,NGPts-I

Alpha = (i-i.) / (NGPts-i.)

CALL Str4 (AINew, Alpha, A, B, DZ)

StrB (i, EdgNum)=AiNew

CONTINUE

ENDIF

StrB(NGPts,EdgNum)=l.

RETUPJ_

END

on the

C

C

C

C

C

C

C

C

C

i0

C

SUBROUTINE

$

$

EdgGPt s (CrvNum, EdgNum, NGPts, XB, YB, ZB, StrB,

x,y,z,s, zx, zy, zz,NDPts,Tensn,

MxBCv s, MxBPt s, MxGS iz )

This subroutine calculates the grid point location along an edge

based on a spline curve fitted through specified nodal points and

given distribution function.

INTEGER CrvNum, EdgNum, NGPts, NDPts(4), i, n

REAL XB(MxGSiz,4), YB(MxGSiz,4), ZB(MxGSiz,4),

$ StrB(MxGSiz,MxBCvs), x(4,MxBPts), y(4,MxBPts),

$ z(4,MxBPts), zx(4,MxBPts), zy(4,MxBPts),

$ zz(4,MxBPts), s(4,MxBPts), Tensn

SRa=S (CrvNum, NDPts (CrvNum))

DO I0 i=l,NGPts

SB=SRa*StrB (i, EdgNum)

CALL Splint (n,s,SB,NDPts,CrvNum, MxBPts)

XB (i, Crvh_um)=SpIVal (s, x, zx, SB, Tensn, n, CrvNum, MxBPt s)

Y5 (i, Crv._um)=Sp!Val (s, y, zy, SB, Tensn, n, CrvNum, MxBPt s)

ZB(i,CrvNum)=SplVal(s, z, zz,SB, Tensn, n,CrvNum, MxBPts)

CONT I_UE

89
RETURN

END

C

C

C

C

C

C

C

C

C

C

C

C

C _-

SUBROUTINE Str4Prm(S0,SI,A,B, DZ)

REAL SO, SI, A, B, DZ, Y, PI

This subroutine calculates the parameters

sided Vinokur stretching function.

A, B, and DZ for the

PI=ACOS (-i.)

A=SQRT (S0/SI)

B=SQRT (S0*SI)

IF (B. GT. 1 .001) THEN

IF (B. LE. 2 .7829681) THEN

Y=B-I

DZ=SQRT (6. *Y) * (i .-0.15*Y+0. 057321429* (Y**2)

$ -0. 024907295* (Y**3) +0. 0077424461* (Y**4)

$ -0.0010794123* (Y**5))

ELSEIF (B. GT. 2. 7829681) THEN

V=LOG (B)

W=I./B - 0.028527431

DZ=V+ (I. +I./V) *LOG (2 .*V) -0. 02041793+0. 24902722"W

$ +i. 9496443* (W*'2)-2. 6294547* (W**3) +8. 56795911* (W**4)

ENDIF

ELSEIF (B.LT. 0 .999) THEN

IF (B.LE. 0. 26938972) THEN

DZ=PI* (I .-B+B**2- (i .+ (PI**2)/6. )* (B**3) +6. 794732* (B**4)

$ -13. 205501* (B** 5) +II. 726095* (B**6))

ELSE

Y=B-I

DZ=SQRT (6. *Y) * (I. +0.15"Y+0. 057321429* (Y**2)

$ +0.048774238* (Y**3) -0.053337753* (Y**4)

$ +0. 075845134* (Y**5))

END!F

ENDIF

RETURI{

END

twoC

C

C

C

C

C

SUBROUTINE Str4(AiNew,Alpha,A,B,DZ)

REAL AiNew, Alpha, A, B, DZ, U, T

This subroutine calculates the value of the two-sided Vinokur

stretching function based on the value of the parameters A, B,

and DZ, and on the value of the "computational" coordinate Alpha.

IF (B. G7. !. 001) THE.<[

U=0 .5+TANH (DZ _:(Aipha-0.5) }/ (2 .*TANH (DZ/2 .))

gO
ELSEIF(B.LT.0. 999)THEN

U=0.5+TAN(DZ*(Alpha-0.5)) / (2. *TAN (DZ/2.) )

ELSE

U=Alpha * (i.+2.* (B-l) * (Alpha-0.5) * (1-Alpha))

ENDIF

T=U/ (A+ (I .-A) *U)

AINew=T

C

C

RETURN

END

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

i00

C

C

C

SUBROUTINE KFctrs(ZoneNo,kS,kl,k2,k3,k4,kXil,kXi2,kEtal,kEta2,

$ kZetal, kZeta2, MxSrfs,MxGSiz)

This subroutine is used to set the k-factors that are to be used.

The value of the k-factors is first set equal to the user specified

values of kXil, kXi2, kEtal, kEta2, kZetal and kZeta2. After that

the user can modify the k-factors for individual grid lines as

desired.

INTEGER ZoneNo

REAL kS(MxSrfs,MxGSiz,MxGSiz), kl(MxSrfs,MxGsiz),

$ k2(MxSrfs,MxGSiz), k3(MxSrfs,MxGSiz), k4(MxSrfs,MxGSiz),

$ kXil, kXi2, kEtal, kEta2, kZetal, kZeta2

Set the starting values of the k-factors:

first, k-factors used in generating interior grid points

DO I00 il=l,MxGSiz

DO I00 i2=l,MxGSiz

kS (i, il, i2) =kEtal

kS (2, il, i2) =kEta2

kS (3, il, i2) =kZetal

kS (4, il, i2) =kZeta2

CONTINUE

then, k-factors used to generate boundary surfaces.

DO 200 il =

kl(l,il

k2 (I, il

k3(l,il

k4 (I, il

kl (2 il

k2 (2 il

k3 (2 il

k4 (2 il

kl (3 i!

k2 (3 il

k3 (3 il

k4 (3 i!

kl (4 il

,MxGSiz

=kXil

=kXi2

=kZetal

=kZeta2

=kXi!

=kXi2

=kZeta!

=kZeta2

=kEta!

=kEta2

=kXi!

=kXi2

=kEtal

gl
2OO

C

C

C

C

C

C

C

C

C

k2 (4, il)=kEta2

k3 (4, il) =kXil

k4 (4, il)=kXi2

CONTINUE

Here, the user can make any desired modification of the Kfactors to improve the grid that he/she is generating. This part

of the subroutine will be case dependent.

RETURN

END

92
REFERENCES

i. Shih, T.I.-P.,et al.: GRID2D/3D - A Computer Program for Generating Grid

Systems in Complex-Shaped Two- and Three-Dimensional Spatial Domains:

Part 1 - Theory and Method. NASA TM-102453, 1990.

2. Bailey, R.T., et al.: GRID2D/3D - A Computer Program for Generating Grid

Systems in Complex-Shaped Two- and Three-Dimensional Spatial Domains:

Part 2 - User's Manual and Program Listing. NASA TM-102454, 1990.

3. Vinokur, M.: On One-Dimensional Stretching Functions for Finite-Difference

Calculations. J. Comput. Phys., vol. 50, May 1983, pp. 215-234. (Also, NASA

CR-3313.)

4. Thompson, J.F.; Warsi, Z.U.A.; and Mastin, C.W.: Numerical Grid Generation.

Elsevier Science Publishing Co., 1985.

93
TABLE 3.1 B Guide To Grid Control Parameters

Parameter

n

o;

i

K¢] K_2

qK-factors)

Type-i

StretchType i

Beml

Beml

Beta2

Range

2 (two boundary method)

4 (four boundary method)

0_<0_oo

O<K<**

1 (grid points)

2 (node points and stretching

function)

0 (no clustering)

1 (clustering near lower

boundary)

2 (clustering near upper

boundary)

3 (symmetric clustering

near both boundaries)

4 (asymmetric clustering

near both boundaries)

For StretchType =1,2 or 3:

1 < Betal < **

For StretchType = 4:

0< Beta <,_

Trial

value

Oto 10

10-3

1.1

2.5

Function

Con_'ols method used to

generatethegrid

Controls the shape

(curvature) of grid lines

Control orthogonality at

boundaries and curvature

of grid lines

Determines whether edge

is defined using grid

points or node points for

splining

Controls the distribution

of grid points along Edge i

Controls amount of

clustering near boundaries

(clustering increases as

Betal _ 1)

Controlamount of

clustering near boundaries

(clustering increases as

Beta _ **)

94
TABLE 3.2 -- Listing of Input File For Generation of Grid System

For Zone 18 of Radial Turbine Coolant Passage

0

0

2

2

5

5

1

4

i0 IL

15 JL

49 KL

I0.000000

20.000000

I0.000000

000000E+00

000000E+00

500000E-03

500000E-03

000000E-03

000000E-03

Type -

11645

11628

11605

11572

11532

11491

11467

11459

11438

11392

11321

11250

11204

11183

11181

11170

11148

11114

11078

11056

11047

11026

10980

10909

10839

10793

10772

10770

10759

10736

10702

10667

10645

10636

10615

10569

10498

10427

Technique

SigmaXi

SigmaEta

SigmaZeta

kXIl

kXI2

kETAI

kETA2

kZETAI

kZETA2

EDGE NO:

01677

01675

01672

01668

01660

01664

01690

01717

01715

01711

01704

01698

01694

01692

01662

01633

01608

01592

01598

01620

01648

01646

01642

01635

01627

01623

01621

01593

01567

01543

01529

01535

01555

01581

01579

01575

01568

01560

1

05842

05843

05844

05845

05843

05853

05878

05902

05904

05907

05911

.05915

05918

05920

05894

O5870

05851

05842

05853

05877

05902

05904

05907

05911

05915

05918

05920

05894

05869

05850

05842

05853

05876

05902

05904

05907

05911

05915

-- 1

---- 2

---- 3

---- 4

---- 5

---- 6

---- 7

---- 8

---- 9

--I0

--ii

--12

--13

--14

--15

--16

--17

--18

--19

--20

--21

--22

--23

--24

--25

--26

--27

--28

--29

--30

--31

--32

--33

--34

--35

--36

--37

--38

95
TABLE 3.2 (continued)

1

10381

10360

10359

10350

10331

10300

10268

10239

10214

10195

10181

Type

11596

11583

11564

11538

11507

11479

11464

11459

11438

11392

11321

11250

11204

11183

11181

11170

11148

11114

11078

11056

11047

11026

10980

10909

10839

10793

10772

10770

10758

10736

10702

10667

10645

10636

10615

10569

10498

10427

10381

10360

01556

01554

01531

01504

01479

01460

01450

01444

01442

01440

01439

- EDGE NO:

01985

01983

01981

01977

01975

01990

02012

02034

02032

02028

02021

02015

02010

O2008

01979

01950

01924

01909

01915

01938

01966

01963

01959

01952

01945

01940

01938

01911

01884

01860

01847

01854

01874

01900

01898

01894

01887

01881

01877

01874

.05918

.05920

.05898

.05875

.05854

.05842

.05838

.05838

.05839

.05841

.05842

2

--39

--40

--41

--42

--43

--44

--45

--46

--47

--48

--49

.05842

05842

05843

05843

05845

05862

05883

05902

05904

05907

05911

05915

05918

05920

05894

05870

05851

05842

05853

05877

05902

05904

05907

05911

05915

05918

05920

05894

05869

05850

05843

05853

05876

05902

05904

05907

05911

05915

05918

05920

-- 1

---- 2

---- 3

---- 4

---- 5

---- 6

---- 7

---- 8

---- 9

--I0

--Ii

--12

--13

--14

--15

-- 6

----' 7

--" 8

--' 9

--20

--21

--22

--23

--24

--25

--26

--27

--28

--29

--30

--31

--32

--33

--34

--35

--36

--37

--38

--39

--40
TABLE 3.2 (continued)

I0

2

20.00

2

.II

.II

3

I.I00

2

20.00

2

.I0

.i0

3

i.I00

1

Ii

II

ii

ii

ii

10358

10345

10315

10274

10234

10199

10170

10147

130

Type

01847

01817

01791

01775

01767

01764

01763

01762

.01762

- EDGE NO: 3

05894

05867

05847

05837

05835

05836

.05838

.05840

.05842

--41

--42

--43

--44

--45

--46

--47

--48

--49

Tension parameter

Number of nodes

645 .01677 .05842

596 .01985 .05842

StretchType

0 Stretching parameter

Type - EDGE NO: 4 .......

Tension parameter

Number of nodes

181 .01439 .05842

130 .01762 .05842

StretchType

BETA

0

Type

737

719

688

642

590

11545

11513

11495

11486

11458

11416

11381

11363

11359

11338

11303

11251

11192

11140

11105

11084

11074

11046

11004

10970

10952

10948

10927

Stretching

- EDGE NO:

01347

01345

01342

01337

01332

01327

01324

01322

01353

01375

01376

01354

01324

01292

01290

01287

01283

01278

01274

01271

01269

01297

01317

01318

01298

01270

01241

01239

parameter BETA

5

05530

05530

05531

05532

05533

05533

05534

05534

05565

05591

05598

05581

05554

05523

05525

05526

05529

05531

05534

05536

05536

05567

05592

05598

05581

05554

05523

05525

-- 1

---- 2

---- 3

---- 4

---- 5

---- 6

---- 7

---- 8

---- 9

--i0

--ii

--12

--13

--14

--15

--16

--17

--18

--19

--20

--21

--22

--23

--24

--25

--26

--27

--28

g7
TABLE 3.2 (continued)

10892

10840

10780

10728

10693

10672

10663

10639

10601

10567

10545

10536

10519

10494

10459

10418

10373

10332

10297

10272

10256

Type

11697

11682

11656

11618

11574

11536

11510

11495

11486

11458

11416

11381

11363

11359

11338

11303

11251

11192

11140

11105

11084

11074

11046

11004

10970

10952

10948

10927

10892

10840

01235

01231

01225

01221

01217

01215

01242

01262

01266

01250

01225

01197

01197

01197

01196

01195

01194

01193

01192

01191

01195

- EDGE NO:

01656

01654

01652

01648

01643

01639

01637

01635

01666

01689

01690

01668

01637

01605

01603

01600

01596

01591

01587

01584

01583

01610

01631

01632

01613

01585

01555

01554

01551

01547

6

05526

05529

05531

05534

05536

05536

05565

05589

05598

05588

05566

05539

05541

05545

05549

05555

05561

05567

05572

05576

05582

--29

--30

--31

--32

--33

--34

--35

--36

--37

--38

--39

--40

--41

--42

--43

--44

--45

--46

--47

--48

--49

05530

05530

05531

05532

05533

05533

05534

05534

05565

05591

05598

05581

05554

05523

05525

05526

05529

05531

05534

05536

05536

05567

05592

05598

05581

05554

05523

05525

05526

05529

-- 1

---- 2

---- 3

---- 4

---- 5

---- 6

---- 7

---- 8

---- 9

--I0

--II

--12

--13

--14

----' 5

-- 6

-- 7

---- 8

-- 9

--20

--21

--22

--23

--24

--25

--26

--27

--28

--29

--30

98
TABLE 3.2 (continued)

10780

10728

10693

10672

10663

10639

10602

10567

10545

10536

10516

10487

10448

10400

10348

10300

10261

10232

10212

Type

01543

01539

01536

01534

01561

01581

01586

01571

01546

01519

01519

01519

01518

01518

01518

01518

01518

01518

01522

- EDGE NO:

.05531

.05534

.05536

.05537

.05565

.05589

.05598

05588

05565

05538

05541

05544

05549

05555

05561

05567

05572

05576

05582

--31

--32

--33

--34

--35

--36

--37

--38

--39

--40

--41

--42

--43

--44

--45

--46

--47

--48

--49

2

20.00 Tension parameter

2 Number of nodes

.11737 .01347 .05530 -- 1

.11697 .01656 .05530 -- 2

3 StretchType

i. I000 Stretching parameter BETA

2 Type - EDGE NO: 8

20.00 Tension parameter

2 Number of nodes

.10256 .01195 .05582 -- 1

.10212 .01522 .05582 -- 2

3 StretchType

I.I000 Stretching parameter BETA

2 Type - EDGE NO: 9

20.00 Tension parameter

2 Number of nodes

.11645 .01677 .05842 -- 1

.11596 .01985 .05842 -- 2

3 StretchType

i.I000 Stretching parameter BETA

2 Type - EDGE NO: I0

20.00 Tension parameter

2 Number of nodes

.11737 .01347 .05530 -- 1

.11697 .01656 .05530 -- 2

3 StretchType

I.i000 Stretching parameter BETA

99
TABLE 3.2 (concluded)

2 Type - EDGE NO: ii

20.00 Tension parameter

4 Number of nodes

.11645 .01677 .05842

.11658 .01628 .05799

.11719 .01407 .05591

.11737 .01347 .05530

3 StretchType

i.I000 Stretching parameter BETA

2 Type - EDGE NO: 12

20.00 Tension parameter

4 Number of nodes

.11596 .01985 .05842

.11610 .01936 .05799

.11678 .01716 .05591

.11697 .01656 .05530

3 StretchType

i. I000 Stretching parameter BETA

2 Type - EDGE NO: 13

20.00 Tension parameter

2 Number of nodes

.10181 .01439 .05842

.10130 .01762 .05842

3 StretchType

I.I000 Stretching parameter BETA

2 Type - EDGE NO: 14 -

20.00 Tension parameter

2 Number of nodes

.10256 .01195 .05582

.10212 .01522 .05582

3 StretchType

i.I000 Stretching parameter BETA

2 Type - EDGE NO: 15

i00.00 Tension parameter

3 Number of nodes

.10181 .01439 .05842

.10253 .01226 .05613

.10256 .01195 .05582

3 StretchType

i.I000 Stretching parameter BETA

2 Type - EDGE NO: 16

I00.00 Tension parameter

3 Number of nodes

.10130 .01762 .05842

.10209 .01552 .05613

.10212 .01522 .05582

3 StretchType

I.I000

-- 1

---- 2

---- 3

---- 4

-- 1

---- 2

---- 3

---- 4

Stretching parameter BETA

-- 1

---- 2

---- 3

-- 1

---- 2

---- 3

i00
Blade

Coolant PassagewithinBlade

Directionof Coolant Flow

Pin Fins

Coolant How

_Axis of Rotation

(a)Cutout view

vider Walls

Coolant Flow

pin Fins

__Tip

Pb

(b) Schematic drawing

--Hub

Figure 1.1 -- Radial turbine coolant passage.

101
n (Technique:

IL

JL

KL

at

cn

o;

k_l

k_2

krll

k_2

k_l

k;2

Type-1

I Information

for

Edge 1

Type-2

Information

for

Edge 2

Type-3

Information

for

Edge 3

n=2 orn =4 )

T_e-m

Information

for

Edge m

StretchType 11

Beta1 Beta 2

StretchType 12

Betal Beta 2

StretchType 15

Betal Beta 2

StretchType 16

Betal Beta 2

m=8 if n=2 (two boundary technique)

m=16 if n--4 (four boundary technique)

Figure 3.1 -- Input-file format for GRID3D-v2.

102
Information for Edge i

if Type-i = 1:

xl Yl zl

X2 Y2 Z2

xrcL y_L z_,a_

NL=IL for i=3,4,9,10,13 and 14

NL=JL for i=11,12,15 and 16

NL=KL for i=1,2,3,4

if Type-i = 2:

NP

Xl Yl Zl

x2 Y2 z2

xha, YNr' zn_,

StretchType i

Betal Beta2

c_ = tension for spline

NP = number of node points

Figure 3.1 (concluded)

103
x

Z

surface 2

(top)

I .I sur face 3

I (back)

l

surface 4

(front )

= $

surface 1

(bottom)

Figure 3.2 -- Edge curve and boundary surface numbering scheme for GRID3D-v2.

104
Zone 4

Zone 3

Zone 2

Zone I

Zone 5

Zone 6

Zone 7

Zone 8

Zone 9

Zone 15

Zone 16

Zone 17

Zone 18

Zone 19

Zone 14

Zone 10

Zone 11

Zone 12

Zone 13

Figure 3.3 --Partitioning of the spatial domain of radial turbine coolant passage into zones

for grid generation.

105
Figure 3.4 m Grid system for zone 18 of radial turbine coolant passage.

106
Figure 3.5 -- Grid system for the whole radial turbine coolant passage (2-D view).

107
Figure 3.6 -- Grid system for the whole radial turbine coolant passage (3-D view).

108
• Report Documentation Page

AdminiltrlIiOn

I. Report No. 2. Government Accession No. 3. Recil:fienrs CI_a_og No.

NASA 'I'M- 103766

4. Title and Subtitle

GRID3D-v2: An Updated Version of the GRID2D/3D Computer

Program for Generating Grid Systems in Complex-Shaped ThreeDimensional Spatial Domains

7. Author(s)

E. Steinthorsson, T.I-P. Shih, and R.J. Roelke

g. PIrforming Organization Name and Address

National Aeronautics and Space Administration

Lewis R_ch Center

Cleveland, Ohio 44135- 3191

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546-0001

5. Report Date

June 1991

6. Pecformin9 Organizalion Code

8. Performing Organization Report No.

E-6028

10. Work Unit No.

535-05-10

I1. _ o_ Grant No.

t3. Type of Report and PerkxI Coveced

Technical Memorandum

14. Sponsonng Agency Code

15. Supplementary Notes

E. Steinthorsson and T. I-P. Shih, Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh,

Pennsylvania 15213; R.J. Roelke, NASA Lewis Research Center. Responsible person, R.J. Roelke, (216) 433-3403.

16. Abstract

In order to generate good quality grid systems for complicated three-dimensional spatial domains, the grid-generation

method used must be able to exert rather precise controls over grid-point distributions. In this report, several techniques

are presented that enhance control of grid-point distribution for a class of algebraic grid-generation methods known as

the two-, four-, and six-boundary methods. These techniques include variable stretching functions from bilinear interpolation, interpolating functions based on tension splines, and normalized "K-factors." The techniques developed in this

study were incorporated into a new versio]l of GRID3D called GRID3D-v2. The usefulness of GRID3D-v2 was demonstrated by using it to generate a three-dimensional grid system in the coolant passage of a radial turbine blade with

serpentine channels and pin fins.

17. KeyWords (Suggested byAuthor(s))

Grid generation

! 19. Secunty Classif. (of the report)

Unclassified

18. Distribution Statement

Unclassified - Unlimited

Subject Category, 61

20. Security Classif. (of this page)

Unclassified

21. No. of pages 22. Price*

112 A06

NASA FORM 1626 OCT 86

"For sale by the National Technical Information Service, Springfield, Virginia 22161