Actual source code: da3.c

petsc-3.4.5 2014-06-29
  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>     /*I   "petscdmda.h"    I*/

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

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

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

 34:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
 35:     PetscViewerGetFormat(viewer, &format);
 36:     if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) {
 37:       DMDALocalInfo info;
 38:       DMDAGetLocalInfo(da,&info);
 39:       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);
 40:       PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D, Z range of indices: %D %D\n",
 41:                                                 info.xs,info.xs+info.xm,info.ys,info.ys+info.ym,info.zs,info.zs+info.zm);
 42: #if !defined(PETSC_USE_COMPLEX)
 43:       if (da->coordinates) {
 44:         PetscInt        last;
 45:         const PetscReal *coors;
 46:         VecGetArrayRead(da->coordinates,&coors);
 47:         VecGetLocalSize(da->coordinates,&last);
 48:         last = last - 3;
 49:         PetscViewerASCIISynchronizedPrintf(viewer,"Lower left corner %G %G %G : Upper right %G %G %G\n",coors[0],coors[1],coors[2],coors[last],coors[last+1],coors[last+2]);
 50:         VecRestoreArrayRead(da->coordinates,&coors);
 51:       }
 52: #endif
 53:       PetscViewerFlush(viewer);
 54:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
 55:     } else {
 56:       DMView_DA_VTK(da,viewer);
 57:     }
 58:   } else if (isdraw) {
 59:     PetscDraw draw;
 60:     PetscReal ymin = -1.0,ymax = (PetscReal)dd->N;
 61:     PetscReal xmin = -1.0,xmax = (PetscReal)((dd->M+2)*dd->P),x,y,ycoord,xcoord;
 62:     PetscInt  k,plane,base,*idx;
 63:     char      node[10];
 64:     PetscBool isnull;

 66:     PetscViewerDrawGetDraw(viewer,0,&draw);
 67:     PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
 68:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
 69:     PetscDrawSynchronizedClear(draw);

 71:     /* first processor draw all node lines */
 72:     if (!rank) {
 73:       for (k=0; k<dd->P; k++) {
 74:         ymin = 0.0; ymax = (PetscReal)(dd->N - 1);
 75:         for (xmin=(PetscReal)(k*(dd->M+1)); xmin<(PetscReal)(dd->M+(k*(dd->M+1))); xmin++) {
 76:           PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
 77:         }

 79:         xmin = (PetscReal)(k*(dd->M+1)); xmax = xmin + (PetscReal)(dd->M - 1);
 80:         for (ymin=0; ymin<(PetscReal)dd->N; ymin++) {
 81:           PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
 82:         }
 83:       }
 84:     }
 85:     PetscDrawSynchronizedFlush(draw);
 86:     PetscDrawPause(draw);

 88:     for (k=0; k<dd->P; k++) {  /*Go through and draw for each plane*/
 89:       if ((k >= dd->zs) && (k < dd->ze)) {
 90:         /* draw my box */
 91:         ymin = dd->ys;
 92:         ymax = dd->ye - 1;
 93:         xmin = dd->xs/dd->w    + (dd->M+1)*k;
 94:         xmax =(dd->xe-1)/dd->w + (dd->M+1)*k;

 96:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
 97:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
 98:         PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
 99:         PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

101:         xmin = dd->xs/dd->w;
102:         xmax =(dd->xe-1)/dd->w;

104:         /* put in numbers*/
105:         base = (dd->base+(dd->xe-dd->xs)*(dd->ye-dd->ys)*(k-dd->zs))/dd->w;

107:         /* Identify which processor owns the box */
108:         sprintf(node,"%d",rank);
109:         PetscDrawString(draw,xmin+(dd->M+1)*k+.2,ymin+.3,PETSC_DRAW_RED,node);

111:         for (y=ymin; y<=ymax; y++) {
112:           for (x=xmin+(dd->M+1)*k; x<=xmax+(dd->M+1)*k; x++) {
113:             sprintf(node,"%d",(int)base++);
114:             PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
115:           }
116:         }

118:       }
119:     }
120:     PetscDrawSynchronizedFlush(draw);
121:     PetscDrawPause(draw);

123:     for (k=0-dd->s; k<dd->P+dd->s; k++) {
124:       /* Go through and draw for each plane */
125:       if ((k >= dd->Zs) && (k < dd->Ze)) {

127:         /* overlay ghost numbers, useful for error checking */
128:         base = (dd->Xe-dd->Xs)*(dd->Ye-dd->Ys)*(k-dd->Zs); idx = dd->idx;
129:         plane=k;
130:         /* Keep z wrap around points on the dradrawg */
131:         if (k<0) plane=dd->P+k;
132:         if (k>=dd->P) plane=k-dd->P;
133:         ymin = dd->Ys; ymax = dd->Ye;
134:         xmin = (dd->M+1)*plane*dd->w;
135:         xmax = (dd->M+1)*plane*dd->w+dd->M*dd->w;
136:         for (y=ymin; y<ymax; y++) {
137:           for (x=xmin+dd->Xs; x<xmin+dd->Xe; x+=dd->w) {
138:             sprintf(node,"%d",(int)(idx[base]/dd->w));
139:             ycoord = y;
140:             /*Keep y wrap around points on drawing */
141:             if (y<0) ycoord = dd->N+y;

143:             if (y>=dd->N) ycoord = y-dd->N;
144:             xcoord = x;   /* Keep x wrap points on drawing */

146:             if (x<xmin) xcoord = xmax - (xmin-x);
147:             if (x>=xmax) xcoord = xmin + (x-xmax);
148:             PetscDrawString(draw,xcoord/dd->w,ycoord,PETSC_DRAW_BLUE,node);
149:             base+=dd->w;
150:           }
151:         }
152:       }
153:     }
154:     PetscDrawSynchronizedFlush(draw);
155:     PetscDrawPause(draw);
156:   } else if (isbinary) {
157:     DMView_DA_Binary(da,viewer);
158: #if defined(PETSC_HAVE_MATLAB_ENGINE)
159:   } else if (ismatlab) {
160:     DMView_DA_Matlab(da,viewer);
161: #endif
162:   }
163:   return(0);
164: }

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


