Actual source code: da3.c


  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_LOAD_BALANCE) {
 36:       PetscInt      i,nmax = 0,nmin = PETSC_MAX_INT,navg = 0,*nz,nzlocal;
 37:       DMDALocalInfo info;
 38:       PetscMPIInt   size;
 39:       MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 40:       DMDAGetLocalInfo(da,&info);
 41:       nzlocal = info.xm*info.ym*info.zm;
 42:       PetscMalloc1(size,&nz);
 43:       MPI_Allgather(&nzlocal,1,MPIU_INT,nz,1,MPIU_INT,PetscObjectComm((PetscObject)da));
 44:       for (i=0; i<(PetscInt)size; i++) {
 45:         nmax = PetscMax(nmax,nz[i]);
 46:         nmin = PetscMin(nmin,nz[i]);
 47:         navg += nz[i];
 48:       }
 49:       PetscFree(nz);
 50:       navg = navg/size;
 51:       PetscViewerASCIIPrintf(viewer,"  Load Balance - Grid Points: Min %D  avg %D  max %D\n",nmin,navg,nmax);
 52:       return(0);
 53:     }
 54:     if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) {
 55:       DMDALocalInfo info;
 56:       DMDAGetLocalInfo(da,&info);
 57:       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);
 58:       PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D, Z range of indices: %D %D\n",
 59:                                                 info.xs,info.xs+info.xm,info.ys,info.ys+info.ym,info.zs,info.zs+info.zm);
 60: #if !defined(PETSC_USE_COMPLEX)
 61:       if (da->coordinates) {
 62:         PetscInt        last;
 63:         const PetscReal *coors;
 64:         VecGetArrayRead(da->coordinates,&coors);
 65:         VecGetLocalSize(da->coordinates,&last);
 66:         last = last - 3;
 67:         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]);
 68:         VecRestoreArrayRead(da->coordinates,&coors);
 69:       }
 70: #endif
 71:       PetscViewerFlush(viewer);
 72:       PetscViewerASCIIPopSynchronized(viewer);
 73:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
 74:       DMView_DA_GLVis(da,viewer);
 75:     } else {
 76:       DMView_DA_VTK(da,viewer);
 77:     }
 78:   } else if (isdraw) {
 79:     PetscDraw      draw;
 80:     PetscReal      ymin = -1.0,ymax = (PetscReal)dd->N;
 81:     PetscReal      xmin = -1.0,xmax = (PetscReal)((dd->M+2)*dd->P),x,y,ycoord,xcoord;
 82:     PetscInt       k,plane,base;
 83:     const PetscInt *idx;
 84:     char           node[10];
 85:     PetscBool      isnull;

 87:     PetscViewerDrawGetDraw(viewer,0,&draw);
 88:     PetscDrawIsNull(draw,&isnull);
 89:     if (isnull) return(0);

 91:     PetscDrawCheckResizedWindow(draw);
 92:     PetscDrawClear(draw);
 93:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);

 95:     PetscDrawCollectiveBegin(draw);
 96:     /* first processor draw all node lines */
 97:     if (!rank) {
 98:       for (k=0; k<dd->P; k++) {
 99:         ymin = 0.0; ymax = (PetscReal)(dd->N - 1);
100:         for (xmin=(PetscReal)(k*(dd->M+1)); xmin<(PetscReal)(dd->M+(k*(dd->M+1))); xmin++) {
101:           PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
102:         }
103:         xmin = (PetscReal)(k*(dd->M+1)); xmax = xmin + (PetscReal)(dd->M - 1);
104:         for (ymin=0; ymin<(PetscReal)dd->N; ymin++) {
105:           PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
106:         }
107:       }
108:     }
109:     PetscDrawCollectiveEnd(draw);
110:     PetscDrawFlush(draw);
111:     PetscDrawPause(draw);

113:     PetscDrawCollectiveBegin(draw);
114:     /*Go through and draw for each plane*/
115:     for (k=0; k<dd->P; k++) {
116:       if ((k >= dd->zs) && (k < dd->ze)) {
117:         /* draw my box */
118:         ymin = dd->ys;
119:         ymax = dd->ye - 1;
120:         xmin = dd->xs/dd->w    + (dd->M+1)*k;
121:         xmax =(dd->xe-1)/dd->w + (dd->M+1)*k;

123:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
124:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
125:         PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
126:         PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

128:         xmin = dd->xs/dd->w;
129:         xmax =(dd->xe-1)/dd->w;

131:         /* identify which processor owns the box */
132:         PetscSNPrintf(node,sizeof(node),"%d",(int)rank);
133:         PetscDrawString(draw,xmin+(dd->M+1)*k+.2,ymin+.3,PETSC_DRAW_RED,node);
134:         /* put in numbers*/
135:         base = (dd->base+(dd->xe-dd->xs)*(dd->ye-dd->ys)*(k-dd->zs))/dd->w;
136:         for (y=ymin; y<=ymax; y++) {
137:           for (x=xmin+(dd->M+1)*k; x<=xmax+(dd->M+1)*k; x++) {
138:             PetscSNPrintf(node,sizeof(node),"%d",(int)base++);
139:             PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
140:           }
141:         }

143:       }
144:     }
145:     PetscDrawCollectiveEnd(draw);
146:     PetscDrawFlush(draw);
147:     PetscDrawPause(draw);

