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:   PetscCheck((x > s) || (bx != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", x, s);
364:   PetscCheck((y > s) || (by != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", y, s);
365:   PetscCheck((z > s) || (bz != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local z-width of domain z %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", z, s);

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

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

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

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

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

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

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

454:   /* Resize all X parameters to reflect w */
455:   s_x = s;
456:   s_y = s;
457:   s_z = s;

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

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

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

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

532:   /* determine who lies on each side of use stored in    n24 n25 n26
533:                                                          n21 n22 n23
534:                                                          n18 n19 n20

536:                                                          n15 n16 n17
537:                                                          n12     n14
538:                                                          n9  n10 n11

540:                                                          n6  n7  n8
541:                                                          n3  n4  n5
542:                                                          n0  n1  n2
543:   */

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

557:   n9  = rank - m - 1;
558:   n10 = rank - m;
559:   n11 = rank - m + 1;
560:   n12 = rank - 1;
561:   n14 = rank + 1;
562:   n15 = rank + m - 1;
563:   n16 = rank + m;
564:   n17 = rank + m + 1;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1068:   PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
1069:   PetscCall(VecScatterCreate(global, from, local, to, &gtol));
1070:   PetscCall(ISDestroy(&to));
1071:   PetscCall(ISDestroy(&from));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1429:   PetscCall(VecDestroy(&local));
1430:   PetscCall(VecDestroy(&global));

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

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

1444:   Collective

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

1467:   Output Parameter:
1468: . da - the resulting distributed array object

1470:   Options Database Keys:
1471: + -dm_view              - Calls `DMView()` at the conclusion of `DMDACreate3d()`
1472: . -da_grid_x <nx>       - number of grid points in x direction
1473: . -da_grid_y <ny>       - number of grid points in y direction
1474: . -da_grid_z <nz>       - number of grid points in z direction
1475: . -da_processors_x <MX> - number of processors in x direction
1476: . -da_processors_y <MY> - number of processors in y direction
1477: . -da_processors_z <MZ> - number of processors in z direction
1478: . -da_bd_x <bx>         - boundary type in x direction
1479: . -da_bd_y <by>         - boundary type in y direction
1480: . -da_bd_z <bz>         - boundary type in x direction
1481: . -da_bd_all <bt>       - boundary type in all directions
1482: . -da_refine_x <rx>     - refinement ratio in x direction
1483: . -da_refine_y <ry>     - refinement ratio in y direction
1484: . -da_refine_z <rz>     - refinement ratio in z directio
1485: - -da_refine <n>        - refine the `DMDA` n times before creating it

1487:   Level: beginner

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

1494:   The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1495:   standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1496:   the standard 27-pt stencil.

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

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

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

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