Actual source code: da3.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  2: /*
  3:    Code for manipulating distributed regular 3d arrays in parallel.
  4:    File created by Peter Mell  7/14/95
  5:  */

  7:  #include <petsc/private/dmdaimpl.h>

  9:  #include <petscdraw.h>
 10: static PetscErrorCode DMView_DA_3d(DM da,PetscViewer viewer)
 11: {
 13:   PetscMPIInt    rank;
 14:   PetscBool      iascii,isdraw,isglvis,isbinary;
 15:   DM_DA          *dd = (DM_DA*)da->data;
 16: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 17:   PetscBool ismatlab;
 18: #endif

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

 23:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 24:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
 25:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
 26:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
 27: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 28:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
 29: #endif
 30:   if (iascii) {
 31:     PetscViewerFormat format;

 33:     PetscViewerASCIIPushSynchronized(viewer);
 34:     PetscViewerGetFormat(viewer, &format);
 35:     if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL && format != PETSC_VIEWER_ASCII_GLVIS) {
 36:       DMDALocalInfo info;
 37:       DMDAGetLocalInfo(da,&info);
 38:       PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D P %D m %D n %D p %D w %D s %D\n",rank,dd->M,dd->N,dd->P,dd->m,dd->n,dd->p,dd->w,dd->s);
 39:       PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D, Z range of indices: %D %D\n",
 40:                                                 info.xs,info.xs+info.xm,info.ys,info.ys+info.ym,info.zs,info.zs+info.zm);
 41: #if !defined(PETSC_USE_COMPLEX)
 42:       if (da->coordinates) {
 43:         PetscInt        last;
 44:         const PetscReal *coors;
 45:         VecGetArrayRead(da->coordinates,&coors);
 46:         VecGetLocalSize(da->coordinates,&last);
 47:         last = last - 3;
 48:         PetscViewerASCIISynchronizedPrintf(viewer,"Lower left corner %g %g %g : Upper right %g %g %g\n",(double)coors[0],(double)coors[1],(double)coors[2],(double)coors[last],(double)coors[last+1],(double)coors[last+2]);
 49:         VecRestoreArrayRead(da->coordinates,&coors);
 50:       }
 51: #endif
 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:     PetscReal      ymin = -1.0,ymax = (PetscReal)dd->N;
 62:     PetscReal      xmin = -1.0,xmax = (PetscReal)((dd->M+2)*dd->P),x,y,ycoord,xcoord;
 63:     PetscInt       k,plane,base;
 64:     const PetscInt *idx;
 65:     char           node[10];
 66:     PetscBool      isnull;

 68:     PetscViewerDrawGetDraw(viewer,0,&draw);
 69:     PetscDrawIsNull(draw,&isnull);
 70:     if (isnull) return(0);

 72:     PetscDrawCheckResizedWindow(draw);
 73:     PetscDrawClear(draw);
 74:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);

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

 94:     PetscDrawCollectiveBegin(draw);
 95:     /*Go through and draw for each plane*/
 96:     for (k=0; k<dd->P; k++) {
 97:       if ((k >= dd->zs) && (k < dd->ze)) {
 98:         /* draw my box */
 99:         ymin = dd->ys;
100:         ymax = dd->ye - 1;
101:         xmin = dd->xs/dd->w    + (dd->M+1)*k;
102:         xmax =(dd->xe-1)/dd->w + (dd->M+1)*k;

104:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
105:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
106:         PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
107:         PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

109:         xmin = dd->xs/dd->w;
110:         xmax =(dd->xe-1)/dd->w;

112:         /* identify which processor owns the box */
113:         PetscSNPrintf(node,sizeof(node),"%d",(int)rank);
114:         PetscDrawString(draw,xmin+(dd->M+1)*k+.2,ymin+.3,PETSC_DRAW_RED,node);
115:         /* put in numbers*/
116:         base = (dd->base+(dd->xe-dd->xs)*(dd->ye-dd->ys)*(k-dd->zs))/dd->w;
117:         for (y=ymin; y<=ymax; y++) {
118:           for (x=xmin+(dd->M+1)*k; x<=xmax+(dd->M+1)*k; x++) {
119:             PetscSNPrintf(node,sizeof(node),"%d",(int)base++);
120:             PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
121:           }
122:         }

124:       }
125:     }
126:     PetscDrawCollectiveEnd(draw);
127:     PetscDrawFlush(draw);
128:     PetscDrawPause(draw);

130:     PetscDrawCollectiveBegin(draw);
131:     for (k=0-dd->s; k<dd->P+dd->s; k++) {
132:       /* Go through and draw for each plane */
133:       if ((k >= dd->Zs) && (k < dd->Ze)) {
134:         /* overlay ghost numbers, useful for error checking */
135:         base = (dd->Xe-dd->Xs)*(dd->Ye-dd->Ys)*(k-dd->Zs)/dd->w;
136:         ISLocalToGlobalMappingGetBlockIndices(da->ltogmap,&idx);
137:         plane=k;
138:         /* Keep z wrap around points on the drawing */
139:         if (k<0) plane=dd->P+k;
140:         if (k>=dd->P) plane=k-dd->P;
141:         ymin = dd->Ys; ymax = dd->Ye;
142:         xmin = (dd->M+1)*plane*dd->w;
143:         xmax = (dd->M+1)*plane*dd->w+dd->M*dd->w;
144:         for (y=ymin; y<ymax; y++) {
145:           for (x=xmin+dd->Xs; x<xmin+dd->Xe; x+=dd->w) {
146:             sprintf(node,"%d",(int)(idx[base]));
147:             ycoord = y;
148:             /*Keep y wrap around points on drawing */
149:             if (y<0) ycoord = dd->N+y;
150:             if (y>=dd->N) ycoord = y-dd->N;
151:             xcoord = x;   /* Keep x wrap points on drawing */
152:             if (x<xmin) xcoord = xmax - (xmin-x);
153:             if (x>=xmax) xcoord = xmin + (x-xmax);
154:             PetscDrawString(draw,xcoord/dd->w,ycoord,PETSC_DRAW_BLUE,node);
155:             base++;
156:           }
157:         }
158:         ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap,&idx);
159:       }
160:     }
161:     PetscDrawCollectiveEnd(draw);
162:     PetscDrawFlush(draw);
163:     PetscDrawPause(draw);
164:     PetscDrawSave(draw);
165:   } else if (isglvis) {
166:     DMView_DA_GLVis(da,viewer);
167:   } else if (isbinary) {
168:     DMView_DA_Binary(da,viewer);
169: #if defined(PETSC_HAVE_MATLAB_ENGINE)
170:   } else if (ismatlab) {
171:     DMView_DA_Matlab(da,viewer);
172: #endif
173:   }
174:   return(0);
175: }