149:     PetscDrawCollectiveBegin(draw);
150:     for (k=0-dd->s; k<dd->P+dd->s; k++) {
151:       /* Go through and draw for each plane */
152:       if ((k >= dd->Zs) && (k < dd->Ze)) {
153:         /* overlay ghost numbers, useful for error checking */
154:         base = (dd->Xe-dd->Xs)*(dd->Ye-dd->Ys)*(k-dd->Zs)/dd->w;
155:         ISLocalToGlobalMappingGetBlockIndices(da->ltogmap,&idx);
156:         plane=k;
157:         /* Keep z wrap around points on the drawing */
158:         if (k<0) plane=dd->P+k;
159:         if (k>=dd->P) plane=k-dd->P;
160:         ymin = dd->Ys; ymax = dd->Ye;
161:         xmin = (dd->M+1)*plane*dd->w;
162:         xmax = (dd->M+1)*plane*dd->w+dd->M*dd->w;
163:         for (y=ymin; y<ymax; y++) {
164:           for (x=xmin+dd->Xs; x<xmin+dd->Xe; x+=dd->w) {
165:             sprintf(node,"%d",(int)(idx[base]));
166:             ycoord = y;
167:             /*Keep y wrap around points on drawing */
168:             if (y<0) ycoord = dd->N+y;
169:             if (y>=dd->N) ycoord = y-dd->N;
170:             xcoord = x;   /* Keep x wrap points on drawing */
171:             if (x<xmin) xcoord = xmax - (xmin-x);
172:             if (x>=xmax) xcoord = xmin + (x-xmax);
173:             PetscDrawString(draw,xcoord/dd->w,ycoord,PETSC_DRAW_BLUE,node);
174:             base++;
175:           }
176:         }
177:         ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap,&idx);
178:       }
179:     }
180:     PetscDrawCollectiveEnd(draw);
181:     PetscDrawFlush(draw);
182:     PetscDrawPause(draw);
183:     PetscDrawSave(draw);
184:   } else if (isglvis) {
185:     DMView_DA_GLVis(da,viewer);
186:   } else if (isbinary) {
187:     DMView_DA_Binary(da,viewer);
188: #if defined(PETSC_HAVE_MATLAB_ENGINE)
189:   } else if (ismatlab) {
190:     DMView_DA_Matlab(da,viewer);
191: #endif
192:   }
193:   return(0);
194: }

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

232:   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");
233:   PetscObjectGetComm((PetscObject) da, &comm);
234: #if !defined(PETSC_USE_64BIT_INDICES)
235:   if (((PetscInt64) M)*((PetscInt64) N)*((PetscInt64) P)*((PetscInt64) dof) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ4(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D by %D (dof) is too large for 32 bit indices",M,N,P,dof);
236: #endif

238:   MPI_Comm_size(comm,&size);
239:   MPI_Comm_rank(comm,&rank);

241:   if (m != PETSC_DECIDE) {
242:     if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
243:     else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
244:   }
245:   if (n != PETSC_DECIDE) {
246:     if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
247:     else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
248:   }
249:   if (p != PETSC_DECIDE) {
250:     if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
251:     else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
252:   }
253:   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);

255:   /* Partition the array among the processors */
256:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
257:     m = size/(n*p);
258:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
259:     n = size/(m*p);
260:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
261:     p = size/(m*n);
262:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
263:     /* try for squarish distribution */
264:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
265:     if (!m) m = 1;
266:     while (m > 0) {
267:       n = size/(m*p);
268:       if (m*n*p == size) break;
269:       m--;
270:     }
271:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
272:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
273:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
274:     /* try for squarish distribution */
275:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
276:     if (!m) m = 1;
277:     while (m > 0) {
278:       p = size/(m*n);
279:       if (m*n*p == size) break;
280:       m--;
281:     }
282:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
283:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
284:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
285:     /* try for squarish distribution */
286:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
287:     if (!n) n = 1;
288:     while (n > 0) {
289:       p = size/(m*n);
290:       if (m*n*p == size) break;
291:       n--;
292:     }
293:     if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
294:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
295:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
296:     /* try for squarish distribution */
297:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
298:     if (!n) n = 1;
299:     while (n > 0) {
300:       pm = size/n;
301:       if (n*pm == size) break;
302:       n--;
303:     }
304:     if (!n) n = 1;
305:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
306:     if (!m) m = 1;
307:     while (m > 0) {
308:       p = size/(m*n);
309:       if (m*n*p == size) break;
310:       m--;
311:     }
312:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
313:   } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

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

