# Thompson, Warsi, Mastin - Numerical Grid Generation

Посмотреть архив целикомПросмотреть файл в отдельном окне: 96802bdecc704dc69482af0cab744ca1.pdf

PREFACE

Numerical grid generation has now become a fairly common tool for use in the

numerical solution of partial differential equations on arbitrarily shaped regions. This is

especially true in computational fluid dynamics, from whence has come much of the impetus

for the development of this technique, but the procedures are equally applicable to all

physical problems that involve field solutions. Numerically generated grids have provided

the key to removing the problem of boundary shape from finite difference methods, and

these grids also can serve for the construction of finite element meshes. With such grids all

numerical algorithms, finite difference or finite element, are implemented on a square grid in

a rectangular computational region regardless of the shape and configuration of the physical

region. (Finite volume methods are effectively a type of conservative finite difference

method on these grids.)

In this text, grid generation and the use thereof in numerical solutions of partial

equations are both discussed. The intent was to provide the necessary basic information,

from both the standpoint of mathematical background and from that of coding

implementation, for numerical solutions of partial differential equations to be constructed on

general regions. Since these numerical solutions are ultimately constructed on a square grid

in a rectangular computational region, any solution algorithm that can treat equations with

variable coefficients is basically applicable, and therefore discussion of specific algorithms is

left to classical texts on the numerical solution of partial differential equations.

The area of numerical grid generation is relatively young in practice, although its roots

in mathematics are old. This somewhat eclectic area involves the engineer’s feel for physical

behavior, the mathematician’s understanding of functional behavior, and a lot of

imagination, with perhaps a little help from Urania. The physics of the problem at hand must

ultimately direct the grid points to congregate so that a functional relationship on these

points can represent the physical solution with sufficient accuracy. The mathematics controls

the points by sensing the gradients in the evolving physical solution, evaluating the accuracy

of the discrete representation of that solution, communicating the needs of the physics to the

points, and by providing mutual communication among the points as they respond to the

physics.

Numerical grid generation can be thought of as a procedure for the orderly distribution

of observers, or sampling stations, over a physical field in such a way that efficient

communication among the observers is possible and that all physical phenomena on the

entire continuous field may be represented with sufficient accuracy by this finite collection

of observations. The structure of an intersecting net of families of coordinate lines allows the

observers to be readily identified in relation to each other, and results in much more simple

coding than would the use of a triangular structure or a random distribution of points. The

grid generation system provides some influence of each observer on the others, so that if one

moves to get into a better position for observation of the solution, its neighbors will follow to

some extent in order to maintain smooth coverage of the field.

Another way to think of the grid is as the structure on which the numerical solution is

built. As the design of the lightest structure requires consideration of the load distribution, so

the most economical distribution of grid points requires that the grid be influenced by both

the geometric configuration and by the physical solution being done thereon. In any case,

since resources are limited in any numerical solution, it is the function of the numerical grid

generation to make the best use of the number of points that are available, and thus to make

the grid points an active part of the numerical solution.

This is a rapidly developing area, being now only about ten years old, and thus is still

in search of new ideas. Therefore no book on the subject at this time could possibly be

considered to be definitive. However, enough material has now accumulated in the literature,

and enough basic concepts have emerged, that a fundamental text is now needed to meet the

needs of the rapidly expanding circle of interest in the area. It is with the knowledge of both

these needs and these limitations that this text has been written. Some of the techniques

discussed will undoubtedly be superceded by better ideas, but the fundamental concepts

should serve for understanding, and hopefully also for some inspiration, of new directions.

The only background assumed of the student is a senior-level understanding of numerical

analysis and partial differential equations. Concepts from differential geometry and tensor

analysis are introduced and explained as needed.

Numerical grid generation draws on various areas of mathematics, and emphasis

throughout is placed on the development of the relations involved, as well as on the

techniques of application. This text is intended to provide the student with the understanding

of both the mathematical background and the application techniques necessary to generate

grids and to develop codes based on numerically generated grids for the numerical solution

of partial differential equations on regions of arbitrary shape.

The writing of this text has been a cooperative effort over the last two years, spurred

on by the institution of a graduate course in numerical grid generation, as well as an annual

short course, at Mississippi State. The students in both of these courses have contributed

significantly in revising the text as it evolved. The last appendix is the result of a class

assignment prepared by Col. Hyun Jin Kim, graduate student in the computational fluid

dynamics program, who also compiled the index. Our colleage, Dr. Helen V.

McConnaughey of Mathematics contributed significantly through continual discussions and

wrote most of Chapter IV.

We are indebted to a large number of former students and fellow researchers around

the world for the development of the ideas that have crystallized into numerical grid

generation. The complete debt can be acknowledged only through mention of the

bibliographies contained in the several surveys cited herein. A list here would either be too

long to note the strongest influences or too short to acknowledge all the significant ones. We

must, however, acknowledge the many long and fruitful discussions with Peter Eiseman of

Columbia University.

0f vital importance is the support that has been provided for the research from which

the developments discussed in this book have emerged, including NASA; the research

offices of the Air Force, Army, and Navy; the National Science Foundation, and various

industrial concerns. The interest and contributions of a number of contract monitors has been

essential over the years. We are especially appreciative of Bud Bobbitt and Jerry South of

NASA Langley Research Center, who provided the initial support for an unknown with an

idea.

Particular debts are owed to W. H. Chu for an idea in the Journal of Computational

Physics in 1971, and to Frank Thames who put the idea into a dissertation.

In the preparation of the text we had the conscientious and untiring efforts of two most

able secretaries, Rita Curry and Susan Triplett, who typed on in good spirits through a year

of numerous revisions and frustrations as the text evolved.

Finally, we were particularly fortunate to have the services of Yeon Seok Chae,

graduate student in the computational fluid dynamics program and illustrator par excellence,

who did all the figures with understanding of the intended meaning as well as artistic

competence. His meticulous efforts were extensions of our thoughts.

Joe F. Thompson

Z. U. A. Warsi

C. Wayne Mastin

Mississippi State, Mississippi

January 1985

I. INTRODUCTION

The numerical solution of partial differential equations requires some discretization of

the field into a collection of points or elemental volumes (cells). The differential equations

are approximated by a set of algebraic equations on this collection, and this system of

algebraic equations is then solved to produce a set of discrete values which approximates the

solution of the partial differential system over the field. The discretization of the field

requires some organization for the solution thereon to be efficient, i.e., it must be possible to

readily identify the points or cells neighboring the computation site. Furthermore, the

discretization must conform to the boundaries of the region in such a way that boundary

conditions can be accurately represented. This organization is provided by a coordinate

system, and the need for alignment with the boundary is reflected in the routine choice of

cartesian coordinates for rectangular regions, cylindrical coordinates for circular regions,

etc., to the extent of the handbook’s resources.

The current interest in numerically-generated, boundary-conforming coordinate

systems arises from this need for organization of the discretization of the field for general

regions, i.e., to provide computationally for arbitrary regions what is available in the

handbook for simple regions. The curvilinear coordinate system covers the field and has

coordinate lines (surfaces) coincident with all boundaries. The distribution of lines should be

smooth, with concentration in regions of strong solution variation, and the system should

ultimately be capable of sensing these variations and dynamically adjusting itself to resolve

them.

A numerically-generated grid is understood here to be the organized set of points

formed by the intersections of the lines of a boundary-conforming curvilinear coordinate

system. The cardinal feature of such a system is that some coordinate line (surface in 3D) is

coincident with each segment of the boundary of the physical region. The use of coordinate

line intersections to define the grid points provides an organizational structure which allows

all computation to be done on a fixed square grid when the partial differential equations of

interest have been transformed so that the curvilinear coordinates replace the cartesian

coordinates as the independent variables.

This grid frees the computational simulation from restriction to certain boundary

shapes and allows general codes to be written in which the boundary shape is specified

simply by input. The boundaries may also be in motion, either as specified externally or in

response to the developing physical solution. Similarly, the coordinate system may adjust to

follow variations developing in the evolving physical solution. In any case, the

numerically-generated grid allows all computation to be done on a fixed square grid in the

computational field which is always rectangular by construction.

In the sections which follow, various configurations for the curvilinear coordinate

system are discussed in Chapter II. In general, the computational field will be made

rectangular, or composed of rectangular sub-regions, and a wide variety of configurations is

possible. Coordinate systems may also be generated separately for sub-regions in the

physical plane and patched together to form a complete system for complex configurations.

The basic transformation relations applicable to the use of general curvilinear coordinate

systems are developed in Chapter III; the construction of numerical solutions of partial

differential equations on those systems is discussed in Chapter IV; and consideration is given

in Chapter V to the evaluation and control of truncation error in the numerical

representations.

Basically, the procedures for the generation of curvilinear coordinate systems are of

two general types: (1) numerical solution of partial differential equations and (2)

construction by algebraic interpolation. In the former, the partial differential system may be

elliptic (Chapter VI), parabolic or hyperbolic (Chapter VII). Included in the elliptic systems

are both the conformal (Chapter X), and the quasi-conformal mappings, the former being

orthogonal. Orthogonal systems (Chapter IX) do not have to be conformal, and may be

generated from hyperbolic systems as well as from elliptic systems. Some procedures

designed to produce coordinates that are nearly orthogonal are also discussed. The algebraic

procedures, discussed in Chapter VIII, include simple normalization of boundary curves,

transfinite interpolation from boundary surfaces, the use of intermediate interpolating

surfaces, and various other related techniques.

Coordinate systems that are orthogonal, or at least nearly orthogonal near the

boundary, make the application of boundary conditions more straightforward. Although

strict orthogonality is not necessary, and conditions involving normal derivatives can

certainly be represented by difference expressions that combine one-sided differences along

the line emerging from the boundary with central expressions along the boundary, the

accuracy deteriorates if the departure from orthogonality is too large. It may also be more

desirable in some cases not to involve adjacent boundary points strongly in the

representation, e.g., on extrapolation boundaries. The implementation of algebraic turbulence

models is more reliable with near-orthogonality at the boundary, since information on local

boundary normals is usually required in such models. The formulation of boundary-layer

equations is also much more straightforward and unambiguous in such systems. Similarly,

algorithms based on the parabolic Navier-Stokes equations require that coordinate lines

approximate the flow streamlines, and the lines normal thereto, especially near solid

boundaries. It is thus better in general, other considerations being equal, for coordinate lines

to be nearly normal to boundaries.

Finally, dynamically-adaptive grids are discussed in Chapter XI. These grids

continually adapt during the course of the solution in order to follow developing gradients in

the physical solution. This topic is at the frontier of numerical grid generation and may well

prove to be one of its most important aspects.

The emphasis throughout is on grids formed by the intersections of coordinate lines of

a curvilinear coordinate system, as opposed to the covering of a field with triangular

elements or a random distribution of points. Neither of these latter collections of points is

suitable for really efficient numerical solutions (although numerical representations can be

constructed on each, of course) because of the cumbersome process of identification of

neighbors of a point and the lack of banded structure in the matrices. Thus the subject of

triangular mesh generators, per se, is not addressed here. (Obviously a triangular mesh can

be produced by construction rectangular mesh diagonals.)

Considerable progress is being made toward the development of the techniques of

numerical grid generation and toward casting them in forms that can be readily applied. A

comprehensive survey of numerical grid generation procedures and applications thereof

through 1981 was given by Thompson, Warsi, and Mastin in Ref. [1], and the conference

proceedings published as Ref. [2] contains a number of expository papers on the area, as

well as current results. Other collections of papers on the area have also appeared (Ref. [3]

and [4]), and a later review through 1983 has been given by Thompson in Ref. [5]. Some

other earlier surveys are noted in Ref. [1]. A later survey by Eiseman is given in Ref. [37].

The present text is meant to be a developmental treatment of the techniques of grid

generation and its applications, not a survey of results, and therefore no attempt is made here

to cite all related references, rather only those needed to illustrate particular points are noted.

The surveys mentioned above should be consulted directly for references to examples of

various applications and related contributions. (Ref [l] gives a short historical development

of the ideas of grid generation.) Other surveys of particular areas of grid generation are cited

later as topics are introduced.

Finally, in regard to implementation, a configuration for the transformed

(computational) field is first established as discussed in Chapter II. The grid is generated

from a generation system constructed as discussed in Chapters VI -- X. (If the grid is to be

adaptive, i.e., coupled with the physical solution done thereon, then the gr1d must be

continually updated as discussed in Chapter XI.) In the construction of the grid, due account

must be taken of the truncation error induced by the grid discussed in Chapter V. The partial

differential equations of the physical problem of interest are transformed according to the

relations given in Chapter III. These transformed equations are then discretized, cf. Chapter

IV, and the resulting set of algebraic equations is solved on the fixed square grid in the

rectangular transformed field.

II. BOUNDARY-CONFORMING COORDINATE SYSTEMS

1. Basic Concepts

To provide a familiar ground from which to view the general development to follow,

consider first a two-dimensional cylindrical coordinate system covering the annular region

between two concentric circles:

Here the curvilinear coordinates (r, ) vary on the intervals [r1,r2] and [0,2 ], respectively.

These curvilinear coordinates are related to the cartesian coordinates (x,y) by the

transformation equations

(1)

The inverse transformation is given by

(2)

Note that one of the curvilinear coordinates, r, is constant on each of the phys1cal

boundaries, while the other coordinate, , varies monotonically over the same range around

each of the boundaries. Note also that the system can be represented as a rectangle on which

the two physical boundaries correspond to the top and bottom sides:

The transformed region, i.e., where the curvilinear coordinates, r and the independent

variables, thus can be thought of as being rectangular, and can be treated as such from a

coding standpoint. These points will be central to what follows.

The curvilinear coordinates (r, ) can be normalized to the interval [0,1] by introducing

the new curvilinear coordinates ( , ), where

(3)

or

(4)

The transformation then may be written

(5a)

(5b)

where now and both vary on the interval [0,1]. This is thus a mapping of the annular

region between the two circles in the physical space onto the unit square in the transformed

space, i.e., each point (x,y) on the annulus corresponds to one, and only one, point ( , ) on

the unit square:

The bottom ( = 0) and top ( = 1) of the square correspond, respectively, to the inner and

outer circles, r = r1, and r = r2. The sides of the square, = 0 and = 1 correspond to = 0

and = 2 , respectively, and hence to the two coincident sides of a branch cut in the

physical space. Therefore, boundary conditions are not to be specified on these sides of the

unit square in the transformed space. Rather these sides are to be considered re-entrant on

each other with points adjacent to one, outside the square, being equivalent to points adjacent

to the other, inside the square.

Conceptually, the physical region can be considered to have been opened at the cut

= 0 and 2 and then deformed into a rectangle to form the transformed region:

Here, point correspondence across the re-entrant boundaries (indicated by the dashed

connecting line) in the transformed region is illustrated by the coincidence of the pair of

circled points. This conceptual device and mode of illustration for the the point

correspondence across re-entrant boundaries will serve later for more general configurations.

These simple concepts extend to more complicated two-dimensional configurations,

the central feature being that one of the curvilinear coordinates is made to be constant on a

boundary curve (as was r above), while the other varies monotonically along that boundary

curve (as does ). The transformation to the rectangle is achieved by making the range and

direction of variation of the varying coordinate the same on each of two opposing boundaries

(as varies from 0 to 2 on each circle above).

The physical space thus transforms to the rectangle shown above regardless of the shape of

the physical region. (It is not necessary to normalize the curvilinear coordinates to the

interval [0,1], and in fact, any normalization can be used. In computational applications the

normalization is more conveniently done to different intervals for each coordinate. The field

in the transformed space is then rectangular, rather than square.) Familiar examples of this

are elliptical coordinates for the region between two confocal ellipses, spherical coordinates

for two spheres, parabolic coordinates for two parabolas, etc.

These, same concepts will be extended later to completely general configurations

involving any number of boundary curves and branch outs. The extension to three

dimensions follows directly, using boundary surfaces instead of curves, i.e., one curvilinear

coordinate will be made constant on a boundary surface, with the other two forming a

two-dimensional coordinate system on the surface.

Returning to the concentric circles, if the functional dependence of on , and/or that of r

on , had been made more general than the simple linear normalizations given by Eq. (4),

the corresponding coordinate lines would have become unequally spaced in the physical

space, while remaining as radial lines and concentric circles:

The transformation, from Eq. (1), is now given by

(6a)

(6b)

In this case the points on the inner and outer circular boundaries are not equally spaced

around the circles in the physical space for equal increments of , although they remain

equally spaced on the top and bottom of the unit square in the transformed space by

construction. The spacing around these circles is determined by the functional dependence of

on , and, since the points are located at equal increments of by construction, this

functional relationship is defined by the placement of these points around the circles. This

point, that the coordinate system in the field is determined from the boundary point

distribution, will be central to the discussion of grid generation to follow. The distribution of

circumferential lines is controlled here by the functional relationship between r and , which

is not related to any boundary point distribution. Thus factors other than the boundary point

distribution may be expected to be involved in grid generation, as well. That the point

distribution on the boundaries may be controlled by direct placement of the points, while the

coordinate line distribution in the field must be controlled by other means will also continue

to appear in the developments to follow.

The one-dimensional functional relationship between and in Eq. (6) requires that

the relative distributions of boundary points around the inner and outer circles be the same.

This restriction can be removed by making a function of , as well as of , while

retaining the periodic nature of the dependence on . In this case the ooordinate lines of

constant will no longer be straight radial lines, although they will continue to connect

corresponding points on the inner and outer circular boundaries. Similarly the

circumferential coordinate lines (lines of constant here) can be made to depart from circles

by making r dependent on both and , but with the restriction that the dependence

vanishes on the inner and outer circular boundaries (where

here).

= 0 and

= 1, respectively,

Obviously certain constraints will have to be placed on the functions ( , ) and r( , ) to

keep the mapping one-to-one. All of these considerations will reappear in the general

developments that follow.

Finally, it should be realized that the intermediate use here of the cylindrical

coordinates (r, ) in defining the transformation between the curvilinear coordinates ( , )

and the cartesian coordinates (x,y) has been only in deference to the familiarity of the

cylindrical coordinates, and such intermediary coordinates will not appear in general. The

generalized statement for the simple configuration under consideration here is as follows:

Find (x,y) and (x,y) in the annular region bounded by the curves x2 + y2 = and x2 + y2

= , subject to the boundary conditions

Specified monotonic variation of over [0,1] on x2 + y2 =

with same sense of direction on each of these two curves.

and on x2 + y2 =

It is the inverse problem that will be treated in fact, however, i.e., find x( , ) and

y( , ) on the unit square in the transformed space (0

1, 0

1), subject to the

boundary conditions

x( ,0) and y( ,0) specified on = 0 such that x2( = 0) + y2 ,0) =

x( ,1) and y( ,1) specified on = 1 such that x2( = 1) + y2( ,1) =

Periodicity in : x(1 + , ) = x( , ) y(1 + , ) = y( , )

The simple form for the transformation given by Eq. (6) is made possible by choosing the

same functional dependence of x and y on on the boundaries, = 0 and = 1. The

