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, &gtol);
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: }