320:   /*
321:      Determine locally owned region
322:      [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
323:   */

325:   if (!lx) {
326:     PetscMalloc1(m, &dd->lx);
327:     lx   = dd->lx;
328:     for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m));
329:   }
330:   x  = lx[rank % m];
331:   xs = 0;
332:   for (i=0; i<(rank%m); i++) xs += lx[i];
333:   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);

335:   if (!ly) {
336:     PetscMalloc1(n, &dd->ly);
337:     ly   = dd->ly;
338:     for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n));
339:   }
340:   y = ly[(rank % (m*n))/m];
341:   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);

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

346:   if (!lz) {
347:     PetscMalloc1(p, &dd->lz);
348:     lz = dd->lz;
349:     for (i=0; i<p; i++) lz[i] = P/p + ((P % p) > (i % p));
350:   }
351:   z = lz[rank/(m*n)];

353:   /* note this is different than x- and y-, as we will handle as an important special
354:    case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z.  This is to deal with 2D problems
355:    in a 3D code.  Additional code for this case is noted with "2d case" comments */
356:   twod = PETSC_FALSE;
357:   if (P == 1) twod = PETSC_TRUE;
358:   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);
359:   zs = 0;
360:   for (i=0; i<(rank/(m*n)); i++) zs += lz[i];
361:   ye = ys + y;
362:   xe = xs + x;
363:   ze = zs + z;

365:   /* determine ghost region (Xs) and region scattered into (IXs)  */
366:   if (xs-s > 0) {
367:     Xs = xs - s; IXs = xs - s;
368:   } else {
369:     if (bx) Xs = xs - s;
370:     else Xs = 0;
371:     IXs = 0;
372:   }
373:   if (xe+s <= M) {
374:     Xe = xe + s; IXe = xe + s;
375:   } else {
376:     if (bx) {
377:       Xs = xs - s; Xe = xe + s;
378:     } else Xe = M;
379:     IXe = M;
380:   }

382:   if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
383:     IXs = xs - s;
384:     IXe = xe + s;
385:     Xs  = xs - s;
386:     Xe  = xe + s;
387:   }

389:   if (ys-s > 0) {
390:     Ys = ys - s; IYs = ys - s;
391:   } else {
392:     if (by) Ys = ys - s;
393:     else Ys = 0;
394:     IYs = 0;
395:   }
396:   if (ye+s <= N) {
397:     Ye = ye + s; IYe = ye + s;
398:   } else {
399:     if (by) Ye = ye + s;
400:     else Ye = N;
401:     IYe = N;
402:   }

404:   if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
405:     IYs = ys - s;
406:     IYe = ye + s;
407:     Ys  = ys - s;
408:     Ye  = ye + s;
409:   }

411:   if (zs-s > 0) {
412:     Zs = zs - s; IZs = zs - s;
413:   } else {
414:     if (bz) Zs = zs - s;
415:     else Zs = 0;
416:     IZs = 0;
417:   }
418:   if (ze+s <= P) {
419:     Ze = ze + s; IZe = ze + s;
420:   } else {
421:     if (bz) Ze = ze + s;
422:     else Ze = P;
423:     IZe = P;
424:   }

426:   if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
427:     IZs = zs - s;
428:     IZe = ze + s;
429:     Zs  = zs - s;
430:     Ze  = ze + s;
431:   }

433:   /* Resize all X parameters to reflect w */
434:   s_x = s;
435:   s_y = s;
436:   s_z = s;

438:   /* determine starting point of each processor */
439:   nn       = x*y*z;
440:   PetscMalloc2(size+1,&bases,size,&ldims);
441:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
442:   bases[0] = 0;
443:   for (i=1; i<=size; i++) bases[i] = ldims[i-1];
444:   for (i=1; i<=size; i++) bases[i] += bases[i-1];
445:   base = bases[rank]*dof;