familiar cylindrical coordinate system is thus a special case of the general grid generation

problem for this simple configuration applicable to the region between two concentric

circles, as is the elliptical coordinate system for two ellipses, etc.

2. Generalization

Generalizing from the above consideration of cylindrical coordinates, the basic idea of

a boundary-conforming curvilinear coordinate system is to have some coordinate line (in 2D,

surface in 3D) coincident with each boundary segment, analogous to the way in which lines

of constant radial coordinate coincide with circles in the cylindrical coordinate system. The

other curvilinear coordinate, analogous to the angular coordinate in the cylindrical system,

will vary along the boundary segment and clearly must do so monotonically, else the same

pair of values of the curvilinear coordinates will occur at two different physical points. (It

should be clear that the curvilinear coordinate that varies along a boundary segment must

have the same direction and range of variation over some opposing segment, e.g., as the

angular variable varies from 0 to 2 over both of two concentric circles in cylindrical

coordinates).

With the values of the curvilinear coordinates thus specified on the boundary, it then

remains to generate values of these coordinates in the field from these boundary values.

There must, or course, be a unique correspondence between the cartesian (or other basis

system) and the curvilinear coordinates, i.e., the mapping of the physical region onto the

transformed region must be one-to-one, so that every point in the physical field corresponds

to one, and only one, point in the transformed field, and vice versa. Coordinate lines of the

same family must not cross, and lines of different families must not cross more than once.

In this chapter a two-dimensional region will be considered in most of the discussions

in the interest of economy of presentation. Generalization to three dimensions will be evident

in most cases and will be mentioned specifically only when necessary. As noted above, the

curvilinear coordinates may be normalized to any intervals, just as the radial and angular

coordinates of the cylindrical coordinate system can be expressed in many different units.

Since the interest of the present discussion is numerical application, it will be generally

convenient to define the increments of all the curvilinear coordinates to be uniformly unity,

and then to normalize these coordinates to the interval [1,N(i)], where N(i) is the total number

of grid points to be used in the i direction. (The three curvilinear coordinates will be

indicated as i, i = 1,2,3, in general. In two dimensions, however, the notation ( , ) will

often be used for the two coordinates 1 and 2 .) The computational field, i.e., the field in

the transformed space, thus will have rectangular boundaries and will be covered by a square

grid. (It will become clear later that the actual values of the increments in the curvilinear

coordinates are immaterial since they do not appear in the final numer1cal expressions.

Therefore no generality is lost in making the grid square and of unit increment in the

transformed field.)

A. Boundary-value Problem -- Physical Region

The generation of the curvilinear coordinate system may be treated as follows: with

the curvilinear coordinates specified on the boundaries, e.g., (x,y) and (x,y) on a

boundary curve (this specification amounting to a constant value for either or on each

segment of , with a specified monotonic variation of the other over the segment), generate

the values, (x,y) and (x,y), in the field bounded by . This is thus a boundary value

problem on the physical field with the curvilinear coordinates ( , ) as the dependent

variables and the cartesian coordinates (x,y) as the independent variables, with boundary

conditions specified on curved boundaries:

(In these discussions, the transformation is assumed to be from cartesian coordinates in the

physical space. The transformation can, however, be from any system of coordinates in the

physical space.)

B. Boundary value Problem - Transformed Region

The problem may be simplified for computation, however, by first transforming so

that the physical cartesian coordinates (x,y) become the dependent variables, with the

curvilinear coordinates ( , ) as the independent variables. Since a constant value of one

curvilinear coordinate, with monotonic variation of the other, has been specified on each

boundary segment, it follows that these boundary segments in the physical field will

correspond to vertical or horizontal lines In the transformed field. Also, since the range of

variation of the curvilinear coordinate varying along a boundary segment has been made the

same over opposing segments, it follows that the transformed field will be composed of

rectangular blocks.

The boundary value problem in the transformed field then involves generating the

values of the physical cartesian coordinates, x( , ) and y( , ), in the transformed field

from the specified boundary values of x( , ) and y( , ) on the rectangular boundary of the

transformed field, the boundary being formed of segments of constant or , i.e., vertical or

horizontal lines. With = constant on a boundary segment, and the increments in taken to

be uniformly unity as discussed above, this boundary value specification is implemented

numerically by distributing the points as desired along the boundary segment and then

assigning the values of the cartesian coordinates of each successive point as boundary values

at the equally spaced boundary points on the bottom (or top) of the transformed field in the

following figure.

Boundary values are not specified on the left and right sides of the transformed field since

these boundaries are re-entrant on each other (analogous to the 0 and 2 lines in the

cylindrical system), as discussed above, and as indicated by the connecting dotted line on the

figure. Points outside one of these re-entrant boundaries are coincident with points at the

same distance inside the other. The problem is thus much more simple in the transformed

field, since the boundaries there are all rectangular, and the computation in the transformed

field thus is on a square grid regardless of the shape of the physical boundaries.

With values of the cartesian coordinates known in the field as functions of the

curvilinear coordinates, the network of intersecting lines formed by contours (surfaces in

3D) on which a curvilinear coordinate is constant, i.e., the curvilinear coordinate system,

provides the needed organization of the discretization with conformation to the physical

boundary. It is also possible to specify intersection angles for the coordinate lines at the

boundaries as well as the point locations.

3. Transformed Region Configurations

As noted above, the generation of the curvilinear coordinate system is done by

devising a scheme for determination of the field values of the cartesian coordinates from

specified values of these coordinates (and/or curvilinear coordinate line intersection angles)

on portions of the boundary of the transformed region. Since the boundary of the

transformed region is comprised of horizontal and vertical line segments, portions of which

correspond to segments of the physical boundary on which a curvilinear coordinate is

specified to be constant, it should be evident that the configuration of the resulting

coordinate system depends on how the boundary correspondence is made, i.e., how the

transformed region is configured.

Some examples of different configurations are given below, from which more

complex configurations can be inferred. In these examples only a minimum number of

coordinate lines are shown in the interest of clarity of presentation tation. In all of these

examples, boundary values of the physical cartesian coordinates (and/or curvilinear

coordinate line intersection angles) are understood to be specified on all boundaries, both

external and internal, of the transformed region except for segments indicated by dotted

lines. These latter segments correspond to branch cuts in the physical space, as is explained

in the examples in which they appear. Such re-entrant boundary segments always occur in

pairs, the members of which are indicated by the dashed connecting lines on each of the

configurations shown. Points outside the field across one segmentof such a pair are

coincident with points inside the field across the other member of the pair. The conceptual

device of opening the physical field at the cuts is used here to help clarify the

correspondence between the physical and transformed fields. In many cases an example of

an actual coordinate system is given as well. References to the use of various configurations

may be found in the surveys given by Ref. [1] and [5], and a number of examples appear in

Ref. [2].

A. Simply-connected Regions

It is natural to define the same curvilinear coordinate to be constant on each member

of a pair of generally opposing boundary segments in the physical plane. Thus, a

simply-connected region formed by four curves is logically treated by transforming to an

empty rectangle:

Similarly, an L-shaped region could remain L-shaped in the transformed region:

Here, for instance, the cartesian coordinates of the desired points on the physical boundary

segment 4-5 are specified as boundary conditions on the vertical line 4-5, in corresponding

order, which forms a portion of the boundary of the transformed region.

The generalization of these ideas to more complicated regions is obvious, the

transformed region being composed of contiguous rectangular blocks. An example follows:

The physical boundary segment on which a single curvilinear coordinate is constant

can have slope discontinuities, however, so that the L-shaped region above could have been

considered to be composed of four segments instead of six, so that the transformed region

becomes a simple rectangle:

Here the cartesian coordinates of the desired points on the physical boundary 5-4-3 are the

specified boundary values from left to right across the top of the transformed region.

Whether or not the boundary slope discontinuity propagates into the field, so that the

coordinate lines in the field exhibit a slope discontinuity as well, depends on how the

coordinate system in the field is generated, as will be discussed later.

It is not necessary that corners on the boundary of the transformed region correspond

to boundary slope discontinuities on the physical boundary and a counter-example follows

next:

In this case, the segment 1-2 on the physical boundary is a line of constant , while the

segment 1-4 is a line of constant . Thus at point 1 we have the following coordinate line

configuration:

The lines through point 1 are as follows:

so that the angle between the two coordinate lines is at point 1, and consequently the

Jacobian of the transformation (the cell area, cf. Chapter III) will vanish at this point. The

coordinate species thus changes on the physical boundary at point 1. (Difference

representations at such special points as this, and others to appear in the following examples,

are discussed in Chapter IV.) Since the species of curvilinear coordinate necessarily changes

at a corner on the transformed region boundary, the identification of a concave corner on the

transformed region boundary with a point on a smooth physical boundary will always result

in a special point of the type illustrated here. (A point of slope discontinuity on the physical

boundary also requires special treatment in difference solutions, since no normal can be

defined thereon. This, however, is inherent in the nature of the physical boundary and is not

related to the construction of the transformed configuration.)

Some slightly more complicated examples of the alternatives introduced above now

follow:

Still another alternative in this case would be to collapse the intrusion 2-3-4-5 to a slit in the

transformed region:

Here the physical cartesian coordinates are specified and are double-valued on the vertical

slit, 2-9-5, in the transformed region. The cartesian coordinates of the desired points on the

physical boundary 2-9 are to be used on the slit in the generation of the grid to the left of the

slit in the transformed region, while those on the physical boundary 5-9 are used for

generation to the right of the slit. Solution values in a numerical solution on such a

coordinate system would also be double-valued on the slit, of course. This

double-valuedness requires extra bookkeeping in the code, since two values of each of the

cartesian coordinates and of the physical solution must be available at the same point in the

transformed region so that difference representations to the left of the slit use the slit values

appropriate to the left side, etc. Difference representations near slits are discussed in Chapter

IV. With the composite grid structure discussed in Section 4, however, this need for

double-valuedness, and the concomitant coding complexity, with the slit configuration can

be avoided.

The point 9 here requires special treatment, since the coordinate line configuration

there is as follows:

The coordinate lines through point 9 are as follows:

Here the slope of the coordinate line on which varies is discontinuous at point 9, and the

line on which varies splits at this point. Such a special point will always occur at the slit

ends with the slit configuration.

B. Multiply-connected Regions

With obstacles in the interior of the field, i.e., with interior boundaries, there are still

more alternative configurations of the transformed region. One possibility is to maintain the

connectivity of the transformed region the same as that of the physical region, as in the

following examples showing two variations of this approach using interior slabs and slits,

respectively, in the transformed region. The slab configuration is as follows:

In coding, points inside the slab in the transformed region are simply skipped in all

computations.

This configuration introduces a special point of the following form at each of the

points corresponding to the slab corners in the transformed field:

The coordinate lines through

point 7 are shown below:

This type of special point, where the coordinate species changes on a smooth line, occurs

when a convex corner in the transformed field is identified with a point on a smooth contour

in the physical field. Both coordinate lines experience slope discontinuities at this point.

The slit configuration is as shown below:

(An obvious varition would be to have the slit vertical.) In this slit configuration, the point 5

and 6 are special points of the form shown on p. 26 characteristic of the slit configuration,

and will require special treatment in difference solutions.

The transformed region could, however, be made simply-connected by introducing a

branch cut in the physical region as illustrated below:

Conceptually this can be viewed as an opening of the field at the out and then a deformation

into a rectangle:

Here the coincident coordinate lines 1-2 and 4-3 form a branch cut, which becomes

re-entrant boundaries on the left and right sides of the transformed region. All derivatives are

continuous across this cut, and points at a horizontal distance outside the right-side boundary

in the transformed region are the same as corresponding points at the same horizontal

distance on the same horizontal line inside the left-side boundary, and vice versa. (In all

discussions of point correspondence across cuts, "distance" means distance in the

transformed region). In coding, the use of a layer of points outside each member of a pair of

re-entrant boundaries in the transformed region holding values corresponding to the

appropriate points inside the other boundary of the pair avoids the need for conditional

choices in difference representations, as discussed in Section 6 of this chapter.

Boundary values are not specified on the cut. (This cut is, of course, analogous to the

coincident 0 and 2 lines in the cylindrical coordinate system discussed above.) At the cut

we have the following coordinate line configuration, as may be seen from the conceptional

deformation to a rectangle:

so that the coordinate species and directions are both continuous across the cut.

This type of configuration is often called an O-type. Another possible configuration is

as shown below, often called a C-type:

Opening the field at the cut we have, conceptually,

with 1-2-3-4 to flatten to the bottom of the rectangle. Here the two members of the pair of

segments forming the branch cut are both on the same side of the transformed region, and

consequently points located at a vertical distance below the segment 1-2, at a horizontal

distance to the left of point 2, coincide with points at the same vertical distance above the

segment 4-3, at the same horizontal distance to the right of point 3. The point 2(3) is a

special point of the type shown on p. 26 for slit configurations.

The coordinate line configuration at the cut in this configuration is as follows:

where it is indicated that varies to the right on the upper side of the cut, but to the left on

the lower side. The direction of variation of also reverses at the cut, so that although the

species and slope of both lines are continuous across the cut, the direction of variation

reverses there.

It is possible to pass onto a different sheet across a branch cut, and discontinuities in

coordinate line species and/or direction occur only when passage is made onto a different

sheet. It is also possible, however, to remain on the same (overlapping) sheet as the cut is

crossed, in which case the species and direction are continuous, and this must be the

interpretation when derivatives are evaluated across the cut, as is discussed in Section 5 to

follow. These concepts are illustrated in the following figure, corresponding to the C-type

configuration given on p. 30:

In the present discussion of configurations, the behavior of the coordinate lines across the cut

will always be described in regard to the passage onto a different sheet, since this is in fact

the case in codes. It is to be understood that complete continuity can always be maintained

by conceptually remaining on the same sheet as the cut is crossed. Much of this complexity

can, however, be avoided with the use of an extra layer of points surrounding the

transformed region as will be discussed in Section 6.

Although in principle any region can be transformed into an empty rectangular block

through the use of branch cuts, the resulting grid point distribution may not necessarily be

reasonable in all of the region. Furthermore, an unreasonable amount of effort may be

required to properly segment the boundary surfaces and to devise an appropriate point

distribution thereon for such a transformation. Some configurations are better treated with a

computational field that has slits or rectangular slabs in it.

Regions of higher connectivity than those shown above are treated in a similar

manner. The level of connectivity may be maintained as in the following illustration:

Here one slit is made horizontal and one vertical just for generality of illustration. Both

could, of course, be of the same orientation. Slabs, rather than slits, could also have been

used. The example has three bodies.

With the transformed region made simply-connected we have, using two branch cuts,

a configuration related to the 0-type shown above for one internal boundary:

The conceptual opening here is as follows:

with segment 2-3-4-5-6-7 opening to the bottom. Here the pairs of segments (1-2,8-7) and

(3-4,6-5) are the branch cuts, which form re-entrant boundaries in the transformed region as

shown. In this case, points outside the right side of the transformed region coincide with

points inside the left side, and vice versa. This cut is of the form described on p. 30, where

both the coordinate species and direction are continuous across the cut. Points below the

bottom segment 3-4, to the left of point 4, coincide with points above the bottom segment

6-5 to the right of point 5. This cut is of the form discussed on p. 31, for which the

coordinate species is continuous across the cut but the direction changes there. There are a

number of other possibilities for placement of the two cuts on the boundary of the

transformed region, of course, some examples of which follow.

It is not necessary to reduce the connectivity of the region completely; rather, a slit or slab

can be used for some of the interior boundaries, while others are placed on the exterior

boundary of the transformed region.

One other possibility in two dimensions is the use of a preliminary analytical

transformation of infinity to a point inside some interior boundary, with the coodinates

resulting therefrom replacing the cartesian coordinates in the physical region. The grid

generation then operates from these transformed coordinates rather than from the cartesian

coordinates. This typically gives a fine grid near the bodies, but may give excessively large

spacing away from the body.

Thus, for example, if points on the two physical boundaries shown below

are transformed according to the complex transformation

z’ = 1/z

where z = x+iy and z’ = x’+iy’, infinity in the x,y system will transform to the origin in the

x’, y’ system, as shown below.

Then with the grid generated numerically from the x’, y’ system the following configuration

results:

References to the use of this approach are made in the survey of Ref. [1]. Somewhat related

to this are various two-dimensional configurations which arise directly from conformal

mapping, cf. Ref. [6] and the survey of Ives on this subject, Ref. [7]. (Conformal mapping is

discussed in Chapter X.)

C. Embedded Regions

In more complicated configurations, one type of coordinate system can be embedded

in another. A simple example of this is shown below, where an 0-type system surrounding an

internal boundary is embedded in a system of a more rectangular form, using what amounts

to a slit configuration.

The conceptual opening of this system is best understood in stages: First considering

only the embedded 0-type system surrounding the interior boundary, we have the region

inside the contour 12-13-6-9 opening as follows:

This then opens to the rectangular central portion of the transformed region shown above,

with the inner boundary contour 8-7-8 collapsing to a slit. The rest of the physical region

then opens as shown below:

These two regions then deform to rectangles and are fitted to the top and bottom of the

rectangle corresponding to the inner system along the contours 12-13 and 9-6 as shown.

Here points at a vertical distance below the segment 11-12 are coincident with points

at the same vertical distance below the segment 10-9 on the same vertical line, and vice

versa, with similar correspondence for the pair of segments 13-14 and 6-5. Points at a

horizontal distance to the left of the segment B-12, at a vertical distance above point 8,

coincide with points at the same horizontal distance to the right of the segment 8-9, at the

same vertical distance below point 8. Similar correspondence holds for the pair 7-13 and 7-6.

Boundary values are specified on the slit 8-7.

The composite system shown on p. 40 can also be represented as a slit configuration in

the transformed region:

with the inner system represented as

and the lower side of the slit considered re-entrant with the left half of the top boundary of

the rectangle corresponding to the inner system, the upper side of the slit being re-entrant

with the right half of this top boundary of the inner region. Now the conceptual opening is as

follows for the inner region:

Difference representations made above the slit thus would use points below the right half of

the top of the inner region in the transformed region, etc. Similarly, representation made

below the left half of the top of the inner region would use points below the slit. The slit is

thus a "black hole" into which coordinate lines from the outer system disappear, to reappear

as part of the inner system. The slit here, matched with the top of the inner system, is then

clearly a branch cut, and passage through the slit onto the inner system is simply passage

onto a different sheet.

Note that the embedded system has its own distinctive species and directions for the

coordinate lines, entirely separate from the outer system. Thus for the inner region the

directions are as follows:

while for the outer region they are as shown below:

Thus at a point on the upper interface, 12-13, between the systems the lines are as follows:

while on the lower interface, 9-6 they are as follows:

Thus both coordinates reverse direction at the lower interface although the species is

continuous, while both the species and directions are continuous across the upper interface.

This again corresponds to passage onto a different sheet, for the interface between the inner

and outer systems, i.e., the segments 12-13 and 9-6, is actually a branch cut.

The points 9(12) and 6(13) here require special notice. For example, at point 9 the

coordinate line configuration is as follows:

The lines through point 9 are as shown below

