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:   PetscFunctionBegin;
 21:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));

 23:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 24:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
 25:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
 26:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
 27: #if defined(PETSC_HAVE_MATLAB)
 28:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
 29: #endif
 30:   if (iascii) {
 31:     PetscViewerFormat format;

 33:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
 34:     PetscCall(PetscViewerGetFormat(viewer, &format));
 35:     if (format == PETSC_VIEWER_LOAD_BALANCE) {
 36:       PetscInt      i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal;
 37:       DMDALocalInfo info;
 38:       PetscMPIInt   size;
 39:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
 40:       PetscCall(DMDAGetLocalInfo(da, &info));
 41:       nzlocal = info.xm * info.ym * info.zm;
 42:       PetscCall(PetscMalloc1(size, &nz));
 43:       PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da)));
 44:       for (i = 0; i < (PetscInt)size; i++) {
 45:         nmax = PetscMax(nmax, nz[i]);
 46:         nmin = PetscMin(nmin, nz[i]);
 47:         navg += nz[i];
 48:       }
 49:       PetscCall(PetscFree(nz));
 50:       navg = navg / size;
 51:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Load Balance - Grid Points: Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax));
 52:       PetscFunctionReturn(PETSC_SUCCESS);
 53:     }
 54:     if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) {
 55:       DMDALocalInfo info;
 56:       PetscCall(DMDAGetLocalInfo(da, &info));
 57:       PetscCall(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));
 58:       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,
 59:                                                    info.xs + info.xm, info.ys, info.ys + info.ym, info.zs, info.zs + info.zm));
 60:       PetscCall(DMGetCoordinates(da, &coordinates));
 61: #if !defined(PETSC_USE_COMPLEX)
 62:       if (coordinates) {
 63:         PetscInt         last;
 64:         const PetscReal *coors;
 65:         PetscCall(VecGetArrayRead(coordinates, &coors));
 66:         PetscCall(VecGetLocalSize(coordinates, &last));
 67:         last = last - 3;
 68:         PetscCall(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]));
 69:         PetscCall(VecRestoreArrayRead(coordinates, &coors));
 70:       }
 71: #endif
 72:       PetscCall(PetscViewerFlush(viewer));
 73:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
 74:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer));
 75:     else PetscCall(DMView_DA_VTK(da, viewer));
 76:   } else if (isdraw) {
 77:     PetscDraw       draw;
 78:     PetscReal       ymin = -1.0, ymax = (PetscReal)dd->N;
 79:     PetscReal       xmin = -1.0, xmax = (PetscReal)((dd->M + 2) * dd->P), x, y, ycoord, xcoord;
 80:     PetscInt        k, plane, base;
 81:     const PetscInt *idx;
 82:     char            node[10];
 83:     PetscBool       isnull;

 85:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
 86:     PetscCall(PetscDrawIsNull(draw, &isnull));
 87:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

 89:     PetscCall(PetscDrawCheckResizedWindow(draw));
 90:     PetscCall(PetscDrawClear(draw));
 91:     PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));

 93:     PetscDrawCollectiveBegin(draw);
 94:     /* first processor draw all node lines */
 95:     if (rank == 0) {
 96:       for (k = 0; k < dd->P; k++) {
 97:         ymin = 0.0;
 98:         ymax = (PetscReal)(dd->N - 1);
 99:         for (xmin = (PetscReal)(k * (dd->M + 1)); xmin < (PetscReal)(dd->M + (k * (dd->M + 1))); xmin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK));
100:         xmin = (PetscReal)(k * (dd->M + 1));
101:         xmax = xmin + (PetscReal)(dd->M - 1);
102:         for (ymin = 0; ymin < (PetscReal)dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
103:       }
104:     }
105:     PetscDrawCollectiveEnd(draw);
106:     PetscCall(PetscDrawFlush(draw));
107:     PetscCall(PetscDrawPause(draw));

109:     PetscDrawCollectiveBegin(draw);
110:     /*Go through and draw for each plane*/
111:     for (k = 0; k < dd->P; k++) {
112:       if ((k >= dd->zs) && (k < dd->ze)) {
113:         /* draw my box */
114:         ymin = dd->ys;
115:         ymax = dd->ye - 1;
116:         xmin = dd->xs / dd->w + (dd->M + 1) * k;
117:         xmax = (dd->xe - 1) / dd->w + (dd->M + 1) * k;

119:         PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
120:         PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
121:         PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
122:         PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED));

124:         xmin = dd->xs / dd->w;
125:         xmax = (dd->xe - 1) / dd->w;

127:         /* identify which processor owns the box */
128:         PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)rank));
129:         PetscCall(PetscDrawString(draw, xmin + (dd->M + 1) * k + .2, ymin + .3, PETSC_DRAW_RED, node));
130:         /* put in numbers*/
131:         base = (dd->base + (dd->xe - dd->xs) * (dd->ye - dd->ys) * (k - dd->zs)) / dd->w;
132:         for (y = ymin; y <= ymax; y++) {
133:           for (x = xmin + (dd->M + 1) * k; x <= xmax + (dd->M + 1) * k; x++) {
134:             PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++));
135:             PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
136:           }
137:         }
138:       }
139:     }
140:     PetscDrawCollectiveEnd(draw);
141:     PetscCall(PetscDrawFlush(draw));
142:     PetscCall(PetscDrawPause(draw));

144:     PetscDrawCollectiveBegin(draw);
145:     for (k = 0 - dd->s; k < dd->P + dd->s; k++) {
146:       /* Go through and draw for each plane */
147:       if ((k >= dd->Zs) && (k < dd->Ze)) {
148:         /* overlay ghost numbers, useful for error checking */
149:         base = (dd->Xe - dd->Xs) * (dd->Ye - dd->Ys) * (k - dd->Zs) / dd->w;
150:         PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx));
151:         plane = k;
152:         /* Keep z wrap around points on the drawing */
153:         if (k < 0) plane = dd->P + k;
154:         if (k >= dd->P) plane = k - dd->P;
155:         ymin = dd->Ys;
156:         ymax = dd->Ye;
157:         xmin = (dd->M + 1) * plane * dd->w;
158:         xmax = (dd->M + 1) * plane * dd->w + dd->M * dd->w;
159:         for (y = ymin; y < ymax; y++) {
160:           for (x = xmin + dd->Xs; x < xmin + dd->Xe; x += dd->w) {
161:             PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, idx[base]));
162:             ycoord = y;
163:             /*Keep y wrap around points on drawing */
164:             if (y < 0) ycoord = dd->N + y;
165:             if (y >= dd->N) ycoord = y - dd->N;
166:             xcoord = x; /* Keep x wrap points on drawing */
167:             if (x < xmin) xcoord = xmax - (xmin - x);
168:             if (x >= xmax) xcoord = xmin + (x - xmax);
169:             PetscCall(PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node));
170:             base++;
171:           }
172:         }
173:         PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
174:       }
175:     }
176:     PetscDrawCollectiveEnd(draw);
177:     PetscCall(PetscDrawFlush(draw));
178:     PetscCall(PetscDrawPause(draw));
179:     PetscCall(PetscDrawSave(draw));
180:   } else if (isglvis) {
181:     PetscCall(DMView_DA_GLVis(da, viewer));
182:   } else if (isbinary) {
183:     PetscCall(DMView_DA_Binary(da, viewer));
184: #if defined(PETSC_HAVE_MATLAB)
185:   } else if (ismatlab) {
186:     PetscCall(DMView_DA_Matlab(da, viewer));
187: #endif
188:   }
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

