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: {
 12:   PetscMPIInt    rank;
 13:   PetscBool      iascii,isdraw,isglvis,isbinary;
 14:   DM_DA          *dd = (DM_DA*)da->data;
 15: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 16:   PetscBool ismatlab;
 17: #endif

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

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

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

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

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

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

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

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

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

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

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

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

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

230:   PetscObjectGetComm((PetscObject) da, &comm);
231: #if !defined(PETSC_USE_64BIT_INDICES)
233: #endif

235:   MPI_Comm_size(comm,&size);
236:   MPI_Comm_rank(comm,&rank);

238:   if (m != PETSC_DECIDE) {
241:   }
242:   if (n != PETSC_DECIDE) {
245:   }
246:   if (p != PETSC_DECIDE) {
249:   }

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


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

322:   if (!lx) {
323:     PetscMalloc1(m, &dd->lx);
324:     lx   = dd->lx;
325:     for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m));
326:   }
327:   x  = lx[rank % m];
328:   xs = 0;
329:   for (i=0; i<(rank%m); i++) xs += lx[i];

332:   if (!ly) {
333:     PetscMalloc1(n, &dd->ly);
334:     ly   = dd->ly;
335:     for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n));
336:   }
337:   y = ly[(rank % (m*n))/m];

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

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

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

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

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

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

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

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

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

430:   /* Resize all X parameters to reflect w */
431:   s_x = s;
432:   s_y = s;
433:   s_z = s;

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

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

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

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

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

508:                                                          n15 n16 n17
509:                                                          n12     n14
510:                                                          n9  n10 n11

512:                                                          n6  n7  n8
513:                                                          n3  n4  n5
514:                                                          n0  n1  n2
515:   */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1359:   VecDestroy(&local);
1360:   VecDestroy(&global);

1362:   dd->gtol      = gtol;
1363:   dd->base      = base;
1364:   da->ops->view = DMView_DA_3d;
1365:   dd->ltol      = NULL;
1366:   dd->ao        = NULL;
1367:   return 0;
1368: }

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

1374:    Collective

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

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

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

1408:    Level: beginner

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

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

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

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

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

1429: @*/
1430: PetscErrorCode  DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M,
1431:                PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[],DM *da)
1432: {
1433:   DMDACreate(comm, da);
1434:   DMSetDimension(*da, 3);
1435:   DMDASetSizes(*da, M, N, P);
1436:   DMDASetNumProcs(*da, m, n, p);
1437:   DMDASetBoundaryType(*da, bx, by, bz);
1438:   DMDASetDof(*da, dof);
1439:   DMDASetStencilType(*da, stencil_type);
1440:   DMDASetStencilWidth(*da, s);
1441:   DMDASetOwnershipRanges(*da, lx, ly, lz);
1442:   return 0;
1443: }