There are thus several changes in species and direction at this point. This type of special

point embodies the form which always occurs with the slit configuration, shown on p. 26,

and occurs here because the embedded region inside the contour 9-6-13-12 is essentially

contained inside a slit defined by the same set of numbers.

The above discussion refers to the slit configuration on p. 41. For the configuration on

p. 40, the lines in the outer region are still as diagrammed on p. 43, but the lines in the inner

region now are as follows:

The coordinate line species and direction given on p. 43 for the upper interface, 12-13, thus

applies here on the entire interface between the two regions.

An alternative treatment of the two special points is to place them inside cells as

shown below:

This results in a six-sided cell surrounding each of these two points which requires special

treatment as discussed in Chapter IV.

Embedded systems can also be constructed in the block configuration:

Here the top of the block, 7-8, in the outer system is re-entrant with the corresponding

segment, 7-8, on a portion of the top of the inner system. The left side of the block, 6-7, and

the bottom of the block, 6-5, are similarly re-entrant with single portions of the top of the

inner system. Finally, the right side of the block, 5-12-8, is re-entrant with two portions, 5-12

and 12-8, of the top of the inner system. Points outside one of these segments in one system

are thus located at corresponding positions inside the other segment of the re-entrant pair in

the other system. The slab sides, matched with the top of the inner system, are thus branch

cuts between the inner and outer systems.

Here the coordinate lines proceed as follows for the outer system:

while those for the inner system are the same as before, as shown on p. 42. This means that

on the left and right sides of the block, i.e., segments 6-7 and 5-8, the line directions are as

follows:

and on the top and bottom, segments 7-8 and 6-5, the directions are as shown below:

There are thus changes in coordinate species and/or direction that are different on each side

of the block.

The point 8 (and points 7,6 and 5) are special points of the following form:

The lines through the point are as shown below:

Here the special points occur in the field instead of on the boundary.

An example of a C-type system embedded in another C-type system is given next:

Here the conceptual opening is as follows: First, considering the system about the upper

body, we have the following configuration:

which, with the body collapsed to a slit, opens to the rectangle in the center of the

transformed region. Next consider the system about the other body:

This opens to a rectangle, with the body flattening to a portion of the bottom, which is fitted

to the first rectangle along the segment 11-13. Finally, the outermost portion opens as

follows:

which opens to a rectangle which is fitted to the first one along the segment 12-14.

Again the embedded region inside the contour 14-12-11-13 can be considered to lie

inside a slit. This contour, which forms the interface between the inner and outer systems, is

actually a branch cut between the two systems, across which there are discontinuities in

coordinate species and directon in the same manner as was discussed above for the previous

embedded system. Points below segment 16-12 coincide with points below segment 17-11 in

this case. Points to the left of segment 15-12, above point 15, are coincident with points to

the right of segment 15-11 below point 15. The slit here is formed of the segments 8-15 and

9-15. The coincident points 11 and 12 here must be taken as a point boundary in the physical

region, i.e., fixed at a specified value. Several special points of the types discussed above are

present here.

An alternative arrangement of the transformed region that corresponds to exactly the

same coordinate system in the physical region is as follows:

Here points below segment 3-4, to the left of point 4, coincide with points above segment

6-5, to the right of point 5. When calculations are made on or above the segment 12-14 on

the larger block, points below this segment coincide with points below the corresponding

segment on the smaller block. Similarly, when calculations are made on or below the

segment 13-11 on the larger block points above this segment coincide with points below the

corresponding segment on the smaller block. Finally, points below the segment 7-8, to the

left of point 8, on the smaller block are coincident with points above the segment 10-9, to the

right of point 9.

This configuration displays explicitly the correspondence of the embedded region

inside the contour 14-12-11-13 to a slit. Conceptually, coordinate lines from the main system

disappear into the slit and emerge into the embedded system. These coordinate lines thus are

continued from the main system onto another sheet representing the embedded system. This

concept of embedded systems, with continuation onto another sheet through a slit adds

considerable flexibility to the grid configurations and is of particular importance with

multiple boundaries and in three dimensions. The composite structure discussed in Section 4

removes much of the coding complexity associated with systems of this type.

D. Other Configurations

Another arrangement of cuts, where the species of coordinate changes on a continuous

line as the cut is crossed, is illustrated below. The transformed region in this case is

composed of three blocks connected by the cuts.

Here points outside one section are coincident with corresponding points inside the adjacent

section.

The coordinate line configuration on the interface on the right side of block A here is

as follows:

This same type of configuration occurs, in different orientations, on each of the interfaces.

These interfaces are branch cuts, so that passage onto the adjacent block amounts to passage

onto another sheet in the same manner discussed above.

As a final configuration for consideration in two dimensions, the following example

shows a case with fewer lines on one side of a slab than on the other. This does not

necessitate the use of different increments of the curvilinear coordinates in the numerical

expressions, because, as has been mentioned, these increments always cancel out anyway.

E. Three-dimensional Regions

All the general concepts illustrated in these examples extend directly to three

dimensions. Interior boundaries in the transformed region can become rectangular solids and

plates, corresponding to the slabs and slits, respectively, illustrated above for two

dimensions. Examples of three-dimensional configurations can be found in the surveys given

by Ref. [8] and [9].

It is also possible to use branch cuts, as illustrated above for two dimensions, to bring

the interior boundaries in the physical region entirely to the exterior boundary of the

transformed region:

Physical space

Computational space

The correspondence between the physical and transformed fields can, however, become

much more complicated in three dimensions, and considerable ingenuity may be required to

visualize this correspondence. For instance, the simple case of polar coordinates corresponds

to a rectangular solid with two opposing sides having the radial coordinate constant thereon,

and two re-entrant sides on which the longitudinal coordinate is constant at 0 and 2 ,

respectively (corresponding to the cut). The remaining two sides correspond to the north and

south polar axes, so that an axis opens to cover an entire side. There is thus a line, i.e., the

axis, in the physical region that corresponds to an entire side in the transformed region.

Three-dimensional grids may be constructed in some cases by simply connecting

corresponding points on two-dimensional grids generated on stacks of planes or curved

surfaces:

It should be noted, however, that this procedure provides no inherent smoothness in

the third direction, except in cases where the stack is formed by an analytical transformation,

such as rotation, translation or scaling, of the two-dimensional systems. An example of such

an analytical transformation of two-dimensional systems is the construciton of a

three-dimensional grid for a curved pipe by rotating and translating (and scaling if the

cross-sectional area of the pipe varies) two-dimensional grids generated for the pipe

cross-section so as to place these transformed two-dimensional grids normal to the pipe axis

at successive locations along the axis:

Another example is the rotation of a two-dimensional grid about an axis to produce an

axi-symmetric grid:

4. Composite Grids

All of the above concepts can be incorporated in a single framework, and the coding

complexity can be greatly reduced, by considering the physical field to be segmented into

sub-regions, bounded by four (six in 3D) generally curved sides, within each of which an

individual coordinate system is generated. The overall coordinate system, covering the entire

physical field, is then formed by joining the sub-systems at the sub-region boundaries. The

degree of continuity with which this juncture is made is a design consideration in regard to

the mode of application intended for the resulting grid.

This segmentation concept is illustrated in the figure below.

The locations of the interfaces between the sub-regions in the physical region are, of course,

arbitrary since these interfaces are not actual boundaries. These interfaces might be fixed,

i.e., the location completely specified just as in the case of actual boundaries, or might be left

to be located by the grid generation procedure. Also the coordinate lines in adjacent

sub-regions might be made to meet at the interface between with complete continuity:

---Interface

with some lesser degree of continuity, e.g., continuous line slope only:

---Interface

or with a discontinuity in slope:

or perhaps not to meet at all:

Naturally, progressively more special treatment at the interface will be required in numerical

applications as more degrees of line continuity at the interface are lost. Procedures for

generating segmented grids with various degrees of interface continuity are discussed later,

and conservative interface conditions are given in Ref. [52], [53].

Now, with regard to placing these concepts in the framework of segmentation, the

sides of an individual subregion (called a "block" hereafter) can be treated as boundaries on

which the coordinate points, and/or the coordinate line intersection angles, are specified, just

as is done for actual boundaries, or a side may be treated as one member of a pair of

re-entrant boundaries, i.e., one side of a branch cut in the physical region across which

complete continuity is established. The other member of the pair may be another side (or

portion thereof) of the same block or may be all (or part of) a side of an adjacent block in the

physical field. Recall that it is not necessary for a coordinate to remain of the same species

across a re-entrant boundary, since the passage is onto a different sheet. This can introduce

some coding complexity, but the treatment is straightforward, and in fact the coding can be

greatly simplified by using an extra layer of points surrounding each block as is discussed in

Section 6.

Some of the general concepts have been embodied in the two-dimensional code

discussed in Ref. [19] and in three recent three-dimensional codes, Ref. [13] Ref. [14], and

[51].

A. Simply-connected Regions

The first L-shaped simply-connected configuration on p. 21 can be interpreted as

being composed of three blocks, with the sides of adjacent blocks forming pairs of re-entrant

boundaries:

or two blocks with a portion of a side of one block re-entrant with an entire side of another

block:

Here, and in the examples to follow, solid lines correspond to physical boundaries, while the

dashed lines correspond to the interfaces between the blocks. The dashed arrows indicate the

linkage between the interfaces. (Obviously, any single block can be broken into any number

of blocks connected by re-entrant boundaries across adjacent sides.) In contrast, the

L-shaped configuration on p. 22 corresponds to the use of a single block. Similarly, the

configuration on p. 24 can be formed with three blocks:

while the first configuration for the same boundary on p. 25 is formed with a single block.

The slit configuration on p. 25 can be formed of three blocks:

or two blocks with only a portion of the adjacent sides of two blocks forming a re-entrant

boundary:

B. Multiply-connected Regions

The configuration with a single cut shown on p. 29 corresponds to the use of a single

block with the left and right sides here being the members of a pair of re-entrant boundaries:

The multiply-connected slab configuration on p. 27 can be broken into four blocks:

Other decompositions should also be immediately conceivable. The slit configuration on p.

28 can be formed with two blocks, again with only portions of adjacent sides serving as

re-entrant boundaries:

or into four blocks, with entire sides as re-entrant boundaries in all cases:

The double-body region on p. 34 opens to a single block as shown there, with portions

of sides as re-entrant boundaries. A five block configuration would use only entire sides as

re-entrant boundaries, however:

There is no real advantage, however, to the five-block system here.

C. Embedded Regions

The segmentation concept is most useful in the construction of embedded coordinate

systems. For instance, the system on p. 40 can be considered to be formed of three blocks as

follows:

Here portions of adjacent sides of the two larger blocks are re-entrant with each other, while

each of the remaining portions of these sides is re-entrant with half of one side of the smaller

block. The left and right sides of the smaller block are re-entrant with each other. This

configuration could also have been constructed with eight blocks:

with only entire sides being involved in re-entrant pairs as shown.

With embedded systems the coordinate species often changes as the re-entrant

boundary is crossed. These systems also show that the blocks need be physically adjacent

only in the physical field, and it is in this sense that "adjacent" is always to be interpreted.

The transformed (computational) field should always be viewed as only a bookkeeping

structure. Various constructions are possible for the configurations on p. 48 and 50, and a

two block structure was actually used on p. 50. A further example follows:

D. Three-dimensional Regions

For general three-dimensional configurations, it is usually very difficult to obtain a

reasonable grid with the entire physical region transformed to a single rectangular block. A

better approach in most cases is to segment the physical region into contiguous sub-regions,

each bounded by six curved surfaces, with each sub-region being transformed into a

rectangular block. An individual grid is generated in each sub-region:

These sub-region grids are patched together to form the overall grads, as in the

two-dimensional cases discussed above. Examples of the use of this segmentation in three

dimensions are found, in particular, in Ref. [11] and [12]. Others are noted in the survey

given by Ref. [9].

As noted above, complete continuity can be achieved at the sub-region interfaces by

noting the correspondence of points exterior to one sub-region with points interior to

another. The necessary bookkeeping can be accomplished, and the coding complexity can be

greatly reduced, by using an auxiliary layer of points just outside each of the six sides of the

computational region, analogous to the procedure mentioned above for two dimensions. A

correspondence is then established in the code between the auxiliary points and the

appropriate points just inside other sub-regions. This approach has recently been

incorporated in an internal region code, Ref. [13], and in two codes for general regions, Ref.

[14] and [51]. This is discussed in more detail in Section 6.

General three-dimensional regions can be built up using sub-regions as follows: First,

point distributions are specified on the edges of a curved surface forming one boundary of a

sub-region:

and a two-dimensional coordinate system is generated on the surface:

When this has been done for all surfaces bounding the sub-region, the three-dimensional

system within the sub-region is generated using the points on the surface grids as boundary

conditions:

In three dimensions it is possible for a line, e.g., a polar axis, in the physical region to

map to an entire side of the computational region as in the illistration below, where the axis

corresponds to the entire left side of the block:

The system illustrated here could be one of several identical blocks joining together to form

a complete system around the axis.

It is illustrated by an exercise that the occurence of a polar axis can be avoided, and

this facilitates the construction of a block structure. Thus a surface grid, having eight

"corners", analogous to the four "corners" on the circle in the 2D grid on p. 23, can be

constructed on the surface of a sphere. This serves much better than a latitude-longitude type

system for joining to adjacent regions. Similarly, the use of the four "corner" system, rather

than a cylindrical system, in a circular pipe allows T-sections and bifurcations to be treated

easily by a composite structure, c.f. Ref. [13].

Generally, grid configurations with polar axis should not be used in composite grid

structures.

E. Overlaid Grids

Another approach to complicated configurations is to overlay coordinate systems of

different types, or those generated for different sub-regions:

Here an appropriate grid is generated to fit each individual component of the configuration,

such that each grid has several lines of overlap with an adjacent grid. Interpolation is then

used in the region of overlap when solutions are done on the composite grid, with iteration

among the various grids. This approach has the advantage of simplicity in the grid

generation, in that the various sub-region grids are only required to overlap, not to fit.

However, there would appear to be problems if regions of strong gradients fall on the

overlap regions. Also the interpolation may have to be constructed differently for different

configurations, so that a general code may be hard to produce. Some applications of such

overlaid grids are noted in Ref. [5].

5. Branch Cuts

As has been noted in the above discussion of transformed field configurations, it is

possible for discontinuities in coordinate species and/or direction to occur at branch cuts, in

the sense of passage onto another sheet. Continuity can be maintained, however, by

conceptually remaining on the same overlapping sheet as the cut is crossed. All derivatives

thus do exist at the cut, but careful attention to difference formulations is necessary to

represent derivatives correctly across the cut. Although the correct representation can be

accomplished directly by surrounding the computational region with an extra layer of points,

as is discussed in Section 6, it is instructive to consider what is required of a correct

representation further here.

A. Point Correspondence

Points on re-entrant boundaries in the transformed region, i.e., on branch cuts in the

physical region, are not special points in the sense used above. Points on re-entrant

boundaries, in fact, differ no more from the other field points than do the points on the 0 and

2 lines in a cylindrical coordinate system. Care must be taken, however, to identify the

interior points coinciding with the extensions from such points beyond the field in the

transformed space. This correspondence was noted above in each of the configurations

shown above, being indicated by the dashed connecting lines joining the two members of a

pair of re-entrant boundaries. There are essentially four types of pairs of re-entrant

boundaries, as illustrated in the following discussion of derivative correspondence. In these

illustrations one exterior point, and its corresponding interior point, are shown for each case.

The converse of the correspondence should be evident in each configuration.

For the configurations involving a change in the coordinate species at the cut, not only

must the coordinate directions be taken into account as the cut is crossed, but also the

coordinate species may need to be interpreted differently from that established across the cut

in order to remain on the same sheet as the cut is crossed. For example, points on an -line

belonging to section A in the figure on p. 52, but located outside the right side of this region,

are coincident with points on a -line of region B at a corresponding distance (in the

transformed region) below the top of this region.

B. Derivative Correspondence

Care must be taken at branch cuts to represent derivatives correctly in relation to the

particular side of the cut on which the derivative is to be used. The existence of branch cuts

indicates that the transformed region is multi-sheeted, and computations must remain on the

same sheet as the cut is crossed. Remaining on the same sheet means continuing the

coordinate lines across the cut coincident with those of the adjacent region, but keeping the

same interpretation of coordinate line species and directions as the cut is crossed, rather than

adopting those of the adjacent region. As noted above, points outside a region across a cut

the transformed space are coincident with points inside the region across the other member

of the pair of re-entrant boundary segments corresponding to the cut in the transformed

space. The positive directions of the curvilinear coordinates to be used at these points inside

the region across the other member of the pair in some cases are the same as the defined

directions there, but in other cases are the opposite directions. As noted above, the

coordinate species may change also.

For cuts located on opposing sides of the transformed region, the proper form is

simply a continuation across the cut. Thus in the configuration on p. 29, with a computation

site on the right side of the transformed region, i.e., on the upper side of the cut in the

physical plane, we have points to the right of the site (below the cut in the physical plane)

coinciding with points to the right of the left side of the transformed region (below the cut in

the physical plane) as noted above. When -derivatives and -derivatives for use outside the

right side of the transformed region are represented inside the left side, the positive

directions of and to be used there are to the right and upward, respectively, as is

illustrated below. (In this and the following figures of the section, the dotted arrows indicate

the proper directions to be used at the interior points coincident with the required exterior

points, i.e., on the same sheet across the cut, while solid arrows indicate the locally

established directions for the coordinate lines, i.e., on a different sheet.)

With the two sides of the cut both located on the same coordinate line, i.e., on the

same side of the transformed region as in the configuration on p. 30, however, the situation

is not as simple as the above. In this case, when the computation site is on the left branch of

the cut in the transformed region (on the lower branch in the physical region), the points

below this boundary in the transformed region coincide with points located above the right

branch of the cut (above the cut in the physical region) at mirror-image positions, as has

been noted earlier. The -derivatives for use at such points below the left branch thus must

be represented at these corresponding points above the right branch. The positive direction

of for purposes of this calculation of derivatives above the right branch, for use below the

left branch, must be taken as downward, not upward. There is a similar reversal in the

interpretation of the positive direction of . This is in accordance with the discussion on p.

31. These interpretations are illustrated below:

In the configuration on p. 40, where two sides of a cut face each other across a void,

there is really no problem of interpretation, since the directions in the configuration are

treated simply as if the void did not exist. This correspondence is as shown below:

In all cases the interpretation of the positive directions of the curvilinear coordinates

must be such as to preserve the direction in the physical region, i.e., on the same sheet, as the

cut is crossed. In the cases where the coordinate species change at the cut, the situation is

even more complicated. Thus on the left side, segment 6-7, of the slab interface between the

inner and outer systems in the embedded configuration on p. 46, where the species changes

across the cut, the correspondence is as follows:

Thus, when a -derivative is needed outside the outer sytem, for use inside the left slab

interface, the positive -direction at the corresponding points inside the inner system must

be taken to coincide with the negative -direction of the inner system. Similarly, an

-derivative would be represented taking the positive -direction to coincide with the

