Actual source code: da2.c

  1: #include <petsc/private/dmdaimpl.h>
  2: #include <petscdraw.h>

  4: static PetscErrorCode DMView_DA_2d(DM da, PetscViewer viewer)
  5: {
  6:   PetscMPIInt rank;
  7:   PetscBool   iascii, isdraw, isglvis, isbinary;
  8:   DM_DA      *dd = (DM_DA *)da->data;
  9: #if defined(PETSC_HAVE_MATLAB)
 10:   PetscBool ismatlab;
 11: #endif

 13:   PetscFunctionBegin;
 14:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));

 16:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 17:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
 18:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
 19:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
 20: #if defined(PETSC_HAVE_MATLAB)
 21:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
 22: #endif
 23:   if (iascii) {
 24:     PetscViewerFormat format;

 26:     PetscCall(PetscViewerGetFormat(viewer, &format));
 27:     if (format == PETSC_VIEWER_LOAD_BALANCE) {
 28:       PetscInt      i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal;
 29:       DMDALocalInfo info;
 30:       PetscMPIInt   size;
 31:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
 32:       PetscCall(DMDAGetLocalInfo(da, &info));
 33:       nzlocal = info.xm * info.ym;
 34:       PetscCall(PetscMalloc1(size, &nz));
 35:       PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da)));
 36:       for (i = 0; i < (PetscInt)size; i++) {
 37:         nmax = PetscMax(nmax, nz[i]);
 38:         nmin = PetscMin(nmin, nz[i]);
 39:         navg += nz[i];
 40:       }
 41:       PetscCall(PetscFree(nz));
 42:       navg = navg / size;
 43:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Load Balance - Grid Points: Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax));
 44:       PetscFunctionReturn(PETSC_SUCCESS);
 45:     }
 46:     if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) {
 47:       DMDALocalInfo info;
 48:       PetscCall(DMDAGetLocalInfo(da, &info));
 49:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
 50:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->m, dd->n, dd->w, dd->s));
 51:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs, info.xs + info.xm, info.ys, info.ys + info.ym));
 52:       PetscCall(PetscViewerFlush(viewer));
 53:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
 54:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer));
 55:     else PetscCall(DMView_DA_VTK(da, viewer));
 56:   } else if (isdraw) {
 57:     PetscDraw       draw;
 58:     double          ymin = -1 * dd->s - 1, ymax = dd->N + dd->s;
 59:     double          xmin = -1 * dd->s - 1, xmax = dd->M + dd->s;
 60:     double          x, y;
 61:     PetscInt        base;
 62:     const PetscInt *idx;
 63:     char            node[10];
 64:     PetscBool       isnull;

 66:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
 67:     PetscCall(PetscDrawIsNull(draw, &isnull));
 68:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

 70:     PetscCall(PetscDrawCheckResizedWindow(draw));
 71:     PetscCall(PetscDrawClear(draw));
 72:     PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));

 74:     PetscDrawCollectiveBegin(draw);
 75:     /* first processor draw all node lines */
 76:     if (rank == 0) {
 77:       ymin = 0.0;
 78:       ymax = dd->N - 1;
 79:       for (xmin = 0; xmin < dd->M; xmin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK));
 80:       xmin = 0.0;
 81:       xmax = dd->M - 1;
 82:       for (ymin = 0; ymin < dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
 83:     }
 84:     PetscDrawCollectiveEnd(draw);
 85:     PetscCall(PetscDrawFlush(draw));
 86:     PetscCall(PetscDrawPause(draw));

 88:     PetscDrawCollectiveBegin(draw);
 89:     /* draw my box */
 90:     xmin = dd->xs / dd->w;
 91:     xmax = (dd->xe - 1) / dd->w;
 92:     ymin = dd->ys;
 93:     ymax = dd->ye - 1;
 94:     PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
 95:     PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
 96:     PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
 97:     PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED));
 98:     /* put in numbers */
 99:     base = (dd->base) / dd->w;
