Actual source code: da2.c


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

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

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

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

 26:     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:       MPI_Comm_size(PetscObjectComm((PetscObject)da), &size);
 32:       DMDAGetLocalInfo(da, &info);
 33:       nzlocal = info.xm * info.ym;
 34:       PetscMalloc1(size, &nz);
 35:       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:       PetscFree(nz);
 42:       navg = navg / size;
 43:       PetscViewerASCIIPrintf(viewer, "  Load Balance - Grid Points: Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax);
 44:       return 0;
 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:       DMDAGetLocalInfo(da, &info);
 49:       PetscViewerASCIIPushSynchronized(viewer);
 50:       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:       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:       PetscViewerFlush(viewer);
 53:       PetscViewerASCIIPopSynchronized(viewer);
 54:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) DMView_DA_GLVis(da, viewer);
 55:     else 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:     PetscViewerDrawGetDraw(viewer, 0, &draw);
 67:     PetscDrawIsNull(draw, &isnull);
 68:     if (isnull) return 0;

 70:     PetscDrawCheckResizedWindow(draw);
 71:     PetscDrawClear(draw);
 72:     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++) 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++) PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK);
 83:     }
 84:     PetscDrawCollectiveEnd(draw);
 85:     PetscDrawFlush(draw);
 86:     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:     PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED);
 95:     PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED);
 96:     PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED);
 97:     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:         PetscSNPrintf(node, sizeof(node), "%d", (int)base++);
103:         PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node);
104:       }
105:     }
106:     PetscDrawCollectiveEnd(draw);
107:     PetscDrawFlush(draw);
108:     PetscDrawPause(draw);

110:     PetscDrawCollectiveBegin(draw);
111:     /* overlay ghost numbers, useful for error checking */
112:     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:           PetscSNPrintf(node, sizeof(node), "%d", (int)(idx[base / dd->w]));
122:           PetscDrawString(draw, x / dd->w, y, PETSC_DRAW_BLUE, node);
123:         }
124:         base++;
125:       }
126:     }
127:     ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx);
128:     PetscDrawCollectiveEnd(draw);
129:     PetscDrawFlush(draw);
130:     PetscDrawPause(draw);
131:     PetscDrawSave(draw);
132:   } else if (isglvis) {
133:     DMView_DA_GLVis(da, viewer);
134:   } else if (isbinary) {
135:     DMView_DA_Binary(da, viewer);
136: #if defined(PETSC_HAVE_MATLAB)
137:   } else if (ismatlab) {
138:     DMView_DA_Matlab(da, viewer);
139: #endif
140:   }
141:   return 0;
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:   (*ctx->func)(0, U, a, ctx->funcctx);
163:   (*ctx->funcisetbase)(U, ctx->funcctx);

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

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

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

181:   VecRestoreArray(U, &ww);
182:   VecRestoreArray(a, &aa);
183:   return 0;
184: }
185: #endif

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

213:   PetscObjectGetComm((PetscObject)da, &comm);
214: #if !defined(PETSC_USE_64BIT_INDICES)
216: #endif

218:   MPI_Comm_size(comm, &size);
219:   MPI_Comm_rank(comm, &rank);

221:   dd->p = 1;
222:   if (m != PETSC_DECIDE) {
225:   }
226:   if (n != PETSC_DECIDE) {
229:   }

231:   if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
232:     if (n != PETSC_DECIDE) {
233:       m = size / n;
234:     } else if (m != PETSC_DECIDE) {
235:       n = size / m;
236:     } else {
237:       /* try for squarish distribution */
238:       m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
239:       if (!m) m = 1;
240:       while (m > 0) {
241:         n = size / m;
242:         if (m * n == size) break;
243:         m--;
244:       }
245:       if (M > N && m < n) {
246:         PetscInt _m = m;
247:         m           = n;
248:         n           = _m;
249:       }
250:     }


257:   /*
258:      Determine locally owned region
259:      xs is the first local node number, x is the number of local nodes
260:   */
261:   if (!lx) {
262:     PetscMalloc1(m, &dd->lx);
263:     lx = dd->lx;
264:     for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i);
265:   }
266:   x  = lx[rank % m];
267:   xs = 0;
268:   for (i = 0; i < (rank % m); i++) xs += lx[i];
269:   if (PetscDefined(USE_DEBUG)) {
270:     left = xs;
271:     for (i = (rank % m); i < m; i++) left += lx[i];
273:   }