192: PetscErrorCode DMSetUp_DA_3D(DM da)
193: {
194:   DM_DA          *dd           = (DM_DA *)da->data;
195:   const PetscInt  M            = dd->M;
196:   const PetscInt  N            = dd->N;
197:   const PetscInt  P            = dd->P;
198:   PetscInt        m            = dd->m;
199:   PetscInt        n            = dd->n;
200:   PetscInt        p            = dd->p;
201:   const PetscInt  dof          = dd->w;
202:   const PetscInt  s            = dd->s;
203:   DMBoundaryType  bx           = dd->bx;
204:   DMBoundaryType  by           = dd->by;
205:   DMBoundaryType  bz           = dd->bz;
206:   DMDAStencilType stencil_type = dd->stencil_type;
207:   PetscInt       *lx           = dd->lx;
208:   PetscInt       *ly           = dd->ly;
209:   PetscInt       *lz           = dd->lz;
210:   MPI_Comm        comm;
211:   PetscMPIInt     rank, size;
212:   PetscInt        xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0;
213:   PetscInt        Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm;
214:   PetscInt        left, right, up, down, bottom, top, i, j, k, *idx, nn;
215:   PetscInt        n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14;
216:   PetscInt        n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26;
217:   PetscInt       *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z;
218:   PetscInt        sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0;
219:   PetscInt        sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0;
220:   PetscInt        sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0;
221:   Vec             local, global;
222:   VecScatter      gtol;
223:   IS              to, from;
224:   PetscBool       twod;

226:   PetscFunctionBegin;
227:   PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && by != DM_BOUNDARY_MIRROR && bz != DM_BOUNDARY_MIRROR), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box stencil");
228:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
229: #if !defined(PETSC_USE_64BIT_INDICES)
230:   PetscCheck(((PetscInt64)M) * ((PetscInt64)N) * ((PetscInt64)P) * ((PetscInt64)dof) <= (PetscInt64)PETSC_MPI_INT_MAX, comm, PETSC_ERR_INT_OVERFLOW, "Mesh of %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " (dof) is too large for 32-bit indices", M, N, P, dof);
231: #endif

233:   PetscCallMPI(MPI_Comm_size(comm, &size));
234:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

236:   if (m != PETSC_DECIDE) {
237:     PetscCheck(m >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m);
238:     PetscCheck(m <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %d", m, size);
239:   }
240:   if (n != PETSC_DECIDE) {
241:     PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %" PetscInt_FMT, n);
242:     PetscCheck(n <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %d", n, size);
243:   }
244:   if (p != PETSC_DECIDE) {
245:     PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Z direction: %" PetscInt_FMT, p);
246:     PetscCheck(p <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Z direction: %" PetscInt_FMT " %d", p, size);
247:   }
248:   PetscCheck(m <= 0 || n <= 0 || p <= 0 || m * n * p == size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m %" PetscInt_FMT " * n %" PetscInt_FMT " * p %" PetscInt_FMT " != size %d", m, n, p, size);

250:   /* Partition the array among the processors */
251:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
252:     m = size / (n * p);
253:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
254:     n = size / (m * p);
255:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
256:     p = size / (m * n);
257:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
258:     /* try for squarish distribution */
259:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N * p)));
260:     if (!m) m = 1;
261:     while (m > 0) {
262:       n = size / (m * p);
263:       if (m * n * p == size) break;
264:       m--;
265:     }
266:     PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad p value: p = %" PetscInt_FMT, p);
267:     if (M > N && m < n) {
268:       PetscInt _m = m;
269:       m           = n;
270:       n           = _m;
271:     }
272:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
273:     /* try for squarish distribution */
274:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
275:     if (!m) m = 1;
276:     while (m > 0) {
277:       p = size / (m * n);
278:       if (m * n * p == size) break;
279:       m--;
280:     }
281:     PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad n value: n = %" PetscInt_FMT, n);
282:     if (M > P && m < p) {
283:       PetscInt _m = m;
284:       m           = p;
285:       p           = _m;
286:     }
287:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
288:     /* try for squarish distribution */
289:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N) * ((PetscReal)size) / ((PetscReal)P * m)));
290:     if (!n) n = 1;
291:     while (n > 0) {
292:       p = size / (m * n);
293:       if (m * n * p == size) break;
294:       n--;
295:     }
296:     PetscCheck(n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad m value: m = %" PetscInt_FMT, n);
297:     if (N > P && n < p) {
298:       PetscInt _n = n;
299:       n           = p;
300:       p           = _n;
301:     }
302:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
303:     /* try for squarish distribution */
304:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.)));
305:     if (!n) n = 1;
306:     while (n > 0) {
307:       pm = size / n;
308:       if (n * pm == size) break;
309:       n--;
310:     }
311:     if (!n) n = 1;
312:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
313:     if (!m) m = 1;
314:     while (m > 0) {
315:       p = size / (m * n);
316:       if (m * n * p == size) break;
317:       m--;
318:     }
319:     if (M > P && m < p) {
320:       PetscInt _m = m;
321:       m           = p;
322:       p           = _m;
323:     }
324:   } else PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition");

326:   PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_PLIB, "Could not find good partition");
327:   PetscCheck(M >= m, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m);
328:   PetscCheck(N >= n, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n);
329:   PetscCheck(P >= p, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in z direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, P, p);

331:   /*
332:      Determine locally owned region
333:      [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
334:   */

336:   if (!lx) {
337:     PetscCall(PetscMalloc1(m, &dd->lx));
338:     lx = dd->lx;
339:     for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m));
340:   }
341:   x  = lx[rank % m];
342:   xs = 0;
343:   for (i = 0; i < (rank % m); i++) xs += lx[i];
344:   PetscCheck(x >= s || (m <= 1 && bx != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, x, s);

346:   if (!ly) {
347:     PetscCall(PetscMalloc1(n, &dd->ly));
348:     ly = dd->ly;
349:     for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n));
350:   }
351:   y = ly[(rank % (m * n)) / m];
352:   PetscCheck(y >= s || (n <= 1 && by != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, y, s);

354:   ys = 0;
355:   for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i];