177: PetscErrorCode  DMSetUp_DA_3D(DM da)
178: {
179:   DM_DA            *dd          = (DM_DA*)da->data;
180:   const PetscInt   M            = dd->M;
181:   const PetscInt   N            = dd->N;
182:   const PetscInt   P            = dd->P;
183:   PetscInt         m            = dd->m;
184:   PetscInt         n            = dd->n;
185:   PetscInt         p            = dd->p;
186:   const PetscInt   dof          = dd->w;
187:   const PetscInt   s            = dd->s;
188:   DMBoundaryType   bx           = dd->bx;
189:   DMBoundaryType   by           = dd->by;
190:   DMBoundaryType   bz           = dd->bz;
191:   DMDAStencilType  stencil_type = dd->stencil_type;
192:   PetscInt         *lx          = dd->lx;
193:   PetscInt         *ly          = dd->ly;
194:   PetscInt         *lz          = dd->lz;
195:   MPI_Comm         comm;
196:   PetscMPIInt      rank,size;
197:   PetscInt         xs = 0,xe,ys = 0,ye,zs = 0,ze,x = 0,y = 0,z = 0;
198:   PetscInt         Xs,Xe,Ys,Ye,Zs,Ze,IXs,IXe,IYs,IYe,IZs,IZe,pm;
199:   PetscInt         left,right,up,down,bottom,top,i,j,k,*idx,nn;
200:   PetscInt         n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n14;
201:   PetscInt         n15,n16,n17,n18,n19,n20,n21,n22,n23,n24,n25,n26;
202:   PetscInt         *bases,*ldims,base,x_t,y_t,z_t,s_t,count,s_x,s_y,s_z;
203:   PetscInt         sn0  = 0,sn1 = 0,sn2 = 0,sn3 = 0,sn5 = 0,sn6 = 0,sn7 = 0;
204:   PetscInt         sn8  = 0,sn9 = 0,sn11 = 0,sn15 = 0,sn24 = 0,sn25 = 0,sn26 = 0;
205:   PetscInt         sn17 = 0,sn18 = 0,sn19 = 0,sn20 = 0,sn21 = 0,sn23 = 0;
206:   Vec              local,global;
207:   VecScatter       gtol;
208:   IS               to,from;
209:   PetscBool        twod;
210:   PetscErrorCode   ierr;


214:   if (stencil_type == DMDA_STENCIL_BOX && (bx == DM_BOUNDARY_MIRROR || by == DM_BOUNDARY_MIRROR || bz == DM_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
215:   if ((bx == DM_BOUNDARY_MIRROR) || (by == DM_BOUNDARY_MIRROR) || (bz == DM_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary not supported yet in 3d");
216:   PetscObjectGetComm((PetscObject) da, &comm);
217: #if !defined(PETSC_USE_64BIT_INDICES)
218:   if (((PetscInt64) M)*((PetscInt64) N)*((PetscInt64) P)*((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);
219: #endif

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

224:   if (m != PETSC_DECIDE) {
225:     if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
226:     else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
227:   }
228:   if (n != PETSC_DECIDE) {
229:     if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
230:     else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
231:   }
232:   if (p != PETSC_DECIDE) {
233:     if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
234:     else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
235:   }
236:   if ((m > 0) && (n > 0) && (p > 0) && (m*n*p != size)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m %D * n %D * p %D != size %d",m,n,p,size);

238:   /* Partition the array among the processors */
239:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
240:     m = size/(n*p);
241:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
242:     n = size/(m*p);
243:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
244:     p = size/(m*n);
245:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
246:     /* try for squarish distribution */
247:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
248:     if (!m) m = 1;
249:     while (m > 0) {
250:       n = size/(m*p);
251:       if (m*n*p == size) break;
252:       m--;
253:     }
254:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
255:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
256:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
257:     /* try for squarish distribution */
258:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
259:     if (!m) m = 1;
260:     while (m > 0) {
261:       p = size/(m*n);
262:       if (m*n*p == size) break;
263:       m--;
264:     }
265:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
266:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
267:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
268:     /* try for squarish distribution */
269:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
270:     if (!n) n = 1;
271:     while (n > 0) {
272:       p = size/(m*n);
273:       if (m*n*p == size) break;
274:       n--;
275:     }
276:     if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
277:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
278:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
279:     /* try for squarish distribution */
280:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
281:     if (!n) n = 1;
282:     while (n > 0) {
283:       pm = size/n;
284:       if (n*pm == size) break;
285:       n--;
286:     }
287:     if (!n) n = 1;
288:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
289:     if (!m) m = 1;
290:     while (m > 0) {
291:       p = size/(m*n);
292:       if (m*n*p == size) break;
293:       m--;
294:     }
295:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
296:   } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

298:   if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_PLIB,"Could not find good partition");
299:   if (M < m) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
300:   if (N < n) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
301:   if (P < p) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in z direction is too fine! %D %D",P,p);

303:   /*
304:      Determine locally owned region
305:      [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
306:   */

308:   if (!lx) {
309:     PetscMalloc1(m, &dd->lx);
310:     lx   = dd->lx;
311:     for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m));
312:   }
313:   x  = lx[rank % m];
314:   xs = 0;
315:   for (i=0; i<(rank%m); i++) xs += lx[i];
316:   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);

318:   if (!ly) {
319:     PetscMalloc1(n, &dd->ly);
320:     ly   = dd->ly;
321:     for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n));
322:   }
323:   y = ly[(rank % (m*n))/m];
324:   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);

326:   ys = 0;
327:   for (i=0; i<(rank % (m*n))/m; i++) ys += ly[i];

329:   if (!lz) {
330:     PetscMalloc1(p, &dd->lz);
331:     lz = dd->lz;
332:     for (i=0; i<p; i++) lz[i] = P/p + ((P % p) > (i % p));
333:   }
334:   z = lz[rank/(m*n)];