positive -direction of the inner system. In an analogous fashion, a -derivative needed

outside the inner system, for use inside the segment 6-7, would be represented at the

corresponding point inside the outer system, i.e., to the left of the left slab side, but with the

positive -direction taken to be the positive -direction of the outer system. An -derivative

would be represented similarly, taking the positive -direction to be the negative -direction

of the outer system.

A -derivative to the left of the right side of the slab in the outer system would be

represented below segment 12-5 or 8-12, as the case may be, but with the positive

-direction taken to be the positive -direction of the inner system. Similarly, an

-derivative would be represented taking the positive -direction to be the negative

direction of the inner system. For a -derivative above the bottom of the slab in the outer

system, the correspondence is to below the segment 5-6 inside the inner system, with the

positive -direction taken to be the negative -direction of the inner system. The

-derivative is represented taking the positive -direction to be the negative -direction of

the inner sytem. Finally, for derivatives below the top of the slab in the outer system, the

correspondence is to below the segment 7-8 inside the inner system, with both the species

and direction of the coordinates unchanged.

The proper interpretation of coordinate species and direction across branch cuts for all

the other configurations discussed above can be inferred directly from these examples. A

conceptual joining of the two members of a pair of re-entrant boundaries in accordance with

the dashed line notation used on the configurations given in this chapter will always show

exactly how to interpret both the coordinate species and directions in order to remain on the

same sheet and thus to maintain continuity in derivative representation across the cut.

Examples of the proper difference representation are given in the following section. The

complexities of this correspondence can be completely avoided, however, by using

surrounding layers around each block in a segmented structure as discussed in the next

section.

6. Implementation

As discussed above, the transformed region is always comprised of contiguous

rectangular blocks by construction. This occurs because of the essential fact that one of the

curvilinear coordinates is defined as constant on each segment of the physical boundary.

Consequently, each segment of the physical boundary corresponds to a plane segment of the

boundary of the transformed region that is parallel to a coordinate plane there. The complete

boundary of the transformed region then is composed of plane segments, all intersecting at

right angles. Although the transformed region may not be a simple six-sided rectangular

solid, it can be broken up into a contiguous collection of such solids, here called blocks.

i cancel from all difference

Now it is noted in Chapter III that the increments

expressions, and that the actual values of the curvilinear coordinates i are immaterial. The

coordinates in the transformed region can thus be considered simple counters identifying the

points on the grid. This being the case, and the transformed region being comprised of a

collection of rectangular blocks, it is convenient to identify the grid points with integer

values of the curvilinear coordinates in each block, and thus to place the cartesian

coordinates of a grid point in ijk, where the subscripts (i,j,k) here indicate position ( i, 2,

3)

in the transformed region. (In coding, a fourth index may be added to identify the

block.) In each block, the curvilinear coordinates are then taken to vary as i = 1,2,...,Ii over

the grid points, where Ii is the number of points in the i-direction. Grid points on a

boundary segment of the transformed region will be placed in ijk with one index fixed.

Now each block has six exterior boundaries, and may also have any number of interior

boundaries (cf. the slab and slit configurations of Section 3), all of which will always be

plane segments intersecting at right angles, although the occur ence of interior boundaries

can be avoided if desired by breaking the block up into a collection of smaller blocks as

discussed in Section 4. The boundary segments in the transformed plane may correspond to

actual segments of the physical boundary, or may correspond to cuts in the physical region.

As discussed in Section 5, these cuts are not physical boundaries, but rather are interfaces

across which the field is re-entrant on itself. A boundary segment in the transformed region

corresponding to such a cut then is an interface across which one block is connected with

complete continuity to another block, or to another side of itself, several examples having

been given above in this chapter.

Depending on the type of grid generation system used (cf., the later chapters), the

cartesian coordinates of the grid points on a physical boundary segment may either be

specified or may be free to move over the boundary in order to satisfy a condition, e.g.,

orthogonality, or the angle at which coordinate lines intersect the boundary.

To set up the configuration of the transformed region, a correspondence is established

between each (exterior or interior) segment of the boundary of the transformed region and

either a segment of the physical boundary or a segment of a cut in the physical region. This

is best illustrated by a series of examples using the configurations of this chapter. The first

step in general is to position points on the physical boundary, or on a cut, which are to

correspond to corners of the transformed region (exterior or interior). As noted in Section 3,

these points do not have to be located at actual corners (slope discontinuities) on the physical

boundary.

For example, considering the two-dimensional simply-connected region on p. 23, four

points on the physical boundary are selected to correspond to the four corners of the empty

rectangle that forms the transformed region here:

Now, considering any one of these four points, one species of curvilinear coordinate will run

from that point to one of the two neighboring corner points, while the other species will run

to the other neighbor:

The corresponding species of coordinates will run to connect opposite pairs of corner points:

Since the curvilinear coordinates are to be assigned integer values at the grid points, i is to

vary from 1 at one corner to a maximum value, Ii, at the next corner, where Ii is the number

of grid points on the boundary segment between these two corners. Thus, proceeding

clockwise from the lower left corner, the cartesian coordinates of the four corner points are

placed in 1,1, 1,J, I,J, and I,1, where I1 = I and I2 = J.

The boundary specification is then completed by positioning I-2 points on the lower and

upper boundary segments of the physical region as desired, and J-2 points on the left and

right segments. The cartesian coordinates of these points on the lower and upper segments

are placed in i,1 and i,J, respectively, for i from 2 to I-1, and those on the left and right

segments are placed, respectively, in

1,j and

I,j for

j from 2 to J-1.

This process of boundary specification can be most easily understood by viewing the

rectangular boundary of the transformed region, with I equally-spaced points along two

opposite sides and J equally-spaced points along the other two sides, conceptually, as being

deformed to fit on the physical boundary. The corners can be located anywhere on the

physical boundary, of course. Here the point distribution on the sides can be conceptually

stretched and compressed to position points as desired along the physical boundary. The

cartesian coordinates of all the selected point locations on the physical boundary are then

placed in as described above.

This conceptual deformation of the rectangular boundary of the transformed region to

fit on the physical boundary serves to quickly illustrate the boundary specification for the

doubly-connected physical field shown on p. 29, which involes a cut. Thus I points are

positioned as desired clockwise around the inner boundary of the physical region from 2 to

3, and I points are positioned as desired, also in clockwise progression, around the outer

boundary from 1 to 4. The cartesian coordinates of these points on the inner boundary are

placed in i,1, and those on the outer boundary in i,J, with i from 1 to I. Note that here the

first and last points must coincide on each boundary, i.e., I,1 = 1,1 and I,J = 1,J. The

left and right sides of the transformed region (i=1 and i= I) are re-entrant boundaries,

corresponding to the cut, and hence values on these boundaries are not set but will be

determined by the generation system. The system must provide that the same value appears

on both of these sides, i.e., I,j = 1,j for all j from 2 to J-1.

The conceptual deformation of the rectangle for a C-type configuration is illustrated

on p. 31. Here, with I1 the number of points on the segments 1-2 and 3-4 (which must have

the same number of points), I2 points are positioned as desired around the inner boundary in

the physical region in a clockwise sense from 2 to 3, and the cartesian coordinates of these

points are placed in i,1 for i from I1 to I1+I2-1.

The first and last of these points must be coincident, i.e. I1,1 = I1 + I2-1,1. Now the

top, and the left and right sides, of the rectangle are deformed here to fit on the outer

boundary of the physical region. (In the illustration given, the two top corners are placed on

the two corners that occur in the physical boundary, a selection that is logical but not

mandatory.) The cartesian coordinates of the J points(positioned as desired on the segment

4-5 of the physical boundary) are placed in I,j, proceeding upward on the physical

boundary from 4 to 5 for j= 1 to J, and those on the segment 1-6 are placed in 1,j, but

proceeding downward on the physical boundary from 1 to 6 for j1 to J. Finally, the cartesian

coordinates of the I selected points on the physical boundary segment 6-5 are placed in i,J,

proceeding clockwise from 6 to 5 for i=1 to I. Since the same number of points must occur

on the top and bottom of the rectangle, we must have I=2(I1-1)+I2. Here the portions of the

lower side of the rectangle, i.e., i from 2 to I1-1, and from I1+I2 to I-1 with j=1, are

re-entrant boundaries corresponding to the cut, and hence no values are to be specified on

these segments. The generation system must make the correspondence i,1 = I-i+1,1 for

i=2 to I1-1 on these segments.

The conceptual deformation of the boundary of the transformed regions also serves for

the slab configuration on p. 27, where the interior rectangle deforms to fit the interior

physical boundary, while the outer rectangle deforms to fit the outer physical boundary. On

the inner boundary, the cartesian coordinates of J2-J1+1 selected points on the segment 5-8

of the physical boundary are placed in I1,j for j from J1 to J2, proceeding upward on the

physical boundary from 5 to 8, where J1 and J2 are the j-indices of the lower and upper

sides, respectively, of the interior rectangle and I1 is the i-index of the left side of this

rectangle. Similarly, J2-J1+1 points are positioned as desired on the segment 6-7 of the

physical boundary and are placed in I2,j, where I2 is the i-index of the right side of the

inner rectangle. Also I2-I1+1 points on the segments 5-6 and 8-7 of the physical boundary

are placed in i,J1 and i,J2, respectively, for i from I1 to I2, proceeding to the right on each

segment. The outer boundary is treated as has been described for an empty rectangle. Here

there will be J1-1 coordinate lines running from left to right below the inner boundary, and

J-J2 lines running above the inner boundary. Similarly, there will be I1-1 lines running

upward to the left of the interior boundary and I-I2 lines to the right. Thus the specifications

of the desired number of coordinate lines running on each side of the inner boundary serves

to determine the indicies I1, I2, J1, and J2. Note that the points inside the slab, i.e., I1 < i <

I2 and J1 < j < J2 are simply excluded from the calculation.

The slit configuration, illustrated on p. 28, can also be treated via the conceptual

deformation, but now with a portion of a line inside the rectangle opening to fit the interior

boundary of the physical region. This requires that provision be made in coding for two

values of the cartesian coordinates to be stored on the slit. If the i-indices of the slit ends, 5

and 6, are I1 and I2, respectively, then the cartesian coordinates of I2-I1+1 points positioned

as desired on the lower portion of the physical interior boundary, again proceeding from 5 to

6, are placed in a one-dimensional array, while the coordinates of the same number of points

selected on the upper portion of the physical interior boundary, again proceeding from 5 to 6,

are placed in another one-dimensional array. The first and last points in one of these arrays

must, of course, coincide with those in the other. Then the generation system must read

values into i,J1 for i from I1 to I2 (J1 being the j-index of the slit) from the former array for

use below the slit, or values from the latter array for use above. (As has been noted, the use

of a composite structure eliminates the need for these two auxiliary arrays.) Note that the

index values I1 and I2 are determined by the number of lines desired to run upward to the

left and right of the interior boundary, respectively, i.e., I1-1 lines on the left and I-I2 on the

right. Similarly, there will be J1-1 lines below the interior boundary, and J-J1 above.

Configurations, such as those illustrated on pp. 24-25, which involve slabs or slits that

intersect the outer boundary are treated similarly, with points inside the slab again being

simply excluded from the calculations. Also multiple slab or slit arrangements are treated by

obvious extensions of the above procedures. Here the indices corresponding to each slab or

slit will be determined by the number of points on the interior boundary segments and the

number of coordinate lines specified to run between the various boundaries. For example, in

the slit configuration shown on p. 33, the ends of the horizontal slit would be at i-indices I1

and I2, where I1-1 lines run vertically to the left of the slit and there are I2-I1+1 points on

the slit. The vertical slit would be at i=I3 where there are I3-I2-1 vertical lines between this

slit and the horizontal slit (and I-I2 lines to the right). Similarly, if the j-indices of the ends of

the vertical slit or J1 and J2, there will be J1-1 horizontal lines below this slit and J-J2 lines

above. With the j-index of the horizontal slit as J3, there will be J3-1 horizontal lines below

this slit and J-J3 above. Provision will now have to be made in coding for two

one-dimensional arrays for each slit to hold the cartesian coordinates of the points on the

segments of the physical interior boundaries corresponding to the two sides of each slit.

Again this coding complexity is avoided in the composite structure.

The use of the conceptual deformation of the rectangle to setup the boundary

configuration for the case with multiple interior boundaries on p. 34 should follow with little

further explanation. Here there must be the same number of points on the pair of segments

2-3 and 6-7, which correspond to the two segments forming the interior boundary on the

right. There must also be the same number of points on the pair, 3-4 and 5-6, corresponding

to the cut connecting the two interior boundaries. Finally the number of points on the outer

boundary must, of course, be the same as that on the bottom boundary. Note also that the

values of the cartesian coordinates placed at 2 must be the same as are placed at 7; those at 3

must be the same as those at 6, and those at 4 the same as at 5. Values are not set on the cuts,

of course, but the generation system must provide that values at points on the segment from

3 to 4 are the same as those on the segment 5-6, but proceeding from 6 to 5. Also values on

the segments 2-1 and 7-8 must be the same, proceeding upward in each case.

Following the conceptual deformation of the rectangular boundaries of the

transformed region and the indexing system illustrated above, it now should be possible to

set up the more complicated configurations such as the embedded regions shown in Section

3C. As noted there, however, the most straightforward and general approach to such more

complicated configurations is to divide the field into contiguous rectangular blocks, each of

which has its own intrinsic set of curvilinear coordinates and hence its own (i,j,k) indexing

system. The necessary correspondence between the individual coordinate systems across the

block interfaces was discussed in some detail in Section 3C. This block structure greatly

simplifies the setup of the configuration. For example, consider the 3-block structure shown

on p. 49 for the physical field shown on p. 48, for which the blocks are as follows:

Here the selected points on the right interior boundary (segment 8-15-9) are placed in i,1 of

the first block, for i from the i-index at 8 to that at 9,proceeding clockwise from 8 to 9 on the

physical boundary. (The difference between these two i-indices here is equal to the number

of points on this interior boundary,less one.) Similary,the selected points on the left interior

boundary (segment 4-5) are placed in i,1 of the second block for i from the i-index at 4 to

that at 5, proceeding clockwise from 4 to 5 on the physical boundary. The selected points on

the outer boundary of the physical region are placed in 1,j of the third block for j from 1 to

J3, in i,J3 for i from 1 to I3, and in I3,j for j from J3 to 1, proceeding from 16 to 1 to 2 to

14 on the physical boundary. Points on the remainder of the physical outer boundary are

placed in 1,j of the second block for j from 1 to J2 and in I2,j for j from J2 to 1,

proceeding from 3 to 17 for the former and from 13 to 6 for the latter, and in

1,j of

the first

block for j from 1 to Jl and in I1,j for j from Jl to 1, proceeding from 7 to 13 for the former

and from 14 to 10 for the latter.

Since the three blocks must fit together we have I3=I2, (I1+1)/2 equal to the difference

in i-indices between 11 and 13 in the second block and to that between 12 and 14 of the third

block. The quantities J1, J2, and J3 determine how many C-type lines occur in each block,

and can be chosen independently. Here the segment 11-13 on the top of the first block

interfaces with the corresponding segment on the top of the second block. The segment

12-14, which forms the remainder of the top of the first block, interfaces with the

corresponding segment on the bottom of the third block. Finally, the segment 12-16, which

forms the remainder of the bottom of the third block, interfaces with segment 11-17, which

forms the remainder of the top of the second block. The segments 3-4 and 6-5 on the bottom

of the second block interface with each other in the order indicated, as do also the segments

7-8 and 10-9 on the bottom of the first block.

In coding, this block structure can be handled by using a fourth index to identify the

block, placing an extra layer around each block, (i=0 and I+1, j=0 and J+1) and providing an

image-point array by which any point of any block can be paired with any point of any other,

or the same, block. Such pairs of points are coincident in the physical region, being on or

across block interfaces, and consequently are to be given the same values of the cartesian

coordinates by the generation system. This imaging extends to the extra layer surrounding

each block, so that appropriate points Inside other blocks can be identified for use in

difference representations on the block interfaces that require points outside the block, (cf.

Section 5).

Interface correspondence then can be established by input by setting the image-point

correspondence on the appropriate block sides, i.e., placing the (i,j,k) indices and block

number of one member of a coincident pair of points in the image-point array at the indices

and block number of the other member of the pair. This correspondence is indicated on the

block diagram on pp. 85-86 by the points enclosed in certain geometric symbols.

Thus, for the 3-block configuration considered above, the indices (I1-i+1, 1) and block

number 1, corresponding to a point on the segment 9-10 of the first block, would be placed

in the image-point array at the point (i,1) on the segment 7-8 of this block, and vice versa. A

similar pairing occurs for points on the segments 3-4 and 5-6 of the second block. The

indices (I2-i+1, J2) and block number 2, corresponding to a point on the segment 11-13 of

the second block would be placed in the image-point array at the point (i,Jl) of the first block

on the segment 13-11 of that block, and vice versa. The indices (I3-I1+i, 1) and block

number 3 (a point on the segment 12-14 of the third block) would be placed in the array at

the point (i,J1) of the first block for a point on the segment 12-14 of that block. Finally, the

indices (i,1) and block number 3 (point on segment 16-12 of the third block) would be

placed in the array at the point (i,J2) of the second block for a point on the segment 17-11 of

that block. The remaining segments all correspond at portions of the physical boundary and

hence do not have image points.

In the same manner the following image correspondence can be set between interior

points and points on the surrounding layers in order to establish difference representations

across the block interfaces: (This correspondence is indicated symbolically on the block

diagram on pp. 85-86 by geometric symbols.):

As noted, all of this information would be input into the image-point array. Then with

values of the cartesian coordinates at the image points on the surrounding layer set equal to

those at the corresponding object point inside one of the blocks, it is possible to use the same

difference representations on the interfaces that are used in the interior.

The discussion given in this chapter should now allow the image-point input to be

constructed for any configuration of interest. As noted, it is not necessary that the coordinate

species remain the same as the interface is crossed. Thus, for instance, a point on the right

side of one block could be paired with one on the bottom of another block. In such a case the

image point of the point (I+1, j) the first block would be the point (j,2) inside the second

block. Similarly the image of the point (i,0) below the second block would be the point (I-1,

i) inside the first block, The correct difference representation across interfaces is thus

automatically established, eliminating the need for the concern with passage onto different

sheets discussed in detail earlier in this chapter.

This greatly simplifies the coding, since with the surrounding layers and the use of the

image points, all of the derivative correspondences are automatic and do not have to be

specified for each configuration. It is only necessary to specify the point correspondence by

input. This construction also allows codes for the numerical solution of partial differential

equations on the grid to be written to operate on rectangular blocks. Then any configuration

can be treated by sweeping over all the blocks. The surrounding layers of points and the

image correspondence provide the proper linkage across the block interfaces. In an implicit

solution the values on the interfaces would have to be updated iteratively in the course of the

solution. The solution for the generation of the grid would similarly keep the interface and

surrounding layer values updated during the course of the iterative solution.

This, of course, maintains completecontinuity across the block interfaces. If complete