100:     for (y = ymin; y <= ymax; y++) {
101:       for (x = xmin; x <= xmax; x++) {
102:         PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++));
103:         PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
104:       }
105:     }
106:     PetscDrawCollectiveEnd(draw);
107:     PetscCall(PetscDrawFlush(draw));
108:     PetscCall(PetscDrawPause(draw));

110:     PetscDrawCollectiveBegin(draw);
111:     /* overlay ghost numbers, useful for error checking */
112:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx));
113:     base = 0;
114:     xmin = dd->Xs;
115:     xmax = dd->Xe;
116:     ymin = dd->Ys;
117:     ymax = dd->Ye;
118:     for (y = ymin; y < ymax; y++) {
119:       for (x = xmin; x < xmax; x++) {
120:         if ((base % dd->w) == 0) {
121:           PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)(idx[base / dd->w])));
122:           PetscCall(PetscDrawString(draw, x / dd->w, y, PETSC_DRAW_BLUE, node));
123:         }
124:         base++;
125:       }
126:     }
127:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
128:     PetscDrawCollectiveEnd(draw);
129:     PetscCall(PetscDrawFlush(draw));
130:     PetscCall(PetscDrawPause(draw));
131:     PetscCall(PetscDrawSave(draw));
132:   } else if (isglvis) {
133:     PetscCall(DMView_DA_GLVis(da, viewer));
134:   } else if (isbinary) {
135:     PetscCall(DMView_DA_Binary(da, viewer));
136: #if defined(PETSC_HAVE_MATLAB)
137:   } else if (ismatlab) {
138:     PetscCall(DMView_DA_Matlab(da, viewer));
139: #endif
140:   }
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: #if defined(new)
145: /*
146:   DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix where local
147:     function lives on a DMDA

149:         y ~= (F(u + ha) - F(u))/h,
150:   where F = nonlinear function, as set by SNESSetFunction()
151:         u = current iterate
152:         h = difference interval
153: */
154: PetscErrorCode DMDAGetDiagonal_MFFD(DM da, Vec U, Vec a)
155: {
156:   PetscScalar   h, *aa, *ww, v;
157:   PetscReal     epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
158:   PetscInt      gI, nI;
159:   MatStencil    stencil;
160:   DMDALocalInfo info;

162:   PetscFunctionBegin;
163:   PetscCall((*ctx->func)(0, U, a, ctx->funcctx));
164:   PetscCall((*ctx->funcisetbase)(U, ctx->funcctx));

166:   PetscCall(VecGetArray(U, &ww));
167:   PetscCall(VecGetArray(a, &aa));

169:   nI = 0;
170:   h  = ww[gI];
171:   if (h == 0.0) h = 1.0;
172:   if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
173:   else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
174:   h *= epsilon;

176:   ww[gI] += h;
177:   PetscCall((*ctx->funci)(i, w, &v, ctx->funcctx));
178:   aa[nI] = (v - aa[nI]) / h;
179:   ww[gI] -= h;
180:   nI++;

182:   PetscCall(VecRestoreArray(U, &ww));
183:   PetscCall(VecRestoreArray(a, &aa));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }
186: #endif