357:   if (!lz) {
358:     PetscCall(PetscMalloc1(p, &dd->lz));
359:     lz = dd->lz;
360:     for (i = 0; i < p; i++) lz[i] = P / p + ((P % p) > (i % p));
361:   }
362:   z = lz[rank / (m * n)];

364:   /* note this is different than x- and y-, as we will handle as an important special
365:    case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z.  This is to deal with 2D problems
366:    in a 3D code.  Additional code for this case is noted with "2d case" comments */
367:   twod = PETSC_FALSE;
368:   if (P == 1) twod = PETSC_TRUE;
369:   else PetscCheck(z >= s || (p <= 1 && bz != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local z-width of domain z %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, z, s);
370:   zs = 0;
371:   for (i = 0; i < (rank / (m * n)); i++) zs += lz[i];
372:   ye = ys + y;
373:   xe = xs + x;
374:   ze = zs + z;

376:   /* determine ghost region (Xs) and region scattered into (IXs)  */
377:   if (xs - s > 0) {
378:     Xs  = xs - s;
379:     IXs = xs - s;
380:   } else {
381:     if (bx) Xs = xs - s;
382:     else Xs = 0;
383:     IXs = 0;
384:   }
385:   if (xe + s <= M) {
386:     Xe  = xe + s;
387:     IXe = xe + s;
388:   } else {
389:     if (bx) {
390:       Xs = xs - s;
391:       Xe = xe + s;
392:     } else Xe = M;
393:     IXe = M;
394:   }

396:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
397:     IXs = xs - s;
398:     IXe = xe + s;
399:     Xs  = xs - s;
400:     Xe  = xe + s;
401:   }

403:   if (ys - s > 0) {
404:     Ys  = ys - s;
405:     IYs = ys - s;
406:   } else {
407:     if (by) Ys = ys - s;
408:     else Ys = 0;
409:     IYs = 0;
410:   }
411:   if (ye + s <= N) {
412:     Ye  = ye + s;
413:     IYe = ye + s;
414:   } else {
415:     if (by) Ye = ye + s;
416:     else Ye = N;
417:     IYe = N;
418:   }

420:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
421:     IYs = ys - s;
422:     IYe = ye + s;
423:     Ys  = ys - s;
424:     Ye  = ye + s;
425:   }

427:   if (zs - s > 0) {
428:     Zs  = zs - s;
429:     IZs = zs - s;
430:   } else {
431:     if (bz) Zs = zs - s;
432:     else Zs = 0;
433:     IZs = 0;
434:   }
435:   if (ze + s <= P) {
436:     Ze  = ze + s;
437:     IZe = ze + s;
438:   } else {
439:     if (bz) Ze = ze + s;
440:     else Ze = P;
441:     IZe = P;
442:   }

444:   if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
445:     IZs = zs - s;
446:     IZe = ze + s;
447:     Zs  = zs - s;
448:     Ze  = ze + s;
449:   }

451:   /* Resize all X parameters to reflect w */
452:   s_x = s;
453:   s_y = s;
454:   s_z = s;

456:   /* determine starting point of each processor */
457:   nn = x * y * z;
458:   PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims));
459:   PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm));
460:   bases[0] = 0;
461:   for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
462:   for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
463:   base = bases[rank] * dof;

465:   /* allocate the base parallel and sequential vectors */
466:   dd->Nlocal = x * y * z * dof;
467:   PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global));
468:   dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof;
469:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local));

471:   /* generate global to local vector scatter and local to global mapping*/

473:   /* global to local must include ghost points within the domain,
474:      but not ghost points outside the domain that aren't periodic */
475:   PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx));
476:   if (stencil_type == DMDA_STENCIL_BOX) {
477:     left   = IXs - Xs;
478:     right  = left + (IXe - IXs);
479:     bottom = IYs - Ys;
480:     top    = bottom + (IYe - IYs);
481:     down   = IZs - Zs;
482:     up     = down + (IZe - IZs);
483:     count  = 0;
484:     for (i = down; i < up; i++) {
485:       for (j = bottom; j < top; j++) {
486:         for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
487:       }
488:     }
489:     PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
490:   } else {
491:     /* This is way ugly! We need to list the funny cross type region */
492:     left   = xs - Xs;
493:     right  = left + x;
494:     bottom = ys - Ys;
495:     top    = bottom + y;
496:     down   = zs - Zs;
497:     up     = down + z;
498:     count  = 0;
499:     /* the bottom chunk */
500:     for (i = (IZs - Zs); i < down; i++) {
501:       for (j = bottom; j < top; j++) {
502:         for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
503:       }
504:     }
505:     /* the middle piece */
506:     for (i = down; i < up; i++) {
507:       /* front */
508:       for (j = (IYs - Ys); j < bottom; j++) {
509:         for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
510:       }
511:       /* middle */
512:       for (j = bottom; j < top; j++) {
513:         for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
514:       }
515:       /* back */
516:       for (j = top; j < top + IYe - ye; j++) {
517:         for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
518:       }
519:     }
520:     /* the top piece */
521:     for (i = up; i < up + IZe - ze; i++) {
522:       for (j = bottom; j < top; j++) {
523:         for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
524:       }
525:     }
526:     PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
527:   }

529:   /* determine who lies on each side of use stored in    n24 n25 n26
530:                                                          n21 n22 n23
531:                                                          n18 n19 n20

533:                                                          n15 n16 n17
534:                                                          n12     n14
535:                                                          n9  n10 n11

537:                                                          n6  n7  n8
538:                                                          n3  n4  n5
539:                                                          n0  n1  n2
540:   */

542:   /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
543:   /* Assume Nodes are Internal to the Cube */
544:   n0 = rank - m * n - m - 1;
545:   n1 = rank - m * n - m;
546:   n2 = rank - m * n - m + 1;
547:   n3 = rank - m * n - 1;
548:   n4 = rank - m * n;
549:   n5 = rank - m * n + 1;
550:   n6 = rank - m * n + m - 1;
551:   n7 = rank - m * n + m;
552:   n8 = rank - m * n + m + 1;

554:   n9  = rank - m - 1;
555:   n10 = rank - m;
556:   n11 = rank - m + 1;
557:   n12 = rank - 1;
558:   n14 = rank + 1;
559:   n15 = rank + m - 1;
560:   n16 = rank + m;
561:   n17 = rank + m + 1;

563:   n18 = rank + m * n - m - 1;
564:   n19 = rank + m * n - m;
565:   n20 = rank + m * n - m + 1;
566:   n21 = rank + m * n - 1;
567:   n22 = rank + m * n;
568:   n23 = rank + m * n + 1;
569:   n24 = rank + m * n + m - 1;
570:   n25 = rank + m * n + m;
571:   n26 = rank + m * n + m + 1;