336:   /* note this is different than x- and y-, as we will handle as an important special
337:    case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z.  This is to deal with 2D problems
338:    in a 3D code.  Additional code for this case is noted with "2d case" comments */
339:   twod = PETSC_FALSE;
340:   if (P == 1) twod = PETSC_TRUE;
341:   else if ((z < s) && ((p > 1) || (bz == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local z-width of domain z %D is smaller than stencil width s %D",z,s);
342:   zs = 0;
343:   for (i=0; i<(rank/(m*n)); i++) zs += lz[i];
344:   ye = ys + y;
345:   xe = xs + x;
346:   ze = zs + z;

348:   /* determine ghost region (Xs) and region scattered into (IXs)  */
349:   if (xs-s > 0) {
350:     Xs = xs - s; IXs = xs - s;
351:   } else {
352:     if (bx) Xs = xs - s;
353:     else Xs = 0;
354:     IXs = 0;
355:   }
356:   if (xe+s <= M) {
357:     Xe = xe + s; IXe = xe + s;
358:   } else {
359:     if (bx) {
360:       Xs = xs - s; Xe = xe + s;
361:     } else Xe = M;
362:     IXe = M;
363:   }

365:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
366:     IXs = xs - s;
367:     IXe = xe + s;
368:     Xs  = xs - s;
369:     Xe  = xe + s;
370:   }

372:   if (ys-s > 0) {
373:     Ys = ys - s; IYs = ys - s;
374:   } else {
375:     if (by) Ys = ys - s;
376:     else Ys = 0;
377:     IYs = 0;
378:   }
379:   if (ye+s <= N) {
380:     Ye = ye + s; IYe = ye + s;
381:   } else {
382:     if (by) Ye = ye + s;
383:     else Ye = N;
384:     IYe = N;
385:   }

387:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
388:     IYs = ys - s;
389:     IYe = ye + s;
390:     Ys  = ys - s;
391:     Ye  = ye + s;
392:   }

394:   if (zs-s > 0) {
395:     Zs = zs - s; IZs = zs - s;
396:   } else {
397:     if (bz) Zs = zs - s;
398:     else Zs = 0;
399:     IZs = 0;
400:   }
401:   if (ze+s <= P) {
402:     Ze = ze + s; IZe = ze + s;
403:   } else {
404:     if (bz) Ze = ze + s;
405:     else Ze = P;
406:     IZe = P;
407:   }

409:   if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
410:     IZs = zs - s;
411:     IZe = ze + s;
412:     Zs  = zs - s;
413:     Ze  = ze + s;
414:   }

416:   /* Resize all X parameters to reflect w */
417:   s_x = s;
418:   s_y = s;
419:   s_z = s;

421:   /* determine starting point of each processor */
422:   nn       = x*y*z;
423:   PetscMalloc2(size+1,&bases,size,&ldims);
424:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
425:   bases[0] = 0;
426:   for (i=1; i<=size; i++) bases[i] = ldims[i-1];
427:   for (i=1; i<=size; i++) bases[i] += bases[i-1];
428:   base = bases[rank]*dof;

430:   /* allocate the base parallel and sequential vectors */
431:   dd->Nlocal = x*y*z*dof;
432:   VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);
433:   dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs)*dof;
434:   VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);

436:   /* generate appropriate vector scatters */
437:   /* local to global inserts non-ghost point region into global */
438:   PetscMalloc1((IXe-IXs)*(IYe-IYs)*(IZe-IZs),&idx);
439:   left   = xs - Xs; right = left + x;
440:   bottom = ys - Ys; top = bottom + y;
441:   down   = zs - Zs; up  = down + z;
442:   count  = 0;
443:   for (i=down; i<up; i++) {
444:     for (j=bottom; j<top; j++) {
445:       for (k=left; k<right; k++) {
446:         idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
447:       }
448:     }
449:   }

451:   /* global to local must include ghost points within the domain,
452:      but not ghost points outside the domain that aren't periodic */
453:   if (stencil_type == DMDA_STENCIL_BOX) {
454:     left   = IXs - Xs; right = left + (IXe-IXs);
455:     bottom = IYs - Ys; top = bottom + (IYe-IYs);
456:     down   = IZs - Zs; up  = down + (IZe-IZs);
457:     count  = 0;
458:     for (i=down; i<up; i++) {
459:       for (j=bottom; j<top; j++) {
460:         for (k=left; k<right; k++) {
461:           idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
462:         }
463:       }
464:     }
465:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
466:   } else {
467:     /* This is way ugly! We need to list the funny cross type region */
468:     left   = xs - Xs; right = left + x;
469:     bottom = ys - Ys; top = bottom + y;
470:     down   = zs - Zs;   up  = down + z;
471:     count  = 0;
472:     /* the bottom chunck */
473:     for (i=(IZs-Zs); i<down; i++) {
474:       for (j=bottom; j<top; j++) {
475:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
476:       }
477:     }
478:     /* the middle piece */
479:     for (i=down; i<up; i++) {
480:       /* front */
481:       for (j=(IYs-Ys); j<bottom; j++) {
482:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
483:       }
484:       /* middle */
485:       for (j=bottom; j<top; j++) {
486:         for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
487:       }
488:       /* back */
489:       for (j=top; j<top+IYe-ye; j++) {
490:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
491:       }
492:     }
493:     /* the top piece */
494:     for (i=up; i<up+IZe-ze; i++) {
495:       for (j=bottom; j<top; j++) {
496:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
497:       }
498:     }
499:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
500:   }

502:   /* determine who lies on each side of use stored in    n24 n25 n26
503:                                                          n21 n22 n23
504:                                                          n18 n19 n20

506:                                                          n15 n16 n17
507:                                                          n12     n14
508:                                                          n9  n10 n11

510:                                                          n6  n7  n8
511:                                                          n3  n4  n5
512:                                                          n0  n1  n2
513:   */

515:   /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
516:   /* Assume Nodes are Internal to the Cube */
517:   n0 = rank - m*n - m - 1;
518:   n1 = rank - m*n - m;
519:   n2 = rank - m*n - m + 1;
520:   n3 = rank - m*n -1;
521:   n4 = rank - m*n;
522:   n5 = rank - m*n + 1;
523:   n6 = rank - m*n + m - 1;
524:   n7 = rank - m*n + m;
525:   n8 = rank - m*n + m + 1;

527:   n9  = rank - m - 1;
528:   n10 = rank - m;
529:   n11 = rank - m + 1;
530:   n12 = rank - 1;
531:   n14 = rank + 1;
532:   n15 = rank + m - 1;
533:   n16 = rank + m;
534:   n17 = rank + m + 1;

536:   n18 = rank + m*n - m - 1;
537:   n19 = rank + m*n - m;
538:   n20 = rank + m*n - m + 1;
539:   n21 = rank + m*n - 1;
540:   n22 = rank + m*n;
541:   n23 = rank + m*n + 1;
542:   n24 = rank + m*n + m - 1;
543:   n25 = rank + m*n + m;
544:   n26 = rank + m*n + m + 1;

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