188: PetscErrorCode DMSetUp_DA_2D(DM da)
189: {
190:   DM_DA          *dd           = (DM_DA *)da->data;
191:   const PetscInt  M            = dd->M;
192:   const PetscInt  N            = dd->N;
193:   PetscInt        m            = dd->m;
194:   PetscInt        n            = dd->n;
195:   const PetscInt  dof          = dd->w;
196:   const PetscInt  s            = dd->s;
197:   DMBoundaryType  bx           = dd->bx;
198:   DMBoundaryType  by           = dd->by;
199:   DMDAStencilType stencil_type = dd->stencil_type;
200:   PetscInt       *lx           = dd->lx;
201:   PetscInt       *ly           = dd->ly;
202:   MPI_Comm        comm;
203:   PetscMPIInt     rank, size;
204:   PetscInt        xs, xe, ys, ye, x, y, Xs, Xe, Ys, Ye, IXs, IXe, IYs, IYe;
205:   PetscInt        up, down, left, right, i, n0, n1, n2, n3, n5, n6, n7, n8, *idx, nn;
206:   PetscInt        xbase, *bases, *ldims, j, x_t, y_t, s_t, base, count;
207:   PetscInt        s_x, s_y; /* s proportionalized to w */
208:   PetscInt        sn0 = 0, sn2 = 0, sn6 = 0, sn8 = 0;
209:   Vec             local, global;
210:   VecScatter      gtol;
211:   IS              to, from;

213:   PetscFunctionBegin;
214:   PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && by != DM_BOUNDARY_MIRROR), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box stencil");
215:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
216: #if !defined(PETSC_USE_64BIT_INDICES)
217:   PetscCheck(((PetscInt64)M) * ((PetscInt64)N) * ((PetscInt64)dof) <= (PetscInt64)PETSC_MPI_INT_MAX, comm, PETSC_ERR_INT_OVERFLOW, "Mesh of %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " (dof) is too large for 32-bit indices", M, N, dof);
218: #endif

220:   PetscCallMPI(MPI_Comm_size(comm, &size));
221:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

223:   dd->p = 1;
224:   if (m != PETSC_DECIDE) {
225:     PetscCheck(m >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m);
226:     PetscCheck(m <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %d", m, size);
227:   }
228:   if (n != PETSC_DECIDE) {
229:     PetscCheck(n >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %" PetscInt_FMT, n);
230:     PetscCheck(n <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %d", n, size);
231:   }

233:   if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
234:     if (n != PETSC_DECIDE) {
235:       m = size / n;
236:     } else if (m != PETSC_DECIDE) {
237:       n = size / m;
238:     } else {
239:       /* try for squarish distribution */
240:       m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
241:       if (!m) m = 1;
242:       while (m > 0) {
243:         n = size / m;
244:         if (m * n == size) break;
245:         m--;
246:       }
247:       if (M > N && m < n) {
248:         PetscInt _m = m;
249:         m           = n;
250:         n           = _m;
251:       }
252:     }
253:     PetscCheck(m * n == size, comm, PETSC_ERR_PLIB, "Unable to create partition, check the size of the communicator and input m and n ");
254:   } else PetscCheck(m * n == size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition");

256:   PetscCheck(M >= m, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m);
257:   PetscCheck(N >= n, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n);

259:   /*
260:      Determine locally owned region
261:      xs is the first local node number, x is the number of local nodes
262:   */
263:   if (!lx) {
264:     PetscCall(PetscMalloc1(m, &dd->lx));
265:     lx = dd->lx;
266:     for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i);
267:   }
268:   x  = lx[rank % m];
269:   xs = 0;
270:   for (i = 0; i < (rank % m); i++) xs += lx[i];
271:   if (PetscDefined(USE_DEBUG)) {
272:     left = xs;
273:     for (i = (rank % m); i < m; i++) left += lx[i];
274:     PetscCheck(left == M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of lx across processors not equal to M: %" PetscInt_FMT " %" PetscInt_FMT, left, M);
275:   }

277:   /*
278:      Determine locally owned region
279:      ys is the first local node number, y is the number of local nodes
280:   */
281:   if (!ly) {
282:     PetscCall(PetscMalloc1(n, &dd->ly));
283:     ly = dd->ly;
284:     for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i);
285:   }
286:   y  = ly[rank / m];
287:   ys = 0;
288:   for (i = 0; i < (rank / m); i++) ys += ly[i];
289:   if (PetscDefined(USE_DEBUG)) {
290:     left = ys;
291:     for (i = (rank / m); i < n; i++) left += ly[i];
292:     PetscCheck(left == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of ly across processors not equal to N: %" PetscInt_FMT " %" PetscInt_FMT, left, N);
293:   }

295:   /*
296:    check if the scatter requires more than one process neighbor or wraps around
297:    the domain more than once
298:   */
299:   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);
300:   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);
301:   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);
302:   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);
303:   xe = xs + x;
304:   ye = ys + y;

