Actual source code: da2.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

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

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

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

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

 28:     PetscViewerGetFormat(viewer, &format);
 29:     if (format == PETSC_VIEWER_LOAD_BALANCE) {
 30:       PetscInt      i,nmax = 0,nmin = PETSC_MAX_INT,navg = 0,*nz,nzlocal;
 31:       DMDALocalInfo info;
 32:       PetscMPIInt   size;
 33:       MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 34:       DMDAGetLocalInfo(da,&info);
 35:       nzlocal = info.xm*info.ym;
 36:       PetscMalloc1(size,&nz);
 37:       MPI_Allgather(&nzlocal,1,MPIU_INT,nz,1,MPIU_INT,PetscObjectComm((PetscObject)da));
 38:       for (i=0; i<(PetscInt)size; i++) {
 39:         nmax = PetscMax(nmax,nz[i]);
 40:         nmin = PetscMin(nmin,nz[i]);
 41:         navg += nz[i];
 42:       }
 43:       PetscFree(nz);
 44:       navg = navg/size;
 45:       PetscViewerASCIIPrintf(viewer,"  Load Balance - Grid Points: Min %D  avg %D  max %D\n",nmin,navg,nmax);
 46:       return(0);
 47:     }
 48:     if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL && format != PETSC_VIEWER_ASCII_GLVIS) {
 49:       DMDALocalInfo info;
 50:       DMDAGetLocalInfo(da,&info);
 51:       PetscViewerASCIIPushSynchronized(viewer);
 52:       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);
 53:       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);
 54:       PetscViewerFlush(viewer);
 55:       PetscViewerASCIIPopSynchronized(viewer);
 56:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
 57:       DMView_DA_GLVis(da,viewer);
 58:     } else {
 59:       DMView_DA_VTK(da,viewer);
 60:     }
 61:   } else if (isdraw) {
 62:     PetscDraw      draw;
 63:     double         ymin = -1*dd->s-1,ymax = dd->N+dd->s;
 64:     double         xmin = -1*dd->s-1,xmax = dd->M+dd->s;
 65:     double         x,y;
 66:     PetscInt       base;
 67:     const PetscInt *idx;
 68:     char           node[10];
 69:     PetscBool      isnull;

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

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

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

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

114:     PetscDrawCollectiveBegin(draw);
115:     /* overlay ghost numbers, useful for error checking */
116:     ISLocalToGlobalMappingGetBlockIndices(da->ltogmap,&idx);
117:     base = 0; xmin = dd->Xs; xmax = dd->Xe; ymin = dd->Ys; 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_ENGINE)
137:   } else if (ismatlab) {
138:     DMView_DA_Matlab(da,viewer);
139: #endif
140:   }
141:   return(0);
142: }


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

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

165:   (*ctx->func)(0,U,a,ctx->funcctx);
166:   (*ctx->funcisetbase)(U,ctx->funcctx);

168:   VecGetArray(U,&ww);
169:   VecGetArray(a,&aa);

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

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

184:   VecRestoreArray(U,&ww);
185:   VecRestoreArray(a,&aa);
186:   return(0);
187: }
188: #endif

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

217:   if (stencil_type == DMDA_STENCIL_BOX && (bx == DM_BOUNDARY_MIRROR || by == DM_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
218:   PetscObjectGetComm((PetscObject)da,&comm);
219: #if !defined(PETSC_USE_64BIT_INDICES)
220:   if (((PetscInt64) M)*((PetscInt64) N)*((PetscInt64) dof) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof);
221: #endif

223:   if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
224:   if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);

226:   MPI_Comm_size(comm,&size);
227:   MPI_Comm_rank(comm,&rank);

229:   dd->p = 1;
230:   if (m != PETSC_DECIDE) {
231:     if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
232:     else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
233:   }
234:   if (n != PETSC_DECIDE) {
235:     if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
236:     else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
237:   }

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

258:   if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
259:   if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);

261:   /*
262:      Determine locally owned region
263:      xs is the first local node number, x is the number of local nodes
264:   */
265:   if (!lx) {
266:     PetscMalloc1(m, &dd->lx);
267:     lx   = dd->lx;
268:     for (i=0; i<m; i++) {
269:       lx[i] = M/m + ((M % m) > i);
270:     }
271:   }
272:   x  = lx[rank % m];
273:   xs = 0;
274:   for (i=0; i<(rank % m); i++) {
275:     xs += lx[i];
276:   }
277: #if defined(PETSC_USE_DEBUG)
278:   left = xs;
279:   for (i=(rank % m); i<m; i++) {
280:     left += lx[i];
281:   }
282:   if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
283: #endif

285:   /*
286:      Determine locally owned region
287:      ys is the first local node number, y is the number of local nodes
288:   */
289:   if (!ly) {
290:     PetscMalloc1(n, &dd->ly);
291:     ly   = dd->ly;
292:     for (i=0; i<n; i++) {
293:       ly[i] = N/n + ((N % n) > i);
294:     }
295:   }
296:   y  = ly[rank/m];
297:   ys = 0;
298:   for (i=0; i<(rank/m); i++) {
299:     ys += ly[i];
300:   }
301: #if defined(PETSC_USE_DEBUG)
302:   left = ys;
303:   for (i=(rank/m); i<n; i++) {
304:     left += ly[i];
305:   }
306:   if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
307: #endif