continuity is not required,

Numerical grid generation has now become a fairly common tool for use in the

numerical solution of partial differential equations on arbitrarily shaped regions. This is

especially true in computational fluid dynamics, from whence has come much of the impetus

for the development of this technique, but the procedures are equally applicable to all

physical problems that involve field solutions. Numerically generated grids have provided

the key to removing the problem of boundary shape from finite difference methods, and

these grids also can serve for the construction of finite element meshes. With such grids all

numerical algorithms, finite difference or finite element, are implemented on a square grid in

a rectangular computational region regardless of the shape and configuration of the physical

region. (Finite volume methods are effectively a type of conservative finite difference

method on these grids.)

In this text, grid generation and the use thereof in numerical solutions of partial

equations are both discussed. The intent was to provide the necessary basic information,

from both the standpoint of mathematical background and from that of coding

implementation, for numerical solutions of partial differential equations to be constructed on

general regions. Since these numerical solutions are ultimately constructed on a square grid

in a rectangular computational region, any solution algorithm that can treat equations with

variable coefficients is basically applicable, and therefore discussion of specific algorithms is

left to classical texts on the numerical solution of partial differential equations.

The area of numerical grid generation is relatively young in practice, although its roots

in mathematics are old. This somewhat eclectic area involves the engineer’s feel for physical

behavior, the mathematician’s understanding of functional behavior, and a lot of

imagination, with perhaps a little help from Urania. The physics of the problem at hand must

ultimately direct the grid points to congregate so that a functional relationship on these

points can represent the physical solution with sufficient accuracy. The mathematics controls

the points by sensing the gradients in the evolving physical solution, evaluating the accuracy

of the discrete representation of that solution, communicating the needs of the physics to the

points, and by providing mutual communication among the points as they respond to the

physics.

Numerical grid generation can be thought of as a procedure for the orderly distribution

of observers, or sampling stations, over a physical field in such a way that efficient

communication among the observers is possible and that all physical phenomena on the

entire continuous field may be represented with sufficient accuracy by this finite collection

of observations. The structure of an intersecting net of families of coordinate lines allows the

observers to be readily identified in relation to each other, and results in much more simple

coding than would the use of a triangular structure or a random distribution of points. The

grid generation system provides some influence of each observer on the others, so that if one

moves to get into a better position for observation of the solution, its neighbors will follow to

some extent in order to maintain smooth coverage of the field.

Another way to think of the grid is as the structure on which the numerical solution is

built. As the design of the lightest structure requires consideration of the load distribution, so

the most economical distribution of grid points requires that the grid be influenced by both

the geometric configuration and by the physical solution being done thereon. In any case,

since resources are limited in any numerical solution, it is the function of the numerical grid

generation to make the best use of the number of points that are available, and thus to make

the grid points an active part of the numerical solution.

This is a rapidly developing area, being now only about ten years old, and thus is still

in search of new ideas. Therefore no book on the subject at this time could possibly be

considered to be definitive. However, enough material has now accumulated in the literature,

and enough basic concepts have emerged, that a fundamental text is now needed to meet the

needs of the rapidly expanding circle of interest in the area. It is with the knowledge of both

these needs and these limitations that this text has been written. Some of the techniques

discussed will undoubtedly be superceded by better ideas, but the fundamental concepts

should serve for understanding, and hopefully also for some inspiration, of new directions.

The only background assumed of the student is a senior-level understanding of numerical

analysis and partial differential equations. Concepts from differential geometry and tensor

analysis are introduced and explained as needed.

Numerical grid generation draws on various areas of mathematics, and emphasis

throughout is placed on the development of the relations involved, as well as on the

techniques of application. This text is intended to provide the student with the understanding

of both the mathematical background and the application techniques necessary to generate

grids and to develop codes based on numerically generated grids for the numerical solution

of partial differential equations on regions of arbitrary shape.

The writing of this text has been a cooperative effort over the last two years, spurred

on by the institution of a graduate course in numerical grid generation, as well as an annual

short course, at Mississippi State. The students in both of these courses have contributed

significantly in revising the text as it evolved. The last appendix is the result of a class

assignment prepared by Col. Hyun Jin Kim, graduate student in the computational fluid

dynamics program, who also compiled the index. Our colleage, Dr. Helen V.

McConnaughey of Mathematics contributed significantly through continual discussions and

wrote most of Chapter IV.

We are indebted to a large number of former students and fellow researchers around

the world for the development of the ideas that have crystallized into numerical grid

generation. The complete debt can be acknowledged only through mention of the

bibliographies contained in the several surveys cited herein. A list here would either be too

long to note the strongest influences or too short to acknowledge all the significant ones. We

must, however, acknowledge the many long and fruitful discussions with Peter Eiseman of

Columbia University.

0f vital importance is the support that has been provided for the research from which

the developments discussed in this book have emerged, including NASA; the research

offices of the Air Force, Army, and Navy; the National Science Foundation, and various

industrial concerns. The interest and contributions of a number of contract monitors has been

essential over the years. We are especially appreciative of Bud Bobbitt and Jerry South of

NASA Langley Research Center, who provided the initial support for an unknown with an

idea.

Particular debts are owed to W. H. Chu for an idea in the Journal of Computational

Physics in 1971, and to Frank Thames who put the idea into a dissertation.

In the preparation of the text we had the conscientious and untiring efforts of two most

able secretaries, Rita Curry and Susan Triplett, who typed on in good spirits through a year

of numerous revisions and frustrations as the text evolved.

Finally, we were particularly fortunate to have the services of Yeon Seok Chae,

graduate student in the computational fluid dynamics program and illustrator par excellence,

who did all the figures with understanding of the intended meaning as well as artistic

competence. His meticulous efforts were extensions of our thoughts.

Joe F. Thompson

Z. U. A. Warsi

C. Wayne Mastin

Mississippi State, Mississippi

January 1985

I. INTRODUCTION

The numerical solution of partial differential equations requires some discretization of

the field into a collection of points or elemental volumes (cells). The differential equations

are approximated by a set of algebraic equations on this collection, and this system of

algebraic equations is then solved to produce a set of discrete values which approximates the

solution of the partial differential system over the field. The discretization of the field

requires some organization for the solution thereon to be efficient, i.e., it must be possible to

readily identify the points or cells neighboring the computation site. Furthermore, the

discretization must conform to the boundaries of the region in such a way that boundary

conditions can be accurately represented. This organization is provided by a coordinate

system, and the need for alignment with the boundary is reflected in the routine choice of

cartesian coordinates for rectangular regions, cylindrical coordinates for circular regions,

etc., to the extent of the handbook’s resources.

The current interest in numerically-generated, boundary-conforming coordinate

systems arises from this need for organization of the discretization of the field for general

regions, i.e., to provide computationally for arbitrary regions what is available in the

handbook for simple regions. The curvilinear coordinate system covers the field and has

coordinate lines (surfaces) coincident with all boundaries. The distribution of lines should be

smooth, with concentration in regions of strong solution variation, and the system should

ultimately be capable of sensing these variations and dynamically adjusting itself to resolve

them.

A numerically-generated grid is understood here to be the organized set of points

formed by the intersections of the lines of a boundary-conforming curvilinear coordinate

system. The cardinal feature of such a system is that some coordinate line (surface in 3D) is

coincident with each segment of the boundary of the physical region. The use of coordinate

line intersections to define the grid points provides an organizational structure which allows

all computation to be done on a fixed square grid when the partial differential equations of

interest have been transformed so that the curvilinear coordinates replace the cartesian

coordinates as the independent variables.

This grid frees the computational simulation from restriction to certain boundary

shapes and allows general codes to be written in which the boundary shape is specified

simply by input. The boundaries may also be in motion, either as specified externally or in

response to the developing physical solution. Similarly, the coordinate system may adjust to

follow variations developing in the evolving physical solution. In any case, the

numerically-generated grid allows all computation to be done on a fixed square grid in the

computational field which is always rectangular by construction.

In the sections which follow, various configurations for the curvilinear coordinate

system are discussed in Chapter II. In general, the computational field will be made

rectangular, or composed of rectangular sub-regions, and a wide variety of configurations is

possible. Coordinate systems may also be generated separately for sub-regions in the

physical plane and patched together to form a complete system for complex configurations.

The basic transformation relations applicable to the use of general curvilinear coordinate

systems are developed in Chapter III; the construction of numerical solutions of partial

differential equations on those systems is discussed in Chapter IV; and consideration is given

in Chapter V to the evaluation and control of truncation error in the numerical

representations.

Basically, the procedures for the generation of curvilinear coordinate systems are of

two general types: (1) numerical solution of partial differential equations and (2)

construction by algebraic interpolation. In the former, the partial differential system may be

elliptic (Chapter VI), parabolic or hyperbolic (Chapter VII). Included in the elliptic systems

are both the conformal (Chapter X), and the quasi-conformal mappings, the former being

orthogonal. Orthogonal systems (Chapter IX) do not have to be conformal, and may be

generated from hyperbolic systems as well as from elliptic systems. Some procedures

designed to produce coordinates that are nearly orthogonal are also discussed. The algebraic

procedures, discussed in Chapter VIII, include simple normalization of boundary curves,

transfinite interpolation from boundary surfaces, the use of intermediate interpolating

surfaces, and various other related techniques.

Coordinate systems that are orthogonal, or at least nearly orthogonal near the

boundary, make the application of boundary conditions more straightforward. Although

strict orthogonality is not necessary, and conditions involving normal derivatives can

certainly be represented by difference expressions that combine one-sided differences along

the line emerging from the boundary with central expressions along the boundary, the

accuracy deteriorates if the departure from orthogonality is too large. It may also be more

desirable in some cases not to involve adjacent boundary points strongly in the

representation, e.g., on extrapolation boundaries. The implementation of algebraic turbulence

models is more reliable with near-orthogonality at the boundary, since information on local

boundary normals is usually required in such models. The formulation of boundary-layer

equations is also much more straightforward and unambiguous in such systems. Similarly,

algorithms based on the parabolic Navier-Stokes equations require that coordinate lines

approximate the flow streamlines, and the lines normal thereto, especially near solid

boundaries. It is thus better in general, other considerations being equal, for coordinate lines

to be nearly normal to boundaries.

Finally, dynamically-adaptive grids are discussed in Chapter XI. These grids

continually adapt during the course of the solution in order to follow developing gradients in

the physical solution. This topic is at the frontier of numerical grid generation and may well

prove to be one of its most important aspects.

The emphasis throughout is on grids formed by the intersections of coordinate lines of

a curvilinear coordinate system, as opposed to the covering of a field with triangular

elements or a random distribution of points. Neither of these latter collections of points is

suitable for really efficient numerical solutions (although numerical representations can be

constructed on each, of course) because of the cumbersome process of identification of

neighbors of a point and the lack of banded structure in the matrices. Thus the subject of

triangular mesh generators, per se, is not addressed here. (Obviously a triangular mesh can

be produced by construction rectangular mesh diagonals.)

Considerable progress is being made toward the development of the techniques of

numerical grid generation and toward casting them in forms that can be readily applied. A

comprehensive survey of numerical grid generation procedures and applications thereof

through 1981 was given by Thompson, Warsi, and Mastin in Ref. [1], and the conference

proceedings published as Ref. [2] contains a number of expository papers on the area, as

well as current results. Other collections of papers on the area have also appeared (Ref. [3]

and [4]), and a later review through 1983 has been given by Thompson in Ref. [5]. Some

other earlier surveys are noted in Ref. [1]. A later survey by Eiseman is given in Ref. [37].

The present text is meant to be a developmental treatment of the techniques of grid

generation and its applications, not a survey of results, and therefore no attempt is made here

to cite all related references, rather only those needed to illustrate particular points are noted.

The surveys mentioned above should be consulted directly for references to examples of

various applications and related contributions. (Ref [l] gives a short historical development

of the ideas of grid generation.) Other surveys of particular areas of grid generation are cited

later as topics are introduced.

Finally, in regard to implementation, a configuration for the transformed

(computational) field is first established as discussed in Chapter II. The grid is generated

from a generation system constructed as discussed in Chapters VI -- X. (If the grid is to be

adaptive, i.e., coupled with the physical solution done thereon, then the gr1d must be

continually updated as discussed in Chapter XI.) In the construction of the grid, due account

must be taken of the truncation error induced by the grid discussed in Chapter V. The partial

differential equations of the physical problem of interest are transformed according to the

relations given in Chapter III. These transformed equations are then discretized, cf. Chapter

IV, and the resulting set of algebraic equations is solved on the fixed square grid in the

rectangular transformed field.

II. BOUNDARY-CONFORMING COORDINATE SYSTEMS

1. Basic Concepts

To provide a familiar ground from which to view the general development to follow,

consider first a two-dimensional cylindrical coordinate system covering the annular region

between two concentric circles:

Here the curvilinear coordinates (r, ) vary on the intervals [r1,r2] and [0,2 ], respectively.

These curvilinear coordinates are related to the cartesian coordinates (x,y) by the

transformation equations

(1)

The inverse transformation is given by

(2)

Note that one of the curvilinear coordinates, r, is constant on each of the phys1cal

boundaries, while the other coordinate, , varies monotonically over the same range around

each of the boundaries. Note also that the system can be represented as a rectangle on which

the two physical boundaries correspond to the top and bottom sides:

The transformed region, i.e., where the curvilinear coordinates, r and the independent

variables, thus can be thought of as being rectangular, and can be treated as such from a

coding standpoint. These points will be central to what follows.

The curvilinear coordinates (r, ) can be normalized to the interval [0,1] by introducing

the new curvilinear coordinates ( , ), where

(3)

or

(4)

The transformation then may be written

(5a)

(5b)

where now and both vary on the interval [0,1]. This is thus a mapping of the annular

region between the two circles in the physical space onto the unit square in the transformed

space, i.e., each point (x,y) on the annulus corresponds to one, and only one, point ( , ) on

the unit square:

The bottom ( = 0) and top ( = 1) of the square correspond, respectively, to the inner and

outer circles, r = r1, and r = r2. The sides of the square, = 0 and = 1 correspond to = 0

and = 2 , respectively, and hence to the two coincident sides of a branch cut in the

physical space. Therefore, boundary conditions are not to be specified on these sides of the

unit square in the transformed space. Rather these sides are to be considered re-entrant on

each other with points adjacent to one, outside the square, being equivalent to points adjacent

to the other, inside the square.

Conceptually, the physical region can be considered to have been opened at the cut

= 0 and 2 and then deformed into a rectangle to form the transformed region:

Here, point correspondence across the re-entrant boundaries (indicated by the dashed

connecting line) in the transformed region is illustrated by the coincidence of the pair of

circled points. This conceptual device and mode of illustration for the the point

correspondence across re-entrant boundaries will serve later for more general configurations.

These simple concepts extend to more complicated two-dimensional configurations,

the central feature being that one of the curvilinear coordinates is made to be constant on a

boundary curve (as was r above), while the other varies monotonically along that boundary

curve (as does ). The transformation to the rectangle is achieved by making the range and

direction of variation of the varying coordinate the same on each of two opposing boundaries

(as varies from 0 to 2 on each circle above).

The physical space thus transforms to the rectangle shown above regardless of the shape of

the physical region. (It is not necessary to normalize the curvilinear coordinates to the

interval [0,1], and in fact, any normalization can be used. In computational applications the

normalization is more conveniently done to different intervals for each coordinate. The field

in the transformed space is then rectangular, rather than square.) Familiar examples of this

are elliptical coordinates for the region between two confocal ellipses, spherical coordinates

for two spheres, parabolic coordinates for two parabolas, etc.

These, same concepts will be extended later to completely general configurations

involving any number of boundary curves and branch outs. The extension to three

dimensions follows directly, using boundary surfaces instead of curves, i.e., one curvilinear

coordinate will be made constant on a boundary surface, with the other two forming a

two-dimensional coordinate system on the surface.

Returning to the concentric circles, if the functional dependence of on , and/or that of r

on , had been made more general than the simple linear normalizations given by Eq. (4),

the corresponding coordinate lines would have become unequally spaced in the physical

space, while remaining as radial lines and concentric circles:

The transformation, from Eq. (1), is now given by

(6a)

(6b)

In this case the points on the inner and outer circular boundaries are not equally spaced

around the circles in the physical space for equal increments of , although they remain

equally spaced on the top and bottom of the unit square in the transformed space by

construction. The spacing around these circles is determined by the functional dependence of

on , and, since the points are located at equal increments of by construction, this

functional relationship is defined by the placement of these points around the circles. This

point, that the coordinate system in the field is determined from the boundary point

distribution, will be central to the discussion of grid generation to follow. The distribution of

circumferential lines is controlled here by the functional relationship between r and , which

is not related to any boundary point distribution. Thus factors other than the boundary point

distribution may be expected to be involved in grid generation, as well. That the point

distribution on the boundaries may be controlled by direct placement of the points, while the

coordinate line distribution in the field must be controlled by other means will also continue

to appear in the developments to follow.

The one-dimensional functional relationship between and in Eq. (6) requires that

the relative distributions of boundary points around the inner and outer circles be the same.

This restriction can be removed by making a function of , as well as of , while

retaining the periodic nature of the dependence on . In this case the ooordinate lines of

constant will no longer be straight radial lines, although they will continue to connect

corresponding points on the inner and outer circular boundaries. Similarly the

circumferential coordinate lines (lines of constant here) can be made to depart from circles

by making r dependent on both and , but with the restriction that the dependence

vanishes on the inner and outer circular boundaries (where

here).

= 0 and

= 1, respectively,

Obviously certain constraints will have to be placed on the functions ( , ) and r( , ) to

keep the mapping one-to-one. All of these considerations will reappear in the general

developments that follow.

Finally, it should be realized that the intermediate use here of the cylindrical

coordinates (r, ) in defining the transformation between the curvilinear coordinates ( , )

and the cartesian coordinates (x,y) has been only in deference to the familiarity of the

cylindrical coordinates, and such intermediary coordinates will not appear in general. The

generalized statement for the simple configuration under consideration here is as follows:

Find (x,y) and (x,y) in the annular region bounded by the curves x2 + y2 = and x2 + y2

= , subject to the boundary conditions

Specified monotonic variation of over [0,1] on x2 + y2 =

with same sense of direction on each of these two curves.

and on x2 + y2 =

It is the inverse problem that will be treated in fact, however, i.e., find x( , ) and

y( , ) on the unit square in the transformed space (0

1, 0

1), subject to the

boundary conditions

x( ,0) and y( ,0) specified on = 0 such that x2( = 0) + y2 ,0) =

x( ,1) and y( ,1) specified on = 1 such that x2( = 1) + y2( ,1) =

Periodicity in : x(1 + , ) = x( , ) y(1 + , ) = y( , )

The simple form for the transformation given by Eq. (6) is made possible by choosing the

same functional dependence of x and y on on the boundaries, = 0 and = 1. The

familiar cylindrical coordinate system is thus a special case of the general grid generation

problem for this simple configuration applicable to the region between two concentric

circles, as is the elliptical coordinate system for two ellipses, etc.

2. Generalization

