Actual source code: da3.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    233:   if (stencil_type == DMDA_STENCIL_BOX && (bx == DM_BOUNDARY_MIRROR || by == DM_BOUNDARY_MIRROR || bz == DM_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
    234:   PetscObjectGetComm((PetscObject) da, &comm);
    235: #if !defined(PETSC_USE_64BIT_INDICES)
    236:   if (((PetscInt64) M)*((PetscInt64) N)*((PetscInt64) P)*((PetscInt64) dof) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ4(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D by %D (dof) is too large for 32 bit indices",M,N,P,dof);
    237: #endif

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

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

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

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

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

    326:   if (!lx) {
    327:     PetscMalloc1(m, &dd->lx);
    328:     lx   = dd->lx;
    329:     for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m));
    330:   }
    331:   x  = lx[rank % m];
    332:   xs = 0;
    333:   for (i=0; i<(rank%m); i++) xs += lx[i];
    334:   if ((x < s) && ((m > 1) || (bx == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);

    336:   if (!ly) {
    337:     PetscMalloc1(n, &dd->ly);
    338:     ly   = dd->ly;
    339:     for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n));
    340:   }
    341:   y = ly[(rank % (m*n))/m];
    342:   if ((y < s) && ((n > 1) || (by == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1379:    Collective

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

    1397:    Output Parameter:
    1398: .  da - the resulting distributed array object

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

    1413:    Level: beginner

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

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

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

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

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

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

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