309:   /*
310:    check if the scatter requires more than one process neighbor or wraps around
311:    the domain more than once
312:   */
313:   if ((x < s) && ((m > 1) || (bx == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);
314:   if ((y < s) && ((n > 1) || (by == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);
315:   xe = xs + x;
316:   ye = ys + y;

318:   /* determine ghost region (Xs) and region scattered into (IXs)  */
319:   if (xs-s > 0) {
320:     Xs = xs - s; IXs = xs - s;
321:   } else {
322:     if (bx) {
323:       Xs = xs - s;
324:     } else {
325:       Xs = 0;
326:     }
327:     IXs = 0;
328:   }
329:   if (xe+s <= M) {
330:     Xe = xe + s; IXe = xe + s;
331:   } else {
332:     if (bx) {
333:       Xs = xs - s; Xe = xe + s;
334:     } else {
335:       Xe = M;
336:     }
337:     IXe = M;
338:   }

340:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
341:     IXs = xs - s;
342:     IXe = xe + s;
343:     Xs  = xs - s;
344:     Xe  = xe + s;
345:   }

347:   if (ys-s > 0) {
348:     Ys = ys - s; IYs = ys - s;
349:   } else {
350:     if (by) {
351:       Ys = ys - s;
352:     } else {
353:       Ys = 0;
354:     }
355:     IYs = 0;
356:   }
357:   if (ye+s <= N) {
358:     Ye = ye + s; IYe = ye + s;
359:   } else {
360:     if (by) {
361:       Ye = ye + s;
362:     } else {
363:       Ye = N;
364:     }
365:     IYe = N;
366:   }

368:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
369:     IYs = ys - s;
370:     IYe = ye + s;
371:     Ys  = ys - s;
372:     Ye  = ye + s;
373:   }

375:   /* stencil length in each direction */
376:   s_x = s;
377:   s_y = s;

379:   /* determine starting point of each processor */
380:   nn       = x*y;
381:   PetscMalloc2(size+1,&bases,size,&ldims);
382:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
383:   bases[0] = 0;
384:   for (i=1; i<=size; i++) {
385:     bases[i] = ldims[i-1];
386:   }
387:   for (i=1; i<=size; i++) {
388:     bases[i] += bases[i-1];
389:   }
390:   base = bases[rank]*dof;

392:   /* allocate the base parallel and sequential vectors */
393:   dd->Nlocal = x*y*dof;
394:   VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);
395:   dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
396:   VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);

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

400:   /* global to local must include ghost points within the domain,
401:      but not ghost points outside the domain that aren't periodic */
402:   PetscMalloc1((IXe-IXs)*(IYe-IYs),&idx);
403:   if (stencil_type == DMDA_STENCIL_BOX) {
404:     left  = IXs - Xs; right = left + (IXe-IXs);
405:     down  = IYs - Ys; up = down + (IYe-IYs);
406:     count = 0;
407:     for (i=down; i<up; i++) {
408:       for (j=left; j<right; j++) {
409:         idx[count++] = j + i*(Xe-Xs);
410:       }
411:     }
412:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);

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


450:   /* determine who lies on each side of us stored in    n6 n7 n8
451:                                                         n3    n5
452:                                                         n0 n1 n2
453:   */

455:   /* Assume the Non-Periodic Case */
456:   n1 = rank - m;
457:   if (rank % m) {
458:     n0 = n1 - 1;
459:   } else {
460:     n0 = -1;
461:   }
462:   if ((rank+1) % m) {
463:     n2 = n1 + 1;
464:     n5 = rank + 1;
465:     n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
466:   } else {
467:     n2 = -1; n5 = -1; n8 = -1;
468:   }
469:   if (rank % m) {
470:     n3 = rank - 1;
471:     n6 = n3 + m; if (n6 >= m*n) n6 = -1;
472:   } else {
473:     n3 = -1; n6 = -1;
474:   }
475:   n7 = rank + m; if (n7 >= m*n) n7 = -1;

477:   if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) {
478:     /* Modify for Periodic Cases */
479:     /* Handle all four corners */
480:     if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
481:     if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
482:     if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
483:     if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;

485:     /* Handle Top and Bottom Sides */
486:     if (n1 < 0) n1 = rank + m * (n-1);
487:     if (n7 < 0) n7 = rank - m * (n-1);
488:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
489:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
490:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
491:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;

493:     /* 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:   } else if (by == DM_BOUNDARY_PERIODIC) {  /* Handle Top and Bottom Sides */
501:     if (n1 < 0) n1 = rank + m * (n-1);
502:     if (n7 < 0) n7 = rank - m * (n-1);
503:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
504:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
505:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
506:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
507:   } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
508:     if (n3 < 0) n3 = rank + (m-1);
509:     if (n5 < 0) n5 = rank - (m-1);
510:     if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
511:     if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
512:     if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
513:     if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
514:   }

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