573:   /* Assume Pieces are on Faces of Cube */

575:   if (xs == 0) { /* First assume not corner or edge */
576:     n0  = rank - 1 - (m * n);
577:     n3  = rank + m - 1 - (m * n);
578:     n6  = rank + 2 * m - 1 - (m * n);
579:     n9  = rank - 1;
580:     n12 = rank + m - 1;
581:     n15 = rank + 2 * m - 1;
582:     n18 = rank - 1 + (m * n);
583:     n21 = rank + m - 1 + (m * n);
584:     n24 = rank + 2 * m - 1 + (m * n);
585:   }

587:   if (xe == M) { /* First assume not corner or edge */
588:     n2  = rank - 2 * m + 1 - (m * n);
589:     n5  = rank - m + 1 - (m * n);
590:     n8  = rank + 1 - (m * n);
591:     n11 = rank - 2 * m + 1;
592:     n14 = rank - m + 1;
593:     n17 = rank + 1;
594:     n20 = rank - 2 * m + 1 + (m * n);
595:     n23 = rank - m + 1 + (m * n);
596:     n26 = rank + 1 + (m * n);
597:   }

599:   if (ys == 0) { /* First assume not corner or edge */
600:     n0  = rank + m * (n - 1) - 1 - (m * n);
601:     n1  = rank + m * (n - 1) - (m * n);
602:     n2  = rank + m * (n - 1) + 1 - (m * n);
603:     n9  = rank + m * (n - 1) - 1;
604:     n10 = rank + m * (n - 1);
605:     n11 = rank + m * (n - 1) + 1;
606:     n18 = rank + m * (n - 1) - 1 + (m * n);
607:     n19 = rank + m * (n - 1) + (m * n);
608:     n20 = rank + m * (n - 1) + 1 + (m * n);
609:   }

611:   if (ye == N) { /* First assume not corner or edge */
612:     n6  = rank - m * (n - 1) - 1 - (m * n);
613:     n7  = rank - m * (n - 1) - (m * n);
614:     n8  = rank - m * (n - 1) + 1 - (m * n);
615:     n15 = rank - m * (n - 1) - 1;
616:     n16 = rank - m * (n - 1);
617:     n17 = rank - m * (n - 1) + 1;
618:     n24 = rank - m * (n - 1) - 1 + (m * n);
619:     n25 = rank - m * (n - 1) + (m * n);
620:     n26 = rank - m * (n - 1) + 1 + (m * n);
621:   }

623:   if (zs == 0) { /* First assume not corner or edge */
624:     n0 = size - (m * n) + rank - m - 1;
625:     n1 = size - (m * n) + rank - m;
626:     n2 = size - (m * n) + rank - m + 1;
627:     n3 = size - (m * n) + rank - 1;
628:     n4 = size - (m * n) + rank;
629:     n5 = size - (m * n) + rank + 1;
630:     n6 = size - (m * n) + rank + m - 1;
631:     n7 = size - (m * n) + rank + m;
632:     n8 = size - (m * n) + rank + m + 1;
633:   }

635:   if (ze == P) { /* First assume not corner or edge */
636:     n18 = (m * n) - (size - rank) - m - 1;
637:     n19 = (m * n) - (size - rank) - m;
638:     n20 = (m * n) - (size - rank) - m + 1;
639:     n21 = (m * n) - (size - rank) - 1;
640:     n22 = (m * n) - (size - rank);
641:     n23 = (m * n) - (size - rank) + 1;
642:     n24 = (m * n) - (size - rank) + m - 1;
643:     n25 = (m * n) - (size - rank) + m;
644:     n26 = (m * n) - (size - rank) + m + 1;
645:   }

647:   if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */
648:     n0 = size - m * n + rank + m - 1 - m;
649:     n3 = size - m * n + rank + m - 1;
650:     n6 = size - m * n + rank + m - 1 + m;
651:   }

653:   if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */
654:     n18 = m * n - (size - rank) + m - 1 - m;
655:     n21 = m * n - (size - rank) + m - 1;
656:     n24 = m * n - (size - rank) + m - 1 + m;
657:   }

659:   if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */
660:     n0  = rank + m * n - 1 - m * n;
661:     n9  = rank + m * n - 1;
662:     n18 = rank + m * n - 1 + m * n;
663:   }

665:   if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */
666:     n6  = rank - m * (n - 1) + m - 1 - m * n;
667:     n15 = rank - m * (n - 1) + m - 1;
668:     n24 = rank - m * (n - 1) + m - 1 + m * n;
669:   }

671:   if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */
672:     n2 = size - (m * n - rank) - (m - 1) - m;
673:     n5 = size - (m * n - rank) - (m - 1);
674:     n8 = size - (m * n - rank) - (m - 1) + m;
675:   }

677:   if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */
678:     n20 = m * n - (size - rank) - (m - 1) - m;
679:     n23 = m * n - (size - rank) - (m - 1);
680:     n26 = m * n - (size - rank) - (m - 1) + m;
681:   }

683:   if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */
684:     n2  = rank + m * (n - 1) - (m - 1) - m * n;
685:     n11 = rank + m * (n - 1) - (m - 1);
686:     n20 = rank + m * (n - 1) - (m - 1) + m * n;
687:   }

689:   if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */
690:     n8  = rank - m * n + 1 - m * n;
691:     n17 = rank - m * n + 1;
692:     n26 = rank - m * n + 1 + m * n;
693:   }

695:   if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */
696:     n0 = size - m + rank - 1;
697:     n1 = size - m + rank;
698:     n2 = size - m + rank + 1;
699:   }

701:   if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */
702:     n18 = m * n - (size - rank) + m * (n - 1) - 1;
703:     n19 = m * n - (size - rank) + m * (n - 1);
704:     n20 = m * n - (size - rank) + m * (n - 1) + 1;
705:   }

707:   if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */
708:     n6 = size - (m * n - rank) - m * (n - 1) - 1;
709:     n7 = size - (m * n - rank) - m * (n - 1);
710:     n8 = size - (m * n - rank) - m * (n - 1) + 1;
711:   }

713:   if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */
714:     n24 = rank - (size - m) - 1;
715:     n25 = rank - (size - m);
716:     n26 = rank - (size - m) + 1;
717:   }

719:   /* Check for Corners */
720:   if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1;
721:   if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1;
722:   if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1);
723:   if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1;
724:   if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m;
725:   if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m;
726:   if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n;
727:   if ((xe == M) && (ye == N) && (ze == P)) n26 = 0;

729:   /* Check for when not X,Y, and Z Periodic */