306:   /* determine ghost region (Xs) and region scattered into (IXs)  */
307:   if (xs - s > 0) {
308:     Xs  = xs - s;
309:     IXs = xs - s;
310:   } else {
311:     if (bx) {
312:       Xs = xs - s;
313:     } else {
314:       Xs = 0;
315:     }
316:     IXs = 0;
317:   }
318:   if (xe + s <= M) {
319:     Xe  = xe + s;
320:     IXe = xe + s;
321:   } else {
322:     if (bx) {
323:       Xs = xs - s;
324:       Xe = xe + s;
325:     } else {
326:       Xe = M;
327:     }
328:     IXe = M;
329:   }

331:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
332:     IXs = xs - s;
333:     IXe = xe + s;
334:     Xs  = xs - s;
335:     Xe  = xe + s;
336:   }

338:   if (ys - s > 0) {
339:     Ys  = ys - s;
340:     IYs = ys - s;
341:   } else {
342:     if (by) {
343:       Ys = ys - s;
344:     } else {
345:       Ys = 0;
346:     }
347:     IYs = 0;
348:   }
349:   if (ye + s <= N) {
350:     Ye  = ye + s;
351:     IYe = ye + s;
352:   } else {
353:     if (by) {
354:       Ye = ye + s;
355:     } else {
356:       Ye = N;
357:     }
358:     IYe = N;
359:   }

361:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
362:     IYs = ys - s;
363:     IYe = ye + s;
364:     Ys  = ys - s;
365:     Ye  = ye + s;
366:   }

368:   /* stencil length in each direction */
369:   s_x = s;
370:   s_y = s;

372:   /* determine starting point of each processor */
373:   nn = x * y;
374:   PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims));
375:   PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm));
376:   bases[0] = 0;
377:   for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
378:   for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
379:   base = bases[rank] * dof;

381:   /* allocate the base parallel and sequential vectors */
382:   dd->Nlocal = x * y * dof;
383:   PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global));
384:   dd->nlocal = (Xe - Xs) * (Ye - Ys) * dof;
385:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local));

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

389:   /* global to local must include ghost points within the domain,
390:      but not ghost points outside the domain that aren't periodic */
391:   PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs), &idx));
392:   if (stencil_type == DMDA_STENCIL_BOX) {
393:     left  = IXs - Xs;
394:     right = left + (IXe - IXs);
395:     down  = IYs - Ys;
396:     up    = down + (IYe - IYs);
397:     count = 0;
398:     for (i = down; i < up; i++) {
399:       for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
400:     }
401:     PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));

403:   } else {
404:     /* must drop into cross shape region */
405:     /*       ---------|
406:             |  top    |
407:          |---         ---| up
408:          |   middle      |
409:          |               |
410:          ----         ---- down
411:             | bottom  |
412:             -----------
413:          Xs xs        xe Xe */
414:     left  = xs - Xs;
415:     right = left + x;
416:     down  = ys - Ys;
417:     up    = down + y;
418:     count = 0;
419:     /* bottom */
420:     for (i = (IYs - Ys); i < down; i++) {
421:       for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
422:     }
423:     /* middle */
424:     for (i = down; i < up; i++) {
425:       for (j = (IXs - Xs); j < (IXe - Xs); j++) idx[count++] = j + i * (Xe - Xs);
426:     }
427:     /* top */
428:     for (i = up; i < up + IYe - ye; i++) {
429:       for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
430:     }
431:     PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
432:   }

434:   /* determine who lies on each side of us stored in    n6 n7 n8
435:                                                         n3    n5
436:                                                         n0 n1 n2
437:   */

439:   /* Assume the Non-Periodic Case */
440:   n1 = rank - m;
441:   if (rank % m) {
442:     n0 = n1 - 1;
443:   } else {
444:     n0 = -1;
445:   }
446:   if ((rank + 1) % m) {
447:     n2 = n1 + 1;
448:     n5 = rank + 1;
449:     n8 = rank + m + 1;
450:     if (n8 >= m * n) n8 = -1;
451:   } else {
452:     n2 = -1;
453:     n5 = -1;
454:     n8 = -1;
455:   }
456:   if (rank % m) {
457:     n3 = rank - 1;
458:     n6 = n3 + m;
459:     if (n6 >= m * n) n6 = -1;
460:   } else {
461:     n3 = -1;
462:     n6 = -1;
463:   }
464:   n7 = rank + m;
465:   if (n7 >= m * n) n7 = -1;