518:   dd->neighbors[0] = n0;
519:   dd->neighbors[1] = n1;
520:   dd->neighbors[2] = n2;
521:   dd->neighbors[3] = n3;
522:   dd->neighbors[4] = rank;
523:   dd->neighbors[5] = n5;
524:   dd->neighbors[6] = n6;
525:   dd->neighbors[7] = n7;
526:   dd->neighbors[8] = n8;

528:   if (stencil_type == DMDA_STENCIL_STAR) {
529:     /* save corner processor numbers */
530:     sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
531:     n0  = n2 = n6 = n8 = -1;
532:   }

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

536:   nn = 0;
537:   xbase = bases[rank];
538:   for (i=1; i<=s_y; i++) {
539:     if (n0 >= 0) { /* left below */
540:       x_t = lx[n0 % m];
541:       y_t = ly[(n0/m)];
542:       s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
543:       for (j=0; j<s_x; j++) idx[nn++] = s_t++;
544:     }

546:     if (n1 >= 0) { /* directly below */
547:       x_t = x;
548:       y_t = ly[(n1/m)];
549:       s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
550:       for (j=0; j<x_t; j++) idx[nn++] = s_t++;
551:     } else if (by == DM_BOUNDARY_MIRROR) {
552:       for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1)  + j;
553:     }

555:     if (n2 >= 0) { /* right below */
556:       x_t = lx[n2 % m];
557:       y_t = ly[(n2/m)];
558:       s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
559:       for (j=0; j<s_x; j++) idx[nn++] = s_t++;
560:     }
561:   }

563:   for (i=0; i<y; i++) {
564:     if (n3 >= 0) { /* directly left */
565:       x_t = lx[n3 % m];
566:       /* y_t = y; */
567:       s_t = bases[n3] + (i+1)*x_t - s_x;
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 + s_x - j;
571:     }

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

575:     if (n5 >= 0) { /* directly right */
576:       x_t = lx[n5 % m];
577:       /* y_t = y; */
578:       s_t = bases[n5] + (i)*x_t;
579:       for (j=0; j<s_x; j++) idx[nn++] = s_t++;
580:     } else if (bx == DM_BOUNDARY_MIRROR) {
581:       for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;
582:     }
583:   }

585:   for (i=1; i<=s_y; i++) {
586:     if (n6 >= 0) { /* left above */
587:       x_t = lx[n6 % m];
588:       /* y_t = ly[(n6/m)]; */
589:       s_t = bases[n6] + (i)*x_t - s_x;
590:       for (j=0; j<s_x; j++) idx[nn++] = s_t++;
591:     }

593:     if (n7 >= 0) { /* directly above */
594:       x_t = x;
595:       /* y_t = ly[(n7/m)]; */
596:       s_t = bases[n7] + (i-1)*x_t;
597:       for (j=0; j<x_t; j++) idx[nn++] = s_t++;
598:     } else if (by == DM_BOUNDARY_MIRROR) {
599:       for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1)  + j;
600:     }

602:     if (n8 >= 0) { /* right above */
603:       x_t = lx[n8 % m];
604:       /* y_t = ly[(n8/m)]; */
605:       s_t = bases[n8] + (i-1)*x_t;
606:       for (j=0; j<s_x; j++) idx[nn++] = s_t++;
607:     }
608:   }

610:   ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);
611:   VecScatterCreate(global,from,local,to,&gtol);
612:   PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
613:   ISDestroy(&to);
614:   ISDestroy(&from);

616:   if (stencil_type == DMDA_STENCIL_STAR) {
617:     n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
618:   }

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

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

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

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

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

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

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

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

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

748:    Collective on MPI_Comm

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

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

769:    Options Database Key:
770: +  -dm_view - Calls DMView() at the conclusion of DMDACreate2d()
771: .  -da_grid_x <nx> - number of grid points in x direction
772: .  -da_grid_y <ny> - number of grid points in y direction
773: .  -da_processors_x <nx> - number of processors in x direction
774: .  -da_processors_y <ny> - number of processors in y direction
775: .  -da_refine_x <rx> - refinement ratio in x direction
776: .  -da_refine_y <ry> - refinement ratio in y direction
777: -  -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: .keywords: distributed array, create, two-dimensional

798: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
799:           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(),
800:           DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()

802: @*/

804: PetscErrorCode  DMDACreate2d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMDAStencilType stencil_type,
805:                           PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
806: {

810:   DMDACreate(comm, da);
811:   DMSetDimension(*da, 2);
812:   DMDASetSizes(*da, M, N, 1);
813:   DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
814:   DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);
815:   DMDASetDof(*da, dof);
816:   DMDASetStencilType(*da, stencil_type);
817:   DMDASetStencilWidth(*da, s);
818:   DMDASetOwnershipRanges(*da, lx, ly, NULL);
819:   return(0);
820: }