731:   /* If not X periodic */
732:   if (bx != DM_BOUNDARY_PERIODIC) {
733:     if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;
734:     if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;
735:   }

737:   /* If not Y periodic */
738:   if (by != DM_BOUNDARY_PERIODIC) {
739:     if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;
740:     if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;
741:   }

743:   /* If not Z periodic */
744:   if (bz != DM_BOUNDARY_PERIODIC) {
745:     if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;
746:     if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;
747:   }

749:   PetscCall(PetscMalloc1(27, &dd->neighbors));

751:   dd->neighbors[0]  = n0;
752:   dd->neighbors[1]  = n1;
753:   dd->neighbors[2]  = n2;
754:   dd->neighbors[3]  = n3;
755:   dd->neighbors[4]  = n4;
756:   dd->neighbors[5]  = n5;
757:   dd->neighbors[6]  = n6;
758:   dd->neighbors[7]  = n7;
759:   dd->neighbors[8]  = n8;
760:   dd->neighbors[9]  = n9;
761:   dd->neighbors[10] = n10;
762:   dd->neighbors[11] = n11;
763:   dd->neighbors[12] = n12;
764:   dd->neighbors[13] = rank;
765:   dd->neighbors[14] = n14;
766:   dd->neighbors[15] = n15;
767:   dd->neighbors[16] = n16;
768:   dd->neighbors[17] = n17;
769:   dd->neighbors[18] = n18;
770:   dd->neighbors[19] = n19;
771:   dd->neighbors[20] = n20;
772:   dd->neighbors[21] = n21;
773:   dd->neighbors[22] = n22;
774:   dd->neighbors[23] = n23;
775:   dd->neighbors[24] = n24;
776:   dd->neighbors[25] = n25;
777:   dd->neighbors[26] = n26;

779:   /* If star stencil then delete the corner neighbors */
780:   if (stencil_type == DMDA_STENCIL_STAR) {
781:     /* save information about corner neighbors */
782:     sn0  = n0;
783:     sn1  = n1;
784:     sn2  = n2;
785:     sn3  = n3;
786:     sn5  = n5;
787:     sn6  = n6;
788:     sn7  = n7;
789:     sn8  = n8;
790:     sn9  = n9;
791:     sn11 = n11;
792:     sn15 = n15;
793:     sn17 = n17;
794:     sn18 = n18;
795:     sn19 = n19;
796:     sn20 = n20;
797:     sn21 = n21;
798:     sn23 = n23;
799:     sn24 = n24;
800:     sn25 = n25;
801:     sn26 = n26;
802:     n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
803:   }

805:   PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx));

807:   nn = 0;
808:   /* Bottom Level */
809:   for (k = 0; k < s_z; k++) {
810:     for (i = 1; i <= s_y; i++) {
811:       if (n0 >= 0) { /* left below */
812:         x_t = lx[n0 % m];
813:         y_t = ly[(n0 % (m * n)) / m];
814:         z_t = lz[n0 / (m * n)];
815:         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;
816:         if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */
817:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
818:       }
819:       if (n1 >= 0) { /* directly below */
820:         x_t = x;
821:         y_t = ly[(n1 % (m * n)) / m];
822:         z_t = lz[n1 / (m * n)];
823:         s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
824:         if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
825:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
826:       }
827:       if (n2 >= 0) { /* right below */
828:         x_t = lx[n2 % m];
829:         y_t = ly[(n2 % (m * n)) / m];
830:         z_t = lz[n2 / (m * n)];
831:         s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
832:         if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
833:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
834:       }
835:     }

837:     for (i = 0; i < y; i++) {
838:       if (n3 >= 0) { /* directly left */
839:         x_t = lx[n3 % m];
840:         y_t = y;
841:         z_t = lz[n3 / (m * n)];
842:         s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
843:         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 */
844:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
845:       }

847:       if (n4 >= 0) { /* middle */
848:         x_t = x;
849:         y_t = y;
850:         z_t = lz[n4 / (m * n)];
851:         s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
852:         if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
853:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
854:       } else if (bz == DM_BOUNDARY_MIRROR) {
855:         for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
856:       }

858:       if (n5 >= 0) { /* directly right */
859:         x_t = lx[n5 % m];
860:         y_t = y;
861:         z_t = lz[n5 / (m * n)];
862:         s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
863:         if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
864:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
865:       }
866:     }

868:     for (i = 1; i <= s_y; i++) {
869:       if (n6 >= 0) { /* left above */
870:         x_t = lx[n6 % m];
871:         y_t = ly[(n6 % (m * n)) / m];
872:         z_t = lz[n6 / (m * n)];
873:         s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
874:         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 */
875:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
876:       }
877:       if (n7 >= 0) { /* directly above */
878:         x_t = x;
879:         y_t = ly[(n7 % (m * n)) / m];
880:         z_t = lz[n7 / (m * n)];
881:         s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
882:         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 */
883:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
884:       }
885:       if (n8 >= 0) { /* right above */
886:         x_t = lx[n8 % m];
887:         y_t = ly[(n8 % (m * n)) / m];
888:         z_t = lz[n8 / (m * n)];
889:         s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
890:         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 */
891:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
892:       }
893:     }
894:   }

896:   /* Middle Level */
897:   for (k = 0; k < z; k++) {
898:     for (i = 1; i <= s_y; i++) {
899:       if (n9 >= 0) { /* left below */
900:         x_t = lx[n9 % m];
901:         y_t = ly[(n9 % (m * n)) / m];
902:         /* z_t = z; */
903:         s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
904:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
905:       }
906:       if (n10 >= 0) { /* directly below */
907:         x_t = x;
908:         y_t = ly[(n10 % (m * n)) / m];
909:         /* z_t = z; */
910:         s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
911:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
912:       } else if (by == DM_BOUNDARY_MIRROR) {
913:         for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
914:       }
915:       if (n11 >= 0) { /* right below */
916:         x_t = lx[n11 % m];
917:         y_t = ly[(n11 % (m * n)) / m];
918:         /* z_t = z; */
919:         s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
920:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
921:       }
922:     }

924:     for (i = 0; i < y; i++) {
925:       if (n12 >= 0) { /* directly left */
926:         x_t = lx[n12 % m];
927:         y_t = y;
928:         /* z_t = z; */
929:         s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
930:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
931:       } else if (bx == DM_BOUNDARY_MIRROR) {
932:         for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
933:       }

935:       /* Interior */
936:       s_t = bases[rank] + i * x + k * x * y;
937:       for (j = 0; j < x; j++) idx[nn++] = s_t++;

939:       if (n14 >= 0) { /* directly right */
940:         x_t = lx[n14 % m];
941:         y_t = y;
942:         /* z_t = z; */
943:         s_t = bases[n14] + i * x_t + k * x_t * y_t;
944:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
945:       } else if (bx == DM_BOUNDARY_MIRROR) {
946:         for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
947:       }
948:     }

950:     for (i = 1; i <= s_y; i++) {
951:       if (n15 >= 0) { /* left above */
952:         x_t = lx[n15 % m];
953:         y_t = ly[(n15 % (m * n)) / m];
954:         /* z_t = z; */
955:         s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
956:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
957:       }
958:       if (n16 >= 0) { /* directly above */
959:         x_t = x;
960:         y_t = ly[(n16 % (m * n)) / m];
961:         /* z_t = z; */
962:         s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
963:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
964:       } else if (by == DM_BOUNDARY_MIRROR) {
965:         for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
966:       }
967:       if (n17 >= 0) { /* right above */
968:         x_t = lx[n17 % m];
969:         y_t = ly[(n17 % (m * n)) / m];
970:         /* z_t = z; */
971:         s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
972:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
973:       }
974:     }
975:   }