447:   /* allocate the base parallel and sequential vectors */
448:   dd->Nlocal = x*y*z*dof;
449:   VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);
450:   dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs)*dof;
451:   VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);

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

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

507:   /* determine who lies on each side of use stored in    n24 n25 n26
508:                                                          n21 n22 n23
509:                                                          n18 n19 n20

511:                                                          n15 n16 n17
512:                                                          n12     n14
513:                                                          n9  n10 n11

515:                                                          n6  n7  n8
516:                                                          n3  n4  n5
517:                                                          n0  n1  n2
518:   */

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

532:   n9  = rank - m - 1;
533:   n10 = rank - m;
534:   n11 = rank - m + 1;
535:   n12 = rank - 1;
536:   n14 = rank + 1;
537:   n15 = rank + m - 1;
538:   n16 = rank + m;
539:   n17 = rank + m + 1;

541:   n18 = rank + m*n - m - 1;
542:   n19 = rank + m*n - m;
543:   n20 = rank + m*n - m + 1;
544:   n21 = rank + m*n - 1;
545:   n22 = rank + m*n;
546:   n23 = rank + m*n + 1;
547:   n24 = rank + m*n + m - 1;
548:   n25 = rank + m*n + m;
549:   n26 = rank + m*n + m + 1;

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

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

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

577:   if (ys==0) { /* First assume not corner or edge */
578:     n0  = rank + m * (n-1) -1 - (m*n);
579:     n1  = rank + m * (n-1)    - (m*n);
580:     n2  = rank + m * (n-1) +1 - (m*n);
581:     n9  = rank + m * (n-1) -1;
582:     n10 = rank + m * (n-1);
583:     n11 = rank + m * (n-1) +1;
584:     n18 = rank + m * (n-1) -1 + (m*n);
585:     n19 = rank + m * (n-1)    + (m*n);
586:     n20 = rank + m * (n-1) +1 + (m*n);
587:   }

589:   if (ye == N) { /* First assume not corner or edge */
590:     n6  = rank - m * (n-1) -1 - (m*n);
591:     n7  = rank - m * (n-1)    - (m*n);
592:     n8  = rank - m * (n-1) +1 - (m*n);
593:     n15 = rank - m * (n-1) -1;
594:     n16 = rank - m * (n-1);
595:     n17 = rank - m * (n-1) +1;
596:     n24 = rank - m * (n-1) -1 + (m*n);
597:     n25 = rank - m * (n-1)    + (m*n);
598:     n26 = rank - m * (n-1) +1 + (m*n);
599:   }

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

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

625:   if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
626:     n0 = size - m*n + rank + m-1 - m;
627:     n3 = size - m*n + rank + m-1;
628:     n6 = size - m*n + rank + m-1 + m;
629:   }

631:   if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
632:     n18 = m*n - (size - rank) + m-1 - m;
633:     n21 = m*n - (size - rank) + m-1;
634:     n24 = m*n - (size - rank) + m-1 + m;
635:   }

637:   if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
638:     n0  = rank + m*n -1 - m*n;
639:     n9  = rank + m*n -1;
640:     n18 = rank + m*n -1 + m*n;
641:   }

643:   if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
644:     n6  = rank - m*(n-1) + m-1 - m*n;
645:     n15 = rank - m*(n-1) + m-1;
646:     n24 = rank - m*(n-1) + m-1 + m*n;
647:   }

649:   if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */
650:     n2 = size - (m*n-rank) - (m-1) - m;
651:     n5 = size - (m*n-rank) - (m-1);
652:     n8 = size - (m*n-rank) - (m-1) + m;
653:   }

655:   if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */
656:     n20 = m*n - (size - rank) - (m-1) - m;
657:     n23 = m*n - (size - rank) - (m-1);
658:     n26 = m*n - (size - rank) - (m-1) + m;
659:   }

661:   if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */
662:     n2  = rank + m*(n-1) - (m-1) - m*n;
663:     n11 = rank + m*(n-1) - (m-1);
664:     n20 = rank + m*(n-1) - (m-1) + m*n;
665:   }

667:   if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */
668:     n8  = rank - m*n +1 - m*n;
669:     n17 = rank - m*n +1;
670:     n26 = rank - m*n +1 + m*n;
671:   }

673:   if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
674:     n0 = size - m + rank -1;
675:     n1 = size - m + rank;
676:     n2 = size - m + rank +1;
677:   }

679:   if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
680:     n18 = m*n - (size - rank) + m*(n-1) -1;
681:     n19 = m*n - (size - rank) + m*(n-1);
682:     n20 = m*n - (size - rank) + m*(n-1) +1;
683:   }