206:   if (stencil_type == DMDA_STENCIL_BOX && (bx == DMDA_BOUNDARY_MIRROR || by == DMDA_BOUNDARY_MIRROR || bz == DMDA_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
207:   PetscObjectGetComm((PetscObject) da, &comm);
208: #if !defined(PETSC_USE_64BIT_INDICES)
209:   if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) P)*((Petsc64bitInt) dof) > (Petsc64bitInt) 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);
210: #endif

212:   MPI_Comm_size(comm,&size);
213:   MPI_Comm_rank(comm,&rank);

215:   if (m != PETSC_DECIDE) {
216:     if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
217:     else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
218:   }
219:   if (n != PETSC_DECIDE) {
220:     if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
221:     else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
222:   }
223:   if (p != PETSC_DECIDE) {
224:     if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
225:     else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
226:   }
227:   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);

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

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

294:   /*
295:      Determine locally owned region
296:      [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
297:   */

299:   if (!lx) {
300:     PetscMalloc(m*sizeof(PetscInt), &dd->lx);
301:     lx   = dd->lx;
302:     for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m));
303:   }
304:   x  = lx[rank % m];
305:   xs = 0;
306:   for (i=0; i<(rank%m); i++) xs += lx[i];
307:   if ((x < s) && ((m > 1) || (bx == DMDA_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);

309:   if (!ly) {
310:     PetscMalloc(n*sizeof(PetscInt), &dd->ly);
311:     ly   = dd->ly;
312:     for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n));
313:   }
314:   y = ly[(rank % (m*n))/m];
315:   if ((y < s) && ((n > 1) || (by == DMDA_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);

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

320:   if (!lz) {
321:     PetscMalloc(p*sizeof(PetscInt), &dd->lz);
322:     lz = dd->lz;
323:     for (i=0; i<p; i++) lz[i] = P/p + ((P % p) > (i % p));
324:   }
325:   z = lz[rank/(m*n)];

327:   /* note this is different than x- and y-, as we will handle as an important special
328:    case when p=P=1 and DMDA_BOUNDARY_PERIODIC and s > z.  This is to deal with 2D problems
329:    in a 3D code.  Additional code for this case is noted with "2d case" comments */
330:   twod = PETSC_FALSE;
331:   if (P == 1) twod = PETSC_TRUE;
332:   else if ((z < s) && ((p > 1) || (bz == DMDA_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);
333:   zs = 0;
334:   for (i=0; i<(rank/(m*n)); i++) zs += lz[i];
335:   ye = ys + y;
336:   xe = xs + x;
337:   ze = zs + z;

339:   /* determine ghost region (Xs) and region scattered into (IXs)  */
340:   if (xs-s > 0) {
341:     Xs = xs - s; IXs = xs - s;
342:   } else {
343:     if (bx) Xs = xs - s;
344:     else Xs = 0;
345:     IXs = 0;
346:   }
347:   if (xe+s <= M) {
348:     Xe = xe + s; IXe = xe + s;
349:   } else {
350:     if (bx) {
351:       Xs = xs - s; Xe = xe + s;
352:     } else Xe = M;
353:     IXe = M;
354:   }

356:   if (bx == DMDA_BOUNDARY_PERIODIC || bx == DMDA_BOUNDARY_MIRROR) {
357:     IXs = xs - s;
358:     IXe = xe + s;
359:     Xs  = xs - s;
360:     Xe  = xe + s;
361:   }

363:   if (ys-s > 0) {
364:     Ys = ys - s; IYs = ys - s;
365:   } else {
366:     if (by) Ys = ys - s;
367:     else Ys = 0;
368:     IYs = 0;
369:   }
370:   if (ye+s <= N) {
371:     Ye = ye + s; IYe = ye + s;
372:   } else {
373:     if (by) Ye = ye + s;
374:     else Ye = N;
375:     IYe = N;
376:   }

378:   if (by == DMDA_BOUNDARY_PERIODIC || by == DMDA_BOUNDARY_MIRROR) {
379:     IYs = ys - s;
380:     IYe = ye + s;
381:     Ys  = ys - s;
382:     Ye  = ye + s;
383:   }

385:   if (zs-s > 0) {
386:     Zs = zs - s; IZs = zs - s;
387:   } else {
388:     if (bz) Zs = zs - s;
389:     else Zs = 0;
390:     IZs = 0;
391:   }
392:   if (ze+s <= P) {
393:     Ze = ze + s; IZe = ze + s;
394:   } else {
395:     if (bz) Ze = ze + s;
396:     else Ze = P;
397:     IZe = P;
398:   }

400:   if (bz == DMDA_BOUNDARY_PERIODIC || bz == DMDA_BOUNDARY_MIRROR) {
401:     IZs = zs - s;
402:     IZe = ze + s;
403:     Zs  = zs - s;
404:     Ze  = ze + s;
405:   }

407:   /* Resize all X parameters to reflect w */
408:   s_x = s;
409:   s_y = s;
410:   s_z = s;

412:   /* determine starting point of each processor */
413:   nn       = x*y*z;
414:   PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);
415:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
416:   bases[0] = 0;
417:   for (i=1; i<=size; i++) bases[i] = ldims[i-1];
418:   for (i=1; i<=size; i++) bases[i] += bases[i-1];
419:   base = bases[rank]*dof;

421:   /* allocate the base parallel and sequential vectors */
422:   dd->Nlocal = x*y*z*dof;
423:   VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);
424:   dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs)*dof;
425:   VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);