977:   /* Upper Level */
978:   for (k = 0; k < s_z; k++) {
979:     for (i = 1; i <= s_y; i++) {
980:       if (n18 >= 0) { /* left below */
981:         x_t = lx[n18 % m];
982:         y_t = ly[(n18 % (m * n)) / m];
983:         /* z_t = lz[n18 / (m*n)]; */
984:         s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
985:         if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */
986:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
987:       }
988:       if (n19 >= 0) { /* directly below */
989:         x_t = x;
990:         y_t = ly[(n19 % (m * n)) / m];
991:         /* z_t = lz[n19 / (m*n)]; */
992:         s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
993:         if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
994:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
995:       }
996:       if (n20 >= 0) { /* right below */
997:         x_t = lx[n20 % m];
998:         y_t = ly[(n20 % (m * n)) / m];
999:         /* z_t = lz[n20 / (m*n)]; */
1000:         s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1001:         if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
1002:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1003:       }
1004:     }

1006:     for (i = 0; i < y; i++) {
1007:       if (n21 >= 0) { /* directly left */
1008:         x_t = lx[n21 % m];
1009:         y_t = y;
1010:         /* z_t = lz[n21 / (m*n)]; */
1011:         s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1012:         if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */
1013:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1014:       }

1016:       if (n22 >= 0) { /* middle */
1017:         x_t = x;
1018:         y_t = y;
1019:         /* z_t = lz[n22 / (m*n)]; */
1020:         s_t = bases[n22] + i * x_t + k * x_t * y_t;
1021:         if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */
1022:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1023:       } else if (bz == DM_BOUNDARY_MIRROR) {
1024:         for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1025:       }

1027:       if (n23 >= 0) { /* directly right */
1028:         x_t = lx[n23 % m];
1029:         y_t = y;
1030:         /* z_t = lz[n23 / (m*n)]; */
1031:         s_t = bases[n23] + i * x_t + k * x_t * y_t;
1032:         if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */
1033:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1034:       }
1035:     }

1037:     for (i = 1; i <= s_y; i++) {
1038:       if (n24 >= 0) { /* left above */
1039:         x_t = lx[n24 % m];
1040:         y_t = ly[(n24 % (m * n)) / m];
1041:         /* z_t = lz[n24 / (m*n)]; */
1042:         s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1043:         if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */
1044:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1045:       }
1046:       if (n25 >= 0) { /* directly above */
1047:         x_t = x;
1048:         y_t = ly[(n25 % (m * n)) / m];
1049:         /* z_t = lz[n25 / (m*n)]; */
1050:         s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1051:         if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */
1052:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1053:       }
1054:       if (n26 >= 0) { /* right above */
1055:         x_t = lx[n26 % m];
1056:         y_t = ly[(n26 % (m * n)) / m];
1057:         /* z_t = lz[n26 / (m*n)]; */
1058:         s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1059:         if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */
1060:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1061:       }
1062:     }
1063:   }

1065:   PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
1066:   PetscCall(VecScatterCreate(global, from, local, to, &gtol));
1067:   PetscCall(ISDestroy(&to));
1068:   PetscCall(ISDestroy(&from));

1070:   if (stencil_type == DMDA_STENCIL_STAR) {
1071:     n0  = sn0;
1072:     n1  = sn1;
1073:     n2  = sn2;
1074:     n3  = sn3;
1075:     n5  = sn5;
1076:     n6  = sn6;
1077:     n7  = sn7;
1078:     n8  = sn8;
1079:     n9  = sn9;
1080:     n11 = sn11;
1081:     n15 = sn15;
1082:     n17 = sn17;
1083:     n18 = sn18;
1084:     n19 = sn19;
1085:     n20 = sn20;
1086:     n21 = sn21;
1087:     n23 = sn23;
1088:     n24 = sn24;
1089:     n25 = sn25;
1090:     n26 = sn26;
1091:   }