548:   if (xs == 0) { /* First assume not corner or edge */
549:     n0  = rank       -1 - (m*n);
550:     n3  = rank + m   -1 - (m*n);
551:     n6  = rank + 2*m -1 - (m*n);
552:     n9  = rank       -1;
553:     n12 = rank + m   -1;
554:     n15 = rank + 2*m -1;
555:     n18 = rank       -1 + (m*n);
556:     n21 = rank + m   -1 + (m*n);
557:     n24 = rank + 2*m -1 + (m*n);
558:   }

560:   if (xe == M) { /* First assume not corner or edge */
561:     n2  = rank -2*m +1 - (m*n);
562:     n5  = rank - m  +1 - (m*n);
563:     n8  = rank      +1 - (m*n);
564:     n11 = rank -2*m +1;
565:     n14 = rank - m  +1;
566:     n17 = rank      +1;
567:     n20 = rank -2*m +1 + (m*n);
568:     n23 = rank - m  +1 + (m*n);
569:     n26 = rank      +1 + (m*n);
570:   }

572:   if (ys==0) { /* First assume not corner or edge */
573:     n0  = rank + m * (n-1) -1 - (m*n);
574:     n1  = rank + m * (n-1)    - (m*n);
575:     n2  = rank + m * (n-1) +1 - (m*n);
576:     n9  = rank + m * (n-1) -1;
577:     n10 = rank + m * (n-1);
578:     n11 = rank + m * (n-1) +1;
579:     n18 = rank + m * (n-1) -1 + (m*n);
580:     n19 = rank + m * (n-1)    + (m*n);
581:     n20 = rank + m * (n-1) +1 + (m*n);
582:   }

584:   if (ye == N) { /* First assume not corner or edge */
585:     n6  = rank - m * (n-1) -1 - (m*n);
586:     n7  = rank - m * (n-1)    - (m*n);
587:     n8  = rank - m * (n-1) +1 - (m*n);
588:     n15 = rank - m * (n-1) -1;
589:     n16 = rank - m * (n-1);
590:     n17 = rank - m * (n-1) +1;
591:     n24 = rank - m * (n-1) -1 + (m*n);
592:     n25 = rank - m * (n-1)    + (m*n);
593:     n26 = rank - m * (n-1) +1 + (m*n);
594:   }

596:   if (zs == 0) { /* First assume not corner or edge */
597:     n0 = size - (m*n) + rank - m - 1;
598:     n1 = size - (m*n) + rank - m;
599:     n2 = size - (m*n) + rank - m + 1;
600:     n3 = size - (m*n) + rank - 1;
601:     n4 = size - (m*n) + rank;
602:     n5 = size - (m*n) + rank + 1;
603:     n6 = size - (m*n) + rank + m - 1;
604:     n7 = size - (m*n) + rank + m;
605:     n8 = size - (m*n) + rank + m + 1;
606:   }

608:   if (ze == P) { /* First assume not corner or edge */
609:     n18 = (m*n) - (size-rank) - m - 1;
610:     n19 = (m*n) - (size-rank) - m;
611:     n20 = (m*n) - (size-rank) - m + 1;
612:     n21 = (m*n) - (size-rank) - 1;
613:     n22 = (m*n) - (size-rank);
614:     n23 = (m*n) - (size-rank) + 1;
615:     n24 = (m*n) - (size-rank) + m - 1;
616:     n25 = (m*n) - (size-rank) + m;
617:     n26 = (m*n) - (size-rank) + m + 1;
618:   }

620:   if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
621:     n0 = size - m*n + rank + m-1 - m;
622:     n3 = size - m*n + rank + m-1;
623:     n6 = size - m*n + rank + m-1 + m;
624:   }

626:   if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
627:     n18 = m*n - (size - rank) + m-1 - m;
628:     n21 = m*n - (size - rank) + m-1;
629:     n24 = m*n - (size - rank) + m-1 + m;
630:   }

632:   if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
633:     n0  = rank + m*n -1 - m*n;
634:     n9  = rank + m*n -1;
635:     n18 = rank + m*n -1 + m*n;
636:   }

638:   if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
639:     n6  = rank - m*(n-1) + m-1 - m*n;
640:     n15 = rank - m*(n-1) + m-1;
641:     n24 = rank - m*(n-1) + m-1 + m*n;
642:   }

644:   if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */
645:     n2 = size - (m*n-rank) - (m-1) - m;
646:     n5 = size - (m*n-rank) - (m-1);
647:     n8 = size - (m*n-rank) - (m-1) + m;
648:   }

650:   if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */
651:     n20 = m*n - (size - rank) - (m-1) - m;
652:     n23 = m*n - (size - rank) - (m-1);
653:     n26 = m*n - (size - rank) - (m-1) + m;
654:   }

656:   if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */
657:     n2  = rank + m*(n-1) - (m-1) - m*n;
658:     n11 = rank + m*(n-1) - (m-1);
659:     n20 = rank + m*(n-1) - (m-1) + m*n;
660:   }

662:   if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */
663:     n8  = rank - m*n +1 - m*n;
664:     n17 = rank - m*n +1;
665:     n26 = rank - m*n +1 + m*n;
666:   }

668:   if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
669:     n0 = size - m + rank -1;
670:     n1 = size - m + rank;
671:     n2 = size - m + rank +1;
672:   }

674:   if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
675:     n18 = m*n - (size - rank) + m*(n-1) -1;
676:     n19 = m*n - (size - rank) + m*(n-1);
677:     n20 = m*n - (size - rank) + m*(n-1) +1;
678:   }

680:   if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
681:     n6 = size - (m*n-rank) - m * (n-1) -1;
682:     n7 = size - (m*n-rank) - m * (n-1);
683:     n8 = size - (m*n-rank) - m * (n-1) +1;
684:   }

686:   if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
687:     n24 = rank - (size-m) -1;
688:     n25 = rank - (size-m);
689:     n26 = rank - (size-m) +1;
690:   }

692:   /* Check for Corners */
693:   if ((xs==0) && (ys==0) && (zs==0)) n0  = size -1;
694:   if ((xs==0) && (ys==0) && (ze==P)) n18 = m*n-1;
695:   if ((xs==0) && (ye==N) && (zs==0)) n6  = (size-1)-m*(n-1);
696:   if ((xs==0) && (ye==N) && (ze==P)) n24 = m-1;
697:   if ((xe==M) && (ys==0) && (zs==0)) n2  = size-m;
698:   if ((xe==M) && (ys==0) && (ze==P)) n20 = m*n-m;
699:   if ((xe==M) && (ye==N) && (zs==0)) n8  = size-m*n;
700:   if ((xe==M) && (ye==N) && (ze==P)) n26 = 0;

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

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

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

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

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

