Actual source code: da3.c
2: /*
3: Code for manipulating distributed regular 3d arrays in parallel.
4: File created by Peter Mell 7/14/95
5: */
7: #include <petsc/private/dmdaimpl.h>
9: #include <petscdraw.h>
10: static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer)
11: {
12: PetscMPIInt rank;
13: PetscBool iascii, isdraw, isglvis, isbinary;
14: DM_DA *dd = (DM_DA *)da->data;
15: Vec coordinates;
16: #if defined(PETSC_HAVE_MATLAB)
17: PetscBool ismatlab;
18: #endif
20: MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank);
22: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
23: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
24: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis);
25: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
26: #if defined(PETSC_HAVE_MATLAB)
27: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab);
28: #endif
29: if (iascii) {
30: PetscViewerFormat format;
32: PetscViewerASCIIPushSynchronized(viewer);
33: PetscViewerGetFormat(viewer, &format);
34: if (format == PETSC_VIEWER_LOAD_BALANCE) {
35: PetscInt i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal;
36: DMDALocalInfo info;
37: PetscMPIInt size;
38: MPI_Comm_size(PetscObjectComm((PetscObject)da), &size);
39: DMDAGetLocalInfo(da, &info);
40: nzlocal = info.xm * info.ym * info.zm;
41: PetscMalloc1(size, &nz);
42: MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da));
43: for (i = 0; i < (PetscInt)size; i++) {
44: nmax = PetscMax(nmax, nz[i]);
45: nmin = PetscMin(nmin, nz[i]);
46: navg += nz[i];
47: }
48: PetscFree(nz);
49: navg = navg / size;
50: PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax);
51: return 0;
52: }
53: if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) {
54: DMDALocalInfo info;
55: DMDAGetLocalInfo(da, &info);
56: PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " P %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " p %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->P, dd->m, dd->n, dd->p, dd->w, dd->s);
57: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Z range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs,
58: info.xs + info.xm, info.ys, info.ys + info.ym, info.zs, info.zs + info.zm));
59: DMGetCoordinates(da, &coordinates);
60: #if !defined(PETSC_USE_COMPLEX)
61: if (coordinates) {
62: PetscInt last;
63: const PetscReal *coors;
64: VecGetArrayRead(coordinates, &coors);
65: VecGetLocalSize(coordinates, &last);
66: last = last - 3;
67: PetscViewerASCIISynchronizedPrintf(viewer, "Lower left corner %g %g %g : Upper right %g %g %g\n", (double)coors[0], (double)coors[1], (double)coors[2], (double)coors[last], (double)coors[last + 1], (double)coors[last + 2]);
68: VecRestoreArrayRead(coordinates, &coors);
69: }
70: #endif
71: PetscViewerFlush(viewer);
72: PetscViewerASCIIPopSynchronized(viewer);
73: } else if (format == PETSC_VIEWER_ASCII_GLVIS) DMView_DA_GLVis(da, viewer);
74: else DMView_DA_VTK(da, viewer);
75: } else if (isdraw) {
76: PetscDraw draw;
77: PetscReal ymin = -1.0, ymax = (PetscReal)dd->N;
78: PetscReal xmin = -1.0, xmax = (PetscReal)((dd->M + 2) * dd->P), x, y, ycoord, xcoord;
79: PetscInt k, plane, base;
80: const PetscInt *idx;
81: char node[10];
82: PetscBool isnull;
84: PetscViewerDrawGetDraw(viewer, 0, &draw);
85: PetscDrawIsNull(draw, &isnull);
86: if (isnull) return 0;
88: PetscDrawCheckResizedWindow(draw);
89: PetscDrawClear(draw);
90: PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax);
92: PetscDrawCollectiveBegin(draw);
93: /* first processor draw all node lines */
94: if (rank == 0) {
95: for (k = 0; k < dd->P; k++) {
96: ymin = 0.0;
97: ymax = (PetscReal)(dd->N - 1);
98: for (xmin = (PetscReal)(k * (dd->M + 1)); xmin < (PetscReal)(dd->M + (k * (dd->M + 1))); xmin++) PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK);
99: xmin = (PetscReal)(k * (dd->M + 1));
100: xmax = xmin + (PetscReal)(dd->M - 1);
101: for (ymin = 0; ymin < (PetscReal)dd->N; ymin++) PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK);
102: }
103: }
104: PetscDrawCollectiveEnd(draw);
105: PetscDrawFlush(draw);
106: PetscDrawPause(draw);
108: PetscDrawCollectiveBegin(draw);
109: /*Go through and draw for each plane*/
110: for (k = 0; k < dd->P; k++) {
111: if ((k >= dd->zs) && (k < dd->ze)) {
112: /* draw my box */
113: ymin = dd->ys;
114: ymax = dd->ye - 1;
115: xmin = dd->xs / dd->w + (dd->M + 1) * k;
116: xmax = (dd->xe - 1) / dd->w + (dd->M + 1) * k;
118: PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED);
119: PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED);
120: PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED);
121: PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED);
123: xmin = dd->xs / dd->w;
124: xmax = (dd->xe - 1) / dd->w;
126: /* identify which processor owns the box */
127: PetscSNPrintf(node, sizeof(node), "%d", (int)rank);
128: PetscDrawString(draw, xmin + (dd->M + 1) * k + .2, ymin + .3, PETSC_DRAW_RED, node);
129: /* put in numbers*/
130: base = (dd->base + (dd->xe - dd->xs) * (dd->ye - dd->ys) * (k - dd->zs)) / dd->w;
131: for (y = ymin; y <= ymax; y++) {
132: for (x = xmin + (dd->M + 1) * k; x <= xmax + (dd->M + 1) * k; x++) {
133: PetscSNPrintf(node, sizeof(node), "%d", (int)base++);
134: PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node);
135: }
136: }
137: }
138: }
139: PetscDrawCollectiveEnd(draw);
140: PetscDrawFlush(draw);
141: PetscDrawPause(draw);
143: PetscDrawCollectiveBegin(draw);
144: for (k = 0 - dd->s; k < dd->P + dd->s; k++) {
145: /* Go through and draw for each plane */
146: if ((k >= dd->Zs) && (k < dd->Ze)) {
147: /* overlay ghost numbers, useful for error checking */
148: base = (dd->Xe - dd->Xs) * (dd->Ye - dd->Ys) * (k - dd->Zs) / dd->w;
149: ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx);
150: plane = k;
151: /* Keep z wrap around points on the drawing */
152: if (k < 0) plane = dd->P + k;
153: if (k >= dd->P) plane = k - dd->P;
154: ymin = dd->Ys;
155: ymax = dd->Ye;
156: xmin = (dd->M + 1) * plane * dd->w;
157: xmax = (dd->M + 1) * plane * dd->w + dd->M * dd->w;
158: for (y = ymin; y < ymax; y++) {
159: for (x = xmin + dd->Xs; x < xmin + dd->Xe; x += dd->w) {
160: sprintf(node, "%d", (int)(idx[base]));
161: ycoord = y;
162: /*Keep y wrap around points on drawing */
163: if (y < 0) ycoord = dd->N + y;
164: if (y >= dd->N) ycoord = y - dd->N;
165: xcoord = x; /* Keep x wrap points on drawing */
166: if (x < xmin) xcoord = xmax - (xmin - x);
167: if (x >= xmax) xcoord = xmin + (x - xmax);
168: PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node);
169: base++;
170: }
171: }
172: ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx);
173: }
174: }
175: PetscDrawCollectiveEnd(draw);
176: PetscDrawFlush(draw);
177: PetscDrawPause(draw);
178: PetscDrawSave(draw);
179: } else if (isglvis) {
180: DMView_DA_GLVis(da, viewer);
181: } else if (isbinary) {
182: DMView_DA_Binary(da, viewer);
183: #if defined(PETSC_HAVE_MATLAB)
184: } else if (ismatlab) {
185: DMView_DA_Matlab(da, viewer);
186: #endif
187: }
188: return 0;
189: }
191: PetscErrorCode DMSetUp_DA_3D(DM da)
192: {
193: DM_DA *dd = (DM_DA *)da->data;
194: const PetscInt M = dd->M;
195: const PetscInt N = dd->N;
196: const PetscInt P = dd->P;
197: PetscInt m = dd->m;
198: PetscInt n = dd->n;
199: PetscInt p = dd->p;
200: const PetscInt dof = dd->w;
201: const PetscInt s = dd->s;
202: DMBoundaryType bx = dd->bx;
203: DMBoundaryType by = dd->by;
204: DMBoundaryType bz = dd->bz;
205: DMDAStencilType stencil_type = dd->stencil_type;
206: PetscInt *lx = dd->lx;
207: PetscInt *ly = dd->ly;
208: PetscInt *lz = dd->lz;
209: MPI_Comm comm;
210: PetscMPIInt rank, size;
211: PetscInt xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0;
212: PetscInt Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm;
213: PetscInt left, right, up, down, bottom, top, i, j, k, *idx, nn;
214: PetscInt n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14;
215: PetscInt n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26;
216: PetscInt *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z;
217: PetscInt sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0;
218: PetscInt sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0;
219: PetscInt sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0;
220: Vec local, global;
221: VecScatter gtol;
222: IS to, from;
223: PetscBool twod;
226: PetscObjectGetComm((PetscObject)da, &comm);
227: #if !defined(PETSC_USE_64BIT_INDICES)
229: #endif
231: MPI_Comm_size(comm, &size);
232: MPI_Comm_rank(comm, &rank);
234: if (m != PETSC_DECIDE) {
237: }
238: if (n != PETSC_DECIDE) {
241: }
242: if (p != PETSC_DECIDE) {
245: }
248: /* Partition the array among the processors */
249: if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
250: m = size / (n * p);
251: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
252: n = size / (m * p);
253: } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
254: p = size / (m * n);
255: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
256: /* try for squarish distribution */
257: m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N * p)));
258: if (!m) m = 1;
259: while (m > 0) {
260: n = size / (m * p);
261: if (m * n * p == size) break;
262: m--;
263: }
265: if (M > N && m < n) {
266: PetscInt _m = m;
267: m = n;
268: n = _m;
269: }
270: } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
271: /* try for squarish distribution */
272: m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
273: if (!m) m = 1;
274: while (m > 0) {
275: p = size / (m * n);
276: if (m * n * p == size) break;
277: m--;
278: }
280: if (M > P && m < p) {
281: PetscInt _m = m;
282: m = p;
283: p = _m;
284: }
285: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
286: /* try for squarish distribution */
287: n = (int)(0.5 + PetscSqrtReal(((PetscReal)N) * ((PetscReal)size) / ((PetscReal)P * m)));
288: if (!n) n = 1;
289: while (n > 0) {
290: p = size / (m * n);
291: if (m * n * p == size) break;
292: n--;
293: }
295: if (N > P && n < p) {
296: PetscInt _n = n;
297: n = p;
298: p = _n;
299: }
300: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
301: /* try for squarish distribution */
302: n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.)));
303: if (!n) n = 1;
304: while (n > 0) {
305: pm = size / n;
306: if (n * pm == size) break;
307: n--;
308: }
309: if (!n) n = 1;
310: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
311: if (!m) m = 1;
312: while (m > 0) {
313: p = size / (m * n);
314: if (m * n * p == size) break;
315: m--;
316: }
317: if (M > P && m < p) {
318: PetscInt _m = m;
319: m = p;
320: p = _m;
321: }
329: /*
330: Determine locally owned region
331: [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
332: */
334: if (!lx) {
335: PetscMalloc1(m, &dd->lx);
336: lx = dd->lx;
337: for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m));
338: }
339: x = lx[rank % m];
340: xs = 0;
341: for (i = 0; i < (rank % m); i++) xs += lx[i];
344: if (!ly) {
345: PetscMalloc1(n, &dd->ly);
346: ly = dd->ly;
347: for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n));
348: }
349: y = ly[(rank % (m * n)) / m];
352: ys = 0;
353: for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i];
355: if (!lz) {
356: PetscMalloc1(p, &dd->lz);
357: lz = dd->lz;
358: for (i = 0; i < p; i++) lz[i] = P / p + ((P % p) > (i % p));
359: }
360: z = lz[rank / (m * n)];
362: /* note this is different than x- and y-, as we will handle as an important special
363: case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems
364: in a 3D code. Additional code for this case is noted with "2d case" comments */
365: twod = PETSC_FALSE;
366: if (P == 1) twod = PETSC_TRUE;
368: zs = 0;
369: for (i = 0; i < (rank / (m * n)); i++) zs += lz[i];
370: ye = ys + y;
371: xe = xs + x;
372: ze = zs + z;
374: /* determine ghost region (Xs) and region scattered into (IXs) */
375: if (xs - s > 0) {
376: Xs = xs - s;
377: IXs = xs - s;
378: } else {
379: if (bx) Xs = xs - s;
380: else Xs = 0;
381: IXs = 0;
382: }
383: if (xe + s <= M) {
384: Xe = xe + s;
385: IXe = xe + s;
386: } else {
387: if (bx) {
388: Xs = xs - s;
389: Xe = xe + s;
390: } else Xe = M;
391: IXe = M;
392: }
394: if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
395: IXs = xs - s;
396: IXe = xe + s;
397: Xs = xs - s;
398: Xe = xe + s;
399: }
401: if (ys - s > 0) {
402: Ys = ys - s;
403: IYs = ys - s;
404: } else {
405: if (by) Ys = ys - s;
406: else Ys = 0;
407: IYs = 0;
408: }
409: if (ye + s <= N) {
410: Ye = ye + s;
411: IYe = ye + s;
412: } else {
413: if (by) Ye = ye + s;
414: else Ye = N;
415: IYe = N;
416: }
418: if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
419: IYs = ys - s;
420: IYe = ye + s;
421: Ys = ys - s;
422: Ye = ye + s;
423: }
425: if (zs - s > 0) {
426: Zs = zs - s;
427: IZs = zs - s;
428: } else {
429: if (bz) Zs = zs - s;
430: else Zs = 0;
431: IZs = 0;
432: }
433: if (ze + s <= P) {
434: Ze = ze + s;
435: IZe = ze + s;
436: } else {
437: if (bz) Ze = ze + s;
438: else Ze = P;
439: IZe = P;
440: }
442: if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
443: IZs = zs - s;
444: IZe = ze + s;
445: Zs = zs - s;
446: Ze = ze + s;
447: }
449: /* Resize all X parameters to reflect w */
450: s_x = s;
451: s_y = s;
452: s_z = s;
454: /* determine starting point of each processor */
455: nn = x * y * z;
456: PetscMalloc2(size + 1, &bases, size, &ldims);
457: MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm);
458: bases[0] = 0;
459: for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
460: for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
461: base = bases[rank] * dof;
463: /* allocate the base parallel and sequential vectors */
464: dd->Nlocal = x * y * z * dof;
465: VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global);
466: dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof;
467: VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local);
469: /* generate global to local vector scatter and local to global mapping*/
471: /* global to local must include ghost points within the domain,
472: but not ghost points outside the domain that aren't periodic */
473: PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx);
474: if (stencil_type == DMDA_STENCIL_BOX) {
475: left = IXs - Xs;
476: right = left + (IXe - IXs);
477: bottom = IYs - Ys;
478: top = bottom + (IYe - IYs);
479: down = IZs - Zs;
480: up = down + (IZe - IZs);
481: count = 0;
482: for (i = down; i < up; i++) {
483: for (j = bottom; j < top; j++) {
484: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
485: }
486: }
487: ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to);
488: } else {
489: /* This is way ugly! We need to list the funny cross type region */
490: left = xs - Xs;
491: right = left + x;
492: bottom = ys - Ys;
493: top = bottom + y;
494: down = zs - Zs;
495: up = down + z;
496: count = 0;
497: /* the bottom chunk */
498: for (i = (IZs - Zs); i < down; i++) {
499: for (j = bottom; j < top; j++) {
500: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
501: }
502: }
503: /* the middle piece */
504: for (i = down; i < up; i++) {
505: /* front */
506: for (j = (IYs - Ys); j < bottom; j++) {
507: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
508: }
509: /* middle */
510: for (j = bottom; j < top; j++) {
511: for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
512: }
513: /* back */
514: for (j = top; j < top + IYe - ye; j++) {
515: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
516: }
517: }
518: /* the top piece */
519: for (i = up; i < up + IZe - ze; i++) {
520: for (j = bottom; j < top; j++) {
521: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
522: }
523: }
524: ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to);
525: }
527: /* determine who lies on each side of use stored in n24 n25 n26
528: n21 n22 n23
529: n18 n19 n20
531: n15 n16 n17
532: n12 n14
533: n9 n10 n11
535: n6 n7 n8
536: n3 n4 n5
537: n0 n1 n2
538: */
540: /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
541: /* Assume Nodes are Internal to the Cube */
542: n0 = rank - m * n - m - 1;
543: n1 = rank - m * n - m;
544: n2 = rank - m * n - m + 1;
545: n3 = rank - m * n - 1;
546: n4 = rank - m * n;
547: n5 = rank - m * n + 1;
548: n6 = rank - m * n + m - 1;
549: n7 = rank - m * n + m;
550: n8 = rank - m * n + m + 1;
552: n9 = rank - m - 1;
553: n10 = rank - m;
554: n11 = rank - m + 1;
555: n12 = rank - 1;
556: n14 = rank + 1;
557: n15 = rank + m - 1;
558: n16 = rank + m;
559: n17 = rank + m + 1;
561: n18 = rank + m * n - m - 1;
562: n19 = rank + m * n - m;
563: n20 = rank + m * n - m + 1;
564: n21 = rank + m * n - 1;
565: n22 = rank + m * n;
566: n23 = rank + m * n + 1;
567: n24 = rank + m * n + m - 1;
568: n25 = rank + m * n + m;
569: n26 = rank + m * n + m + 1;
571: /* Assume Pieces are on Faces of Cube */
573: if (xs == 0) { /* First assume not corner or edge */
574: n0 = rank - 1 - (m * n);
575: n3 = rank + m - 1 - (m * n);
576: n6 = rank + 2 * m - 1 - (m * n);
577: n9 = rank - 1;
578: n12 = rank + m - 1;
579: n15 = rank + 2 * m - 1;
580: n18 = rank - 1 + (m * n);
581: n21 = rank + m - 1 + (m * n);
582: n24 = rank + 2 * m - 1 + (m * n);
583: }
585: if (xe == M) { /* First assume not corner or edge */
586: n2 = rank - 2 * m + 1 - (m * n);
587: n5 = rank - m + 1 - (m * n);
588: n8 = rank + 1 - (m * n);
589: n11 = rank - 2 * m + 1;
590: n14 = rank - m + 1;
591: n17 = rank + 1;
592: n20 = rank - 2 * m + 1 + (m * n);
593: n23 = rank - m + 1 + (m * n);
594: n26 = rank + 1 + (m * n);
595: }
597: if (ys == 0) { /* First assume not corner or edge */
598: n0 = rank + m * (n - 1) - 1 - (m * n);
599: n1 = rank + m * (n - 1) - (m * n);
600: n2 = rank + m * (n - 1) + 1 - (m * n);
601: n9 = rank + m * (n - 1) - 1;
602: n10 = rank + m * (n - 1);
603: n11 = rank + m * (n - 1) + 1;
604: n18 = rank + m * (n - 1) - 1 + (m * n);
605: n19 = rank + m * (n - 1) + (m * n);
606: n20 = rank + m * (n - 1) + 1 + (m * n);
607: }
609: if (ye == N) { /* First assume not corner or edge */
610: n6 = rank - m * (n - 1) - 1 - (m * n);
611: n7 = rank - m * (n - 1) - (m * n);
612: n8 = rank - m * (n - 1) + 1 - (m * n);
613: n15 = rank - m * (n - 1) - 1;
614: n16 = rank - m * (n - 1);
615: n17 = rank - m * (n - 1) + 1;
616: n24 = rank - m * (n - 1) - 1 + (m * n);
617: n25 = rank - m * (n - 1) + (m * n);
618: n26 = rank - m * (n - 1) + 1 + (m * n);
619: }
621: if (zs == 0) { /* First assume not corner or edge */
622: n0 = size - (m * n) + rank - m - 1;
623: n1 = size - (m * n) + rank - m;
624: n2 = size - (m * n) + rank - m + 1;
625: n3 = size - (m * n) + rank - 1;
626: n4 = size - (m * n) + rank;
627: n5 = size - (m * n) + rank + 1;
628: n6 = size - (m * n) + rank + m - 1;
629: n7 = size - (m * n) + rank + m;
630: n8 = size - (m * n) + rank + m + 1;
631: }
633: if (ze == P) { /* First assume not corner or edge */
634: n18 = (m * n) - (size - rank) - m - 1;
635: n19 = (m * n) - (size - rank) - m;
636: n20 = (m * n) - (size - rank) - m + 1;
637: n21 = (m * n) - (size - rank) - 1;
638: n22 = (m * n) - (size - rank);
639: n23 = (m * n) - (size - rank) + 1;
640: n24 = (m * n) - (size - rank) + m - 1;
641: n25 = (m * n) - (size - rank) + m;
642: n26 = (m * n) - (size - rank) + m + 1;
643: }
645: if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */
646: n0 = size - m * n + rank + m - 1 - m;
647: n3 = size - m * n + rank + m - 1;
648: n6 = size - m * n + rank + m - 1 + m;
649: }
651: if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */
652: n18 = m * n - (size - rank) + m - 1 - m;
653: n21 = m * n - (size - rank) + m - 1;
654: n24 = m * n - (size - rank) + m - 1 + m;
655: }
657: if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */
658: n0 = rank + m * n - 1 - m * n;
659: n9 = rank + m * n - 1;
660: n18 = rank + m * n - 1 + m * n;
661: }
663: if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */
664: n6 = rank - m * (n - 1) + m - 1 - m * n;
665: n15 = rank - m * (n - 1) + m - 1;
666: n24 = rank - m * (n - 1) + m - 1 + m * n;
667: }
669: if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */
670: n2 = size - (m * n - rank) - (m - 1) - m;
671: n5 = size - (m * n - rank) - (m - 1);
672: n8 = size - (m * n - rank) - (m - 1) + m;
673: }
675: if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */
676: n20 = m * n - (size - rank) - (m - 1) - m;
677: n23 = m * n - (size - rank) - (m - 1);
678: n26 = m * n - (size - rank) - (m - 1) + m;
679: }
681: if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */
682: n2 = rank + m * (n - 1) - (m - 1) - m * n;
683: n11 = rank + m * (n - 1) - (m - 1);
684: n20 = rank + m * (n - 1) - (m - 1) + m * n;
685: }
687: if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */
688: n8 = rank - m * n + 1 - m * n;
689: n17 = rank - m * n + 1;
690: n26 = rank - m * n + 1 + m * n;
691: }
693: if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */
694: n0 = size - m + rank - 1;
695: n1 = size - m + rank;
696: n2 = size - m + rank + 1;
697: }
699: if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */
700: n18 = m * n - (size - rank) + m * (n - 1) - 1;
701: n19 = m * n - (size - rank) + m * (n - 1);
702: n20 = m * n - (size - rank) + m * (n - 1) + 1;
703: }
705: if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */
706: n6 = size - (m * n - rank) - m * (n - 1) - 1;
707: n7 = size - (m * n - rank) - m * (n - 1);
708: n8 = size - (m * n - rank) - m * (n - 1) + 1;
709: }
711: if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */
712: n24 = rank - (size - m) - 1;
713: n25 = rank - (size - m);
714: n26 = rank - (size - m) + 1;
715: }
717: /* Check for Corners */
718: if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1;
719: if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1;
720: if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1);
721: if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1;
722: if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m;
723: if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m;
724: if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n;
725: if ((xe == M) && (ye == N) && (ze == P)) n26 = 0;
727: /* Check for when not X,Y, and Z Periodic */
729: /* If not X periodic */
730: if (bx != DM_BOUNDARY_PERIODIC) {
731: if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;
732: if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;
733: }
735: /* If not Y periodic */
736: if (by != DM_BOUNDARY_PERIODIC) {
737: if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;
738: if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;
739: }
741: /* If not Z periodic */
742: if (bz != DM_BOUNDARY_PERIODIC) {
743: if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;
744: if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;
745: }
747: PetscMalloc1(27, &dd->neighbors);
749: dd->neighbors[0] = n0;
750: dd->neighbors[1] = n1;
751: dd->neighbors[2] = n2;
752: dd->neighbors[3] = n3;
753: dd->neighbors[4] = n4;
754: dd->neighbors[5] = n5;
755: dd->neighbors[6] = n6;
756: dd->neighbors[7] = n7;
757: dd->neighbors[8] = n8;
758: dd->neighbors[9] = n9;
759: dd->neighbors[10] = n10;
760: dd->neighbors[11] = n11;
761: dd->neighbors[12] = n12;
762: dd->neighbors[13] = rank;
763: dd->neighbors[14] = n14;
764: dd->neighbors[15] = n15;
765: dd->neighbors[16] = n16;
766: dd->neighbors[17] = n17;
767: dd->neighbors[18] = n18;
768: dd->neighbors[19] = n19;
769: dd->neighbors[20] = n20;
770: dd->neighbors[21] = n21;
771: dd->neighbors[22] = n22;
772: dd->neighbors[23] = n23;
773: dd->neighbors[24] = n24;
774: dd->neighbors[25] = n25;
775: dd->neighbors[26] = n26;
777: /* If star stencil then delete the corner neighbors */
778: if (stencil_type == DMDA_STENCIL_STAR) {
779: /* save information about corner neighbors */
780: sn0 = n0;
781: sn1 = n1;
782: sn2 = n2;
783: sn3 = n3;
784: sn5 = n5;
785: sn6 = n6;
786: sn7 = n7;
787: sn8 = n8;
788: sn9 = n9;
789: sn11 = n11;
790: sn15 = n15;
791: sn17 = n17;
792: sn18 = n18;
793: sn19 = n19;
794: sn20 = n20;
795: sn21 = n21;
796: sn23 = n23;
797: sn24 = n24;
798: sn25 = n25;
799: sn26 = n26;
800: n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
801: }
803: PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx);
805: nn = 0;
806: /* Bottom Level */
807: for (k = 0; k < s_z; k++) {
808: for (i = 1; i <= s_y; i++) {
809: if (n0 >= 0) { /* left below */
810: x_t = lx[n0 % m];
811: y_t = ly[(n0 % (m * n)) / m];
812: z_t = lz[n0 / (m * n)];
813: s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x - (s_z - k - 1) * x_t * y_t;
814: if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */
815: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
816: }
817: if (n1 >= 0) { /* directly below */
818: x_t = x;
819: y_t = ly[(n1 % (m * n)) / m];
820: z_t = lz[n1 / (m * n)];
821: s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
822: if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
823: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
824: }
825: if (n2 >= 0) { /* right below */
826: x_t = lx[n2 % m];
827: y_t = ly[(n2 % (m * n)) / m];
828: z_t = lz[n2 / (m * n)];
829: s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
830: if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
831: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
832: }
833: }
835: for (i = 0; i < y; i++) {
836: if (n3 >= 0) { /* directly left */
837: x_t = lx[n3 % m];
838: y_t = y;
839: z_t = lz[n3 / (m * n)];
840: s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
841: if (twod && (s_t < 0)) s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - x_t * y_t; /* 2D case */
842: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
843: }
845: if (n4 >= 0) { /* middle */
846: x_t = x;
847: y_t = y;
848: z_t = lz[n4 / (m * n)];
849: s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
850: if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
851: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
852: } else if (bz == DM_BOUNDARY_MIRROR) {
853: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
854: }
856: if (n5 >= 0) { /* directly right */
857: x_t = lx[n5 % m];
858: y_t = y;
859: z_t = lz[n5 / (m * n)];
860: s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
861: if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
862: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
863: }
864: }
866: for (i = 1; i <= s_y; i++) {
867: if (n6 >= 0) { /* left above */
868: x_t = lx[n6 % m];
869: y_t = ly[(n6 % (m * n)) / m];
870: z_t = lz[n6 / (m * n)];
871: s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
872: if (twod && (s_t < 0)) s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - x_t * y_t; /* 2D case */
873: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
874: }
875: if (n7 >= 0) { /* directly above */
876: x_t = x;
877: y_t = ly[(n7 % (m * n)) / m];
878: z_t = lz[n7 / (m * n)];
879: s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
880: if (twod && (s_t < 0)) s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
881: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
882: }
883: if (n8 >= 0) { /* right above */
884: x_t = lx[n8 % m];
885: y_t = ly[(n8 % (m * n)) / m];
886: z_t = lz[n8 / (m * n)];
887: s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
888: if (twod && (s_t < 0)) s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
889: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
890: }
891: }
892: }
894: /* Middle Level */
895: for (k = 0; k < z; k++) {
896: for (i = 1; i <= s_y; i++) {
897: if (n9 >= 0) { /* left below */
898: x_t = lx[n9 % m];
899: y_t = ly[(n9 % (m * n)) / m];
900: /* z_t = z; */
901: s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
902: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
903: }
904: if (n10 >= 0) { /* directly below */
905: x_t = x;
906: y_t = ly[(n10 % (m * n)) / m];
907: /* z_t = z; */
908: s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
909: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
910: } else if (by == DM_BOUNDARY_MIRROR) {
911: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
912: }
913: if (n11 >= 0) { /* right below */
914: x_t = lx[n11 % m];
915: y_t = ly[(n11 % (m * n)) / m];
916: /* z_t = z; */
917: s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
918: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
919: }
920: }
922: for (i = 0; i < y; i++) {
923: if (n12 >= 0) { /* directly left */
924: x_t = lx[n12 % m];
925: y_t = y;
926: /* z_t = z; */
927: s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
928: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
929: } else if (bx == DM_BOUNDARY_MIRROR) {
930: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
931: }
933: /* Interior */
934: s_t = bases[rank] + i * x + k * x * y;
935: for (j = 0; j < x; j++) idx[nn++] = s_t++;
937: if (n14 >= 0) { /* directly right */
938: x_t = lx[n14 % m];
939: y_t = y;
940: /* z_t = z; */
941: s_t = bases[n14] + i * x_t + k * x_t * y_t;
942: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
943: } else if (bx == DM_BOUNDARY_MIRROR) {
944: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
945: }
946: }
948: for (i = 1; i <= s_y; i++) {
949: if (n15 >= 0) { /* left above */
950: x_t = lx[n15 % m];
951: y_t = ly[(n15 % (m * n)) / m];
952: /* z_t = z; */
953: s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
954: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
955: }
956: if (n16 >= 0) { /* directly above */
957: x_t = x;
958: y_t = ly[(n16 % (m * n)) / m];
959: /* z_t = z; */
960: s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
961: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
962: } else if (by == DM_BOUNDARY_MIRROR) {
963: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
964: }
965: if (n17 >= 0) { /* right above */
966: x_t = lx[n17 % m];
967: y_t = ly[(n17 % (m * n)) / m];
968: /* z_t = z; */
969: s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
970: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
971: }
972: }
973: }
975: /* Upper Level */
976: for (k = 0; k < s_z; k++) {
977: for (i = 1; i <= s_y; i++) {
978: if (n18 >= 0) { /* left below */
979: x_t = lx[n18 % m];
980: y_t = ly[(n18 % (m * n)) / m];
981: /* z_t = lz[n18 / (m*n)]; */
982: s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
983: if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */
984: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
985: }
986: if (n19 >= 0) { /* directly below */
987: x_t = x;
988: y_t = ly[(n19 % (m * n)) / m];
989: /* z_t = lz[n19 / (m*n)]; */
990: s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
991: if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
992: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
993: }
994: if (n20 >= 0) { /* right below */
995: x_t = lx[n20 % m];
996: y_t = ly[(n20 % (m * n)) / m];
997: /* z_t = lz[n20 / (m*n)]; */
998: s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
999: if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
1000: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1001: }
1002: }
1004: for (i = 0; i < y; i++) {
1005: if (n21 >= 0) { /* directly left */
1006: x_t = lx[n21 % m];
1007: y_t = y;
1008: /* z_t = lz[n21 / (m*n)]; */
1009: s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1010: if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */
1011: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1012: }
1014: if (n22 >= 0) { /* middle */
1015: x_t = x;
1016: y_t = y;
1017: /* z_t = lz[n22 / (m*n)]; */
1018: s_t = bases[n22] + i * x_t + k * x_t * y_t;
1019: if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */
1020: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1021: } else if (bz == DM_BOUNDARY_MIRROR) {
1022: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1023: }
1025: if (n23 >= 0) { /* directly right */
1026: x_t = lx[n23 % m];
1027: y_t = y;
1028: /* z_t = lz[n23 / (m*n)]; */
1029: s_t = bases[n23] + i * x_t + k * x_t * y_t;
1030: if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */
1031: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1032: }
1033: }
1035: for (i = 1; i <= s_y; i++) {
1036: if (n24 >= 0) { /* left above */
1037: x_t = lx[n24 % m];
1038: y_t = ly[(n24 % (m * n)) / m];
1039: /* z_t = lz[n24 / (m*n)]; */
1040: s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1041: if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */
1042: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1043: }
1044: if (n25 >= 0) { /* directly above */
1045: x_t = x;
1046: y_t = ly[(n25 % (m * n)) / m];
1047: /* z_t = lz[n25 / (m*n)]; */
1048: s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1049: if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */
1050: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1051: }
1052: if (n26 >= 0) { /* right above */
1053: x_t = lx[n26 % m];
1054: y_t = ly[(n26 % (m * n)) / m];
1055: /* z_t = lz[n26 / (m*n)]; */
1056: s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1057: if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */
1058: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1059: }
1060: }
1061: }
1063: ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from);
1064: VecScatterCreate(global, from, local, to, >ol);
1065: ISDestroy(&to);
1066: ISDestroy(&from);
1068: if (stencil_type == DMDA_STENCIL_STAR) {
1069: n0 = sn0;
1070: n1 = sn1;
1071: n2 = sn2;
1072: n3 = sn3;
1073: n5 = sn5;
1074: n6 = sn6;
1075: n7 = sn7;
1076: n8 = sn8;
1077: n9 = sn9;
1078: n11 = sn11;
1079: n15 = sn15;
1080: n17 = sn17;
1081: n18 = sn18;
1082: n19 = sn19;
1083: n20 = sn20;
1084: n21 = sn21;
1085: n23 = sn23;
1086: n24 = sn24;
1087: n25 = sn25;
1088: n26 = sn26;
1089: }
1091: if (((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz))) {
1092: /*
1093: Recompute the local to global mappings, this time keeping the
1094: information about the cross corner processor numbers.
1095: */
1096: nn = 0;
1097: /* Bottom Level */
1098: for (k = 0; k < s_z; k++) {
1099: for (i = 1; i <= s_y; i++) {
1100: if (n0 >= 0) { /* left below */
1101: x_t = lx[n0 % m];
1102: y_t = ly[(n0 % (m * n)) / m];
1103: z_t = lz[n0 / (m * n)];
1104: s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x - (s_z - k - 1) * x_t * y_t;
1105: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1106: } else if (Xs - xs < 0 && Ys - ys < 0 && Zs - zs < 0) {
1107: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1108: }
1109: if (n1 >= 0) { /* directly below */
1110: x_t = x;
1111: y_t = ly[(n1 % (m * n)) / m];
1112: z_t = lz[n1 / (m * n)];
1113: s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1114: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1115: } else if (Ys - ys < 0 && Zs - zs < 0) {
1116: for (j = 0; j < x; j++) idx[nn++] = -1;
1117: }
1118: if (n2 >= 0) { /* right below */
1119: x_t = lx[n2 % m];
1120: y_t = ly[(n2 % (m * n)) / m];
1121: z_t = lz[n2 / (m * n)];
1122: s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1123: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1124: } else if (xe - Xe < 0 && Ys - ys < 0 && Zs - zs < 0) {
1125: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1126: }
1127: }
1129: for (i = 0; i < y; i++) {
1130: if (n3 >= 0) { /* directly left */
1131: x_t = lx[n3 % m];
1132: y_t = y;
1133: z_t = lz[n3 / (m * n)];
1134: s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1135: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1136: } else if (Xs - xs < 0 && Zs - zs < 0) {
1137: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1138: }
1140: if (n4 >= 0) { /* middle */
1141: x_t = x;
1142: y_t = y;
1143: z_t = lz[n4 / (m * n)];
1144: s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1145: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1146: } else if (Zs - zs < 0) {
1147: if (bz == DM_BOUNDARY_MIRROR) {
1148: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
1149: } else {
1150: for (j = 0; j < x; j++) idx[nn++] = -1;
1151: }
1152: }
1154: if (n5 >= 0) { /* directly right */
1155: x_t = lx[n5 % m];
1156: y_t = y;
1157: z_t = lz[n5 / (m * n)];
1158: s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1159: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1160: } else if (xe - Xe < 0 && Zs - zs < 0) {
1161: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1162: }
1163: }
1165: for (i = 1; i <= s_y; i++) {
1166: if (n6 >= 0) { /* left above */
1167: x_t = lx[n6 % m];
1168: y_t = ly[(n6 % (m * n)) / m];
1169: z_t = lz[n6 / (m * n)];
1170: s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1171: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1172: } else if (Xs - xs < 0 && ye - Ye < 0 && Zs - zs < 0) {
1173: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1174: }
1175: if (n7 >= 0) { /* directly above */
1176: x_t = x;
1177: y_t = ly[(n7 % (m * n)) / m];
1178: z_t = lz[n7 / (m * n)];
1179: s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1180: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1181: } else if (ye - Ye < 0 && Zs - zs < 0) {
1182: for (j = 0; j < x; j++) idx[nn++] = -1;
1183: }
1184: if (n8 >= 0) { /* right above */
1185: x_t = lx[n8 % m];
1186: y_t = ly[(n8 % (m * n)) / m];
1187: z_t = lz[n8 / (m * n)];
1188: s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1189: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1190: } else if (xe - Xe < 0 && ye - Ye < 0 && Zs - zs < 0) {
1191: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1192: }
1193: }
1194: }
1196: /* Middle Level */
1197: for (k = 0; k < z; k++) {
1198: for (i = 1; i <= s_y; i++) {
1199: if (n9 >= 0) { /* left below */
1200: x_t = lx[n9 % m];
1201: y_t = ly[(n9 % (m * n)) / m];
1202: /* z_t = z; */
1203: s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1204: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1205: } else if (Xs - xs < 0 && Ys - ys < 0) {
1206: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1207: }
1208: if (n10 >= 0) { /* directly below */
1209: x_t = x;
1210: y_t = ly[(n10 % (m * n)) / m];
1211: /* z_t = z; */
1212: s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1213: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1214: } else if (Ys - ys < 0) {
1215: if (by == DM_BOUNDARY_MIRROR) {
1216: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
1217: } else {
1218: for (j = 0; j < x; j++) idx[nn++] = -1;
1219: }
1220: }
1221: if (n11 >= 0) { /* right below */
1222: x_t = lx[n11 % m];
1223: y_t = ly[(n11 % (m * n)) / m];
1224: /* z_t = z; */
1225: s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1226: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1227: } else if (xe - Xe < 0 && Ys - ys < 0) {
1228: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1229: }
1230: }
1232: for (i = 0; i < y; i++) {
1233: if (n12 >= 0) { /* directly left */
1234: x_t = lx[n12 % m];
1235: y_t = y;
1236: /* z_t = z; */
1237: s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
1238: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1239: } else if (Xs - xs < 0) {
1240: if (bx == DM_BOUNDARY_MIRROR) {
1241: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
1242: } else {
1243: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1244: }
1245: }
1247: /* Interior */
1248: s_t = bases[rank] + i * x + k * x * y;
1249: for (j = 0; j < x; j++) idx[nn++] = s_t++;
1251: if (n14 >= 0) { /* directly right */
1252: x_t = lx[n14 % m];
1253: y_t = y;
1254: /* z_t = z; */
1255: s_t = bases[n14] + i * x_t + k * x_t * y_t;
1256: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1257: } else if (xe - Xe < 0) {
1258: if (bx == DM_BOUNDARY_MIRROR) {
1259: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
1260: } else {
1261: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1262: }
1263: }
1264: }
1266: for (i = 1; i <= s_y; i++) {
1267: if (n15 >= 0) { /* left above */
1268: x_t = lx[n15 % m];
1269: y_t = ly[(n15 % (m * n)) / m];
1270: /* z_t = z; */
1271: s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
1272: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1273: } else if (Xs - xs < 0 && ye - Ye < 0) {
1274: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1275: }
1276: if (n16 >= 0) { /* directly above */
1277: x_t = x;
1278: y_t = ly[(n16 % (m * n)) / m];
1279: /* z_t = z; */
1280: s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
1281: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1282: } else if (ye - Ye < 0) {
1283: if (by == DM_BOUNDARY_MIRROR) {
1284: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
1285: } else {
1286: for (j = 0; j < x; j++) idx[nn++] = -1;
1287: }
1288: }
1289: if (n17 >= 0) { /* right above */
1290: x_t = lx[n17 % m];
1291: y_t = ly[(n17 % (m * n)) / m];
1292: /* z_t = z; */
1293: s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
1294: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1295: } else if (xe - Xe < 0 && ye - Ye < 0) {
1296: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1297: }
1298: }
1299: }
1301: /* Upper Level */
1302: for (k = 0; k < s_z; k++) {
1303: for (i = 1; i <= s_y; i++) {
1304: if (n18 >= 0) { /* left below */
1305: x_t = lx[n18 % m];
1306: y_t = ly[(n18 % (m * n)) / m];
1307: /* z_t = lz[n18 / (m*n)]; */
1308: s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1309: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1310: } else if (Xs - xs < 0 && Ys - ys < 0 && ze - Ze < 0) {
1311: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1312: }
1313: if (n19 >= 0) { /* directly below */
1314: x_t = x;
1315: y_t = ly[(n19 % (m * n)) / m];
1316: /* z_t = lz[n19 / (m*n)]; */
1317: s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1318: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1319: } else if (Ys - ys < 0 && ze - Ze < 0) {
1320: for (j = 0; j < x; j++) idx[nn++] = -1;
1321: }
1322: if (n20 >= 0) { /* right below */
1323: x_t = lx[n20 % m];
1324: y_t = ly[(n20 % (m * n)) / m];
1325: /* z_t = lz[n20 / (m*n)]; */
1326: s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1327: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1328: } else if (xe - Xe < 0 && Ys - ys < 0 && ze - Ze < 0) {
1329: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1330: }
1331: }
1333: for (i = 0; i < y; i++) {
1334: if (n21 >= 0) { /* directly left */
1335: x_t = lx[n21 % m];
1336: y_t = y;
1337: /* z_t = lz[n21 / (m*n)]; */
1338: s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1339: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1340: } else if (Xs - xs < 0 && ze - Ze < 0) {
1341: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1342: }
1344: if (n22 >= 0) { /* middle */
1345: x_t = x;
1346: y_t = y;
1347: /* z_t = lz[n22 / (m*n)]; */
1348: s_t = bases[n22] + i * x_t + k * x_t * y_t;
1349: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1350: } else if (ze - Ze < 0) {
1351: if (bz == DM_BOUNDARY_MIRROR) {
1352: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1353: } else {
1354: for (j = 0; j < x; j++) idx[nn++] = -1;
1355: }
1356: }
1358: if (n23 >= 0) { /* directly right */
1359: x_t = lx[n23 % m];
1360: y_t = y;
1361: /* z_t = lz[n23 / (m*n)]; */
1362: s_t = bases[n23] + i * x_t + k * x_t * y_t;
1363: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1364: } else if (xe - Xe < 0 && ze - Ze < 0) {
1365: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1366: }
1367: }
1369: for (i = 1; i <= s_y; i++) {
1370: if (n24 >= 0) { /* left above */
1371: x_t = lx[n24 % m];
1372: y_t = ly[(n24 % (m * n)) / m];
1373: /* z_t = lz[n24 / (m*n)]; */
1374: s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1375: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1376: } else if (Xs - xs < 0 && ye - Ye < 0 && ze - Ze < 0) {
1377: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1378: }
1379: if (n25 >= 0) { /* directly above */
1380: x_t = x;
1381: y_t = ly[(n25 % (m * n)) / m];
1382: /* z_t = lz[n25 / (m*n)]; */
1383: s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1384: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1385: } else if (ye - Ye < 0 && ze - Ze < 0) {
1386: for (j = 0; j < x; j++) idx[nn++] = -1;
1387: }
1388: if (n26 >= 0) { /* right above */
1389: x_t = lx[n26 % m];
1390: y_t = ly[(n26 % (m * n)) / m];
1391: /* z_t = lz[n26 / (m*n)]; */
1392: s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1393: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1394: } else if (xe - Xe < 0 && ye - Ye < 0 && ze - Ze < 0) {
1395: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1396: }
1397: }
1398: }
1399: }
1400: /*
1401: Set the local to global ordering in the global vector, this allows use
1402: of VecSetValuesLocal().
1403: */
1404: ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap);
1406: PetscFree2(bases, ldims);
1407: dd->m = m;
1408: dd->n = n;
1409: dd->p = p;
1410: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1411: dd->xs = xs * dof;
1412: dd->xe = xe * dof;
1413: dd->ys = ys;
1414: dd->ye = ye;
1415: dd->zs = zs;
1416: dd->ze = ze;
1417: dd->Xs = Xs * dof;
1418: dd->Xe = Xe * dof;
1419: dd->Ys = Ys;
1420: dd->Ye = Ye;
1421: dd->Zs = Zs;
1422: dd->Ze = Ze;
1424: VecDestroy(&local);
1425: VecDestroy(&global);
1427: dd->gtol = gtol;
1428: dd->base = base;
1429: da->ops->view = DMView_DA_3d;
1430: dd->ltol = NULL;
1431: dd->ao = NULL;
1432: return 0;
1433: }
1435: /*@C
1436: DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1437: regular array data that is distributed across some processors.
1439: Collective
1441: Input Parameters:
1442: + comm - MPI communicator
1443: . bx,by,bz - type of ghost nodes the array have.
1444: Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1445: . stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`)
1446: . M,N,P - global dimension in each direction of the array
1447: . m,n,p - corresponding number of processors in each dimension
1448: (or `PETSC_DECIDE` to have calculated)
1449: . dof - number of degrees of freedom per node
1450: . s - stencil width
1451: - lx, ly, lz - arrays containing the number of nodes in each cell along
1452: the x, y, and z coordinates, or NULL. If non-null, these
1453: must be of length as m,n,p and the corresponding
1454: m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
1455: the ly[] must N, sum of the lz[] must be P
1457: Output Parameter:
1458: . da - the resulting distributed array object
1460: Options Database Keys:
1461: + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()`
1462: . -da_grid_x <nx> - number of grid points in x direction
1463: . -da_grid_y <ny> - number of grid points in y direction
1464: . -da_grid_z <nz> - number of grid points in z direction
1465: . -da_processors_x <MX> - number of processors in x direction
1466: . -da_processors_y <MY> - number of processors in y direction
1467: . -da_processors_z <MZ> - number of processors in z direction
1468: . -da_refine_x <rx> - refinement ratio in x direction
1469: . -da_refine_y <ry> - refinement ratio in y direction
1470: . -da_refine_z <rz>- refinement ratio in z directio
1471: - -da_refine <n> - refine the DMDA n times before creating it
1473: Level: beginner
1475: Notes:
1476: The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1477: standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1478: the standard 27-pt stencil.
1480: The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
1481: The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
1482: and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed.
1484: You must call `DMSetUp()` after this call before using this `DM`.
1486: If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1487: but before `DMSetUp()`.
1489: .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
1490: `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
1491: `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
1492: `DMStagCreate3d()`
1493: @*/
1494: PetscErrorCode DMDACreate3d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt P, PetscInt m, PetscInt n, PetscInt p, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[], DM *da)
1495: {
1496: DMDACreate(comm, da);
1497: DMSetDimension(*da, 3);
1498: DMDASetSizes(*da, M, N, P);
1499: DMDASetNumProcs(*da, m, n, p);
1500: DMDASetBoundaryType(*da, bx, by, bz);
1501: DMDASetDof(*da, dof);
1502: DMDASetStencilType(*da, stencil_type);
1503: DMDASetStencilWidth(*da, s);
1504: DMDASetOwnershipRanges(*da, lx, ly, lz);
1505: return 0;
1506: }