467:   if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) {
468:     /* Modify for Periodic Cases */
469:     /* Handle all four corners */
470:     if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m - 1;
471:     if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
472:     if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size - m;
473:     if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size - 1;

475:     /* Handle Top and Bottom Sides */
476:     if (n1 < 0) n1 = rank + m * (n - 1);
477:     if (n7 < 0) n7 = rank - m * (n - 1);
478:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
479:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
480:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
481:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;

483:     /* Handle Left and Right Sides */
484:     if (n3 < 0) n3 = rank + (m - 1);
485:     if (n5 < 0) n5 = rank - (m - 1);
486:     if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
487:     if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
488:     if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
489:     if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
490:   } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
491:     if (n1 < 0) n1 = rank + m * (n - 1);
492:     if (n7 < 0) n7 = rank - m * (n - 1);
493:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
494:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
495:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
496:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;
497:   } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
498:     if (n3 < 0) n3 = rank + (m - 1);
499:     if (n5 < 0) n5 = rank - (m - 1);
500:     if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
501:     if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
502:     if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
503:     if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
504:   }

506:   PetscCall(PetscMalloc1(9, &dd->neighbors));

508:   dd->neighbors[0] = n0;
509:   dd->neighbors[1] = n1;
510:   dd->neighbors[2] = n2;
511:   dd->neighbors[3] = n3;
512:   dd->neighbors[4] = rank;
513:   dd->neighbors[5] = n5;
514:   dd->neighbors[6] = n6;
515:   dd->neighbors[7] = n7;
516:   dd->neighbors[8] = n8;

518:   if (stencil_type == DMDA_STENCIL_STAR) {
519:     /* save corner processor numbers */
520:     sn0 = n0;
521:     sn2 = n2;
522:     sn6 = n6;
523:     sn8 = n8;
524:     n0 = n2 = n6 = n8 = -1;
525:   }

527:   PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx));

529:   nn    = 0;
530:   xbase = bases[rank];
531:   for (i = 1; i <= s_y; i++) {
532:     if (n0 >= 0) { /* left below */
533:       x_t = lx[n0 % m];
534:       y_t = ly[n0 / m];
535:       s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
536:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
537:     }

539:     if (n1 >= 0) { /* directly below */
540:       x_t = x;
541:       y_t = ly[n1 / m];
542:       s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
543:       for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
544:     } else if (by == DM_BOUNDARY_MIRROR) {
545:       for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
546:     }

548:     if (n2 >= 0) { /* right below */
549:       x_t = lx[n2 % m];
550:       y_t = ly[n2 / m];
551:       s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
552:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
553:     }
554:   }

556:   for (i = 0; i < y; i++) {
557:     if (n3 >= 0) { /* directly left */
558:       x_t = lx[n3 % m];
559:       /* y_t = y; */
560:       s_t = bases[n3] + (i + 1) * x_t - s_x;
561:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
562:     } else if (bx == DM_BOUNDARY_MIRROR) {
563:       for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
564:     }

566:     for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */

568:     if (n5 >= 0) { /* directly right */
569:       x_t = lx[n5 % m];
570:       /* y_t = y; */
571:       s_t = bases[n5] + (i)*x_t;
572:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
573:     } else if (bx == DM_BOUNDARY_MIRROR) {
574:       for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
575:     }
576:   }

578:   for (i = 1; i <= s_y; i++) {
579:     if (n6 >= 0) { /* left above */
580:       x_t = lx[n6 % m];
581:       /* y_t = ly[n6 / m]; */
582:       s_t = bases[n6] + (i)*x_t - s_x;
583:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
584:     }

586:     if (n7 >= 0) { /* directly above */
587:       x_t = x;
588:       /* y_t = ly[n7 / m]; */
589:       s_t = bases[n7] + (i - 1) * x_t;
590:       for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
591:     } else if (by == DM_BOUNDARY_MIRROR) {
592:       for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
593:     }

595:     if (n8 >= 0) { /* right above */
596:       x_t = lx[n8 % m];
597:       /* y_t = ly[n8 / m]; */
598:       s_t = bases[n8] + (i - 1) * x_t;
599:       for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
600:     }
601:   }