427:   /* generate appropriate vector scatters */
428:   /* local to global inserts non-ghost point region into global */
429:   VecGetOwnershipRange(global,&start,&end);
430:   ISCreateStride(comm,x*y*z*dof,start,1,&to);

432:   PetscMalloc(x*y*z*sizeof(PetscInt),&idx);
433:   left   = xs - Xs; right = left + x;
434:   bottom = ys - Ys; top = bottom + y;
435:   down   = zs - Zs; up  = down + z;
436:   count  = 0;
437:   for (i=down; i<up; i++) {
438:     for (j=bottom; j<top; j++) {
439:       for (k=left; k<right; k++) {
440:         idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
441:       }
442:     }
443:   }

445:   ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);
446:   VecScatterCreate(local,from,global,to,&ltog);
447:   PetscLogObjectParent(da,ltog);
448:   ISDestroy(&from);
449:   ISDestroy(&to);

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:     count = (IXe-IXs)*(IYe-IYs)*(IZe-IZs);
455:     PetscMalloc(count*sizeof(PetscInt),&idx);

457:     left   = IXs - Xs; right = left + (IXe-IXs);
458:     bottom = IYs - Ys; top = bottom + (IYe-IYs);
459:     down   = IZs - Zs; up  = down + (IZe-IZs);
460:     count  = 0;
461:     for (i=down; i<up; i++) {
462:       for (j=bottom; j<top; j++) {
463:         for (k=left; k<right; k++) {
464:           idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
465:         }
466:       }
467:     }
468:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);