724:   dd->neighbors[0]  = n0;
725:   dd->neighbors[1]  = n1;
726:   dd->neighbors[2]  = n2;
727:   dd->neighbors[3]  = n3;
728:   dd->neighbors[4]  = n4;
729:   dd->neighbors[5]  = n5;
730:   dd->neighbors[6]  = n6;
731:   dd->neighbors[7]  = n7;
732:   dd->neighbors[8]  = n8;
733:   dd->neighbors[9]  = n9;
734:   dd->neighbors[10] = n10;
735:   dd->neighbors[11] = n11;
736:   dd->neighbors[12] = n12;
737:   dd->neighbors[13] = rank;
738:   dd->neighbors[14] = n14;
739:   dd->neighbors[15] = n15;
740:   dd->neighbors[16] = n16;
741:   dd->neighbors[17] = n17;
742:   dd->neighbors[18] = n18;
743:   dd->neighbors[19] = n19;
744:   dd->neighbors[20] = n20;
745:   dd->neighbors[21] = n21;
746:   dd->neighbors[22] = n22;
747:   dd->neighbors[23] = n23;
748:   dd->neighbors[24] = n24;
749:   dd->neighbors[25] = n25;
750:   dd->neighbors[26] = n26;

752:   /* If star stencil then delete the corner neighbors */
753:   if (stencil_type == DMDA_STENCIL_STAR) {
754:     /* save information about corner neighbors */
755:     sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7;
756:     sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18;
757:     sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25;
758:     sn26 = n26;
759:     n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
760:   }

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

764:   nn = 0;
765:   /* Bottom Level */
766:   for (k=0; k<s_z; k++) {
767:     for (i=1; i<=s_y; i++) {
768:       if (n0 >= 0) { /* left below */
769:         x_t = lx[n0 % m];
770:         y_t = ly[(n0 % (m*n))/m];
771:         z_t = lz[n0 / (m*n)];
772:         s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x - (s_z-k-1)*x_t*y_t;
773:         if (twod && (s_t < 0)) s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x; /* 2D case */
774:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
775:       }
776:       if (n1 >= 0) { /* directly below */
777:         x_t = x;
778:         y_t = ly[(n1 % (m*n))/m];
779:         z_t = lz[n1 / (m*n)];
780:         s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
781:         if (twod && (s_t < 0)) s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */
782:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
783:       }
784:       if (n2 >= 0) { /* right below */
785:         x_t = lx[n2 % m];
786:         y_t = ly[(n2 % (m*n))/m];
787:         z_t = lz[n2 / (m*n)];
788:         s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
789:         if (twod && (s_t < 0)) s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */
790:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
791:       }
792:     }

794:     for (i=0; i<y; i++) {
795:       if (n3 >= 0) { /* directly left */
796:         x_t = lx[n3 % m];
797:         y_t = y;
798:         z_t = lz[n3 / (m*n)];
799:         s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
800:         if (twod && (s_t < 0)) s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - x_t*y_t; /* 2D case */
801:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
802:       }

804:       if (n4 >= 0) { /* middle */
805:         x_t = x;
806:         y_t = y;
807:         z_t = lz[n4 / (m*n)];
808:         s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
809:         if (twod && (s_t < 0)) s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */
810:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
811:       } else if (bz == DM_BOUNDARY_MIRROR) {
812:         for (j=0; j<x; j++) idx[nn++] = 0;
813:       }

815:       if (n5 >= 0) { /* directly right */
816:         x_t = lx[n5 % m];
817:         y_t = y;
818:         z_t = lz[n5 / (m*n)];
819:         s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
820:         if (twod && (s_t < 0)) s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */
821:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
822:       }
823:     }

825:     for (i=1; i<=s_y; i++) {
826:       if (n6 >= 0) { /* left above */
827:         x_t = lx[n6 % m];
828:         y_t = ly[(n6 % (m*n))/m];
829:         z_t = lz[n6 / (m*n)];
830:         s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
831:         if (twod && (s_t < 0)) s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - x_t*y_t; /* 2D case */
832:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
833:       }
834:       if (n7 >= 0) { /* directly above */
835:         x_t = x;
836:         y_t = ly[(n7 % (m*n))/m];
837:         z_t = lz[n7 / (m*n)];
838:         s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
839:         if (twod && (s_t < 0)) s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */
840:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
841:       }
842:       if (n8 >= 0) { /* right above */
843:         x_t = lx[n8 % m];
844:         y_t = ly[(n8 % (m*n))/m];
845:         z_t = lz[n8 / (m*n)];
846:         s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
847:         if (twod && (s_t < 0)) s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */
848:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
849:       }
850:     }
851:   }

853:   /* Middle Level */
854:   for (k=0; k<z; k++) {
855:     for (i=1; i<=s_y; i++) {
856:       if (n9 >= 0) { /* left below */
857:         x_t = lx[n9 % m];
858:         y_t = ly[(n9 % (m*n))/m];
859:         /* z_t = z; */
860:         s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
861:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
862:       }
863:       if (n10 >= 0) { /* directly below */
864:         x_t = x;
865:         y_t = ly[(n10 % (m*n))/m];
866:         /* z_t = z; */
867:         s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
868:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
869:       }  else if (by == DM_BOUNDARY_MIRROR) {
870:         for (j=0; j<x; j++) idx[nn++] = 0;
871:       }
872:       if (n11 >= 0) { /* right below */
873:         x_t = lx[n11 % m];
874:         y_t = ly[(n11 % (m*n))/m];
875:         /* z_t = z; */
876:         s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
877:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
878:       }
879:     }

881:     for (i=0; i<y; i++) {
882:       if (n12 >= 0) { /* directly left */
883:         x_t = lx[n12 % m];
884:         y_t = y;
885:         /* z_t = z; */
886:         s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
887:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
888:       }  else if (bx == DM_BOUNDARY_MIRROR) {
889:         for (j=0; j<s_x; j++) idx[nn++] = 0;
890:       }

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

896:       if (n14 >= 0) { /* directly right */
897:         x_t = lx[n14 % m];
898:         y_t = y;
899:         /* z_t = z; */
900:         s_t = bases[n14] + i*x_t + k*x_t*y_t;
901:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
902:       } else if (bx == DM_BOUNDARY_MIRROR) {
903:         for (j=0; j<s_x; j++) idx[nn++] = 0;
904:       }
905:     }

907:     for (i=1; i<=s_y; i++) {
908:       if (n15 >= 0) { /* left above */
909:         x_t = lx[n15 % m];
910:         y_t = ly[(n15 % (m*n))/m];
911:         /* z_t = z; */
912:         s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
913:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
914:       }
915:       if (n16 >= 0) { /* directly above */
916:         x_t = x;
917:         y_t = ly[(n16 % (m*n))/m];
918:         /* z_t = z; */
919:         s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
920:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
921:       } else if (by == DM_BOUNDARY_MIRROR) {
922:         for (j=0; j<x; j++) idx[nn++] = 0;
923:       }
924:       if (n17 >= 0) { /* right above */
925:         x_t = lx[n17 % m];
926:         y_t = ly[(n17 % (m*n))/m];
927:         /* z_t = z; */
928:         s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
929:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
930:       }
931:     }
932:   }