603:   PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
604:   PetscCall(VecScatterCreate(global, from, local, to, &gtol));
605:   PetscCall(ISDestroy(&to));
606:   PetscCall(ISDestroy(&from));

608:   if (stencil_type == DMDA_STENCIL_STAR) {
609:     n0 = sn0;
610:     n2 = sn2;
611:     n6 = sn6;
612:     n8 = sn8;
613:   }

615:   if ((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC)) {
616:     /*
617:         Recompute the local to global mappings, this time keeping the
618:       information about the cross corner processor numbers and any ghosted
619:       but not periodic indices.
620:     */
621:     nn    = 0;
622:     xbase = bases[rank];
623:     for (i = 1; i <= s_y; i++) {
624:       if (n0 >= 0) { /* left below */
625:         x_t = lx[n0 % m];
626:         y_t = ly[n0 / m];
627:         s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
628:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
629:       } else if (xs - Xs > 0 && ys - Ys > 0) {
630:         for (j = 0; j < s_x; j++) idx[nn++] = -1;
631:       }
632:       if (n1 >= 0) { /* directly below */
633:         x_t = x;
634:         y_t = ly[n1 / m];
635:         s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
636:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
637:       } else if (ys - Ys > 0) {
638:         if (by == DM_BOUNDARY_MIRROR) {
639:           for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
640:         } else {
641:           for (j = 0; j < x; j++) idx[nn++] = -1;
642:         }
643:       }
644:       if (n2 >= 0) { /* right below */
645:         x_t = lx[n2 % m];
646:         y_t = ly[n2 / m];
647:         s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
648:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
649:       } else if (Xe - xe > 0 && ys - Ys > 0) {
650:         for (j = 0; j < s_x; j++) idx[nn++] = -1;
651:       }
652:     }

654:     for (i = 0; i < y; i++) {
655:       if (n3 >= 0) { /* directly left */
656:         x_t = lx[n3 % m];
657:         /* y_t = y; */
658:         s_t = bases[n3] + (i + 1) * x_t - s_x;
659:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
660:       } else if (xs - Xs > 0) {
661:         if (bx == DM_BOUNDARY_MIRROR) {
662:           for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
663:         } else {
664:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
665:         }
666:       }

668:       for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */

670:       if (n5 >= 0) { /* directly right */
671:         x_t = lx[n5 % m];
672:         /* y_t = y; */
673:         s_t = bases[n5] + (i)*x_t;
674:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
675:       } else if (Xe - xe > 0) {
676:         if (bx == DM_BOUNDARY_MIRROR) {
677:           for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
678:         } else {
679:           for (j = 0; j < s_x; j++) idx[nn++] = -1;
680:         }
681:       }
682:     }

684:     for (i = 1; i <= s_y; i++) {
685:       if (n6 >= 0) { /* left above */
686:         x_t = lx[n6 % m];
687:         /* y_t = ly[n6 / m]; */
688:         s_t = bases[n6] + (i)*x_t - s_x;
689:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
690:       } else if (xs - Xs > 0 && Ye - ye > 0) {
691:         for (j = 0; j < s_x; j++) idx[nn++] = -1;
692:       }
693:       if (n7 >= 0) { /* directly above */
694:         x_t = x;
695:         /* y_t = ly[n7 / m]; */
696:         s_t = bases[n7] + (i - 1) * x_t;
697:         for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
698:       } else if (Ye - ye > 0) {
699:         if (by == DM_BOUNDARY_MIRROR) {
700:           for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
701:         } else {
702:           for (j = 0; j < x; j++) idx[nn++] = -1;
703:         }
704:       }
705:       if (n8 >= 0) { /* right above */
706:         x_t = lx[n8 % m];
707:         /* y_t = ly[n8 / m]; */
708:         s_t = bases[n8] + (i - 1) * x_t;
709:         for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
710:       } else if (Xe - xe > 0 && Ye - ye > 0) {
711:         for (j = 0; j < s_x; j++) idx[nn++] = -1;
712:       }
713:     }
714:   }
715:   /*
716:      Set the local to global ordering in the global vector, this allows use
717:      of VecSetValuesLocal().
718:   */
719:   PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap));