Generalizing from the above consideration of cylindrical coordinates, the basic idea of

a boundary-conforming curvilinear coordinate system is to have some coordinate line (in 2D,

surface in 3D) coincident with each boundary segment, analogous to the way in which lines

of constant radial coordinate coincide with circles in the cylindrical coordinate system. The

other curvilinear coordinate, analogous to the angular coordinate in the cylindrical system,

will vary along the boundary segment and clearly must do so monotonically, else the same

pair of values of the curvilinear coordinates will occur at two different physical points. (It

should be clear that the curvilinear coordinate that varies along a boundary segment must

have the same direction and range of variation over some opposing segment, e.g., as the

angular variable varies from 0 to 2 over both of two concentric circles in cylindrical

coordinates).

With the values of the curvilinear coordinates thus specified on the boundary, it then

remains to generate values of these coordinates in the field from these boundary values.

There must, or course, be a unique correspondence between the cartesian (or other basis

system) and the curvilinear coordinates, i.e., the mapping of the physical region onto the

transformed region must be one-to-one, so that every point in the physical field corresponds

to one, and only one, point in the transformed field, and vice versa. Coordinate lines of the

same family must not cross, and lines of different families must not cross more than once.

In this chapter a two-dimensional region will be considered in most of the discussions

in the interest of economy of presentation. Generalization to three dimensions will be evident

in most cases and will be mentioned specifically only when necessary. As noted above, the

curvilinear coordinates may be normalized to any intervals, just as the radial and angular

coordinates of the cylindrical coordinate system can be expressed in many different units.

Since the interest of the present discussion is numerical application, it will be generally

convenient to define the increments of all the curvilinear coordinates to be uniformly unity,

and then to normalize these coordinates to the interval [1,N(i)], where N(i) is the total number

of grid points to be used in the i direction. (The three curvilinear coordinates will be

indicated as i, i = 1,2,3, in general. In two dimensions, however, the notation ( , ) will

often be used for the two coordinates 1 and 2 .) The computational field, i.e., the field in

the transformed space, thus will have rectangular boundaries and will be covered by a square

grid. (It will become clear later that the actual values of the increments in the curvilinear

coordinates are immaterial since they do not appear in the final numer1cal expressions.

Therefore no generality is lost in making the grid square and of unit increment in the

transformed field.)

A. Boundary-value Problem -- Physical Region

The generation of the curvilinear coordinate system may be treated as follows: with

the curvilinear coordinates specified on the boundaries, e.g., (x,y) and (x,y) on a

boundary curve (this specification amounting to a constant value for either or on each

segment of , with a specified monotonic variation of the other over the segment), generate

the values, (x,y) and (x,y), in the field bounded by . This is thus a boundary value

problem on the physical field with the curvilinear coordinates ( , ) as the dependent

variables and the cartesian coordinates (x,y) as the independent variables, with boundary

conditions specified on curved boundaries:

(In these discussions, the transformation is assumed to be from cartesian coordinates in the

physical space. The transformation can, however, be from any system of coordinates in the

physical space.)

B. Boundary value Problem - Transformed Region

The problem may be simplified for computation, however, by first transforming so

that the physical cartesian coordinates (x,y) become the dependent variables, with the

curvilinear coordinates ( , ) as the independent variables. Since a constant value of one

curvilinear coordinate, with monotonic variation of the other, has been specified on each

boundary segment, it follows that these boundary segments in the physical field will

correspond to vertical or horizontal lines In the transformed field. Also, since the range of

variation of the curvilinear coordinate varying along a boundary segment has been made the

same over opposing segments, it follows that the transformed field will be composed of

rectangular blocks.

The boundary value problem in the transformed field then involves generating the

values of the physical cartesian coordinates, x( , ) and y( , ), in the transformed field

from the specified boundary values of x( , ) and y( , ) on the rectangular boundary of the

transformed field, the boundary being formed of segments of constant or , i.e., vertical or

horizontal lines. With = constant on a boundary segment, and the increments in taken to

be uniformly unity as discussed above, this boundary value specification is implemented

numerically by distributing the points as desired along the boundary segment and then

assigning the values of the cartesian coordinates of each successive point as boundary values

at the equally spaced boundary points on the bottom (or top) of the transformed field in the

following figure.

Boundary values are not specified on the left and right sides of the transformed field since

these boundaries are re-entrant on each other (analogous to the 0 and 2 lines in the

cylindrical system), as discussed above, and as indicated by the connecting dotted line on the

figure. Points outside one of these re-entrant boundaries are coincident with points at the

same distance inside the other. The problem is thus much more simple in the transformed

field, since the boundaries there are all rectangular, and the computation in the transformed

field thus is on a square grid regardless of the shape of the physical boundaries.

With values of the cartesian coordinates known in the field as functions of the

curvilinear coordinates, the network of intersecting lines formed by contours (surfaces in

3D) on which a curvilinear coordinate is constant, i.e., the curvilinear coordinate system,

provides the needed organization of the discretization with conformation to the physical

boundary. It is also possible to specify intersection angles for the coordinate lines at the

boundaries as well as the point locations.

3. Transformed Region Configurations

As noted above, the generation of the curvilinear coordinate system is done by

devising a scheme for determination of the field values of the cartesian coordinates from

specified values of these coordinates (and/or curvilinear coordinate line intersection angles)

on portions of the boundary of the transformed region. Since the boundary of the

transformed region is comprised of horizontal and vertical line segments, portions of which

correspond to segments of the physical boundary on which a curvilinear coordinate is

specified to be constant, it should be evident that the configuration of the resulting

coordinate system depends on how the boundary correspondence is made, i.e., how the

transformed region is configured.

Some examples of different configurations are given below, from which more

complex configurations can be inferred. In these examples only a minimum number of

coordinate lines are shown in the interest of clarity of presentation tation. In all of these

examples, boundary values of the physical cartesian coordinates (and/or curvilinear

coordinate line intersection angles) are understood to be specified on all boundaries, both

external and internal, of the transformed region except for segments indicated by dotted

lines. These latter segments correspond to branch cuts in the physical space, as is explained

in the examples in which they appear. Such re-entrant boundary segments always occur in

pairs, the members of which are indicated by the dashed connecting lines on each of the

configurations shown. Points outside the field across one segmentof such a pair are

coincident with points inside the field across the other member of the pair. The conceptual

device of opening the physical field at the cuts is used here to help clarify the

correspondence between the physical and transformed fields. In many cases an example of

an actual coordinate system is given as well. References to the use of various configurations

may be found in the surveys given by Ref. [1] and [5], and a number of examples appear in

Ref. [2].

A. Simply-connected Regions

It is natural to define the same curvilinear coordinate to be constant on each member

of a pair of generally opposing boundary segments in the physical plane. Thus, a

simply-connected region formed by four curves is logically treated by transforming to an

empty rectangle:

Similarly, an L-shaped region could remain L-shaped in the transformed region:

Here, for instance, the cartesian coordinates of the desired points on the physical boundary

segment 4-5 are specified as boundary conditions on the vertical line 4-5, in corresponding

order, which forms a portion of the boundary of the transformed region.

The generalization of these ideas to more complicated regions is obvious, the

transformed region being composed of contiguous rectangular blocks. An example follows:

The physical boundary segment on which a single curvilinear coordinate is constant

can have slope discontinuities, however, so that the L-shaped region above could have been

considered to be composed of four segments instead of six, so that the transformed region

becomes a simple rectangle:

Here the cartesian coordinates of the desired points on the physical boundary 5-4-3 are the

specified boundary values from left to right across the top of the transformed region.

Whether or not the boundary slope discontinuity propagates into the field, so that the

coordinate lines in the field exhibit a slope discontinuity as well, depends on how the

coordinate system in the field is generated, as will be discussed later.

It is not necessary that corners on the boundary of the transformed region correspond

to boundary slope discontinuities on the physical boundary and a counter-example follows

next:

In this case, the segment 1-2 on the physical boundary is a line of constant , while the

segment 1-4 is a line of constant . Thus at point 1 we have the following coordinate line

configuration:

The lines through point 1 are as follows:

so that the angle between the two coordinate lines is at point 1, and consequently the

Jacobian of the transformation (the cell area, cf. Chapter III) will vanish at this point. The

coordinate species thus changes on the physical boundary at point 1. (Difference

representations at such special points as this, and others to appear in the following examples,

are discussed in Chapter IV.) Since the species of curvilinear coordinate necessarily changes

at a corner on the transformed region boundary, the identification of a concave corner on the

transformed region boundary with a point on a smooth physical boundary will always result

in a special point of the type illustrated here. (A point of slope discontinuity on the physical

boundary also requires special treatment in difference solutions, since no normal can be

defined thereon. This, however, is inherent in the nature of the physical boundary and is not

related to the construction of the transformed configuration.)

Some slightly more complicated examples of the alternatives introduced above now

follow:

Still another alternative in this case would be to collapse the intrusion 2-3-4-5 to a slit in the

transformed region:

Here the physical cartesian coordinates are specified and are double-valued on the vertical

slit, 2-9-5, in the transformed region. The cartesian coordinates of the desired points on the

physical boundary 2-9 are to be used on the slit in the generation of the grid to the left of the

slit in the transformed region, while those on the physical boundary 5-9 are used for

generation to the right of the slit. Solution values in a numerical solution on such a

coordinate system would also be double-valued on the slit, of course. This

double-valuedness requires extra bookkeeping in the code, since two values of each of the

cartesian coordinates and of the physical solution must be available at the same point in the

transformed region so that difference representations to the left of the slit use the slit values

appropriate to the left side, etc. Difference representations near slits are discussed in Chapter

IV. With the composite grid structure discussed in Section 4, however, this need for

double-valuedness, and the concomitant coding complexity, with the slit configuration can

be avoided.

The point 9 here requires special treatment, since the coordinate line configuration

there is as follows:

The coordinate lines through point 9 are as follows:

Here the slope of the coordinate line on which varies is discontinuous at point 9, and the

line on which varies splits at this point. Such a special point will always occur at the slit

ends with the slit configuration.

B. Multiply-connected Regions

With obstacles in the interior of the field, i.e., with interior boundaries, there are still

more alternative configurations of the transformed region. One possibility is to maintain the

connectivity of the transformed region the same as that of the physical region, as in the

following examples showing two variations of this approach using interior slabs and slits,

respectively, in the transformed region. The slab configuration is as follows:

In coding, points inside the slab in the transformed region are simply skipped in all

computations.

This configuration introduces a special point of the following form at each of the

points corresponding to the slab corners in the transformed field:

The coordinate lines through

point 7 are shown below:

This type of special point, where the coordinate species changes on a smooth line, occurs

when a convex corner in the transformed field is identified with a point on a smooth contour

in the physical field. Both coordinate lines experience slope discontinuities at this point.

The slit configuration is as shown below:

(An obvious varition would be to have the slit vertical.) In this slit configuration, the point 5

and 6 are special points of the form shown on p. 26 characteristic of the slit configuration,

and will require special treatment in difference solutions.

The transformed region could, however, be made simply-connected by introducing a

branch cut in the physical region as illustrated below:

Conceptually this can be viewed as an opening of the field at the out and then a deformation

into a rectangle:

Here the coincident coordinate lines 1-2 and 4-3 form a branch cut, which becomes

re-entrant boundaries on the left and right sides of the transformed region. All derivatives are

continuous across this cut, and points at a horizontal distance outside the right-side boundary

in the transformed region are the same as corresponding points at the same horizontal

distance on the same horizontal line inside the left-side boundary, and vice versa. (In all

discussions of point correspondence across cuts, "distance" means distance in the

transformed region). In coding, the use of a layer of points outside each member of a pair of

re-entrant boundaries in the transformed region holding values corresponding to the

appropriate points inside the other boundary of the pair avoids the need for conditional

choices in difference representations, as discussed in Section 6 of this chapter.

Boundary values are not specified on the cut. (This cut is, of course, analogous to the

coincident 0 and 2 lines in the cylindrical coordinate system discussed above.) At the cut

we have the following coordinate line configuration, as may be seen from the conceptional

deformation to a rectangle:

so that the coordinate species and directions are both continuous across the cut.

This type of configuration is often called an O-type. Another possible configuration is

as shown below, often called a C-type:

Opening the field at the cut we have, conceptually,

with 1-2-3-4 to flatten to the bottom of the rectangle. Here the two members of the pair of

segments forming the branch cut are both on the same side of the transformed region, and

consequently points located at a vertical distance below the segment 1-2, at a horizontal

distance to the left of point 2, coincide with points at the same vertical distance above the

segment 4-3, at the same horizontal distance to the right of point 3. The point 2(3) is a

special point of the type shown on p. 26 for slit configurations.

The coordinate line configuration at the cut in this configuration is as follows:

where it is indicated that varies to the right on the upper side of the cut, but to the left on

the lower side. The direction of variation of also reverses at the cut, so that although the

species and slope of both lines are continuous across the cut, the direction of variation

reverses there.

It is possible to pass onto a different sheet across a branch cut, and discontinuities in

coordinate line species and/or direction occur only when passage is made onto a different

sheet. It is also possible, however, to remain on the same (overlapping) sheet as the cut is

crossed, in which case the species and direction are continuous, and this must be the

interpretation when derivatives are evaluated across the cut, as is discussed in Section 5 to

follow. These concepts are illustrated in the following figure, corresponding to the C-type

configuration given on p. 30:

In the present discussion of configurations, the behavior of the coordinate lines across the cut

will always be described in regard to the passage onto a different sheet, since this is in fact

the case in codes. It is to be understood that complete continuity can always be maintained

by conceptually remaining on the same sheet as the cut is crossed. Much of this complexity

can, however, be avoided with the use of an extra layer of points surrounding the

transformed region as will be discussed in Section 6.

Although in principle any region can be transformed into an empty rectangular block

through the use of branch cuts, the resulting grid point distribution may not necessarily be

reasonable in all of the region. Furthermore, an unreasonable amount of effort may be

required to properly segment the boundary surfaces and to devise an appropriate point

distribution thereon for such a transformation. Some configurations are better treated with a

computational field that has slits or rectangular slabs in it.

Regions of higher connectivity than those shown above are treated in a similar

manner. The level of connectivity may be maintained as in the following illustration:

Here one slit is made horizontal and one vertical just for generality of illustration. Both

could, of course, be of the same orientation. Slabs, rather than slits, could also have been

used. The example has three bodies.

With the transformed region made simply-connected we have, using two branch cuts,

a configuration related to the 0-type shown above for one internal boundary:

The conceptual opening here is as follows:

with segment 2-3-4-5-6-7 opening to the bottom. Here the pairs of segments (1-2,8-7) and

(3-4,6-5) are the branch cuts, which form re-entrant boundaries in the transformed region as

shown. In this case, points outside the right side of the transformed region coincide with

points inside the left side, and vice versa. This cut is of the form described on p. 30, where

both the coordinate species and direction are continuous across the cut. Points below the

bottom segment 3-4, to the left of point 4, coincide with points above the bottom segment

6-5 to the right of point 5. This cut is of the form discussed on p. 31, for which the

coordinate species is continuous across the cut but the direction changes there. There are a

number of other possibilities for placement of the two cuts on the boundary of the

transformed region, of course, some examples of which follow.

It is not necessary to reduce the connectivity of the region completely; rather, a slit or slab

can be used for some of the interior boundaries, while others are placed on the exterior

boundary of the transformed region.

One other possibility in two dimensions is the use of a preliminary analytical

transformation of infinity to a point inside some interior boundary, with the coodinates

resulting therefrom replacing the cartesian coordinates in the physical region. The grid

generation then operates from these transformed coordinates rather than from the cartesian

coordinates. This typically gives a fine grid near the bodies, but may give excessively large

spacing away from the body.

Thus, for example, if points on the two physical boundaries shown below

are transformed according to the complex transformation

z’ = 1/z

where z = x+iy and z’ = x’+iy’, infinity in the x,y system will transform to the origin in the

x’, y’ system, as shown below.

Then with the grid generated numerically from the x’, y’ system the following configuration

results:

References to the use of this approach are made in the survey of Ref. [1]. Somewhat related

to this are various two-dimensional configurations which arise directly from conformal

mapping, cf. Ref. [6] and the survey of Ives on this subject, Ref. [7]. (Conformal mapping is

discussed in Chapter X.)

C. Embedded Regions

In more complicated configurations, one type of coordinate system can be embedded

in another. A simple example of this is shown below, where an 0-type system surrounding an

internal boundary is embedded in a system of a more rectangular form, using what amounts

to a slit configuration.

The conceptual opening of this system is best understood in stages: First considering

only the embedded 0-type system surrounding the interior boundary, we have the region

inside the contour 12-13-6-9 opening as follows:

This then opens to the rectangular central portion of the transformed region shown above,

with the inner boundary contour 8-7-8 collapsing to a slit. The rest of the physical region

then opens as shown below:

These two regions then deform to rectangles and are fitted to the top and bottom of the

rectangle corresponding to the inner system along the contours 12-13 and 9-6 as shown.

Here points at a vertical distance below the segment 11-12 are coincident with points

at the same vertical distance below the segment 10-9 on the same vertical line, and vice

versa, with similar correspondence for the pair of segments 13-14 and 6-5. Points at a

horizontal distance to the left of the segment B-12, at a vertical distance above point 8,

coincide with points at the same horizontal distance to the right of the segment 8-9, at the

same vertical distance below point 8. Similar correspondence holds for the pair 7-13 and 7-6.

Boundary values are specified on the slit 8-7.

The composite system shown on p. 40 can also be represented as a slit configuration in

the transformed region:

with the inner system represented as

and the lower side of the slit considered re-entrant with the left half of the top boundary of

the rectangle corresponding to the inner system, the upper side of the slit being re-entrant

with the right half of this top boundary of the inner region. Now the conceptual opening is as

follows for the inner region:

Difference representations made above the slit thus would use points below the right half of

the top of the inner region in the transformed region, etc. Similarly, representation made

below the left half of the top of the inner region would use points below the slit. The slit is

thus a "black hole" into which coordinate lines from the outer system disappear, to reappear

as part of the inner system. The slit here, matched with the top of the inner system, is then

clearly a branch cut, and passage through the slit onto the inner system is simply passage

onto a different sheet.

Note that the embedded system has its own distinctive species and directions for the

coordinate lines, entirely separate from the outer system. Thus for the inner region the

directions are as follows:

while for the outer region they are as shown below:

Thus at a point on the upper interface, 12-13, between the systems the lines are as follows:

while on the lower interface, 9-6 they are as follows:

Thus both coordinates reverse direction at the lower interface although the species is

continuous, while both the species and directions are continuous across the upper interface.

This again corresponds to passage onto a different sheet, for the interface between the inner

and outer systems, i.e., the segments 12-13 and 9-6, is actually a branch cut.

The points 9(12) and 6(13) here require special notice. For example, at point 9 the

coordinate line configuration is as follows:

The lines through point 9 are as shown below

There are thus several changes in species and direction at this point. This type of special

point embodies the form which always occurs with the slit configuration, shown on p. 26,

and occurs here because the embedded region inside the contour 9-6-13-12 is essentially

