23 int i,
j,
k, idim, ngh, ileft;
26 double xiL, xiR, dxi, dvol;
41 for (idim = 0; idim < 3; idim++) {
44 (GXYZ + idim)->
dV =
ARRAY_1D (GXYZ[idim].np_tot,
double);
45 (GXYZ + idim)->
r_1 =
ARRAY_1D (GXYZ[idim].np_tot,
double);
46 (GXYZ + idim)->
ct =
ARRAY_1D (GXYZ[idim].np_tot,
double);
66 for (i = 0; i <= iend; i++) {
73 #if GEOMETRY == CARTESIAN
75 if (i == 0) GG->
A[-1] = 1.0;
78 #elif GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
80 if (i == 0) GG->
A[-1] = fabs(xl);
81 GG->
dV[
i] = fabs(x)*
dx;
82 GG->
xgc[
i] = x + dx*dx/(12.0*
x);
84 #elif GEOMETRY == SPHERICAL
86 if (i == 0) GG->
A[-1] = xl*
xl;
88 GG->
dV[
i] = fabs(xr*xr*xr - xl*xl*xl)/3.0;
89 GG->
xgc[
i] = x + 2.0*x*dx*dx/(12.0*x*x + dx*
dx);
99 for (j = 0; j <= jend; j++) {
105 #if GEOMETRY != SPHERICAL
107 if (j == 0) GG->
A[-1] = 1.0;
112 GG->
A[
j] = fabs(sin(xr));
113 if (j == 0) GG->
A[-1] = fabs(sin(xl));
114 GG->
dV[
j] = fabs(cos(xl) - cos(xr));
116 GG->
xgc[
j] = (sin(xr) - sin(xl) + xl*cos(xl) - xr*cos(xr));
117 GG->
xgc[
j] /= cos(xl) - cos(xr);
118 GG->
ct[
j] = 1.0/tan(x);
127 for (k = 0; k <= kend; k++) {
134 if (k == 0) GG->
A[-1] = 1.0;
145 for (i = 0; i < GXYZ[idim].
np_tot; i++) {
146 GXYZ[idim].
inv_dx[
i] = 1.0/(GXYZ[idim].
dx[
i]);
149 for (i = 0; i < GXYZ[idim].
np_tot-1; i++) {
150 GXYZ[idim].
inv_dxi[
i] = 2.0/(GXYZ[idim].
dx[
i] + GXYZ[idim].
dx[i+1]);
165 return (grid[0].
dx[i]);
178 #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
179 return (grid[1].
dx[j]);
180 #elif GEOMETRY == POLAR || GEOMETRY == SPHERICAL
181 return (fabs(grid[0].
xgc[i])*grid[1].
dx[j]);
195 #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
196 return (grid[2].
dx[k]);
197 #elif GEOMETRY == CYLINDRICAL
198 return (fabs(grid[0].
xgc[i])*grid[2].
dx[k]);
199 #elif GEOMETRY == SPHERICAL
200 return (fabs(grid[0].
xgc[i]*sin(grid[1].
xgc[j]))*grid[2].
dx[k]);
219 #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
223 #elif GEOMETRY == POLAR
230 static double *inv_dl;
232 if (inv_dl == NULL) {
244 #elif GEOMETRY == SPHERICAL
248 static double *inv_dl2, *inv_dl3;
250 if (inv_dl2 == NULL) {
double Length_1(int i, int j, int k, Grid *grid)
void MakeGeometry(Grid *GXYZ)
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
int lend
Local end index for the local array.
double Length_3(int i, int j, int k, Grid *grid)
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
double * GetInverse_dl(const Grid *grid)
int g_dir
Specifies the current sweep or direction of integration.
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
double Length_2(int i, int j, int k, Grid *grid)
int nghost
Number of ghost zones.
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
double * r_1
Geometrical factor 1/r.
#define ARRAY_1D(nx, type)
double * ct
Geometrical factor cot(theta).
int np_tot
Total number of points in the local domain (boundaries included).
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
double * A
Right interface area, A[i] = .