1093:   if (((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz))) {
1094:     /*
1095:         Recompute the local to global mappings, this time keeping the
1096:       information about the cross corner processor numbers.
1097:     */
1098:     nn = 0;
1099:     /* Bottom Level */
1100:     for (k = 0; k < s_z; k++) {
1101:       for (i = 1; i <= s_y; i++) {
1102:         if (n0 >= 0) { /* left below */
1103:           x_t = lx[n0 % m];
1104:           y_t = ly[(n0 % (m * n)) / m];
1105:           z_t = lz[n0 / (m * n)];
1106:           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;
1107:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1108:         } else if (Xs - xs < 0 && Ys - ys < 0 && Zs - zs < 0) {
1109:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1110:         }
1111:         if (n1 >= 0) { /* directly below */
1112:           x_t = x;
1113:           y_t = ly[(n1 % (m * n)) / m];
1114:           z_t = lz[n1 / (m * n)];
1115:           s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1116:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1117:         } else if (Ys - ys < 0 && Zs - zs < 0) {
1118:           for (j = 0; j < x; j++) idx[nn++] = -1;
1119:         }
1120:         if (n2 >= 0) { /* right below */
1121:           x_t = lx[n2 % m];
1122:           y_t = ly[(n2 % (m * n)) / m];
1123:           z_t = lz[n2 / (m * n)];
1124:           s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1125:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1126:         } else if (xe - Xe < 0 && Ys - ys < 0 && Zs - zs < 0) {
1127:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1128:         }
1129:       }

1131:       for (i = 0; i < y; i++) {
1132:         if (n3 >= 0) { /* directly left */
1133:           x_t = lx[n3 % m];
1134:           y_t = y;
1135:           z_t = lz[n3 / (m * n)];
1136:           s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1137:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1138:         } else if (Xs - xs < 0 && Zs - zs < 0) {
1139:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1140:         }

1142:         if (n4 >= 0) { /* middle */
1143:           x_t = x;
1144:           y_t = y;
1145:           z_t = lz[n4 / (m * n)];
1146:           s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1147:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1148:         } else if (Zs - zs < 0) {
1149:           if (bz == DM_BOUNDARY_MIRROR) {
1150:             for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
1151:           } else {
1152:             for (j = 0; j < x; j++) idx[nn++] = -1;
1153:           }
1154:         }

1156:         if (n5 >= 0) { /* directly right */
1157:           x_t = lx[n5 % m];
1158:           y_t = y;
1159:           z_t = lz[n5 / (m * n)];
1160:           s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1161:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1162:         } else if (xe - Xe < 0 && Zs - zs < 0) {
1163:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1164:         }
1165:       }

1167:       for (i = 1; i <= s_y; i++) {
1168:         if (n6 >= 0) { /* left above */
1169:           x_t = lx[n6 % m];
1170:           y_t = ly[(n6 % (m * n)) / m];
1171:           z_t = lz[n6 / (m * n)];
1172:           s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1173:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1174:         } else if (Xs - xs < 0 && ye - Ye < 0 && Zs - zs < 0) {
1175:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1176:         }
1177:         if (n7 >= 0) { /* directly above */
1178:           x_t = x;
1179:           y_t = ly[(n7 % (m * n)) / m];
1180:           z_t = lz[n7 / (m * n)];
1181:           s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1182:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1183:         } else if (ye - Ye < 0 && Zs - zs < 0) {
1184:           for (j = 0; j < x; j++) idx[nn++] = -1;
1185:         }
1186:         if (n8 >= 0) { /* right above */
1187:           x_t = lx[n8 % m];
1188:           y_t = ly[(n8 % (m * n)) / m];
1189:           z_t = lz[n8 / (m * n)];
1190:           s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1191:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1192:         } else if (xe - Xe < 0 && ye - Ye < 0 && Zs - zs < 0) {
1193:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1194:         }
1195:       }
1196:     }

1198:     /* Middle Level */
1199:     for (k = 0; k < z; k++) {
1200:       for (i = 1; i <= s_y; i++) {
1201:         if (n9 >= 0) { /* left below */
1202:           x_t = lx[n9 % m];
1203:           y_t = ly[(n9 % (m * n)) / m];
1204:           /* z_t = z; */
1205:           s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1206:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1207:         } else if (Xs - xs < 0 && Ys - ys < 0) {
1208:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1209:         }
1210:         if (n10 >= 0) { /* directly below */
1211:           x_t = x;
1212:           y_t = ly[(n10 % (m * n)) / m];
1213:           /* z_t = z; */
1214:           s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1215:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1216:         } else if (Ys - ys < 0) {
1217:           if (by == DM_BOUNDARY_MIRROR) {
1218:             for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
1219:           } else {
1220:             for (j = 0; j < x; j++) idx[nn++] = -1;
1221:           }
1222:         }
1223:         if (n11 >= 0) { /* right below */
1224:           x_t = lx[n11 % m];
1225:           y_t = ly[(n11 % (m * n)) / m];
1226:           /* z_t = z; */
1227:           s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1228:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1229:         } else if (xe - Xe < 0 && Ys - ys < 0) {
1230:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1231:         }
1232:       }

1234:       for (i = 0; i < y; i++) {
1235:         if (n12 >= 0) { /* directly left */
1236:           x_t = lx[n12 % m];
1237:           y_t = y;
1238:           /* z_t = z; */
1239:           s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
1240:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1241:         } else if (Xs - xs < 0) {
1242:           if (bx == DM_BOUNDARY_MIRROR) {
1243:             for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
1244:           } else {
1245:             for (j = 0; j < s_x; j++) idx[nn++] = -1;
1246:           }
1247:         }

1249:         /* Interior */
1250:         s_t = bases[rank] + i * x + k * x * y;
1251:         for (j = 0; j < x; j++) idx[nn++] = s_t++;

1253:         if (n14 >= 0) { /* directly right */
1254:           x_t = lx[n14 % m];
1255:           y_t = y;
1256:           /* z_t = z; */
1257:           s_t = bases[n14] + i * x_t + k * x_t * y_t;
1258:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1259:         } else if (xe - Xe < 0) {
1260:           if (bx == DM_BOUNDARY_MIRROR) {
1261:             for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
1262:           } else {
1263:             for (j = 0; j < s_x; j++) idx[nn++] = -1;
1264:           }
1265:         }
1266:       }

1268:       for (i = 1; i <= s_y; i++) {
1269:         if (n15 >= 0) { /* left above */
1270:           x_t = lx[n15 % m];
1271:           y_t = ly[(n15 % (m * n)) / m];
1272:           /* z_t = z; */
1273:           s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
1274:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1275:         } else if (Xs - xs < 0 && ye - Ye < 0) {
1276:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1277:         }
1278:         if (n16 >= 0) { /* directly above */
1279:           x_t = x;
1280:           y_t = ly[(n16 % (m * n)) / m];
1281:           /* z_t = z; */
1282:           s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
1283:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1284:         } else if (ye - Ye < 0) {
1285:           if (by == DM_BOUNDARY_MIRROR) {
1286:             for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
1287:           } else {
1288:             for (j = 0; j < x; j++) idx[nn++] = -1;
1289:           }
1290:         }
1291:         if (n17 >= 0) { /* right above */
1292:           x_t = lx[n17 % m];
1293:           y_t = ly[(n17 % (m * n)) / m];
1294:           /* z_t = z; */
1295:           s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
1296:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1297:         } else if (xe - Xe < 0 && ye - Ye < 0) {
1298:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1299:         }
1300:       }
1301:     }

1303:     /* Upper Level */
1304:     for (k = 0; k < s_z; k++) {
1305:       for (i = 1; i <= s_y; i++) {
1306:         if (n18 >= 0) { /* left below */
1307:           x_t = lx[n18 % m];
1308:           y_t = ly[(n18 % (m * n)) / m];
1309:           /* z_t = lz[n18 / (m*n)]; */
1310:           s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1311:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1312:         } else if (Xs - xs < 0 && Ys - ys < 0 && ze - Ze < 0) {
1313:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1314:         }
1315:         if (n19 >= 0) { /* directly below */
1316:           x_t = x;
1317:           y_t = ly[(n19 % (m * n)) / m];
1318:           /* z_t = lz[n19 / (m*n)]; */
1319:           s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1320:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1321:         } else if (Ys - ys < 0 && ze - Ze < 0) {
1322:           for (j = 0; j < x; j++) idx[nn++] = -1;
1323:         }
1324:         if (n20 >= 0) { /* right below */
1325:           x_t = lx[n20 % m];
1326:           y_t = ly[(n20 % (m * n)) / m];
1327:           /* z_t = lz[n20 / (m*n)]; */
1328:           s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1329:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1330:         } else if (xe - Xe < 0 && Ys - ys < 0 && ze - Ze < 0) {
1331:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1332:         }
1333:       }

1335:       for (i = 0; i < y; i++) {
1336:         if (n21 >= 0) { /* directly left */
1337:           x_t = lx[n21 % m];
1338:           y_t = y;
1339:           /* z_t = lz[n21 / (m*n)]; */
1340:           s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1341:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1342:         } else if (Xs - xs < 0 && ze - Ze < 0) {
1343:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1344:         }

1346:         if (n22 >= 0) { /* middle */
1347:           x_t = x;
1348:           y_t = y;
1349:           /* z_t = lz[n22 / (m*n)]; */
1350:           s_t = bases[n22] + i * x_t + k * x_t * y_t;
1351:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1352:         } else if (ze - Ze < 0) {
1353:           if (bz == DM_BOUNDARY_MIRROR) {
1354:             for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1355:           } else {
1356:             for (j = 0; j < x; j++) idx[nn++] = -1;
1357:           }
1358:         }

1360:         if (n23 >= 0) { /* directly right */
1361:           x_t = lx[n23 % m];
1362:           y_t = y;
1363:           /* z_t = lz[n23 / (m*n)]; */
1364:           s_t = bases[n23] + i * x_t + k * x_t * y_t;
1365:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1366:         } else if (xe - Xe < 0 && ze - Ze < 0) {
1367:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1368:         }
1369:       }

1371:       for (i = 1; i <= s_y; i++) {
1372:         if (n24 >= 0) { /* left above */
1373:           x_t = lx[n24 % m];
1374:           y_t = ly[(n24 % (m * n)) / m];
1375:           /* z_t = lz[n24 / (m*n)]; */
1376:           s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1377:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1378:         } else if (Xs - xs < 0 && ye - Ye < 0 && ze - Ze < 0) {
1379:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1380:         }
1381:         if (n25 >= 0) { /* directly above */
1382:           x_t = x;
1383:           y_t = ly[(n25 % (m * n)) / m];
1384:           /* z_t = lz[n25 / (m*n)]; */
1385:           s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1386:           for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1387:         } else if (ye - Ye < 0 && ze - Ze < 0) {
1388:           for (j = 0; j < x; j++) idx[nn++] = -1;
1389:         }
1390:         if (n26 >= 0) { /* right above */
1391:           x_t = lx[n26 % m];
1392:           y_t = ly[(n26 % (m * n)) / m];
1393:           /* z_t = lz[n26 / (m*n)]; */
1394:           s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1395:           for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1396:         } else if (xe - Xe < 0 && ye - Ye < 0 && ze - Ze < 0) {
1397:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
1398:         }
1399:       }
1400:     }
1401:   }
1402:   /*
1403:      Set the local to global ordering in the global vector, this allows use
1404:      of VecSetValuesLocal().
1405:   */
1406:   PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap));