934:   /* Upper Level */
935:   for (k=0; k<s_z; k++) {
936:     for (i=1; i<=s_y; i++) {
937:       if (n18 >= 0) { /* left below */
938:         x_t = lx[n18 % m];
939:         y_t = ly[(n18 % (m*n))/m];
940:         /* z_t = lz[n18 / (m*n)]; */
941:         s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
942:         if (twod && (s_t >= M*N*P)) s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t; /* 2d case */
943:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
944:       }
945:       if (n19 >= 0) { /* directly below */
946:         x_t = x;
947:         y_t = ly[(n19 % (m*n))/m];
948:         /* z_t = lz[n19 / (m*n)]; */
949:         s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
950:         if (twod && (s_t >= M*N*P)) s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */
951:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
952:       }
953:       if (n20 >= 0) { /* right below */
954:         x_t = lx[n20 % m];
955:         y_t = ly[(n20 % (m*n))/m];
956:         /* z_t = lz[n20 / (m*n)]; */
957:         s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
958:         if (twod && (s_t >= M*N*P)) s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */
959:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
960:       }
961:     }

963:     for (i=0; i<y; i++) {
964:       if (n21 >= 0) { /* directly left */
965:         x_t = lx[n21 % m];
966:         y_t = y;
967:         /* z_t = lz[n21 / (m*n)]; */
968:         s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
969:         if (twod && (s_t >= M*N*P)) s_t = bases[n21] + (i+1)*x_t - s_x;  /* 2d case */
970:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
971:       }

973:       if (n22 >= 0) { /* middle */
974:         x_t = x;
975:         y_t = y;
976:         /* z_t = lz[n22 / (m*n)]; */
977:         s_t = bases[n22] + i*x_t + k*x_t*y_t;
978:         if (twod && (s_t >= M*N*P)) s_t = bases[n22] + i*x_t; /* 2d case */
979:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
980:       } else if (bz == DM_BOUNDARY_MIRROR) {
981:         for (j=0; j<x; j++) idx[nn++] = 0;
982:       }

984:       if (n23 >= 0) { /* directly right */
985:         x_t = lx[n23 % m];
986:         y_t = y;
987:         /* z_t = lz[n23 / (m*n)]; */
988:         s_t = bases[n23] + i*x_t + k*x_t*y_t;
989:         if (twod && (s_t >= M*N*P)) s_t = bases[n23] + i*x_t; /* 2d case */
990:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
991:       }
992:     }

994:     for (i=1; i<=s_y; i++) {
995:       if (n24 >= 0) { /* left above */
996:         x_t = lx[n24 % m];
997:         y_t = ly[(n24 % (m*n))/m];
998:         /* z_t = lz[n24 / (m*n)]; */
999:         s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
1000:         if (twod && (s_t >= M*N*P)) s_t = bases[n24] + i*x_t - s_x; /* 2d case */
1001:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1002:       }
1003:       if (n25 >= 0) { /* directly above */
1004:         x_t = x;
1005:         y_t = ly[(n25 % (m*n))/m];
1006:         /* z_t = lz[n25 / (m*n)]; */
1007:         s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1008:         if (twod && (s_t >= M*N*P)) s_t = bases[n25] + (i-1)*x_t; /* 2d case */
1009:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1010:       }
1011:       if (n26 >= 0) { /* right above */
1012:         x_t = lx[n26 % m];
1013:         y_t = ly[(n26 % (m*n))/m];
1014:         /* z_t = lz[n26 / (m*n)]; */
1015:         s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1016:         if (twod && (s_t >= M*N*P)) s_t = bases[n26] + (i-1)*x_t; /* 2d case */
1017:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1018:       }
1019:     }
1020:   }

1022:   ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);
1023:   VecScatterCreate(global,from,local,to,&gtol);
1024:   PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
1025:   ISDestroy(&to);
1026:   ISDestroy(&from);

1028:   if (stencil_type == DMDA_STENCIL_STAR) {
1029:     n0  = sn0;  n1  = sn1;  n2  = sn2;  n3  = sn3;  n5  = sn5;  n6  = sn6; n7 = sn7;
1030:     n8  = sn8;  n9  = sn9;  n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18;
1031:     n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25;
1032:     n26 = sn26;
1033:   }

1035:   if (((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz))) {
1036:     /*
1037:         Recompute the local to global mappings, this time keeping the
1038:       information about the cross corner processor numbers.
1039:     */
1040:     nn = 0;
1041:     /* Bottom Level */
1042:     for (k=0; k<s_z; k++) {
1043:       for (i=1; i<=s_y; i++) {
1044:         if (n0 >= 0) { /* left below */
1045:           x_t = lx[n0 % m];
1046:           y_t = ly[(n0 % (m*n))/m];
1047:           z_t = lz[n0 / (m*n)];
1048:           s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x - (s_z-k-1)*x_t*y_t;
1049:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1050:         } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) {
1051:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1052:         }
1053:         if (n1 >= 0) { /* directly below */
1054:           x_t = x;
1055:           y_t = ly[(n1 % (m*n))/m];
1056:           z_t = lz[n1 / (m*n)];
1057:           s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1058:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1059:         } else if (Ys-ys < 0 && Zs-zs < 0) {
1060:           for (j=0; j<x; j++) idx[nn++] = -1;
1061:         }
1062:         if (n2 >= 0) { /* right below */
1063:           x_t = lx[n2 % m];
1064:           y_t = ly[(n2 % (m*n))/m];
1065:           z_t = lz[n2 / (m*n)];
1066:           s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1067:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1068:         } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) {
1069:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1070:         }
1071:       }