275:   /*
276:      Determine locally owned region
277:      ys is the first local node number, y is the number of local nodes
278:   */
279:   if (!ly) {
280:     PetscMalloc1(n, &dd->ly);
281:     ly = dd->ly;
282:     for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i);
283:   }
284:   y  = ly[rank / m];
285:   ys = 0;
286:   for (i = 0; i < (rank / m); i++) ys += ly[i];
287:   if (PetscDefined(USE_DEBUG)) {
288:     left = ys;
289:     for (i = (rank / m); i < n; i++) left += ly[i];
291:   }

293:   /*
294:    check if the scatter requires more than one process neighbor or wraps around
295:    the domain more than once
296:   */
299:   xe = xs + x;
300:   ye = ys + y;

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

327:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
328:     IXs = xs - s;
329:     IXe = xe + s;
330:     Xs  = xs - s;
331:     Xe  = xe + s;
332:   }

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

357:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
358:     IYs = ys - s;
359:     IYe = ye + s;
360:     Ys  = ys - s;
361:     Ye  = ye + s;
362:   }

364:   /* stencil length in each direction */
365:   s_x = s;
366:   s_y = s;

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

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

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

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

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

430:   /* determine who lies on each side of us stored in    n6 n7 n8
431:                                                         n3    n5
432:                                                         n0 n1 n2
433:   */

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

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

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

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

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

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

514:   if (stencil_type == DMDA_STENCIL_STAR) {
515:     /* save corner processor numbers */
516:     sn0 = n0;
517:     sn2 = n2;
518:     sn6 = n6;
519:     sn8 = n8;
520:     n0 = n2 = n6 = n8 = -1;
521:   }

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

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

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

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

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

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

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

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

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

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

599:   ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from);
600:   VecScatterCreate(global, from, local, to, &gtol);
601:   ISDestroy(&to);
602:   ISDestroy(&from);

604:   if (stencil_type == DMDA_STENCIL_STAR) {
605:     n0 = sn0;
606:     n2 = sn2;
607:     n6 = sn6;
608:     n8 = sn8;
609:   }

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

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

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

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

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

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

734:   VecDestroy(&local);
735:   VecDestroy(&global);

737:   dd->gtol      = gtol;
738:   dd->base      = base;
739:   da->ops->view = DMView_DA_2d;
740:   dd->ltol      = NULL;
741:   dd->ao        = NULL;
742:   return 0;
743: }

745: /*@C
746:    DMDACreate2d -  Creates an object that will manage the communication of  two-dimensional
747:    regular array data that is distributed across some processors.

749:    Collective

751:    Input Parameters:
752: +  comm - MPI communicator
753: .  bx,by - type of ghost nodes the array have.
754:          Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
755: .  stencil_type - stencil type.  Use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
756: .  M,N - global dimension in each direction of the array
757: .  m,n - corresponding number of processors in each dimension
758:          (or `PETSC_DECIDE` to have calculated)
759: .  dof - number of degrees of freedom per node
760: .  s - stencil width
761: -  lx, ly - arrays containing the number of nodes in each cell along
762:            the x and y coordinates, or NULL. If non-null, these
763:            must be of length as m and n, and the corresponding
764:            m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
765:            must be M, and the sum of the ly[] entries must be N.

767:    Output Parameter:
768: .  da - the resulting distributed array object

770:    Options Database Keys:
771: +  -dm_view - Calls `DMView()` at the conclusion of `DMDACreate2d()`
772: .  -da_grid_x <nx> - number of grid points in x direction
773: .  -da_grid_y <ny> - number of grid points in y direction
774: .  -da_processors_x <nx> - number of processors in x direction
775: .  -da_processors_y <ny> - number of processors in y direction
776: .  -da_refine_x <rx> - refinement ratio in x direction
777: .  -da_refine_y <ry> - refinement ratio in y direction
778: -  -da_refine <n> - refine the DMDA n times before creating

780:    Level: beginner

782:    Notes:
783:    The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
784:    standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
785:    the standard 9-pt stencil.

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

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

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

796: .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
797:           `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
798:           `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
799:           `DMStagCreate2d()`
800: @*/

802: 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)
803: {
804:   DMDACreate(comm, da);
805:   DMSetDimension(*da, 2);
806:   DMDASetSizes(*da, M, N, 1);
807:   DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
808:   DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);
809:   DMDASetDof(*da, dof);
810:   DMDASetStencilType(*da, stencil_type);
811:   DMDASetStencilWidth(*da, s);
812:   DMDASetOwnershipRanges(*da, lx, ly, NULL);
813:   return 0;
814: }