1408:   PetscCall(PetscFree2(bases, ldims));
1409:   dd->m = m;
1410:   dd->n = n;
1411:   dd->p = p;
1412:   /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1413:   dd->xs = xs * dof;
1414:   dd->xe = xe * dof;
1415:   dd->ys = ys;
1416:   dd->ye = ye;
1417:   dd->zs = zs;
1418:   dd->ze = ze;
1419:   dd->Xs = Xs * dof;
1420:   dd->Xe = Xe * dof;
1421:   dd->Ys = Ys;
1422:   dd->Ye = Ye;
1423:   dd->Zs = Zs;
1424:   dd->Ze = Ze;

1426:   PetscCall(VecDestroy(&local));
1427:   PetscCall(VecDestroy(&global));

1429:   dd->gtol      = gtol;
1430:   dd->base      = base;
1431:   da->ops->view = DMView_DA_3d;
1432:   dd->ltol      = NULL;
1433:   dd->ao        = NULL;
1434:   PetscFunctionReturn(PETSC_SUCCESS);
1435: }

1437: /*@C
1438:    DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1439:    regular array data that is distributed across some processors.

1441:    Collective

1443:    Input Parameters:
1444: +  comm - MPI communicator
1445: .  bx,by,bz - type of ghost nodes the array have.
1446:          Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1447: .  stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`)
1448: .  M,N,P - global dimension in each direction of the array
1449: .  m,n,p - corresponding number of processors in each dimension
1450:            (or `PETSC_DECIDE` to have calculated)
1451: .  dof - number of degrees of freedom per node
1452: .  s - stencil width
1453: -  lx, ly, lz - arrays containing the number of nodes in each cell along
1454:           the x, y, and z coordinates, or NULL. If non-null, these
1455:           must be of length as m,n,p and the corresponding
1456:           m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
1457:           the ly[] must N, sum of the lz[] must be P

1459:    Output Parameter:
1460: .  da - the resulting distributed array object

1462:    Options Database Keys:
1463: +  -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()`
1464: .  -da_grid_x <nx> - number of grid points in x direction
1465: .  -da_grid_y <ny> - number of grid points in y direction
1466: .  -da_grid_z <nz> - number of grid points in z direction
1467: .  -da_processors_x <MX> - number of processors in x direction
1468: .  -da_processors_y <MY> - number of processors in y direction
1469: .  -da_processors_z <MZ> - number of processors in z direction
1470: .  -da_refine_x <rx> - refinement ratio in x direction
1471: .  -da_refine_y <ry> - refinement ratio in y direction
1472: .  -da_refine_z <rz>- refinement ratio in z directio
1473: -  -da_refine <n> - refine the DMDA n times before creating it

1475:    Level: beginner

1477:    Notes:
1478:    The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1479:    standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1480:    the standard 27-pt stencil.

1482:    The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
1483:    The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
1484:    and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed.

1486:    You must call `DMSetUp()` after this call before using this `DM`.

1488:    If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1489:    but before `DMSetUp()`.

1491: .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
1492:           `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
1493:           `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
1494:           `DMStagCreate3d()`
1495: @*/
1496: 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)
1497: {
1498:   PetscFunctionBegin;
1499:   PetscCall(DMDACreate(comm, da));
1500:   PetscCall(DMSetDimension(*da, 3));
1501:   PetscCall(DMDASetSizes(*da, M, N, P));
1502:   PetscCall(DMDASetNumProcs(*da, m, n, p));
1503:   PetscCall(DMDASetBoundaryType(*da, bx, by, bz));
1504:   PetscCall(DMDASetDof(*da, dof));
1505:   PetscCall(DMDASetStencilType(*da, stencil_type));
1506:   PetscCall(DMDASetStencilWidth(*da, s));
1507:   PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz));
1508:   PetscFunctionReturn(PETSC_SUCCESS);
1509: }