Actual source code: da3.c

  1: /*
  2:    Code for manipulating distributed regular 3d arrays in parallel.
  3:    File created by Peter Mell  7/14/95
  4:  */

  6: #include <petsc/private/dmdaimpl.h>

  8: #include <petscdraw.h>
  9: static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer)
 10: {
 11:   PetscMPIInt rank;
 12:   PetscBool   iascii, isdraw, isglvis, isbinary;
 13:   DM_DA      *dd = (DM_DA *)da->data;
 14:   Vec         coordinates;
 15: #if defined(PETSC_HAVE_MATLAB)
 16:   PetscBool ismatlab;
 17: #endif

 19:   PetscFunctionBegin;
 20:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));

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

 32:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
 33:     PetscCall(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:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
 39:       PetscCall(DMDAGetLocalInfo(da, &info));
 40:       nzlocal = info.xm * info.ym * info.zm;
 41:       PetscCall(PetscMalloc1(size, &nz));
 42:       PetscCallMPI(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:       PetscCall(PetscFree(nz));
 49:       navg = navg / size;
 50:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Load Balance - Grid Points: Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax));
 51:       PetscFunctionReturn(PETSC_SUCCESS);
 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:       PetscCall(DMDAGetLocalInfo(da, &info));
 56:       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));
 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:       PetscCall(DMGetCoordinates(da, &coordinates));
 60: #if !defined(PETSC_USE_COMPLEX)
 61:       if (coordinates) {
 62:         PetscInt         last;
 63:         const PetscReal *coors;
 64:         PetscCall(VecGetArrayRead(coordinates, &coors));
 65:         PetscCall(VecGetLocalSize(coordinates, &last));
 66:         last = last - 3;
 67:         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]));
 68:         PetscCall(VecRestoreArrayRead(coordinates, &coors));
 69:       }
 70: #endif
 71:       PetscCall(PetscViewerFlush(viewer));
 72:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
 73:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer));
 74:     else PetscCall(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:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
 85:     PetscCall(PetscDrawIsNull(draw, &isnull));
 86:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

 88:     PetscCall(PetscDrawCheckResizedWindow(draw));
 89:     PetscCall(PetscDrawClear(draw));
 90:     PetscCall(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++) PetscCall(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++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
102:       }
103:     }
104:     PetscDrawCollectiveEnd(draw);
105:     PetscCall(PetscDrawFlush(draw));
106:     PetscCall(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:         PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
119:         PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
120:         PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
121:         PetscCall(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:         PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)rank));
128:         PetscCall(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:             PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++));
134:             PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
135:           }
136:         }
137:       }
138:     }
139:     PetscDrawCollectiveEnd(draw);
140:     PetscCall(PetscDrawFlush(draw));
141:     PetscCall(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:         PetscCall(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:             PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, 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:             PetscCall(PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node));
169:             base++;
170:           }
171:         }
172:         PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
173:       }
174:     }
175:     PetscDrawCollectiveEnd(draw);
176:     PetscCall(PetscDrawFlush(draw));
177:     PetscCall(PetscDrawPause(draw));
178:     PetscCall(PetscDrawSave(draw));
179:   } else if (isglvis) {
180:     PetscCall(DMView_DA_GLVis(da, viewer));
181:   } else if (isbinary) {
182:     PetscCall(DMView_DA_Binary(da, viewer));
183: #if defined(PETSC_HAVE_MATLAB)
184:   } else if (ismatlab) {
185:     PetscCall(DMView_DA_Matlab(da, viewer));
186: #endif
187:   }
188:   PetscFunctionReturn(PETSC_SUCCESS);
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;

225:   PetscFunctionBegin;
226:   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");
227:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
228: #if !defined(PETSC_USE_64BIT_INDICES)
229:   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);
230: #endif

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

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

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

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

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

335:   if (!lx) {
336:     PetscCall(PetscMalloc1(m, &dd->lx));
337:     lx = dd->lx;
338:     for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m));
339:   }
340:   x  = lx[rank % m];
341:   xs = 0;
342:   for (i = 0; i < (rank % m); i++) xs += lx[i];
343:   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);

345:   if (!ly) {
346:     PetscCall(PetscMalloc1(n, &dd->ly));
347:     ly = dd->ly;
348:     for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n));
349:   }
350:   y = ly[(rank % (m * n)) / m];
351:   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);

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

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

363:   /* note this is different than x- and y-, as we will handle as an important special
364:    case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z.  This is to deal with 2D problems
365:    in a 3D code.  Additional code for this case is noted with "2d case" comments */
366:   twod = PETSC_FALSE;
367:   if (P == 1) twod = PETSC_TRUE;
368:   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);
369:   zs = 0;
370:   for (i = 0; i < (rank / (m * n)); i++) zs += lz[i];
371:   ye = ys + y;
372:   xe = xs + x;
373:   ze = zs + z;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1436: /*@C
1437:   DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1438:   regular array data that is distributed across one or more MPI processes.

1440:   Collective

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

1463:   Output Parameter:
1464: . da - the resulting distributed array object

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

1479:   Level: beginner

1481:   Notes:
1482:   If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the
1483:   corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`,
1484:   sum of the `ly` must `N`, sum of the `lz` must be `P`.

1486:   The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1487:   standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1488:   the standard 27-pt stencil.

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

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

1496:   To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1497:   but before `DMSetUp()`.

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