1073:       for (i=0; i<y; i++) {
1074:         if (n3 >= 0) { /* directly left */
1075:           x_t = lx[n3 % m];
1076:           y_t = y;
1077:           z_t = lz[n3 / (m*n)];
1078:           s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1079:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1080:         } else if (Xs-xs < 0 && Zs-zs < 0) {
1081:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1082:         }

1084:         if (n4 >= 0) { /* middle */
1085:           x_t = x;
1086:           y_t = y;
1087:           z_t = lz[n4 / (m*n)];
1088:           s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1089:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1090:         } else if (Zs-zs < 0) {
1091:           if (bz == DM_BOUNDARY_MIRROR) {
1092:             for (j=0; j<x; j++) idx[nn++] = 0;
1093:           } else {
1094:             for (j=0; j<x; j++) idx[nn++] = -1;
1095:           }
1096:         }

1098:         if (n5 >= 0) { /* directly right */
1099:           x_t = lx[n5 % m];
1100:           y_t = y;
1101:           z_t = lz[n5 / (m*n)];
1102:           s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1103:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1104:         } else if (xe-Xe < 0 && Zs-zs < 0) {
1105:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1106:         }
1107:       }

1109:       for (i=1; i<=s_y; i++) {
1110:         if (n6 >= 0) { /* left above */
1111:           x_t = lx[n6 % m];
1112:           y_t = ly[(n6 % (m*n))/m];
1113:           z_t = lz[n6 / (m*n)];
1114:           s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1115:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1116:         } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) {
1117:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1118:         }
1119:         if (n7 >= 0) { /* directly above */
1120:           x_t = x;
1121:           y_t = ly[(n7 % (m*n))/m];
1122:           z_t = lz[n7 / (m*n)];
1123:           s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1124:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1125:         } else if (ye-Ye < 0 && Zs-zs < 0) {
1126:           for (j=0; j<x; j++) idx[nn++] = -1;
1127:         }
1128:         if (n8 >= 0) { /* right above */
1129:           x_t = lx[n8 % m];
1130:           y_t = ly[(n8 % (m*n))/m];
1131:           z_t = lz[n8 / (m*n)];
1132:           s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1133:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1134:         } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) {
1135:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1136:         }
1137:       }
1138:     }

1140:     /* Middle Level */
1141:     for (k=0; k<z; k++) {
1142:       for (i=1; i<=s_y; i++) {
1143:         if (n9 >= 0) { /* left below */
1144:           x_t = lx[n9 % m];
1145:           y_t = ly[(n9 % (m*n))/m];
1146:           /* z_t = z; */
1147:           s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1148:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1149:         } else if (Xs-xs < 0 && Ys-ys < 0) {
1150:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1151:         }
1152:         if (n10 >= 0) { /* directly below */
1153:           x_t = x;
1154:           y_t = ly[(n10 % (m*n))/m];
1155:           /* z_t = z; */
1156:           s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1157:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1158:         } else if (Ys-ys < 0) {
1159:           if (by == DM_BOUNDARY_MIRROR) {
1160:             for (j=0; j<x; j++) idx[nn++] = -1;
1161:           } else {
1162:             for (j=0; j<x; j++) idx[nn++] = -1;
1163:           }
1164:         }
1165:         if (n11 >= 0) { /* right below */
1166:           x_t = lx[n11 % m];
1167:           y_t = ly[(n11 % (m*n))/m];
1168:           /* z_t = z; */
1169:           s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1170:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1171:         } else if (xe-Xe < 0 && Ys-ys < 0) {
1172:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1173:         }
1174:       }

1176:       for (i=0; i<y; i++) {
1177:         if (n12 >= 0) { /* directly left */
1178:           x_t = lx[n12 % m];
1179:           y_t = y;
1180:           /* z_t = z; */
1181:           s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
1182:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1183:         } else if (Xs-xs < 0) {
1184:           if (bx == DM_BOUNDARY_MIRROR) {
1185:             for (j=0; j<s_x; j++) idx[nn++] = 0;
1186:           } else {
1187:             for (j=0; j<s_x; j++) idx[nn++] = -1;
1188:           }
1189:         }

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

1195:         if (n14 >= 0) { /* directly right */
1196:           x_t = lx[n14 % m];
1197:           y_t = y;
1198:           /* z_t = z; */
1199:           s_t = bases[n14] + i*x_t + k*x_t*y_t;
1200:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1201:         } else if (xe-Xe < 0) {
1202:           if (bx == DM_BOUNDARY_MIRROR) {
1203:             for (j=0; j<s_x; j++) idx[nn++] = 0;
1204:           } else {
1205:             for (j=0; j<s_x; j++) idx[nn++] = -1;
1206:           }
1207:         }
1208:       }

1210:       for (i=1; i<=s_y; i++) {
1211:         if (n15 >= 0) { /* left above */
1212:           x_t = lx[n15 % m];
1213:           y_t = ly[(n15 % (m*n))/m];
1214:           /* z_t = z; */
1215:           s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
1216:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1217:         } else if (Xs-xs < 0 && ye-Ye < 0) {
1218:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1219:         }
1220:         if (n16 >= 0) { /* directly above */
1221:           x_t = x;
1222:           y_t = ly[(n16 % (m*n))/m];
1223:           /* z_t = z; */
1224:           s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
1225:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1226:         } else if (ye-Ye < 0) {
1227:           if (by == DM_BOUNDARY_MIRROR) {
1228:             for (j=0; j<x; j++) idx[nn++] = 0;
1229:           } else {
1230:             for (j=0; j<x; j++) idx[nn++] = -1;
1231:           }
1232:         }
1233:         if (n17 >= 0) { /* right above */
1234:           x_t = lx[n17 % m];
1235:           y_t = ly[(n17 % (m*n))/m];
1236:           /* z_t = z; */
1237:           s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
1238:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1239:         } else if (xe-Xe < 0 && ye-Ye < 0) {
1240:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1241:         }
1242:       }
1243:     }

