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_ENGINE)
 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_ENGINE)
 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 %D  avg %D  max %D\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 %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);
 51:       PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\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) {
 55:       DMView_DA_GLVis(da,viewer);
 56:     } else {
 57:       DMView_DA_VTK(da,viewer);
 58:     }
 59:   } else if (isdraw) {
 60:     PetscDraw      draw;
 61:     double         ymin = -1*dd->s-1,ymax = dd->N+dd->s;
 62:     double         xmin = -1*dd->s-1,xmax = dd->M+dd->s;
 63:     double         x,y;
 64:     PetscInt       base;
 65:     const PetscInt *idx;
 66:     char           node[10];
 67:     PetscBool      isnull;

 70:     PetscViewerDrawGetDraw(viewer,0,&draw);
 71:     PetscDrawIsNull(draw,&isnull);
 72:     if (isnull) return 0;

 74:     PetscDrawCheckResizedWindow(draw);
 75:     PetscDrawClear(draw);
 76:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);

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

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

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

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

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

161:   (*ctx->func)(0,U,a,ctx->funcctx);
162:   (*ctx->funcisetbase)(U,ctx->funcctx);

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

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

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

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

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

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

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

220:   dd->p = 1;
221:   if (m != PETSC_DECIDE) {
224:   }
225:   if (n != PETSC_DECIDE) {
228:   }

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


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

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

300:   /*
301:    check if the scatter requires more than one process neighbor or wraps around
302:    the domain more than once
303:   */
306:   xe = xs + x;
307:   ye = ys + y;

309:   /* determine ghost region (Xs) and region scattered into (IXs)  */
310:   if (xs-s > 0) {
311:     Xs = xs - s; IXs = xs - s;
312:   } else {
313:     if (bx) {
314:       Xs = xs - s;
315:     } else {
316:       Xs = 0;
317:     }
318:     IXs = 0;
319:   }
320:   if (xe+s <= M) {
321:     Xe = xe + s; IXe = xe + s;
322:   } else {
323:     if (bx) {
324:       Xs = xs - s; 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; IYs = ys - s;
340:   } else {
341:     if (by) {
342:       Ys = ys - s;
343:     } else {
344:       Ys = 0;
345:     }
346:     IYs = 0;
347:   }
348:   if (ye+s <= N) {
349:     Ye = ye + s; IYe = ye + s;
350:   } else {
351:     if (by) {
352:       Ye = ye + s;
353:     } else {
354:       Ye = N;
355:     }
356:     IYe = N;
357:   }

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

366:   /* stencil length in each direction */
367:   s_x = s;
368:   s_y = s;

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

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

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

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

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

440:   /* determine who lies on each side of us stored in    n6 n7 n8
441:                                                         n3    n5
442:                                                         n0 n1 n2
443:   */

445:   /* Assume the Non-Periodic Case */
446:   n1 = rank - m;
447:   if (rank % m) {
448:     n0 = n1 - 1;
449:   } else {
450:     n0 = -1;
451:   }
452:   if ((rank+1) % m) {
453:     n2 = n1 + 1;
454:     n5 = rank + 1;
455:     n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
456:   } else {
457:     n2 = -1; n5 = -1; n8 = -1;
458:   }
459:   if (rank % m) {
460:     n3 = rank - 1;
461:     n6 = n3 + m; if (n6 >= m*n) n6 = -1;
462:   } else {
463:     n3 = -1; n6 = -1;
464:   }
465:   n7 = rank + m; 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:   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; sn2 = n2; sn6 = n6; sn8 = n8;
521:     n0  = n2 = n6 = n8 = -1;
522:   }

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

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

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

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

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

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

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

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

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

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

600:   ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);
601:   VecScatterCreate(global,from,local,to,&gtol);
602:   PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
603:   ISDestroy(&to);
604:   ISDestroy(&from);

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

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

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

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

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

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

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

723:   VecDestroy(&local);
724:   VecDestroy(&global);

726:   dd->gtol      = gtol;
727:   dd->base      = base;
728:   da->ops->view = DMView_DA_2d;
729:   dd->ltol      = NULL;
730:   dd->ao        = NULL;
731:   return 0;
732: }

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

738:    Collective

740:    Input Parameters:
741: +  comm - MPI communicator
742: .  bx,by - type of ghost nodes the array have.
743:          Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC.
744: .  stencil_type - stencil type.  Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
745: .  M,N - global dimension in each direction of the array
746: .  m,n - corresponding number of processors in each dimension
747:          (or PETSC_DECIDE to have calculated)
748: .  dof - number of degrees of freedom per node
749: .  s - stencil width
750: -  lx, ly - arrays containing the number of nodes in each cell along
751:            the x and y coordinates, or NULL. If non-null, these
752:            must be of length as m and n, and the corresponding
753:            m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
754:            must be M, and the sum of the ly[] entries must be N.

756:    Output Parameter:
757: .  da - the resulting distributed array object

759:    Options Database Key:
760: +  -dm_view - Calls DMView() at the conclusion of DMDACreate2d()
761: .  -da_grid_x <nx> - number of grid points in x direction
762: .  -da_grid_y <ny> - number of grid points in y direction
763: .  -da_processors_x <nx> - number of processors in x direction
764: .  -da_processors_y <ny> - number of processors in y direction
765: .  -da_refine_x <rx> - refinement ratio in x direction
766: .  -da_refine_y <ry> - refinement ratio in y direction
767: -  -da_refine <n> - refine the DMDA n times before creating

769:    Level: beginner

771:    Notes:
772:    The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
773:    standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
774:    the standard 9-pt stencil.

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

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

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

785: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
786:           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(),
787:           DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges(),
788:           DMStagCreate2d()

790: @*/

792: PetscErrorCode  DMDACreate2d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMDAStencilType stencil_type,
793:                              PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
794: {
795:   DMDACreate(comm, da);
796:   DMSetDimension(*da, 2);
797:   DMDASetSizes(*da, M, N, 1);
798:   DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
799:   DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);
800:   DMDASetDof(*da, dof);
801:   DMDASetStencilType(*da, stencil_type);
802:   DMDASetStencilWidth(*da, s);
803:   DMDASetOwnershipRanges(*da, lx, ly, NULL);
804:   return 0;
805: }