470:   } else {
471:     /* This is way ugly! We need to list the funny cross type region */
472:     count = ((ys-IYs) + (IYe-ye))*x*z + ((xs-IXs) + (IXe-xe))*y*z + ((zs-IZs) + (IZe-ze))*x*y + x*y*z;
473:     PetscMalloc(count*sizeof(PetscInt),&idx);

475:     left   = xs - Xs; right = left + x;
476:     bottom = ys - Ys; top = bottom + y;
477:     down   = zs - Zs;   up  = down + z;
478:     count  = 0;
479:     /* the bottom chunck */
480:     for (i=(IZs-Zs); i<down; i++) {
481:       for (j=bottom; j<top; j++) {
482:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
483:       }
484:     }
485:     /* the middle piece */
486:     for (i=down; i<up; i++) {
487:       /* front */
488:       for (j=(IYs-Ys); j<bottom; j++) {
489:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
490:       }
491:       /* middle */
492:       for (j=bottom; j<top; j++) {
493:         for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
494:       }
495:       /* back */
496:       for (j=top; j<top+IYe-ye; j++) {
497:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
498:       }
499:     }
500:     /* the top piece */
501:     for (i=up; i<up+IZe-ze; i++) {
502:       for (j=bottom; j<top; j++) {
503:         for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
504:       }
505:     }
506:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
507:   }

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

513:                                                          n15 n16 n17
514:                                                          n12     n14
515:                                                          n9  n10 n11

517:                                                          n6  n7  n8
518:                                                          n3  n4  n5
519:                                                          n0  n1  n2
520:   */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

729:   PetscMalloc(27*sizeof(PetscInt),&dd->neighbors);

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

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

769:   PetscMalloc((Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt),&idx);
770:   PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt));

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

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

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

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

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

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

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

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

904:       if (n14 >= 0) { /* directly right */
905:         x_t = lx[n14 % m];
906:         y_t = y;
907:         /* z_t = z; */
908:         s_t = bases[n14] + i*x_t + k*x_t*y_t;
909:         for (j=0; j<s_x; j++) idx[nn++] = s_t++;
910:       } else if (bx == DMDA_BOUNDARY_MIRROR) {
911:         for (j=0; j<s_x; j++) idx[nn++] = 0;
912:       }
913:     }

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

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

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

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

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

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