1245:     /* Upper Level */
1246:     for (k=0; k<s_z; k++) {
1247:       for (i=1; i<=s_y; i++) {
1248:         if (n18 >= 0) { /* left below */
1249:           x_t = lx[n18 % m];
1250:           y_t = ly[(n18 % (m*n))/m];
1251:           /* z_t = lz[n18 / (m*n)]; */
1252:           s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1253:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1254:         } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) {
1255:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1256:         }
1257:         if (n19 >= 0) { /* directly below */
1258:           x_t = x;
1259:           y_t = ly[(n19 % (m*n))/m];
1260:           /* z_t = lz[n19 / (m*n)]; */
1261:           s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1262:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1263:         } else if (Ys-ys < 0 && ze-Ze < 0) {
1264:           for (j=0; j<x; j++) idx[nn++] = -1;
1265:         }
1266:         if (n20 >= 0) { /* right below */
1267:           x_t = lx[n20 % m];
1268:           y_t = ly[(n20 % (m*n))/m];
1269:           /* z_t = lz[n20 / (m*n)]; */
1270:           s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1271:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1272:         } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) {
1273:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1274:         }
1275:       }

1277:       for (i=0; i<y; i++) {
1278:         if (n21 >= 0) { /* directly left */
1279:           x_t = lx[n21 % m];
1280:           y_t = y;
1281:           /* z_t = lz[n21 / (m*n)]; */
1282:           s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
1283:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1284:         } else if (Xs-xs < 0 && ze-Ze < 0) {
1285:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1286:         }

1288:         if (n22 >= 0) { /* middle */
1289:           x_t = x;
1290:           y_t = y;
1291:           /* z_t = lz[n22 / (m*n)]; */
1292:           s_t = bases[n22] + i*x_t + k*x_t*y_t;
1293:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1294:         } else if (ze-Ze < 0) {
1295:           if (bz == DM_BOUNDARY_MIRROR) {
1296:             for (j=0; j<x; j++) idx[nn++] = 0;
1297:           } else {
1298:             for (j=0; j<x; j++) idx[nn++] = -1;
1299:           }
1300:         }

1302:         if (n23 >= 0) { /* directly right */
1303:           x_t = lx[n23 % m];
1304:           y_t = y;
1305:           /* z_t = lz[n23 / (m*n)]; */
1306:           s_t = bases[n23] + i*x_t + k*x_t*y_t;
1307:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1308:         } else if (xe-Xe < 0 && ze-Ze < 0) {
1309:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1310:         }
1311:       }

1313:       for (i=1; i<=s_y; i++) {
1314:         if (n24 >= 0) { /* left above */
1315:           x_t = lx[n24 % m];
1316:           y_t = ly[(n24 % (m*n))/m];
1317:           /* z_t = lz[n24 / (m*n)]; */
1318:           s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
1319:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1320:         } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) {
1321:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1322:         }
1323:         if (n25 >= 0) { /* directly above */
1324:           x_t = x;
1325:           y_t = ly[(n25 % (m*n))/m];
1326:           /* z_t = lz[n25 / (m*n)]; */
1327:           s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1328:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1329:         } else if (ye-Ye < 0 && ze-Ze < 0) {
1330:           for (j=0; j<x; j++) idx[nn++] = -1;
1331:         }
1332:         if (n26 >= 0) { /* right above */
1333:           x_t = lx[n26 % m];
1334:           y_t = ly[(n26 % (m*n))/m];
1335:           /* z_t = lz[n26 / (m*n)]; */
1336:           s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1337:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1338:         } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) {
1339:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1340:         }
1341:       }
1342:     }
1343:   }
1344:   /*
1345:      Set the local to global ordering in the global vector, this allows use
1346:      of VecSetValuesLocal().
1347:   */
1348:   ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);
1349:   PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);

1351:   PetscFree2(bases,ldims);
1352:   dd->m = m;  dd->n  = n;  dd->p  = p;
1353:   /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1354:   dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze;
1355:   dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze;

1357:   VecDestroy(&local);
1358:   VecDestroy(&global);

1360:   dd->gtol      = gtol;
1361:   dd->base      = base;
1362:   da->ops->view = DMView_DA_3d;
1363:   dd->ltol      = NULL;
1364:   dd->ao        = NULL;
1365:   return(0);
1366: }


1369: /*@C
1370:    DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1371:    regular array data that is distributed across some processors.

1373:    Collective on MPI_Comm

1375:    Input Parameters:
1376: +  comm - MPI communicator
1377: .  bx,by,bz - type of ghost nodes the array have.
1378:          Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC.
1379: .  stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX)
1380: .  M,N,P - global dimension in each direction of the array 
1381: .  m,n,p - corresponding number of processors in each dimension
1382:            (or PETSC_DECIDE to have calculated)
1383: .  dof - number of degrees of freedom per node
1384: .  s - stencil width
1385: -  lx, ly, lz - arrays containing the number of nodes in each cell along
1386:           the x, y, and z coordinates, or NULL. If non-null, these
1387:           must be of length as m,n,p and the corresponding
1388:           m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
1389:           the ly[] must N, sum of the lz[] must be P

1391:    Output Parameter:
1392: .  da - the resulting distributed array object

1394:    Options Database Key:
1395: +  -dm_view - Calls DMView() at the conclusion of DMDACreate3d()
1396: .  -da_grid_x <nx> - number of grid points in x direction, if M < 0
1397: .  -da_grid_y <ny> - number of grid points in y direction, if N < 0
1398: .  -da_grid_z <nz> - number of grid points in z direction, if P < 0
1399: .  -da_processors_x <MX> - number of processors in x direction
1400: .  -da_processors_y <MY> - number of processors in y direction
1401: .  -da_processors_z <MZ> - number of processors in z direction
1402: .  -da_refine_x <rx> - refinement ratio in x direction
1403: .  -da_refine_y <ry> - refinement ratio in y direction
1404: .  -da_refine_z <rz>- refinement ratio in z directio
1405: -  -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0

1407:    Level: beginner

1409:    Notes:
1410:    The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
1411:    standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
1412:    the standard 27-pt stencil.

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

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

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

1423: .keywords: distributed array, create, three-dimensional

1425: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
1426:           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(),
1427:           DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()

1429: @*/
1430: PetscErrorCode  DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M,
1431:                PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[],DM *da)
1432: {

1436:   DMDACreate(comm, da);
1437:   DMSetDimension(*da, 3);
1438:   DMDASetSizes(*da, M, N, P);
1439:   DMDASetNumProcs(*da, m, n, p);
1440:   DMDASetBoundaryType(*da, bx, by, bz);
1441:   DMDASetDof(*da, dof);
1442:   DMDASetStencilType(*da, stencil_type);
1443:   DMDASetStencilWidth(*da, s);
1444:   DMDASetOwnershipRanges(*da, lx, ly, lz);
1445:   return(0);
1446: }