contained inside a slit defined by the same set of numbers.

The above discussion refers to the slit configuration on p. 41. For the configuration on

p. 40, the lines in the outer region are still as diagrammed on p. 43, but the lines in the inner

region now are as follows:

The coordinate line species and direction given on p. 43 for the upper interface, 12-13, thus

applies here on the entire interface between the two regions.

An alternative treatment of the two special points is to place them inside cells as

shown below:

This results in a six-sided cell surrounding each of these two points which requires special

treatment as discussed in Chapter IV.

Embedded systems can also be constructed in the block configuration:

Here the top of the block, 7-8, in the outer system is re-entrant with the corresponding

segment, 7-8, on a portion of the top of the inner system. The left side of the block, 6-7, and

the bottom of the block, 6-5, are similarly re-entrant with single portions of the top of the

inner system. Finally, the right side of the block, 5-12-8, is re-entrant with two portions, 5-12

and 12-8, of the top of the inner system. Points outside one of these segments in one system

are thus located at corresponding positions inside the other segment of the re-entrant pair in

the other system. The slab sides, matched with the top of the inner system, are thus branch

cuts between the inner and outer systems.

Here the coordinate lines proceed as follows for the outer system:

while those for the inner system are the same as before, as shown on p. 42. This means that

on the left and right sides of the block, i.e., segments 6-7 and 5-8, the line directions are as

follows:

and on the top and bottom, segments 7-8 and 6-5, the directions are as shown below:

There are thus changes in coordinate species and/or direction that are different on each side

of the block.

The point 8 (and points 7,6 and 5) are special points of the following form:

The lines through the point are as shown below:

Here the special points occur in the field instead of on the boundary.

An example of a C-type system embedded in another C-type system is given next:

Here the conceptual opening is as follows: First, considering the system about the upper

body, we have the following configuration:

which, with the body collapsed to a slit, opens to the rectangle in the center of the

transformed region. Next consider the system about the other body:

This opens to a rectangle, with the body flattening to a portion of the bottom, which is fitted

to the first rectangle along the segment 11-13. Finally, the outermost portion opens as

follows:

which opens to a rectangle which is fitted to the first one along the segment 12-14.

Again the embedded region inside the contour 14-12-11-13 can be considered to lie

inside a slit. This contour, which forms the interface between the inner and outer systems, is

actually a branch cut between the two systems, across which there are discontinuities in

coordinate species and directon in the same manner as was discussed above for the previous

embedded system. Points below segment 16-12 coincide with points below segment 17-11 in

this case. Points to the left of segment 15-12, above point 15, are coincident with points to

the right of segment 15-11 below point 15. The slit here is formed of the segments 8-15 and

9-15. The coincident points 11 and 12 here must be taken as a point boundary in the physical

region, i.e., fixed at a specified value. Several special points of the types discussed above are

present here.

An alternative arrangement of the transformed region that corresponds to exactly the

same coordinate system in the physical region is as follows:

Here points below segment 3-4, to the left of point 4, coincide with points above segment

6-5, to the right of point 5. When calculations are made on or above the segment 12-14 on

the larger block, points below this segment coincide with points below the corresponding

segment on the smaller block. Similarly, when calculations are made on or below the

segment 13-11 on the larger block points above this segment coincide with points below the

corresponding segment on the smaller block. Finally, points below the segment 7-8, to the

left of point 8, on the smaller block are coincident with points above the segment 10-9, to the

right of point 9.

This configuration displays explicitly the correspondence of the embedded region

inside the contour 14-12-11-13 to a slit. Conceptually, coordinate lines from the main system

disappear into the slit and emerge into the embedded system. These coordinate lines thus are

continued from the main system onto another sheet representing the embedded system. This

concept of embedded systems, with continuation onto another sheet through a slit adds

considerable flexibility to the grid configurations and is of particular importance with

multiple boundaries and in three dimensions. The composite structure discussed in Section 4

removes much of the coding complexity associated with systems of this type.

D. Other Configurations

Another arrangement of cuts, where the species of coordinate changes on a continuous

line as the cut is crossed, is illustrated below. The transformed region in this case is

composed of three blocks connected by the cuts.

Here points outside one section are coincident with corresponding points inside the adjacent

section.

The coordinate line configuration on the interface on the right side of block A here is

as follows:

This same type of configuration occurs, in different orientations, on each of the interfaces.

These interfaces are branch cuts, so that passage onto the adjacent block amounts to passage

onto another sheet in the same manner discussed above.

As a final configuration for consideration in two dimensions, the following example

shows a case with fewer lines on one side of a slab than on the other. This does not

necessitate the use of different increments of the curvilinear coordinates in the numerical

expressions, because, as has been mentioned, these increments always cancel out anyway.

E. Three-dimensional Regions

All the general concepts illustrated in these examples extend directly to three

dimensions. Interior boundaries in the transformed region can become rectangular solids and

plates, corresponding to the slabs and slits, respectively, illustrated above for two

dimensions. Examples of three-dimensional configurations can be found in the surveys given

by Ref. [8] and [9].

It is also possible to use branch cuts, as illustrated above for two dimensions, to bring

the interior boundaries in the physical region entirely to the exterior boundary of the

transformed region:

Physical space

Computational space

The correspondence between the physical and transformed fields can, however, become

much more complicated in three dimensions, and considerable ingenuity may be required to

visualize this correspondence. For instance, the simple case of polar coordinates corresponds

to a rectangular solid with two opposing sides having the radial coordinate constant thereon,

and two re-entrant sides on which the longitudinal coordinate is constant at 0 and 2 ,

respectively (corresponding to the cut). The remaining two sides correspond to the north and

south polar axes, so that an axis opens to cover an entire side. There is thus a line, i.e., the

axis, in the physical region that corresponds to an entire side in the transformed region.

Three-dimensional grids may be constructed in some cases by simply connecting

corresponding points on two-dimensional grids generated on stacks of planes or curved

surfaces:

It should be noted, however, that this procedure provides no inherent smoothness in

the third direction, except in cases where the stack is formed by an analytical transformation,

such as rotation, translation or scaling, of the two-dimensional systems. An example of such

an analytical transformation of two-dimensional systems is the construciton of a

three-dimensional grid for a curved pipe by rotating and translating (and scaling if the

cross-sectional area of the pipe varies) two-dimensional grids generated for the pipe

cross-section so as to place these transformed two-dimensional grids normal to the pipe axis

at successive locations along the axis:

Another example is the rotation of a two-dimensional grid about an axis to produce an

axi-symmetric grid:

4. Composite Grids

All of the above concepts can be incorporated in a single framework, and the coding

complexity can be greatly reduced, by considering the physical field to be segmented into

sub-regions, bounded by four (six in 3D) generally curved sides, within each of which an

individual coordinate system is generated. The overall coordinate system, covering the entire

physical field, is then formed by joining the sub-systems at the sub-region boundaries. The

degree of continuity with which this juncture is made is a design consideration in regard to

the mode of application intended for the resulting grid.

This segmentation concept is illustrated in the figure below.

The locations of the interfaces between the sub-regions in the physical region are, of course,

arbitrary since these interfaces are not actual boundaries. These interfaces might be fixed,

i.e., the location completely specified just as in the case of actual boundaries, or might be left

to be located by the grid generation procedure. Also the coordinate lines in adjacent

sub-regions might be made to meet at the interface between with complete continuity:

---Interface

with some lesser degree of continuity, e.g., continuous line slope only:

---Interface

or with a discontinuity in slope:

or perhaps not to meet at all:

Naturally, progressively more special treatment at the interface will be required in numerical

applications as more degrees of line continuity at the interface are lost. Procedures for

generating segmented grids with various degrees of interface continuity are discussed later,

and conservative interface conditions are given in Ref. [52], [53].

Now, with regard to placing these concepts in the framework of segmentation, the

sides of an individual subregion (called a "block" hereafter) can be treated as boundaries on

which the coordinate points, and/or the coordinate line intersection angles, are specified, just

as is done for actual boundaries, or a side may be treated as one member of a pair of

re-entrant boundaries, i.e., one side of a branch cut in the physical region across which

complete continuity is established. The other member of the pair may be another side (or

portion thereof) of the same block or may be all (or part of) a side of an adjacent block in the

physical field. Recall that it is not necessary for a coordinate to remain of the same species

across a re-entrant boundary, since the passage is onto a different sheet. This can introduce

some coding complexity, but the treatment is straightforward, and in fact the coding can be

greatly simplified by using an extra layer of points surrounding each block as is discussed in

Section 6.

Some of the general concepts have been embodied in the two-dimensional code

discussed in Ref. [19] and in three recent three-dimensional codes, Ref. [13] Ref. [14], and

[51].

A. Simply-connected Regions

The first L-shaped simply-connected configuration on p. 21 can be interpreted as

being composed of three blocks, with the sides of adjacent blocks forming pairs of re-entrant

boundaries:

or two blocks with a portion of a side of one block re-entrant with an entire side of another

block:

Here, and in the examples to follow, solid lines correspond to physical boundaries, while the

dashed lines correspond to the interfaces between the blocks. The dashed arrows indicate the

linkage between the interfaces. (Obviously, any single block can be broken into any number

of blocks connected by re-entrant boundaries across adjacent sides.) In contrast, the

L-shaped configuration on p. 22 corresponds to the use of a single block. Similarly, the

configuration on p. 24 can be formed with three blocks:

while the first configuration for the same boundary on p. 25 is formed with a single block.

The slit configuration on p. 25 can be formed of three blocks:

or two blocks with only a portion of the adjacent sides of two blocks forming a re-entrant

boundary:

B. Multiply-connected Regions

The configuration with a single cut shown on p. 29 corresponds to the use of a single

block with the left and right sides here being the members of a pair of re-entrant boundaries:

The multiply-connected slab configuration on p. 27 can be broken into four blocks:

Other decompositions should also be immediately conceivable. The slit configuration on p.

28 can be formed with two blocks, again with only portions of adjacent sides serving as

re-entrant boundaries:

or into four blocks, with entire sides as re-entrant boundaries in all cases:

The double-body region on p. 34 opens to a single block as shown there, with portions

of sides as re-entrant boundaries. A five block configuration would use only entire sides as

re-entrant boundaries, however:

There is no real advantage, however, to the five-block system here.

C. Embedded Regions

The segmentation concept is most useful in the construction of embedded coordinate

systems. For instance, the system on p. 40 can be considered to be formed of three blocks as

follows:

Here portions of adjacent sides of the two larger blocks are re-entrant with each other, while

each of the remaining portions of these sides is re-entrant with half of one side of the smaller

block. The left and right sides of the smaller block are re-entrant with each other. This

configuration could also have been constructed with eight blocks:

with only entire sides being involved in re-entrant pairs as shown.

With embedded systems the coordinate species often changes as the re-entrant

boundary is crossed. These systems also show that the blocks need be physically adjacent

only in the physical field, and it is in this sense that "adjacent" is always to be interpreted.

The transformed (computational) field should always be viewed as only a bookkeeping

structure. Various constructions are possible for the configurations on p. 48 and 50, and a

two block structure was actually used on p. 50. A further example follows:

D. Three-dimensional Regions

For general three-dimensional configurations, it is usually very difficult to obtain a

reasonable grid with the entire physical region transformed to a single rectangular block. A

better approach in most cases is to segment the physical region into contiguous sub-regions,

each bounded by six curved surfaces, with each sub-region being transformed into a

rectangular block. An individual grid is generated in each sub-region:

These sub-region grids are patched together to form the overall grads, as in the

two-dimensional cases discussed above. Examples of the use of this segmentation in three

dimensions are found, in particular, in Ref. [11] and [12]. Others are noted in the survey

given by Ref. [9].

As noted above, complete continuity can be achieved at the sub-region interfaces by

noting the correspondence of points exterior to one sub-region with points interior to

another. The necessary bookkeeping can be accomplished, and the coding complexity can be

greatly reduced, by using an auxiliary layer of points just outside each of the six sides of the

computational region, analogous to the procedure mentioned above for two dimensions. A

correspondence is then established in the code between the auxiliary points and the

appropriate points just inside other sub-regions. This approach has recently been

incorporated in an internal region code, Ref. [13], and in two codes for general regions, Ref.

[14] and [51]. This is discussed in more detail in Section 6.

General three-dimensional regions can be built up using sub-regions as follows: First,

point distributions are specified on the edges of a curved surface forming one boundary of a

sub-region:

and a two-dimensional coordinate system is generated on the surface:

When this has been done for all surfaces bounding the sub-region, the three-dimensional

system within the sub-region is generated using the points on the surface grids as boundary

conditions:

In three dimensions it is possible for a line, e.g., a polar axis, in the physical region to

map to an entire side of the computational region as in the illistration below, where the axis

corresponds to the entire left side of the block:

The system illustrated here could be one of several identical blocks joining together to form

a complete system around the axis.

It is illustrated by an exercise that the occurence of a polar axis can be avoided, and

this facilitates the construction of a block structure. Thus a surface grid, having eight

"corners", analogous to the four "corners" on the circle in the 2D grid on p. 23, can be

constructed on the surface of a sphere. This serves much better than a latitude-longitude type

system for joining to adjacent regions. Similarly, the use of the four "corner" system, rather

than a cylindrical system, in a circular pipe allows T-sections and bifurcations to be treated

easily by a composite structure, c.f. Ref. [13].

Generally, grid configurations with polar axis should not be used in composite grid

structures.

E. Overlaid Grids

Another approach to complicated configurations is to overlay coordinate systems of

different types, or those generated for different sub-regions:

Here an appropriate grid is generated to fit each individual component of the configuration,

such that each grid has several lines of overlap with an adjacent grid. Interpolation is then

used in the region of overlap when solutions are done on the composite grid, with iteration

among the various grids. This approach has the advantage of simplicity in the grid

generation, in that the various sub-region grids are only required to overlap, not to fit.

However, there would appear to be problems if regions of strong gradients fall on the

overlap regions. Also the interpolation may have to be constructed differently for different

configurations, so that a general code may be hard to produce. Some applications of such

overlaid grids are noted in Ref. [5].

5. Branch Cuts

As has been noted in the above discussion of transformed field configurations, it is

possible for discontinuities in coordinate species and/or direction to occur at branch cuts, in

the sense of passage onto another sheet. Continuity can be maintained, however, by

conceptually remaining on the same overlapping sheet as the cut is crossed. All derivatives

thus do exist at the cut, but careful attention to difference formulations is necessary to

represent derivatives correctly across the cut. Although the correct representation can be

accomplished directly by surrounding the computational region with an extra layer of points,

as is discussed in Section 6, it is instructive to consider what is required of a correct

representation further here.

A. Point Correspondence

Points on re-entrant boundaries in the transformed region, i.e., on branch cuts in the

physical region, are not special points in the sense used above. Points on re-entrant

boundaries, in fact, differ no more from the other field points than do the points on the 0 and

2 lines in a cylindrical coordinate system. Care must be taken, however, to identify the

interior points coinciding with the extensions from such points beyond the field in the

transformed space. This correspondence was noted above in each of the configurations

shown above, being indicated by the dashed connecting lines joining the two members of a

pair of re-entrant boundaries. There are essentially four types of pairs of re-entrant

boundaries, as illustrated in the following discussion of derivative correspondence. In these

illustrations one exterior point, and its corresponding interior point, are shown for each case.

The converse of the correspondence should be evident in each configuration.

For the configurations involving a change in the coordinate species at the cut, not only

must the coordinate directions be taken into account as the cut is crossed, but also the

coordinate species may need to be interpreted differently from that established across the cut

in order to remain on the same sheet as the cut is crossed. For example, points on an -line

belonging to section A in the figure on p. 52, but located outside the right side of this region,

are coincident with points on a -line of region B at a corresponding distance (in the

transformed region) below the top of this region.

B. Derivative Correspondence

Care must be taken at branch cuts to represent derivatives correctly in relation to the

particular side of the cut on which the derivative is to be used. The existence of branch cuts

indicates that the transformed region is multi-sheeted, and computations must remain on the

same sheet as the cut is crossed. Remaining on the same sheet means continuing the

coordinate lines across the cut coincident with those of the adjacent region, but keeping the

same interpretation of coordinate line species and directions as the cut is crossed, rather than

adopting those of the adjacent region. As noted above, points outside a region across a cut

the transformed space are coincident with points inside the region across the other member

of the pair of re-entrant boundary segments corresponding to the cut in the transformed

space. The positive directions of the curvilinear coordinates to be used at these points inside

the region across the other member of the pair in some cases are the same as the defined

directions there, but in other cases are the opposite directions. As noted above, the

coordinate species may change also.

For cuts located on opposing sides of the transformed region, the proper form is

simply a continuation across the cut. Thus in the configuration on p. 29, with a computation

site on the right side of the transformed region, i.e., on the upper side of the cut in the

physical plane, we have points to the right of the site (below the cut in the physical plane)

coinciding with points to the right of the left side of the transformed region (below the cut in

the physical plane) as noted above. When -derivatives and -derivatives for use outside the

right side of the transformed region are represented inside the left side, the positive

directions of and to be used there are to the right and upward, respectively, as is

illustrated below. (In this and the following figures of the section, the dotted arrows indicate

the proper directions to be used at the interior points coincident with the required exterior

points, i.e., on the same sheet across the cut, while solid arrows indicate the locally

established directions for the coordinate lines, i.e., on a different sheet.)

With the two sides of the cut both located on the same coordinate line, i.e., on the

same side of the transformed region as in the configuration on p. 30, however, the situation

is not as simple as the above. In this case, when the computation site is on the left branch of

the cut in the transformed region (on the lower branch in the physical region), the points

below this boundary in the transformed region coincide with points located above the right

branch of the cut (above the cut in the physical region) at mirror-image positions, as has

been noted earlier. The -derivatives for use at such points below the left branch thus must

be represented at these corresponding points above the right branch. The positive direction

of for purposes of this calculation of derivatives above the right branch, for use below the

left branch, must be taken as downward, not upward. There is a similar reversal in the

interpretation of the positive direction of . This is in accordance with the discussion on p.

31. These interpretations are illustrated below:

In the configuration on p. 40, where two sides of a cut face each other across a void,

there is really no problem of interpretation, since the directions in the configuration are

treated simply as if the void did not exist. This correspondence is as shown below:

In all cases the interpretation of the positive directions of the curvilinear coordinates

must be such as to preserve the direction in the physical region, i.e., on the same sheet, as the

cut is crossed. In the cases where the coordinate species change at the cut, the situation is

even more complicated. Thus on the left side, segment 6-7, of the slab interface between the

inner and outer systems in the embedded configuration on p. 46, where the species changes

across the cut, the correspondence is as follows:

Thus, when a -derivative is needed outside the outer sytem, for use inside the left slab

interface, the positive -direction at the corresponding points inside the inner system must

be taken to coincide with the negative -direction of the inner system. Similarly, an

-derivative would be represented taking the positive -direction to coincide with the

positive -direction of the inner system. In an analogous fashion, a -derivative needed

outside the inner system, for use inside the segment 6-7, would be represented at the