721:   PetscCall(PetscFree2(bases, ldims));
722:   dd->m = m;
723:   dd->n = n;
724:   /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
725:   dd->xs = xs * dof;
726:   dd->xe = xe * dof;
727:   dd->ys = ys;
728:   dd->ye = ye;
729:   dd->zs = 0;
730:   dd->ze = 1;
731:   dd->Xs = Xs * dof;
732:   dd->Xe = Xe * dof;
733:   dd->Ys = Ys;
734:   dd->Ye = Ye;
735:   dd->Zs = 0;
736:   dd->Ze = 1;

738:   PetscCall(VecDestroy(&local));
739:   PetscCall(VecDestroy(&global));

741:   dd->gtol      = gtol;
742:   dd->base      = base;
743:   da->ops->view = DMView_DA_2d;
744:   dd->ltol      = NULL;
745:   dd->ao        = NULL;
746:   PetscFunctionReturn(PETSC_SUCCESS);
747: }

749: /*@C
750:   DMDACreate2d -  Creates an object that will manage the communication of two-dimensional
751:   regular array data that is distributed across one or more MPI processes.

753:   Collective

755:   Input Parameters:
756: + comm         - MPI communicator
757: . bx           - type of ghost nodes the x array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
758: . by           - type of ghost nodes the y array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
759: . stencil_type - stencil type.  Use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
760: . M            - global dimension in x direction of the array
761: . N            - global dimension in y direction of the array
762: . m            - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated)
763: . n            - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated)
764: . dof          - number of degrees of freedom per node
765: . s            - stencil width
766: . lx           - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`.
767: - ly           - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`.

769:   Output Parameter:
770: . da - the resulting distributed array object

772:   Options Database Keys:
773: + -dm_view              - Calls `DMView()` at the conclusion of `DMDACreate2d()`
774: . -da_grid_x <nx>       - number of grid points in x direction
775: . -da_grid_y <ny>       - number of grid points in y direction
776: . -da_processors_x <nx> - number of processors in x direction
777: . -da_processors_y <ny> - number of processors in y direction
778: . -da_bd_x <bx>         - boundary type in x direction
779: . -da_bd_y <by>         - boundary type in y direction
780: . -da_bd_all <bt>       - boundary type in all directions
781: . -da_refine_x <rx>     - refinement ratio in x direction
782: . -da_refine_y <ry>     - refinement ratio in y direction
783: - -da_refine <n>        - refine the `DMDA` n times before creating

785:   Level: beginner

787:   Notes:
788:   If `lx` or `ly` are non-null, these must be of length as `m` and `n`, and the corresponding
789:   `m` and `n` cannot be `PETSC_DECIDE`. The sum of the `lx` entries must be `M`, and the sum of
790:   the `ly` entries must be `N`.

792:   The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
793:   standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
794:   the standard 9-pt stencil.

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

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

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

805: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
806:           `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
807:           `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
808:           `DMStagCreate2d()`, `DMBoundaryType`
809: @*/
810: PetscErrorCode DMDACreate2d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt m, PetscInt n, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], DM *da)
811: {
812:   PetscFunctionBegin;
813:   PetscCall(DMDACreate(comm, da));
814:   PetscCall(DMSetDimension(*da, 2));
815:   PetscCall(DMDASetSizes(*da, M, N, 1));
816:   PetscCall(DMDASetNumProcs(*da, m, n, PETSC_DECIDE));
817:   PetscCall(DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE));
818:   PetscCall(DMDASetDof(*da, dof));
819:   PetscCall(DMDASetStencilType(*da, stencil_type));
820:   PetscCall(DMDASetStencilWidth(*da, s));
821:   PetscCall(DMDASetOwnershipRanges(*da, lx, ly, NULL));
822:   PetscFunctionReturn(PETSC_SUCCESS);
823: }