1030:   ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);
1031:   VecScatterCreate(global,from,local,to,&gtol);
1032:   PetscLogObjectParent(da,gtol);
1033:   ISDestroy(&to);
1034:   ISDestroy(&from);

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

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

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

1095:         if (n4 >= 0) { /* middle */
1096:           x_t = x;
1097:           y_t = y;
1098:           z_t = lz[n4 / (m*n)];
1099:           s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1100:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1101:         } else if (Zs-zs < 0) {
1102:           if (bz == DMDA_BOUNDARY_MIRROR) {
1103:             for (j=0; j<x; j++) idx[nn++] = 0;
1104:           } else {
1105:             for (j=0; j<x; j++) idx[nn++] = -1;
1106:           }
1107:         }

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

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

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

1187:       for (i=0; i<y; i++) {
1188:         if (n12 >= 0) { /* directly left */
1189:           x_t = lx[n12 % m];
1190:           y_t = y;
1191:           /* z_t = z; */
1192:           s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
1193:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1194:         } else if (Xs-xs < 0) {
1195:           if (bx == DMDA_BOUNDARY_MIRROR) {
1196:             for (j=0; j<s_x; j++) idx[nn++] = 0;
1197:           } else {
1198:             for (j=0; j<s_x; j++) idx[nn++] = -1;
1199:           }
1200:         }

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

1206:         if (n14 >= 0) { /* directly right */
1207:           x_t = lx[n14 % m];
1208:           y_t = y;
1209:           /* z_t = z; */
1210:           s_t = bases[n14] + i*x_t + k*x_t*y_t;
1211:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1212:         } else if (xe-Xe < 0) {
1213:           if (bx == DMDA_BOUNDARY_MIRROR) {
1214:             for (j=0; j<s_x; j++) idx[nn++] = 0;
1215:           } else {
1216:             for (j=0; j<s_x; j++) idx[nn++] = -1;
1217:           }
1218:         }
1219:       }

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

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

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

1299:         if (n22 >= 0) { /* middle */
1300:           x_t = x;
1301:           y_t = y;
1302:           /* z_t = lz[n22 / (m*n)]; */
1303:           s_t = bases[n22] + i*x_t + k*x_t*y_t;
1304:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1305:         } else if (ze-Ze < 0) {
1306:           if (bz == DMDA_BOUNDARY_MIRROR) {
1307:             for (j=0; j<x; j++) idx[nn++] = 0;
1308:           } else {
1309:             for (j=0; j<x; j++) idx[nn++] = -1;
1310:           }
1311:         }

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

1324:       for (i=1; i<=s_y; i++) {
1325:         if (n24 >= 0) { /* left above */
1326:           x_t = lx[n24 % m];
1327:           y_t = ly[(n24 % (m*n))/m];
1328:           /* z_t = lz[n24 / (m*n)]; */
1329:           s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
1330:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1331:         } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) {
1332:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1333:         }
1334:         if (n25 >= 0) { /* directly above */
1335:           x_t = x;
1336:           y_t = ly[(n25 % (m*n))/m];
1337:           /* z_t = lz[n25 / (m*n)]; */
1338:           s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1339:           for (j=0; j<x_t; j++) idx[nn++] = s_t++;
1340:         } else if (ye-Ye < 0 && ze-Ze < 0) {
1341:           for (j=0; j<x; j++) idx[nn++] = -1;
1342:         }
1343:         if (n26 >= 0) { /* right above */
1344:           x_t = lx[n26 % m];
1345:           y_t = ly[(n26 % (m*n))/m];
1346:           /* z_t = lz[n26 / (m*n)]; */
1347:           s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1348:           for (j=0; j<s_x; j++) idx[nn++] = s_t++;
1349:         } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) {
1350:           for (j=0; j<s_x; j++) idx[nn++] = -1;
1351:         }
1352:       }
1353:     }
1354:   }
1355:   /*
1356:      Set the local to global ordering in the global vector, this allows use
1357:      of VecSetValuesLocal().
1358:   */
1359:   ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,&ltogis);
1360:   PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);
1361:   PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));
1362:   ISGetIndices(ltogis, &idx_full);
1363:   PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));
1364:   ISRestoreIndices(ltogis, &idx_full);
1365:   ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);
1366:   PetscLogObjectParent(da,da->ltogmap);
1367:   ISDestroy(&ltogis);
1368:   ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);
1369:   PetscLogObjectParent(da,da->ltogmap);

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