corresponding point inside the outer system, i.e., to the left of the left slab side, but with the

positive -direction taken to be the positive -direction of the outer system. An -derivative

would be represented similarly, taking the positive -direction to be the negative -direction

of the outer system.

A -derivative to the left of the right side of the slab in the outer system would be

represented below segment 12-5 or 8-12, as the case may be, but with the positive

-direction taken to be the positive -direction of the inner system. Similarly, an

-derivative would be represented taking the positive -direction to be the negative

direction of the inner system. For a -derivative above the bottom of the slab in the outer

system, the correspondence is to below the segment 5-6 inside the inner system, with the

positive -direction taken to be the negative -direction of the inner system. The

-derivative is represented taking the positive -direction to be the negative -direction of

the inner sytem. Finally, for derivatives below the top of the slab in the outer system, the

correspondence is to below the segment 7-8 inside the inner system, with both the species

and direction of the coordinates unchanged.

The proper interpretation of coordinate species and direction across branch cuts for all

the other configurations discussed above can be inferred directly from these examples. A

conceptual joining of the two members of a pair of re-entrant boundaries in accordance with

the dashed line notation used on the configurations given in this chapter will always show

exactly how to interpret both the coordinate species and directions in order to remain on the

same sheet and thus to maintain continuity in derivative representation across the cut.

Examples of the proper difference representation are given in the following section. The

complexities of this correspondence can be completely avoided, however, by using

surrounding layers around each block in a segmented structure as discussed in the next

section.

6. Implementation

As discussed above, the transformed region is always comprised of contiguous

rectangular blocks by construction. This occurs because of the essential fact that one of the

curvilinear coordinates is defined as constant on each segment of the physical boundary.

Consequently, each segment of the physical boundary corresponds to a plane segment of the

boundary of the transformed region that is parallel to a coordinate plane there. The complete

boundary of the transformed region then is composed of plane segments, all intersecting at

right angles. Although the transformed region may not be a simple six-sided rectangular

solid, it can be broken up into a contiguous collection of such solids, here called blocks.

i cancel from all difference

Now it is noted in Chapter III that the increments

expressions, and that the actual values of the curvilinear coordinates i are immaterial. The

coordinates in the transformed region can thus be considered simple counters identifying the

points on the grid. This being the case, and the transformed region being comprised of a

collection of rectangular blocks, it is convenient to identify the grid points with integer

values of the curvilinear coordinates in each block, and thus to place the cartesian

coordinates of a grid point in ijk, where the subscripts (i,j,k) here indicate position ( i, 2,

3)

in the transformed region. (In coding, a fourth index may be added to identify the

block.) In each block, the curvilinear coordinates are then taken to vary as i = 1,2,...,Ii over

the grid points, where Ii is the number of points in the i-direction. Grid points on a

boundary segment of the transformed region will be placed in ijk with one index fixed.

Now each block has six exterior boundaries, and may also have any number of interior

boundaries (cf. the slab and slit configurations of Section 3), all of which will always be

plane segments intersecting at right angles, although the occur ence of interior boundaries

can be avoided if desired by breaking the block up into a collection of smaller blocks as

discussed in Section 4. The boundary segments in the transformed plane may correspond to

actual segments of the physical boundary, or may correspond to cuts in the physical region.

As discussed in Section 5, these cuts are not physical boundaries, but rather are interfaces

across which the field is re-entrant on itself. A boundary segment in the transformed region

corresponding to such a cut then is an interface across which one block is connected with

complete continuity to another block, or to another side of itself, several examples having

been given above in this chapter.

Depending on the type of grid generation system used (cf., the later chapters), the

cartesian coordinates of the grid points on a physical boundary segment may either be

specified or may be free to move over the boundary in order to satisfy a condition, e.g.,

orthogonality, or the angle at which coordinate lines intersect the boundary.

To set up the configuration of the transformed region, a correspondence is established

between each (exterior or interior) segment of the boundary of the transformed region and

either a segment of the physical boundary or a segment of a cut in the physical region. This

is best illustrated by a series of examples using the configurations of this chapter. The first

step in general is to position points on the physical boundary, or on a cut, which are to

correspond to corners of the transformed region (exterior or interior). As noted in Section 3,

these points do not have to be located at actual corners (slope discontinuities) on the physical

boundary.

For example, considering the two-dimensional simply-connected region on p. 23, four

points on the physical boundary are selected to correspond to the four corners of the empty

rectangle that forms the transformed region here:

Now, considering any one of these four points, one species of curvilinear coordinate will run

from that point to one of the two neighboring corner points, while the other species will run

to the other neighbor:

The corresponding species of coordinates will run to connect opposite pairs of corner points:

Since the curvilinear coordinates are to be assigned integer values at the grid points, i is to

vary from 1 at one corner to a maximum value, Ii, at the next corner, where Ii is the number

of grid points on the boundary segment between these two corners. Thus, proceeding

clockwise from the lower left corner, the cartesian coordinates of the four corner points are

placed in 1,1, 1,J, I,J, and I,1, where I1 = I and I2 = J.

The boundary specification is then completed by positioning I-2 points on the lower and

upper boundary segments of the physical region as desired, and J-2 points on the left and

right segments. The cartesian coordinates of these points on the lower and upper segments

are placed in i,1 and i,J, respectively, for i from 2 to I-1, and those on the left and right

segments are placed, respectively, in

1,j and

I,j for

j from 2 to J-1.

This process of boundary specification can be most easily understood by viewing the

rectangular boundary of the transformed region, with I equally-spaced points along two

opposite sides and J equally-spaced points along the other two sides, conceptually, as being

deformed to fit on the physical boundary. The corners can be located anywhere on the

physical boundary, of course. Here the point distribution on the sides can be conceptually

stretched and compressed to position points as desired along the physical boundary. The

cartesian coordinates of all the selected point locations on the physical boundary are then

placed in as described above.

This conceptual deformation of the rectangular boundary of the transformed region to

fit on the physical boundary serves to quickly illustrate the boundary specification for the

doubly-connected physical field shown on p. 29, which involes a cut. Thus I points are

positioned as desired clockwise around the inner boundary of the physical region from 2 to

3, and I points are positioned as desired, also in clockwise progression, around the outer

boundary from 1 to 4. The cartesian coordinates of these points on the inner boundary are

placed in i,1, and those on the outer boundary in i,J, with i from 1 to I. Note that here the

first and last points must coincide on each boundary, i.e., I,1 = 1,1 and I,J = 1,J. The

left and right sides of the transformed region (i=1 and i= I) are re-entrant boundaries,

corresponding to the cut, and hence values on these boundaries are not set but will be

determined by the generation system. The system must provide that the same value appears

on both of these sides, i.e., I,j = 1,j for all j from 2 to J-1.

The conceptual deformation of the rectangle for a C-type configuration is illustrated

on p. 31. Here, with I1 the number of points on the segments 1-2 and 3-4 (which must have

the same number of points), I2 points are positioned as desired around the inner boundary in

the physical region in a clockwise sense from 2 to 3, and the cartesian coordinates of these

points are placed in i,1 for i from I1 to I1+I2-1.

The first and last of these points must be coincident, i.e. I1,1 = I1 + I2-1,1. Now the

top, and the left and right sides, of the rectangle are deformed here to fit on the outer

boundary of the physical region. (In the illustration given, the two top corners are placed on

the two corners that occur in the physical boundary, a selection that is logical but not

mandatory.) The cartesian coordinates of the J points(positioned as desired on the segment

4-5 of the physical boundary) are placed in I,j, proceeding upward on the physical

boundary from 4 to 5 for j= 1 to J, and those on the segment 1-6 are placed in 1,j, but

proceeding downward on the physical boundary from 1 to 6 for j1 to J. Finally, the cartesian

coordinates of the I selected points on the physical boundary segment 6-5 are placed in i,J,

proceeding clockwise from 6 to 5 for i=1 to I. Since the same number of points must occur

on the top and bottom of the rectangle, we must have I=2(I1-1)+I2. Here the portions of the

lower side of the rectangle, i.e., i from 2 to I1-1, and from I1+I2 to I-1 with j=1, are

re-entrant boundaries corresponding to the cut, and hence no values are to be specified on

these segments. The generation system must make the correspondence i,1 = I-i+1,1 for

i=2 to I1-1 on these segments.

The conceptual deformation of the boundary of the transformed regions also serves for

the slab configuration on p. 27, where the interior rectangle deforms to fit the interior

physical boundary, while the outer rectangle deforms to fit the outer physical boundary. On

the inner boundary, the cartesian coordinates of J2-J1+1 selected points on the segment 5-8

of the physical boundary are placed in I1,j for j from J1 to J2, proceeding upward on the

physical boundary from 5 to 8, where J1 and J2 are the j-indices of the lower and upper

sides, respectively, of the interior rectangle and I1 is the i-index of the left side of this

rectangle. Similarly, J2-J1+1 points are positioned as desired on the segment 6-7 of the

physical boundary and are placed in I2,j, where I2 is the i-index of the right side of the

inner rectangle. Also I2-I1+1 points on the segments 5-6 and 8-7 of the physical boundary

are placed in i,J1 and i,J2, respectively, for i from I1 to I2, proceeding to the right on each

segment. The outer boundary is treated as has been described for an empty rectangle. Here

there will be J1-1 coordinate lines running from left to right below the inner boundary, and

J-J2 lines running above the inner boundary. Similarly, there will be I1-1 lines running

upward to the left of the interior boundary and I-I2 lines to the right. Thus the specifications

of the desired number of coordinate lines running on each side of the inner boundary serves

to determine the indicies I1, I2, J1, and J2. Note that the points inside the slab, i.e., I1 < i <

I2 and J1 < j < J2 are simply excluded from the calculation.

The slit configuration, illustrated on p. 28, can also be treated via the conceptual

deformation, but now with a portion of a line inside the rectangle opening to fit the interior

boundary of the physical region. This requires that provision be made in coding for two

values of the cartesian coordinates to be stored on the slit. If the i-indices of the slit ends, 5

and 6, are I1 and I2, respectively, then the cartesian coordinates of I2-I1+1 points positioned

as desired on the lower portion of the physical interior boundary, again proceeding from 5 to

6, are placed in a one-dimensional array, while the coordinates of the same number of points

selected on the upper portion of the physical interior boundary, again proceeding from 5 to 6,

are placed in another one-dimensional array. The first and last points in one of these arrays

must, of course, coincide with those in the other. Then the generation system must read

values into i,J1 for i from I1 to I2 (J1 being the j-index of the slit) from the former array for

use below the slit, or values from the latter array for use above. (As has been noted, the use

of a composite structure eliminates the need for these two auxiliary arrays.) Note that the

index values I1 and I2 are determined by the number of lines desired to run upward to the

left and right of the interior boundary, respectively, i.e., I1-1 lines on the left and I-I2 on the

right. Similarly, there will be J1-1 lines below the interior boundary, and J-J1 above.

Configurations, such as those illustrated on pp. 24-25, which involve slabs or slits that

intersect the outer boundary are treated similarly, with points inside the slab again being

simply excluded from the calculations. Also multiple slab or slit arrangements are treated by

obvious extensions of the above procedures. Here the indices corresponding to each slab or

slit will be determined by the number of points on the interior boundary segments and the

number of coordinate lines specified to run between the various boundaries. For example, in

the slit configuration shown on p. 33, the ends of the horizontal slit would be at i-indices I1

and I2, where I1-1 lines run vertically to the left of the slit and there are I2-I1+1 points on

the slit. The vertical slit would be at i=I3 where there are I3-I2-1 vertical lines between this

slit and the horizontal slit (and I-I2 lines to the right). Similarly, if the j-indices of the ends of

the vertical slit or J1 and J2, there will be J1-1 horizontal lines below this slit and J-J2 lines

above. With the j-index of the horizontal slit as J3, there will be J3-1 horizontal lines below

this slit and J-J3 above. Provision will now have to be made in coding for two

one-dimensional arrays for each slit to hold the cartesian coordinates of the points on the

segments of the physical interior boundaries corresponding to the two sides of each slit.

Again this coding complexity is avoided in the composite structure.

The use of the conceptual deformation of the rectangle to setup the boundary

configuration for the case with multiple interior boundaries on p. 34 should follow with little

further explanation. Here there must be the same number of points on the pair of segments

2-3 and 6-7, which correspond to the two segments forming the interior boundary on the

right. There must also be the same number of points on the pair, 3-4 and 5-6, corresponding

to the cut connecting the two interior boundaries. Finally the number of points on the outer

boundary must, of course, be the same as that on the bottom boundary. Note also that the

values of the cartesian coordinates placed at 2 must be the same as are placed at 7; those at 3

must be the same as those at 6, and those at 4 the same as at 5. Values are not set on the cuts,

of course, but the generation system must provide that values at points on the segment from

3 to 4 are the same as those on the segment 5-6, but proceeding from 6 to 5. Also values on

the segments 2-1 and 7-8 must be the same, proceeding upward in each case.

Following the conceptual deformation of the rectangular boundaries of the

transformed region and the indexing system illustrated above, it now should be possible to

set up the more complicated configurations such as the embedded regions shown in Section

3C. As noted there, however, the most straightforward and general approach to such more

complicated configurations is to divide the field into contiguous rectangular blocks, each of

which has its own intrinsic set of curvilinear coordinates and hence its own (i,j,k) indexing

system. The necessary correspondence between the individual coordinate systems across the

block interfaces was discussed in some detail in Section 3C. This block structure greatly

simplifies the setup of the configuration. For example, consider the 3-block structure shown

on p. 49 for the physical field shown on p. 48, for which the blocks are as follows:

Here the selected points on the right interior boundary (segment 8-15-9) are placed in i,1 of

the first block, for i from the i-index at 8 to that at 9,proceeding clockwise from 8 to 9 on the

physical boundary. (The difference between these two i-indices here is equal to the number

of points on this interior boundary,less one.) Similary,the selected points on the left interior

boundary (segment 4-5) are placed in i,1 of the second block for i from the i-index at 4 to

that at 5, proceeding clockwise from 4 to 5 on the physical boundary. The selected points on

the outer boundary of the physical region are placed in 1,j of the third block for j from 1 to

J3, in i,J3 for i from 1 to I3, and in I3,j for j from J3 to 1, proceeding from 16 to 1 to 2 to

14 on the physical boundary. Points on the remainder of the physical outer boundary are

placed in 1,j of the second block for j from 1 to J2 and in I2,j for j from J2 to 1,

proceeding from 3 to 17 for the former and from 13 to 6 for the latter, and in

1,j of

the first

block for j from 1 to Jl and in I1,j for j from Jl to 1, proceeding from 7 to 13 for the former

and from 14 to 10 for the latter.

Since the three blocks must fit together we have I3=I2, (I1+1)/2 equal to the difference

in i-indices between 11 and 13 in the second block and to that between 12 and 14 of the third

block. The quantities J1, J2, and J3 determine how many C-type lines occur in each block,

and can be chosen independently. Here the segment 11-13 on the top of the first block

interfaces with the corresponding segment on the top of the second block. The segment

12-14, which forms the remainder of the top of the first block, interfaces with the

corresponding segment on the bottom of the third block. Finally, the segment 12-16, which

forms the remainder of the bottom of the third block, interfaces with segment 11-17, which

forms the remainder of the top of the second block. The segments 3-4 and 6-5 on the bottom

of the second block interface with each other in the order indicated, as do also the segments

7-8 and 10-9 on the bottom of the first block.

In coding, this block structure can be handled by using a fourth index to identify the

block, placing an extra layer around each block, (i=0 and I+1, j=0 and J+1) and providing an

image-point array by which any point of any block can be paired with any point of any other,

or the same, block. Such pairs of points are coincident in the physical region, being on or

across block interfaces, and consequently are to be given the same values of the cartesian

coordinates by the generation system. This imaging extends to the extra layer surrounding

each block, so that appropriate points Inside other blocks can be identified for use in

difference representations on the block interfaces that require points outside the block, (cf.

Section 5).

Interface correspondence then can be established by input by setting the image-point

correspondence on the appropriate block sides, i.e., placing the (i,j,k) indices and block

number of one member of a coincident pair of points in the image-point array at the indices

and block number of the other member of the pair. This correspondence is indicated on the

block diagram on pp. 85-86 by the points enclosed in certain geometric symbols.

Thus, for the 3-block configuration considered above, the indices (I1-i+1, 1) and block

number 1, corresponding to a point on the segment 9-10 of the first block, would be placed

in the image-point array at the point (i,1) on the segment 7-8 of this block, and vice versa. A

similar pairing occurs for points on the segments 3-4 and 5-6 of the second block. The

indices (I2-i+1, J2) and block number 2, corresponding to a point on the segment 11-13 of

the second block would be placed in the image-point array at the point (i,Jl) of the first block

on the segment 13-11 of that block, and vice versa. The indices (I3-I1+i, 1) and block

number 3 (a point on the segment 12-14 of the third block) would be placed in the array at

the point (i,J1) of the first block for a point on the segment 12-14 of that block. Finally, the

indices (i,1) and block number 3 (point on segment 16-12 of the third block) would be

placed in the array at the point (i,J2) of the second block for a point on the segment 17-11 of

that block. The remaining segments all correspond at portions of the physical boundary and

hence do not have image points.

In the same manner the following image correspondence can be set between interior

points and points on the surrounding layers in order to establish difference representations

across the block interfaces: (This correspondence is indicated symbolically on the block

diagram on pp. 85-86 by geometric symbols.):

As noted, all of this information would be input into the image-point array. Then with

values of the cartesian coordinates at the image points on the surrounding layer set equal to

those at the corresponding object point inside one of the blocks, it is possible to use the same

difference representations on the interfaces that are used in the interior.

The discussion given in this chapter should now allow the image-point input to be

constructed for any configuration of interest. As noted, it is not necessary that the coordinate

species remain the same as the interface is crossed. Thus, for instance, a point on the right

side of one block could be paired with one on the bottom of another block. In such a case the

image point of the point (I+1, j) the first block would be the point (j,2) inside the second

block. Similarly the image of the point (i,0) below the second block would be the point (I-1,

i) inside the first block, The correct difference representation across interfaces is thus

automatically established, eliminating the need for the concern with passage onto different

sheets discussed in detail earlier in this chapter.

This greatly simplifies the coding, since with the surrounding layers and the use of the

image points, all of the derivative correspondences are automatic and do not have to be

specified for each configuration. It is only necessary to specify the point correspondence by

input. This construction also allows codes for the numerical solution of partial differential

equations on the grid to be written to operate on rectangular blocks. Then any configuration

can be treated by sweeping over all the blocks. The surrounding layers of points and the

image correspondence provide the proper linkage across the block interfaces. In an implicit

solution the values on the interfaces would have to be updated iteratively in the course of the

solution. The solution for the generation of the grid would similarly keep the interface and

surrounding layer values updated during the course of the iterative solution.

This, of course, maintains completecontinuity across the block interfaces. If complete

continuity is not required,

Чтобы не видеть никакую рекламу на сайте, нужно стать

**VIP**-пользователем.

Это можно сделать совершенно бесплатно. Читайте подробности тут.