685:   if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
686:     n6 = size - (m*n-rank) - m * (n-1) -1;
687:     n7 = size - (m*n-rank) - m * (n-1);
688:     n8 = size - (m*n-rank) - m * (n-1) +1;
689:   }

691:   if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
692:     n24 = rank - (size-m) -1;
693:     n25 = rank - (size-m);
694:     n26 = rank - (size-m) +1;
695:   }

697:   /* Check for Corners */
698:   if ((xs==0) && (ys==0) && (zs==0)) n0  = size -1;
699:   if ((xs==0) && (ys==0) && (ze==P)) n18 = m*n-1;
700:   if ((xs==0) && (ye==N) && (zs==0)) n6  = (size-1)-m*(n-1);
701:   if ((xs==0) && (ye==N) && (ze==P)) n24 = m-1;
702:   if ((xe==M) && (ys==0) && (zs==0)) n2  = size-m;
703:   if ((xe==M) && (ys==0) && (ze==P)) n20 = m*n-m;
704:   if ((xe==M) && (ye==N) && (zs==0)) n8  = size-m*n;
705:   if ((xe==M) && (ye==N) && (ze==P)) n26 = 0;

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

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

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

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

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

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

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

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

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

799:     for (i=0; i<y; i++) {
800:       if (n3 >= 0) { /* directly left */
801:         x_t = lx[n3 % m];
802:         y_t = y;
803:         z_t = lz[n3 / (m*n)];
804:         s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
805:         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 */
806:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
807:       }

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

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

830:     for (i=1; i<=s_y; i++) {
831:       if (n6 >= 0) { /* left above */
832:         x_t = lx[n6 % m];
833:         y_t = ly[(n6 % (m*n))/m];
834:         z_t = lz[n6 / (m*n)];
835:         s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
836:         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 */
837:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
838:       }
839:       if (n7 >= 0) { /* directly above */
840:         x_t = x;
841:         y_t = ly[(n7 % (m*n))/m];
842:         z_t = lz[n7 / (m*n)];
843:         s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
844:         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 */
845:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
846:       }
847:       if (n8 >= 0) { /* right above */
848:         x_t = lx[n8 % m];
849:         y_t = ly[(n8 % (m*n))/m];
850:         z_t = lz[n8 / (m*n)];
851:         s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
852:         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 */
853:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
854:       }
855:     }
856:   }

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

886:     for (i=0; i<y; i++) {
887:       if (n12 >= 0) { /* directly left */
888:         x_t = lx[n12 % m];
889:         y_t = y;
890:         /* z_t = z; */
891:         s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
892:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
893:       }  else if (bx == DM_BOUNDARY_MIRROR) {
894:         for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k*x*y + i*x;
895:       }

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

901:       if (n14 >= 0) { /* directly right */
902:         x_t = lx[n14 % m];
903:         y_t = y;
904:         /* z_t = z; */
905:         s_t = bases[n14] + i*x_t + k*x_t*y_t;
906:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
907:       } else if (bx == DM_BOUNDARY_MIRROR) {
908:         for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k*x*y + i*x;
909:       }
910:     }

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

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

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

978:       if (n22 >= 0) { /* middle */
979:         x_t = x;
980:         y_t = y;
981:         /* z_t = lz[n22 / (m*n)]; */
982:         s_t = bases[n22] + i*x_t + k*x_t*y_t;
983:         if (twod && (s_t >= M*N*P)) s_t = bases[n22] + i*x_t; /* 2d case */
984:         for (j=0; j<x_t; j++) idx[nn++] = s_t++;
985:       } else if (bz == DM_BOUNDARY_MIRROR) {
986:         for (j=0; j<x; j++) idx[nn++] = bases[rank] + j + (z-k-1)*x*y + i*x;
987:       }

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

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

1027:   ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);
1028:   VecScatterCreate(global,from,local,to,&gtol);
1029:   PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
1030:   ISDestroy(&to);
1031:   ISDestroy(&from);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1362:   VecDestroy(&local);
1363:   VecDestroy(&global);

1365:   dd->gtol      = gtol;
1366:   dd->base      = base;
1367:   da->ops->view = DMView_DA_3d;
1368:   dd->ltol      = NULL;
1369:   dd->ao        = NULL;
1370:   return(0);
1371: }

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

1377:    Collective

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

1395:    Output Parameter:
1396: .  da - the resulting distributed array object

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

1411:    Level: beginner

1413:    Notes:
1414:    The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
1415:    standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
1416:    the standard 27-pt stencil.

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

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

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

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

1432: @*/
1433: PetscErrorCode  DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M,
1434:                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)
1435: {

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