1377:   VecDestroy(&local);
1378:   VecDestroy(&global);

1380:   dd->gtol      = gtol;
1381:   dd->ltog      = ltog;
1382:   dd->idx       = idx_cpy;
1383:   dd->Nl        = nn*dof;
1384:   dd->base      = base;
1385:   da->ops->view = DMView_DA_3d;
1386:   dd->ltol      = NULL;
1387:   dd->ao        = NULL;
1388:   return(0);
1389: }


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

1398:    Collective on MPI_Comm

1400:    Input Parameters:
1401: +  comm - MPI communicator
1402: .  bx,by,bz - type of ghost nodes the array have.
1403:          Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC.
1404: .  stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX)
1405: .  M,N,P - global dimension in each direction of the array (use -M, -N, and or -P to indicate that it may be set to a different value
1406:             from the command line with -da_grid_x <M> -da_grid_y <N> -da_grid_z <P>)
1407: .  m,n,p - corresponding number of processors in each dimension
1408:            (or PETSC_DECIDE to have calculated)
1409: .  dof - number of degrees of freedom per node
1410: .  s - stencil width
1411: -  lx, ly, lz - arrays containing the number of nodes in each cell along
1412:           the x, y, and z coordinates, or NULL. If non-null, these
1413:           must be of length as m,n,p and the corresponding
1414:           m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
1415:           the ly[] must N, sum of the lz[] must be P

1417:    Output Parameter:
1418: .  da - the resulting distributed array object

1420:    Options Database Key:
1421: +  -dm_view - Calls DMView() at the conclusion of DMDACreate3d()
1422: .  -da_grid_x <nx> - number of grid points in x direction, if M < 0
1423: .  -da_grid_y <ny> - number of grid points in y direction, if N < 0
1424: .  -da_grid_z <nz> - number of grid points in z direction, if P < 0
1425: .  -da_processors_x <MX> - number of processors in x direction
1426: .  -da_processors_y <MY> - number of processors in y direction
1427: .  -da_processors_z <MZ> - number of processors in z direction
1428: .  -da_refine_x <rx> - refinement ratio in x direction
1429: .  -da_refine_y <ry> - refinement ratio in y direction
1430: .  -da_refine_z <rz>- refinement ratio in z directio
1431: -  -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0

1433:    Level: beginner

1435:    Notes:
1436:    The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
1437:    standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
1438:    the standard 27-pt stencil.

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

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

1446: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
1447:           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(),
1448:           DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()

1450: @*/
1451: PetscErrorCode  DMDACreate3d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz,DMDAStencilType stencil_type,PetscInt M,
1452:                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)
1453: {

1457:   DMDACreate(comm, da);
1458:   DMDASetDim(*da, 3);
1459:   DMDASetSizes(*da, M, N, P);
1460:   DMDASetNumProcs(*da, m, n, p);
1461:   DMDASetBoundaryType(*da, bx, by, bz);
1462:   DMDASetDof(*da, dof);
1463:   DMDASetStencilType(*da, stencil_type);
1464:   DMDASetStencilWidth(*da, s);
1465:   DMDASetOwnershipRanges(*da, lx, ly, lz);
1466:   /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */
1467:   DMSetFromOptions(*da);
1468:   DMSetUp(*da);
1469:   DMViewFromOptions(*da,NULL,"-dm_view");
1470:   return(0);
1471: }