Actual source code: da3.c

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

     7:  #include src/dm/da/daimpl.h

     11: PetscErrorCode DAView_3d(DA da,PetscViewer viewer)
     12: {
     14:   PetscMPIInt    rank;
     15:   PetscTruth     iascii,isdraw;

     18:   MPI_Comm_rank(da->comm,&rank);

     20:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     21:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
     22:   if (iascii) {
     23:     PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D P %D m %D n %D p %D w %D s %D\n",
     24:                rank,da->M,da->N,da->P,da->m,da->n,da->p,da->w,da->s);
     25:     PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D, Z range of indices: %D %D\n",
     26:                da->xs,da->xe,da->ys,da->ye,da->zs,da->ze);
     27: #if !defined(PETSC_USE_COMPLEX)
     28:     if (da->coordinates) {
     29:       PetscInt  last;
     30:       PetscReal *coors;
     31:       VecGetArray(da->coordinates,&coors);
     32:       VecGetLocalSize(da->coordinates,&last);
     33:       last = last - 3;
     34:       PetscViewerASCIISynchronizedPrintf(viewer,"Lower left corner %g %g %g : Upper right %g %g %g\n",
     35:                coors[0],coors[1],coors[2],coors[last],coors[last+1],coors[last+2]);
     36:       VecRestoreArray(da->coordinates,&coors);
     37:     }
     38: #endif
     39:     PetscViewerFlush(viewer);
     40:   } else if (isdraw) {
     41:     PetscDraw       draw;
     42:     PetscReal     ymin = -1.0,ymax = (PetscReal)da->N;
     43:     PetscReal     xmin = -1.0,xmax = (PetscReal)((da->M+2)*da->P),x,y,ycoord,xcoord;
     44:     PetscInt        k,plane,base,*idx;
     45:     char       node[10];
     46:     PetscTruth isnull;

     48:     PetscViewerDrawGetDraw(viewer,0,&draw);
     49:     PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
     50:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
     51:     PetscDrawSynchronizedClear(draw);

     53:     /* first processor draw all node lines */
     54:     if (!rank) {
     55:       for (k=0; k<da->P; k++) {
     56:         ymin = 0.0; ymax = (PetscReal)(da->N - 1);
     57:         for (xmin=(PetscReal)(k*(da->M+1)); xmin<(PetscReal)(da->M+(k*(da->M+1))); xmin++) {
     58:           PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
     59:         }
     60:
     61:         xmin = (PetscReal)(k*(da->M+1)); xmax = xmin + (PetscReal)(da->M - 1);
     62:         for (ymin=0; ymin<(PetscReal)da->N; ymin++) {
     63:           PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
     64:         }
     65:       }
     66:     }
     67:     PetscDrawSynchronizedFlush(draw);
     68:     PetscDrawPause(draw);

     70:     for (k=0; k<da->P; k++) {  /*Go through and draw for each plane*/
     71:       if ((k >= da->zs) && (k < da->ze)) {
     72:         /* draw my box */
     73:         ymin = da->ys;
     74:         ymax = da->ye - 1;
     75:         xmin = da->xs/da->w    + (da->M+1)*k;
     76:         xmax =(da->xe-1)/da->w + (da->M+1)*k;

     78:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
     79:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
     80:         PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
     81:         PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

     83:         xmin = da->xs/da->w;
     84:         xmax =(da->xe-1)/da->w;

     86:         /* put in numbers*/
     87:         base = (da->base+(da->xe-da->xs)*(da->ye-da->ys)*(k-da->zs))/da->w;

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

     93:         for (y=ymin; y<=ymax; y++) {
     94:           for (x=xmin+(da->M+1)*k; x<=xmax+(da->M+1)*k; x++) {
     95:             sprintf(node,"%d",(int)base++);
     96:             PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
     97:           }
     98:         }
     99:
    100:       }
    101:     }
    102:     PetscDrawSynchronizedFlush(draw);
    103:     PetscDrawPause(draw);

    105:     for (k=0-da->s; k<da->P+da->s; k++) {
    106:       /* Go through and draw for each plane */
    107:       if ((k >= da->Zs) && (k < da->Ze)) {
    108:
    109:         /* overlay ghost numbers, useful for error checking */
    110:         base = (da->Xe-da->Xs)*(da->Ye-da->Ys)*(k-da->Zs); idx = da->idx;
    111:         plane=k;
    112:         /* Keep z wrap around points on the dradrawg */
    113:         if (k<0)    { plane=da->P+k; }
    114:         if (k>=da->P) { plane=k-da->P; }
    115:         ymin = da->Ys; ymax = da->Ye;
    116:         xmin = (da->M+1)*plane*da->w;
    117:         xmax = (da->M+1)*plane*da->w+da->M*da->w;
    118:         for (y=ymin; y<ymax; y++) {
    119:           for (x=xmin+da->Xs; x<xmin+da->Xe; x+=da->w) {
    120:             sprintf(node,"%d",(int)(idx[base]/da->w));
    121:             ycoord = y;
    122:             /*Keep y wrap around points on drawing */
    123:             if (y<0)      { ycoord = da->N+y; }

    125:             if (y>=da->N) { ycoord = y-da->N; }
    126:             xcoord = x;   /* Keep x wrap points on drawing */

    128:             if (x<xmin)  { xcoord = xmax - (xmin-x); }
    129:             if (x>=xmax) { xcoord = xmin + (x-xmax); }
    130:             PetscDrawString(draw,xcoord/da->w,ycoord,PETSC_DRAW_BLUE,node);
    131:             base+=da->w;
    132:           }
    133:         }
    134:       }
    135:     }
    136:     PetscDrawSynchronizedFlush(draw);
    137:     PetscDrawPause(draw);
    138:   } else {
    139:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for DA 3d",((PetscObject)viewer)->type_name);
    140:   }
    141:   return(0);
    142: }

    144: EXTERN PetscErrorCode DAPublish_Petsc(PetscObject);

    148: /*@C
    149:    DACreate3d - Creates an object that will manage the communication of three-dimensional
    150:    regular array data that is distributed across some processors.

    152:    Collective on MPI_Comm

    154:    Input Parameters:
    155: +  comm - MPI communicator
    156: .  wrap - type of periodicity the array should have, if any.  Use one
    157:           of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, DA_XYPERIODIC, DA_XYZPERIODIC, DA_XZPERIODIC, or DA_YZPERIODIC.
    158: .  stencil_type - Type of stencil (DA_STENCIL_STAR or DA_STENCIL_BOX)
    159: .  M,N,P - global dimension in each direction of the array (use -M, -N, and or -P to indicate that it may be set to a different value
    160:             from the command line with -da_grid_x <M> -da_grid_y <N> -da_grid_z <P>)
    161: .  m,n,p - corresponding number of processors in each dimension
    162:            (or PETSC_DECIDE to have calculated)
    163: .  dof - number of degrees of freedom per node
    164: .  lx, ly, lz - arrays containing the number of nodes in each cell along
    165:           the x, y, and z coordinates, or PETSC_NULL. If non-null, these
    166:           must be of length as m,n,p and the corresponding
    167:           m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
    168:           the ly[] must n, sum of the lz[] must be P
    169: -  s - stencil width

    171:    Output Parameter:
    172: .  inra - the resulting distributed array object

    174:    Options Database Key:
    175: +  -da_view - Calls DAView() at the conclusion of DACreate3d()
    176: .  -da_grid_x <nx> - number of grid points in x direction, if M < 0
    177: .  -da_grid_y <ny> - number of grid points in y direction, if N < 0
    178: .  -da_grid_z <nz> - number of grid points in z direction, if P < 0
    179: .  -da_refine_x - refinement ratio in x direction
    180: .  -da_refine_y - refinement ratio in y direction
    181: -  -da_refine_y - refinement ratio in z direction

    183:    Level: beginner

    185:    Notes:
    186:    The stencil type DA_STENCIL_STAR with width 1 corresponds to the
    187:    standard 7-pt stencil, while DA_STENCIL_BOX with width 1 denotes
    188:    the standard 27-pt stencil.

    190:    The array data itself is NOT stored in the DA, it is stored in Vec objects;
    191:    The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
    192:    and DACreateLocalVector() and calls to VecDuplicate() if more are needed.

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

    196: .seealso: DADestroy(), DAView(), DACreate1d(), DACreate2d(), DAGlobalToLocalBegin(), DAGetRefinementFactor(),
    197:           DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(), DASetRefinementFactor(),
    198:           DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()

    200: @*/
    201: PetscErrorCode PETSCDM_DLLEXPORT DACreate3d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type,PetscInt M,
    202:                PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof,PetscInt s,PetscInt *lx,PetscInt *ly,PetscInt *lz,DA *inra)
    203: {
    205:   PetscMPIInt    rank,size;
    206:   PetscInt       xs = 0,xe,ys = 0,ye,zs = 0,ze,x = 0,y = 0,z = 0,Xs,Xe,Ys,Ye,Zs,Ze,start,end,pm;
    207:   PetscInt       left,up,down,bottom,top,i,j,k,*idx,nn,*flx = 0,*fly = 0,*flz = 0;
    208:   PetscInt       n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n14;
    209:   PetscInt       n15,n16,n17,n18,n19,n20,n21,n22,n23,n24,n25,n26;
    210:   PetscInt       *bases,*ldims,x_t,y_t,z_t,s_t,base,count,s_x,s_y,s_z;
    211:   PetscInt       tM = M,tN = N,tP = P;
    212:   PetscInt       sn0 = 0,sn1 = 0,sn2 = 0,sn3 = 0,sn5 = 0,sn6 = 0,sn7 = 0;
    213:   PetscInt       sn8 = 0,sn9 = 0,sn11 = 0,sn15 = 0,sn24 = 0,sn25 = 0,sn26 = 0;
    214:   PetscInt       sn17 = 0,sn18 = 0,sn19 = 0,sn20 = 0,sn21 = 0,sn23 = 0,refine_x = 2, refine_y = 2, refine_z = 2;
    215:   PetscTruth     flg1;
    216:   DA             da;
    217:   Vec            local,global;
    218:   VecScatter     ltog,gtol;
    219:   IS             to,from;

    223:   *inra = 0;
    224: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    225:   DMInitializePackage(PETSC_NULL);
    226: #endif

    228:   if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
    229:   if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);

    231:   PetscOptionsBegin(comm,PETSC_NULL,"3d DA Options","DA");
    232:     if (M < 0){
    233:       tM   = -M;
    234:       PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate3d",tM,&tM,PETSC_NULL);
    235:     }
    236:     if (N < 0){
    237:       tN   = -N;
    238:       PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DACreate3d",tN,&tN,PETSC_NULL);
    239:     }
    240:     if (P < 0){
    241:       tP   = -P;
    242:       PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DACreate3d",tP,&tP,PETSC_NULL);
    243:     }
    244:     PetscOptionsInt("-da_processors_x","Number of processors in x direction","DACreate3d",m,&m,PETSC_NULL);
    245:     PetscOptionsInt("-da_processors_y","Number of processors in y direction","DACreate3d",n,&n,PETSC_NULL);
    246:     PetscOptionsInt("-da_processors_z","Number of processors in z direction","DACreate3d",p,&p,PETSC_NULL);
    247:     PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DASetRefinementFactor",refine_x,&refine_x,PETSC_NULL);
    248:     PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DASetRefinementFactor",refine_y,&refine_y,PETSC_NULL);
    249:     PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DASetRefinementFactor",refine_z,&refine_z,PETSC_NULL);
    250:   PetscOptionsEnd();
    251:   M = tM; N = tN; P = tP;

    253:   PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
    254:   da->bops->publish           = DAPublish_Petsc;
    255:   da->ops->createglobalvector = DACreateGlobalVector;
    256:   da->ops->getinterpolation   = DAGetInterpolation;
    257:   da->ops->getcoloring        = DAGetColoring;
    258:   da->ops->getmatrix          = DAGetMatrix;
    259:   da->ops->refine             = DARefine;

    261:   PetscLogObjectMemory(da,sizeof(struct _p_DA));
    262:   da->dim        = 3;
    263:   da->interptype = DA_Q1;
    264:   da->refine_x   = refine_x;
    265:   da->refine_y   = refine_y;
    266:   da->refine_z   = refine_z;
    267:   PetscMalloc(dof*sizeof(char*),&da->fieldname);
    268:   PetscMemzero(da->fieldname,dof*sizeof(char*));

    270:   MPI_Comm_size(comm,&size);
    271:   MPI_Comm_rank(comm,&rank);

    273:   if (m != PETSC_DECIDE) {
    274:     if (m < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);}
    275:     else if (m > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);}
    276:   }
    277:   if (n != PETSC_DECIDE) {
    278:     if (n < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);}
    279:     else if (n > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);}
    280:   }
    281:   if (p != PETSC_DECIDE) {
    282:     if (p < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);}
    283:     else if (p > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);}
    284:   }

    286:   /* Partition the array among the processors */
    287:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
    288:     m = size/(n*p);
    289:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
    290:     n = size/(m*p);
    291:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
    292:     p = size/(m*n);
    293:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
    294:     /* try for squarish distribution */
    295:     m = (int)(0.5 + sqrt(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
    296:     if (!m) m = 1;
    297:     while (m > 0) {
    298:       n = size/(m*p);
    299:       if (m*n*p == size) break;
    300:       m--;
    301:     }
    302:     if (!m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
    303:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
    304:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
    305:     /* try for squarish distribution */
    306:     m = (int)(0.5 + sqrt(((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) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
    314:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
    315:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
    316:     /* try for squarish distribution */
    317:     n = (int)(0.5 + sqrt(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
    318:     if (!n) n = 1;
    319:     while (n > 0) {
    320:       p = size/(m*n);
    321:       if (m*n*p == size) break;
    322:       n--;
    323:     }
    324:     if (!n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
    325:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
    326:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
    327:     /* try for squarish distribution */
    328:     n = (PetscInt)(0.5 + pow(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),1./3.));
    329:     if (!n) n = 1;
    330:     while (n > 0) {
    331:       pm = size/n;
    332:       if (n*pm == size) break;
    333:       n--;
    334:     }
    335:     if (!n) n = 1;
    336:     m = (PetscInt)(0.5 + sqrt(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
    337:     if (!m) m = 1;
    338:     while (m > 0) {
    339:       p = size/(m*n);
    340:       if (m*n*p == size) break;
    341:       m--;
    342:     }
    343:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
    344:   } else if (m*n*p != size) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

    346:   if (m*n*p != size) SETERRQ(PETSC_ERR_PLIB,"Could not find good partition");
    347:   if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
    348:   if (N < n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
    349:   if (P < p) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in z direction is too fine! %D %D",P,p);

    351:   /*
    352:      Determine locally owned region
    353:      [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
    354:   */
    355:   if (!lx) { /* user decided distribution */
    356:     PetscMalloc(m*sizeof(PetscInt),&lx);
    357:     flx = lx;
    358:     for (i=0; i<m; i++) {
    359:       lx[i] = M/m + ((M % m) > (i % m));
    360:     }
    361:   }
    362:   x  = lx[rank % m];
    363:   xs = 0;
    364:   for (i=0; i<(rank%m); i++) { xs += lx[i];}
    365:   if (m > 1 && x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column width is too thin for stencil! %D %D",x,s);

    367:   if (!ly) { /* user decided distribution */
    368:     PetscMalloc(n*sizeof(PetscInt),&ly);
    369:     fly = ly;
    370:     for (i=0; i<n; i++) {
    371:       ly[i] = N/n + ((N % n) > (i % n));
    372:     }
    373:   }
    374:   y  = ly[(rank % (m*n))/m];
    375:   if (n > 1 && y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row width is too thin for stencil! %D %D",y,s);
    376:   ys = 0;
    377:   for (i=0; i<(rank % (m*n))/m; i++) { ys += ly[i];}

    379:   if (!lz) { /* user decided distribution */
    380:     PetscMalloc(p*sizeof(PetscInt),&lz);
    381:     flz = lz;
    382:     for (i=0; i<p; i++) {
    383:       lz[i] = P/p + ((P % p) > (i % p));
    384:     }
    385:   }
    386:   z  = lz[rank/(m*n)];
    387:   if (p > 1 && z < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Plane width is too thin for stencil! %D %D",z,s);
    388:   zs = 0;
    389:   for (i=0; i<(rank/(m*n)); i++) { zs += lz[i];}
    390:   ye = ys + y;
    391:   xe = xs + x;
    392:   ze = zs + z;

    394:   /* determine ghost region */
    395:   /* Assume No Periodicity */
    396:   if (xs-s > 0) Xs = xs - s; else Xs = 0;
    397:   if (ys-s > 0) Ys = ys - s; else Ys = 0;
    398:   if (zs-s > 0) Zs = zs - s; else Zs = 0;
    399:   if (xe+s <= M) Xe = xe + s; else Xe = M;
    400:   if (ye+s <= N) Ye = ye + s; else Ye = N;
    401:   if (ze+s <= P) Ze = ze + s; else Ze = P;

    403:   /* X Periodic */
    404:   if (DAXPeriodic(wrap)){
    405:     Xs = xs - s;
    406:     Xe = xe + s;
    407:   }

    409:   /* Y Periodic */
    410:   if (DAYPeriodic(wrap)){
    411:     Ys = ys - s;
    412:     Ye = ye + s;
    413:   }

    415:   /* Z Periodic */
    416:   if (DAZPeriodic(wrap)){
    417:     Zs = zs - s;
    418:     Ze = ze + s;
    419:   }

    421:   /* Resize all X parameters to reflect w */
    422:   x   *= dof;
    423:   xs  *= dof;
    424:   xe  *= dof;
    425:   Xs  *= dof;
    426:   Xe  *= dof;
    427:   s_x  = s*dof;
    428:   s_y  = s;
    429:   s_z  = s;

    431:   /* determine starting point of each processor */
    432:   nn       = x*y*z;
    433:   PetscMalloc((2*size+1)*sizeof(PetscInt),&bases);
    434:   ldims    = (PetscInt*)(bases+size+1);
    435:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
    436:   bases[0] = 0;
    437:   for (i=1; i<=size; i++) {
    438:     bases[i] = ldims[i-1];
    439:   }
    440:   for (i=1; i<=size; i++) {
    441:     bases[i] += bases[i-1];
    442:   }

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

    452:   /* generate appropriate vector scatters */
    453:   /* local to global inserts non-ghost point region into global */
    454:   VecGetOwnershipRange(global,&start,&end);
    455:   ISCreateStride(comm,x*y*z,start,1,&to);

    457:   left   = xs - Xs;
    458:   bottom = ys - Ys; top = bottom + y;
    459:   down   = zs - Zs; up  = down + z;
    460:   count  = x*(top-bottom)*(up-down);
    461:   PetscMalloc(count*sizeof(PetscInt)/dof,&idx);
    462:   count  = 0;
    463:   for (i=down; i<up; i++) {
    464:     for (j=bottom; j<top; j++) {
    465:       for (k=0; k<x; k += dof) {
    466:         idx[count++] = (left+j*(Xe-Xs))+i*(Xe-Xs)*(Ye-Ys) + k;
    467:       }
    468:     }
    469:   }
    470:   ISCreateBlock(comm,dof,count,idx,&from);
    471:   PetscFree(idx);

    473:   VecScatterCreate(local,from,global,to,&ltog);
    474:   PetscLogObjectParent(da,to);
    475:   PetscLogObjectParent(da,from);
    476:   PetscLogObjectParent(da,ltog);
    477:   ISDestroy(from);
    478:   ISDestroy(to);

    480:   /* global to local must include ghost points */
    481:   if (stencil_type == DA_STENCIL_BOX) {
    482:     ISCreateStride(comm,(Xe-Xs)*(Ye-Ys)*(Ze-Zs),0,1,&to);
    483:   } else {
    484:     /* This is way ugly! We need to list the funny cross type region */
    485:     /* the bottom chunck */
    486:     left   = xs - Xs;
    487:     bottom = ys - Ys; top = bottom + y;
    488:     down   = zs - Zs;   up  = down + z;
    489:     count  = down*(top-bottom)*x + (up-down)*(bottom*x  + (top-bottom)*(Xe-Xs) + (Ye-Ys-top)*x) + (Ze-Zs-up)*(top-bottom)*x;
    490:     PetscMalloc(count*sizeof(PetscInt)/dof,&idx);
    491:     count  = 0;
    492:     for (i=0; i<down; i++) {
    493:       for (j=bottom; j<top; j++) {
    494:         for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
    495:       }
    496:     }
    497:     /* the middle piece */
    498:     for (i=down; i<up; i++) {
    499:       /* front */
    500:       for (j=0; j<bottom; j++) {
    501:         for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
    502:       }
    503:       /* middle */
    504:       for (j=bottom; j<top; j++) {
    505:         for (k=0; k<Xe-Xs; k += dof) idx[count++] = j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
    506:       }
    507:       /* back */
    508:       for (j=top; j<Ye-Ys; j++) {
    509:         for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
    510:       }
    511:     }
    512:     /* the top piece */
    513:     for (i=up; i<Ze-Zs; i++) {
    514:       for (j=bottom; j<top; j++) {
    515:         for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
    516:       }
    517:     }
    518:     ISCreateBlock(comm,dof,count,idx,&to);
    519:     PetscFree(idx);
    520:   }

    522:   /* determine who lies on each side of use stored in    n24 n25 n26
    523:                                                          n21 n22 n23
    524:                                                          n18 n19 n20

    526:                                                          n15 n16 n17
    527:                                                          n12     n14
    528:                                                          n9  n10 n11

    530:                                                          n6  n7  n8
    531:                                                          n3  n4  n5
    532:                                                          n0  n1  n2
    533:   */
    534:
    535:   /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
    536:
    537:   /* Assume Nodes are Internal to the Cube */
    538:
    539:   n0  = rank - m*n - m - 1;
    540:   n1  = rank - m*n - m;
    541:   n2  = rank - m*n - m + 1;
    542:   n3  = rank - m*n -1;
    543:   n4  = rank - m*n;
    544:   n5  = rank - m*n + 1;
    545:   n6  = rank - m*n + m - 1;
    546:   n7  = rank - m*n + m;
    547:   n8  = rank - m*n + m + 1;

    549:   n9  = rank - m - 1;
    550:   n10 = rank - m;
    551:   n11 = rank - m + 1;
    552:   n12 = rank - 1;
    553:   n14 = rank + 1;
    554:   n15 = rank + m - 1;
    555:   n16 = rank + m;
    556:   n17 = rank + m + 1;

    558:   n18 = rank + m*n - m - 1;
    559:   n19 = rank + m*n - m;
    560:   n20 = rank + m*n - m + 1;
    561:   n21 = rank + m*n - 1;
    562:   n22 = rank + m*n;
    563:   n23 = rank + m*n + 1;
    564:   n24 = rank + m*n + m - 1;
    565:   n25 = rank + m*n + m;
    566:   n26 = rank + m*n + m + 1;

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

    570:   if (xs == 0) { /* First assume not corner or edge */
    571:     n0  = rank       -1 - (m*n);
    572:     n3  = rank + m   -1 - (m*n);
    573:     n6  = rank + 2*m -1 - (m*n);
    574:     n9  = rank       -1;
    575:     n12 = rank + m   -1;
    576:     n15 = rank + 2*m -1;
    577:     n18 = rank       -1 + (m*n);
    578:     n21 = rank + m   -1 + (m*n);
    579:     n24 = rank + 2*m -1 + (m*n);
    580:    }

    582:   if (xe == M*dof) { /* First assume not corner or edge */
    583:     n2  = rank -2*m +1 - (m*n);
    584:     n5  = rank - m  +1 - (m*n);
    585:     n8  = rank      +1 - (m*n);
    586:     n11 = rank -2*m +1;
    587:     n14 = rank - m  +1;
    588:     n17 = rank      +1;
    589:     n20 = rank -2*m +1 + (m*n);
    590:     n23 = rank - m  +1 + (m*n);
    591:     n26 = rank      +1 + (m*n);
    592:   }

    594:   if (ys==0) { /* First assume not corner or edge */
    595:     n0  = rank + m * (n-1) -1 - (m*n);
    596:     n1  = rank + m * (n-1)    - (m*n);
    597:     n2  = rank + m * (n-1) +1 - (m*n);
    598:     n9  = rank + m * (n-1) -1;
    599:     n10 = rank + m * (n-1);
    600:     n11 = rank + m * (n-1) +1;
    601:     n18 = rank + m * (n-1) -1 + (m*n);
    602:     n19 = rank + m * (n-1)    + (m*n);
    603:     n20 = rank + m * (n-1) +1 + (m*n);
    604:   }

    606:   if (ye == N) { /* First assume not corner or edge */
    607:     n6  = rank - m * (n-1) -1 - (m*n);
    608:     n7  = rank - m * (n-1)    - (m*n);
    609:     n8  = rank - m * (n-1) +1 - (m*n);
    610:     n15 = rank - m * (n-1) -1;
    611:     n16 = rank - m * (n-1);
    612:     n17 = rank - m * (n-1) +1;
    613:     n24 = rank - m * (n-1) -1 + (m*n);
    614:     n25 = rank - m * (n-1)    + (m*n);
    615:     n26 = rank - m * (n-1) +1 + (m*n);
    616:   }
    617:
    618:   if (zs == 0) { /* First assume not corner or edge */
    619:     n0 = size - (m*n) + rank - m - 1;
    620:     n1 = size - (m*n) + rank - m;
    621:     n2 = size - (m*n) + rank - m + 1;
    622:     n3 = size - (m*n) + rank - 1;
    623:     n4 = size - (m*n) + rank;
    624:     n5 = size - (m*n) + rank + 1;
    625:     n6 = size - (m*n) + rank + m - 1;
    626:     n7 = size - (m*n) + rank + m ;
    627:     n8 = size - (m*n) + rank + m + 1;
    628:   }

    630:   if (ze == P) { /* First assume not corner or edge */
    631:     n18 = (m*n) - (size-rank) - m - 1;
    632:     n19 = (m*n) - (size-rank) - m;
    633:     n20 = (m*n) - (size-rank) - m + 1;
    634:     n21 = (m*n) - (size-rank) - 1;
    635:     n22 = (m*n) - (size-rank);
    636:     n23 = (m*n) - (size-rank) + 1;
    637:     n24 = (m*n) - (size-rank) + m - 1;
    638:     n25 = (m*n) - (size-rank) + m;
    639:     n26 = (m*n) - (size-rank) + m + 1;
    640:   }

    642:   if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
    643:     n0 = size - m*n + rank + m-1 - m;
    644:     n3 = size - m*n + rank + m-1;
    645:     n6 = size - m*n + rank + m-1 + m;
    646:   }
    647:
    648:   if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
    649:     n18 = m*n - (size - rank) + m-1 - m;
    650:     n21 = m*n - (size - rank) + m-1;
    651:     n24 = m*n - (size - rank) + m-1 + m;
    652:   }

    654:   if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
    655:     n0  = rank + m*n -1 - m*n;
    656:     n9  = rank + m*n -1;
    657:     n18 = rank + m*n -1 + m*n;
    658:   }

    660:   if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
    661:     n6  = rank - m*(n-1) + m-1 - m*n;
    662:     n15 = rank - m*(n-1) + m-1;
    663:     n24 = rank - m*(n-1) + m-1 + m*n;
    664:   }

    666:   if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
    667:     n2 = size - (m*n-rank) - (m-1) - m;
    668:     n5 = size - (m*n-rank) - (m-1);
    669:     n8 = size - (m*n-rank) - (m-1) + m;
    670:   }

    672:   if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
    673:     n20 = m*n - (size - rank) - (m-1) - m;
    674:     n23 = m*n - (size - rank) - (m-1);
    675:     n26 = m*n - (size - rank) - (m-1) + m;
    676:   }

    678:   if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
    679:     n2  = rank + m*(n-1) - (m-1) - m*n;
    680:     n11 = rank + m*(n-1) - (m-1);
    681:     n20 = rank + m*(n-1) - (m-1) + m*n;
    682:   }

    684:   if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
    685:     n8  = rank - m*n +1 - m*n;
    686:     n17 = rank - m*n +1;
    687:     n26 = rank - m*n +1 + m*n;
    688:   }

    690:   if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
    691:     n0 = size - m + rank -1;
    692:     n1 = size - m + rank;
    693:     n2 = size - m + rank +1;
    694:   }

    696:   if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
    697:     n18 = m*n - (size - rank) + m*(n-1) -1;
    698:     n19 = m*n - (size - rank) + m*(n-1);
    699:     n20 = m*n - (size - rank) + m*(n-1) +1;
    700:   }

    702:   if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
    703:     n6 = size - (m*n-rank) - m * (n-1) -1;
    704:     n7 = size - (m*n-rank) - m * (n-1);
    705:     n8 = size - (m*n-rank) - m * (n-1) +1;
    706:   }

    708:   if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
    709:     n24 = rank - (size-m) -1;
    710:     n25 = rank - (size-m);
    711:     n26 = rank - (size-m) +1;
    712:   }

    714:   /* Check for Corners */
    715:   if ((xs==0)   && (ys==0) && (zs==0)) { n0  = size -1;}
    716:   if ((xs==0)   && (ys==0) && (ze==P)) { n18 = m*n-1;}
    717:   if ((xs==0)   && (ye==N) && (zs==0)) { n6  = (size-1)-m*(n-1);}
    718:   if ((xs==0)   && (ye==N) && (ze==P)) { n24 = m-1;}
    719:   if ((xe==M*dof) && (ys==0) && (zs==0)) { n2  = size-m;}
    720:   if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
    721:   if ((xe==M*dof) && (ye==N) && (zs==0)) { n8  = size-m*n;}
    722:   if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}

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

    726:   /* If not X periodic */
    727:   if ((wrap != DA_XPERIODIC)  && (wrap != DA_XYPERIODIC) &&
    728:      (wrap != DA_XZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
    729:     if (xs==0)   {n0  = n3  = n6  = n9  = n12 = n15 = n18 = n21 = n24 = -2;}
    730:     if (xe==M*dof) {n2  = n5  = n8  = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
    731:   }

    733:   /* If not Y periodic */
    734:   if ((wrap != DA_YPERIODIC)  && (wrap != DA_XYPERIODIC) &&
    735:       (wrap != DA_YZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
    736:     if (ys==0)   {n0  = n1  = n2  = n9  = n10 = n11 = n18 = n19 = n20 = -2;}
    737:     if (ye==N)   {n6  = n7  = n8  = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
    738:   }

    740:   /* If not Z periodic */
    741:   if ((wrap != DA_ZPERIODIC)  && (wrap != DA_XZPERIODIC) &&
    742:       (wrap != DA_YZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
    743:     if (zs==0)   {n0  = n1  = n2  = n3  = n4  = n5  = n6  = n7  = n8  = -2;}
    744:     if (ze==P)   {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
    745:   }

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

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

    762:   nn = 0;

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

    790:     for (i=0; i<y; i++) {
    791:       if (n3 >= 0) { /* directly left */
    792:         x_t = lx[n3 % m]*dof;
    793:         y_t = y;
    794:         z_t = lz[n3 / (m*n)];
    795:         s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    796:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    797:       }

    799:       if (n4 >= 0) { /* middle */
    800:         x_t = x;
    801:         y_t = y;
    802:         z_t = lz[n4 / (m*n)];
    803:         s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    804:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    805:       }

    807:       if (n5 >= 0) { /* directly right */
    808:         x_t = lx[n5 % m]*dof;
    809:         y_t = y;
    810:         z_t = lz[n5 / (m*n)];
    811:         s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    812:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    813:       }
    814:     }

    816:     for (i=1; i<=s_y; i++) {
    817:       if (n6 >= 0) { /* left above */
    818:         x_t = lx[n6 % m]*dof;
    819:         y_t = ly[(n6 % (m*n))/m];
    820:         z_t = lz[n6 / (m*n)];
    821:         s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    822:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    823:       }
    824:       if (n7 >= 0) { /* directly above */
    825:         x_t = x;
    826:         y_t = ly[(n7 % (m*n))/m];
    827:         z_t = lz[n7 / (m*n)];
    828:         s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    829:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    830:       }
    831:       if (n8 >= 0) { /* right above */
    832:         x_t = lx[n8 % m]*dof;
    833:         y_t = ly[(n8 % (m*n))/m];
    834:         z_t = lz[n8 / (m*n)];
    835:         s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    836:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    837:       }
    838:     }
    839:   }

    841:   /* Middle Level */
    842:   for (k=0; k<z; k++) {
    843:     for (i=1; i<=s_y; i++) {
    844:       if (n9 >= 0) { /* left below */
    845:         x_t = lx[n9 % m]*dof;
    846:         y_t = ly[(n9 % (m*n))/m];
    847:         /* z_t = z; */
    848:         s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    849:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    850:       }
    851:       if (n10 >= 0) { /* directly below */
    852:         x_t = x;
    853:         y_t = ly[(n10 % (m*n))/m];
    854:         /* z_t = z; */
    855:         s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    856:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    857:       }
    858:       if (n11 >= 0) { /* right below */
    859:         x_t = lx[n11 % m]*dof;
    860:         y_t = ly[(n11 % (m*n))/m];
    861:         /* z_t = z; */
    862:         s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    863:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    864:       }
    865:     }

    867:     for (i=0; i<y; i++) {
    868:       if (n12 >= 0) { /* directly left */
    869:         x_t = lx[n12 % m]*dof;
    870:         y_t = y;
    871:         /* z_t = z; */
    872:         s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
    873:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    874:       }

    876:       /* Interior */
    877:       s_t = bases[rank] + i*x + k*x*y;
    878:       for (j=0; j<x; j++) { idx[nn++] = s_t++;}

    880:       if (n14 >= 0) { /* directly right */
    881:         x_t = lx[n14 % m]*dof;
    882:         y_t = y;
    883:         /* z_t = z; */
    884:         s_t = bases[n14] + i*x_t + k*x_t*y_t;
    885:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    886:       }
    887:     }

    889:     for (i=1; i<=s_y; i++) {
    890:       if (n15 >= 0) { /* left above */
    891:         x_t = lx[n15 % m]*dof;
    892:         y_t = ly[(n15 % (m*n))/m];
    893:         /* z_t = z; */
    894:         s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
    895:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    896:       }
    897:       if (n16 >= 0) { /* directly above */
    898:         x_t = x;
    899:         y_t = ly[(n16 % (m*n))/m];
    900:         /* z_t = z; */
    901:         s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
    902:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    903:       }
    904:       if (n17 >= 0) { /* right above */
    905:         x_t = lx[n17 % m]*dof;
    906:         y_t = ly[(n17 % (m*n))/m];
    907:         /* z_t = z; */
    908:         s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
    909:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    910:       }
    911:     }
    912:   }
    913:
    914:   /* Upper Level */
    915:   for (k=0; k<s_z; k++) {
    916:     for (i=1; i<=s_y; i++) {
    917:       if (n18 >= 0) { /* left below */
    918:         x_t = lx[n18 % m]*dof;
    919:         y_t = ly[(n18 % (m*n))/m];
    920:         /* z_t = lz[n18 / (m*n)]; */
    921:         s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    922:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    923:       }
    924:       if (n19 >= 0) { /* directly below */
    925:         x_t = x;
    926:         y_t = ly[(n19 % (m*n))/m];
    927:         /* z_t = lz[n19 / (m*n)]; */
    928:         s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    929:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    930:       }
    931:       if (n20 >= 0) { /* right below */
    932:         x_t = lx[n20 % m]*dof;
    933:         y_t = ly[(n20 % (m*n))/m];
    934:         /* z_t = lz[n20 / (m*n)]; */
    935:         s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    936:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    937:       }
    938:     }

    940:     for (i=0; i<y; i++) {
    941:       if (n21 >= 0) { /* directly left */
    942:         x_t = lx[n21 % m]*dof;
    943:         y_t = y;
    944:         /* z_t = lz[n21 / (m*n)]; */
    945:         s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
    946:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    947:       }

    949:       if (n22 >= 0) { /* middle */
    950:         x_t = x;
    951:         y_t = y;
    952:         /* z_t = lz[n22 / (m*n)]; */
    953:         s_t = bases[n22] + i*x_t + k*x_t*y_t;
    954:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    955:       }

    957:       if (n23 >= 0) { /* directly right */
    958:         x_t = lx[n23 % m]*dof;
    959:         y_t = y;
    960:         /* z_t = lz[n23 / (m*n)]; */
    961:         s_t = bases[n23] + i*x_t + k*x_t*y_t;
    962:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    963:       }
    964:     }

    966:     for (i=1; i<=s_y; i++) {
    967:       if (n24 >= 0) { /* left above */
    968:         x_t = lx[n24 % m]*dof;
    969:         y_t = ly[(n24 % (m*n))/m];
    970:         /* z_t = lz[n24 / (m*n)]; */
    971:         s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
    972:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    973:       }
    974:       if (n25 >= 0) { /* directly above */
    975:         x_t = x;
    976:         y_t = ly[(n25 % (m*n))/m];
    977:         /* z_t = lz[n25 / (m*n)]; */
    978:         s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
    979:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    980:       }
    981:       if (n26 >= 0) { /* right above */
    982:         x_t = lx[n26 % m]*dof;
    983:         y_t = ly[(n26 % (m*n))/m];
    984:         /* z_t = lz[n26 / (m*n)]; */
    985:         s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
    986:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    987:       }
    988:     }
    989:   }
    990:   base = bases[rank];
    991:   {
    992:     PetscInt nnn = nn/dof,*iidx;
    993:     PetscMalloc(nnn*sizeof(PetscInt),&iidx);
    994:     for (i=0; i<nnn; i++) {
    995:       iidx[i] = idx[dof*i];
    996:     }
    997:     ISCreateBlock(comm,dof,nnn,iidx,&from);
    998:     PetscFree(iidx);
    999:   }
    1000:   VecScatterCreate(global,from,local,to,&gtol);
    1001:   PetscLogObjectParent(da,gtol);
    1002:   PetscLogObjectParent(da,to);
    1003:   PetscLogObjectParent(da,from);
    1004:   ISDestroy(to);
    1005:   ISDestroy(from);
    1006:   da->stencil_type = stencil_type;
    1007:   da->M  = M;  da->N  = N; da->P = P;
    1008:   da->m  = m;  da->n  = n; da->p = p;
    1009:   da->w  = dof;  da->s  = s;
    1010:   da->xs = xs; da->xe = xe; da->ys = ys; da->ye = ye; da->zs = zs; da->ze = ze;
    1011:   da->Xs = Xs; da->Xe = Xe; da->Ys = Ys; da->Ye = Ye; da->Zs = Zs; da->Ze = Ze;

    1013:   VecDestroy(local);
    1014:   VecDestroy(global);

    1016:   if (stencil_type == DA_STENCIL_STAR) {
    1017:     /*
    1018:         Recompute the local to global mappings, this time keeping the
    1019:       information about the cross corner processor numbers.
    1020:     */
    1021:     n0  = sn0;  n1  = sn1;  n2  = sn2;  n3  = sn3;  n5  = sn5;  n6  = sn6; n7 = sn7;
    1022:     n8  = sn8;  n9  = sn9;  n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18;
    1023:     n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25;
    1024:     n26 = sn26;

    1026:     nn = 0;

    1028:     /* Bottom Level */
    1029:     for (k=0; k<s_z; k++) {
    1030:       for (i=1; i<=s_y; i++) {
    1031:         if (n0 >= 0) { /* left below */
    1032:           x_t = lx[n0 % m]*dof;
    1033:           y_t = ly[(n0 % (m*n))/m];
    1034:           z_t = lz[n0 / (m*n)];
    1035:           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;
    1036:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1037:         }
    1038:         if (n1 >= 0) { /* directly below */
    1039:           x_t = x;
    1040:           y_t = ly[(n1 % (m*n))/m];
    1041:           z_t = lz[n1 / (m*n)];
    1042:           s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
    1043:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1044:         }
    1045:         if (n2 >= 0) { /* right below */
    1046:           x_t = lx[n2 % m]*dof;
    1047:           y_t = ly[(n2 % (m*n))/m];
    1048:           z_t = lz[n2 / (m*n)];
    1049:           s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
    1050:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1051:         }
    1052:       }

    1054:       for (i=0; i<y; i++) {
    1055:         if (n3 >= 0) { /* directly left */
    1056:           x_t = lx[n3 % m]*dof;
    1057:           y_t = y;
    1058:           z_t = lz[n3 / (m*n)];
    1059:           s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1060:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1061:         }

    1063:         if (n4 >= 0) { /* middle */
    1064:           x_t = x;
    1065:           y_t = y;
    1066:           z_t = lz[n4 / (m*n)];
    1067:           s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1068:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1069:         }

    1071:         if (n5 >= 0) { /* directly right */
    1072:           x_t = lx[n5 % m]*dof;
    1073:           y_t = y;
    1074:           z_t = lz[n5 / (m*n)];
    1075:           s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1076:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1077:         }
    1078:       }

    1080:       for (i=1; i<=s_y; i++) {
    1081:         if (n6 >= 0) { /* left above */
    1082:           x_t = lx[n6 % m]*dof;
    1083:           y_t = ly[(n6 % (m*n))/m];
    1084:           z_t = lz[n6 / (m*n)];
    1085:           s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1086:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1087:         }
    1088:         if (n7 >= 0) { /* directly above */
    1089:           x_t = x;
    1090:           y_t = ly[(n7 % (m*n))/m];
    1091:           z_t = lz[n7 / (m*n)];
    1092:           s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1093:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1094:         }
    1095:         if (n8 >= 0) { /* right above */
    1096:           x_t = lx[n8 % m]*dof;
    1097:           y_t = ly[(n8 % (m*n))/m];
    1098:           z_t = lz[n8 / (m*n)];
    1099:           s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1100:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1101:         }
    1102:       }
    1103:     }

    1105:     /* Middle Level */
    1106:     for (k=0; k<z; k++) {
    1107:       for (i=1; i<=s_y; i++) {
    1108:         if (n9 >= 0) { /* left below */
    1109:           x_t = lx[n9 % m]*dof;
    1110:           y_t = ly[(n9 % (m*n))/m];
    1111:           /* z_t = z; */
    1112:           s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    1113:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1114:         }
    1115:         if (n10 >= 0) { /* directly below */
    1116:           x_t = x;
    1117:           y_t = ly[(n10 % (m*n))/m];
    1118:           /* z_t = z; */
    1119:           s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1120:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1121:         }
    1122:         if (n11 >= 0) { /* right below */
    1123:           x_t = lx[n11 % m]*dof;
    1124:           y_t = ly[(n11 % (m*n))/m];
    1125:           /* z_t = z; */
    1126:           s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1127:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1128:         }
    1129:       }

    1131:       for (i=0; i<y; i++) {
    1132:         if (n12 >= 0) { /* directly left */
    1133:           x_t = lx[n12 % m]*dof;
    1134:           y_t = y;
    1135:           /* z_t = z; */
    1136:           s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
    1137:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1138:         }

    1140:         /* Interior */
    1141:         s_t = bases[rank] + i*x + k*x*y;
    1142:         for (j=0; j<x; j++) { idx[nn++] = s_t++;}

    1144:         if (n14 >= 0) { /* directly right */
    1145:           x_t = lx[n14 % m]*dof;
    1146:           y_t = y;
    1147:           /* z_t = z; */
    1148:           s_t = bases[n14] + i*x_t + k*x_t*y_t;
    1149:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1150:         }
    1151:       }

    1153:       for (i=1; i<=s_y; i++) {
    1154:         if (n15 >= 0) { /* left above */
    1155:           x_t = lx[n15 % m]*dof;
    1156:           y_t = ly[(n15 % (m*n))/m];
    1157:           /* z_t = z; */
    1158:           s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
    1159:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1160:         }
    1161:         if (n16 >= 0) { /* directly above */
    1162:           x_t = x;
    1163:           y_t = ly[(n16 % (m*n))/m];
    1164:           /* z_t = z; */
    1165:           s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
    1166:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1167:         }
    1168:         if (n17 >= 0) { /* right above */
    1169:           x_t = lx[n17 % m]*dof;
    1170:           y_t = ly[(n17 % (m*n))/m];
    1171:           /* z_t = z; */
    1172:           s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
    1173:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1174:         }
    1175:       }
    1176:     }
    1177:
    1178:     /* Upper Level */
    1179:     for (k=0; k<s_z; k++) {
    1180:       for (i=1; i<=s_y; i++) {
    1181:         if (n18 >= 0) { /* left below */
    1182:           x_t = lx[n18 % m]*dof;
    1183:           y_t = ly[(n18 % (m*n))/m];
    1184:           /* z_t = lz[n18 / (m*n)]; */
    1185:           s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    1186:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1187:         }
    1188:         if (n19 >= 0) { /* directly below */
    1189:           x_t = x;
    1190:           y_t = ly[(n19 % (m*n))/m];
    1191:           /* z_t = lz[n19 / (m*n)]; */
    1192:           s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1193:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1194:         }
    1195:         if (n20 >= 0) { /* right below */
    1196:           x_t = lx[n20 % m]*dof;
    1197:           y_t = ly[(n20 % (m*n))/m];
    1198:           /* z_t = lz[n20 / (m*n)]; */
    1199:           s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1200:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1201:         }
    1202:       }

    1204:       for (i=0; i<y; i++) {
    1205:         if (n21 >= 0) { /* directly left */
    1206:           x_t = lx[n21 % m]*dof;
    1207:           y_t = y;
    1208:           /* z_t = lz[n21 / (m*n)]; */
    1209:           s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
    1210:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1211:         }

    1213:         if (n22 >= 0) { /* middle */
    1214:           x_t = x;
    1215:           y_t = y;
    1216:           /* z_t = lz[n22 / (m*n)]; */
    1217:           s_t = bases[n22] + i*x_t + k*x_t*y_t;
    1218:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1219:         }

    1221:         if (n23 >= 0) { /* directly right */
    1222:           x_t = lx[n23 % m]*dof;
    1223:           y_t = y;
    1224:           /* z_t = lz[n23 / (m*n)]; */
    1225:           s_t = bases[n23] + i*x_t + k*x_t*y_t;
    1226:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1227:         }
    1228:       }

    1230:       for (i=1; i<=s_y; i++) {
    1231:         if (n24 >= 0) { /* left above */
    1232:           x_t = lx[n24 % m]*dof;
    1233:           y_t = ly[(n24 % (m*n))/m];
    1234:           /* z_t = lz[n24 / (m*n)]; */
    1235:           s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
    1236:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1237:         }
    1238:         if (n25 >= 0) { /* directly above */
    1239:           x_t = x;
    1240:           y_t = ly[(n25 % (m*n))/m];
    1241:           /* z_t = lz[n25 / (m*n)]; */
    1242:           s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
    1243:           for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1244:         }
    1245:         if (n26 >= 0) { /* right above */
    1246:           x_t = lx[n26 % m]*dof;
    1247:           y_t = ly[(n26 % (m*n))/m];
    1248:           /* z_t = lz[n26 / (m*n)]; */
    1249:           s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
    1250:           for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1251:         }
    1252:       }
    1253:     }
    1254:   }
    1255:   /* redo idx to include "missing" ghost points */
    1256:   /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
    1257:
    1258:   /* Assume Nodes are Internal to the Cube */
    1259:
    1260:   n0  = rank - m*n - m - 1;
    1261:   n1  = rank - m*n - m;
    1262:   n2  = rank - m*n - m + 1;
    1263:   n3  = rank - m*n -1;
    1264:   n4  = rank - m*n;
    1265:   n5  = rank - m*n + 1;
    1266:   n6  = rank - m*n + m - 1;
    1267:   n7  = rank - m*n + m;
    1268:   n8  = rank - m*n + m + 1;

    1270:   n9  = rank - m - 1;
    1271:   n10 = rank - m;
    1272:   n11 = rank - m + 1;
    1273:   n12 = rank - 1;
    1274:   n14 = rank + 1;
    1275:   n15 = rank + m - 1;
    1276:   n16 = rank + m;
    1277:   n17 = rank + m + 1;

    1279:   n18 = rank + m*n - m - 1;
    1280:   n19 = rank + m*n - m;
    1281:   n20 = rank + m*n - m + 1;
    1282:   n21 = rank + m*n - 1;
    1283:   n22 = rank + m*n;
    1284:   n23 = rank + m*n + 1;
    1285:   n24 = rank + m*n + m - 1;
    1286:   n25 = rank + m*n + m;
    1287:   n26 = rank + m*n + m + 1;

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

    1291:   if (xs == 0) { /* First assume not corner or edge */
    1292:     n0  = rank       -1 - (m*n);
    1293:     n3  = rank + m   -1 - (m*n);
    1294:     n6  = rank + 2*m -1 - (m*n);
    1295:     n9  = rank       -1;
    1296:     n12 = rank + m   -1;
    1297:     n15 = rank + 2*m -1;
    1298:     n18 = rank       -1 + (m*n);
    1299:     n21 = rank + m   -1 + (m*n);
    1300:     n24 = rank + 2*m -1 + (m*n);
    1301:    }

    1303:   if (xe == M*dof) { /* First assume not corner or edge */
    1304:     n2  = rank -2*m +1 - (m*n);
    1305:     n5  = rank - m  +1 - (m*n);
    1306:     n8  = rank      +1 - (m*n);
    1307:     n11 = rank -2*m +1;
    1308:     n14 = rank - m  +1;
    1309:     n17 = rank      +1;
    1310:     n20 = rank -2*m +1 + (m*n);
    1311:     n23 = rank - m  +1 + (m*n);
    1312:     n26 = rank      +1 + (m*n);
    1313:   }

    1315:   if (ys==0) { /* First assume not corner or edge */
    1316:     n0  = rank + m * (n-1) -1 - (m*n);
    1317:     n1  = rank + m * (n-1)    - (m*n);
    1318:     n2  = rank + m * (n-1) +1 - (m*n);
    1319:     n9  = rank + m * (n-1) -1;
    1320:     n10 = rank + m * (n-1);
    1321:     n11 = rank + m * (n-1) +1;
    1322:     n18 = rank + m * (n-1) -1 + (m*n);
    1323:     n19 = rank + m * (n-1)    + (m*n);
    1324:     n20 = rank + m * (n-1) +1 + (m*n);
    1325:   }

    1327:   if (ye == N) { /* First assume not corner or edge */
    1328:     n6  = rank - m * (n-1) -1 - (m*n);
    1329:     n7  = rank - m * (n-1)    - (m*n);
    1330:     n8  = rank - m * (n-1) +1 - (m*n);
    1331:     n15 = rank - m * (n-1) -1;
    1332:     n16 = rank - m * (n-1);
    1333:     n17 = rank - m * (n-1) +1;
    1334:     n24 = rank - m * (n-1) -1 + (m*n);
    1335:     n25 = rank - m * (n-1)    + (m*n);
    1336:     n26 = rank - m * (n-1) +1 + (m*n);
    1337:   }
    1338:
    1339:   if (zs == 0) { /* First assume not corner or edge */
    1340:     n0 = size - (m*n) + rank - m - 1;
    1341:     n1 = size - (m*n) + rank - m;
    1342:     n2 = size - (m*n) + rank - m + 1;
    1343:     n3 = size - (m*n) + rank - 1;
    1344:     n4 = size - (m*n) + rank;
    1345:     n5 = size - (m*n) + rank + 1;
    1346:     n6 = size - (m*n) + rank + m - 1;
    1347:     n7 = size - (m*n) + rank + m ;
    1348:     n8 = size - (m*n) + rank + m + 1;
    1349:   }

    1351:   if (ze == P) { /* First assume not corner or edge */
    1352:     n18 = (m*n) - (size-rank) - m - 1;
    1353:     n19 = (m*n) - (size-rank) - m;
    1354:     n20 = (m*n) - (size-rank) - m + 1;
    1355:     n21 = (m*n) - (size-rank) - 1;
    1356:     n22 = (m*n) - (size-rank);
    1357:     n23 = (m*n) - (size-rank) + 1;
    1358:     n24 = (m*n) - (size-rank) + m - 1;
    1359:     n25 = (m*n) - (size-rank) + m;
    1360:     n26 = (m*n) - (size-rank) + m + 1;
    1361:   }

    1363:   if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
    1364:     n0 = size - m*n + rank + m-1 - m;
    1365:     n3 = size - m*n + rank + m-1;
    1366:     n6 = size - m*n + rank + m-1 + m;
    1367:   }
    1368:
    1369:   if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
    1370:     n18 = m*n - (size - rank) + m-1 - m;
    1371:     n21 = m*n - (size - rank) + m-1;
    1372:     n24 = m*n - (size - rank) + m-1 + m;
    1373:   }

    1375:   if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
    1376:     n0  = rank + m*n -1 - m*n;
    1377:     n9  = rank + m*n -1;
    1378:     n18 = rank + m*n -1 + m*n;
    1379:   }

    1381:   if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
    1382:     n6  = rank - m*(n-1) + m-1 - m*n;
    1383:     n15 = rank - m*(n-1) + m-1;
    1384:     n24 = rank - m*(n-1) + m-1 + m*n;
    1385:   }

    1387:   if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
    1388:     n2 = size - (m*n-rank) - (m-1) - m;
    1389:     n5 = size - (m*n-rank) - (m-1);
    1390:     n8 = size - (m*n-rank) - (m-1) + m;
    1391:   }

    1393:   if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
    1394:     n20 = m*n - (size - rank) - (m-1) - m;
    1395:     n23 = m*n - (size - rank) - (m-1);
    1396:     n26 = m*n - (size - rank) - (m-1) + m;
    1397:   }

    1399:   if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
    1400:     n2  = rank + m*(n-1) - (m-1) - m*n;
    1401:     n11 = rank + m*(n-1) - (m-1);
    1402:     n20 = rank + m*(n-1) - (m-1) + m*n;
    1403:   }

    1405:   if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
    1406:     n8  = rank - m*n +1 - m*n;
    1407:     n17 = rank - m*n +1;
    1408:     n26 = rank - m*n +1 + m*n;
    1409:   }

    1411:   if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
    1412:     n0 = size - m + rank -1;
    1413:     n1 = size - m + rank;
    1414:     n2 = size - m + rank +1;
    1415:   }

    1417:   if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
    1418:     n18 = m*n - (size - rank) + m*(n-1) -1;
    1419:     n19 = m*n - (size - rank) + m*(n-1);
    1420:     n20 = m*n - (size - rank) + m*(n-1) +1;
    1421:   }

    1423:   if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
    1424:     n6 = size - (m*n-rank) - m * (n-1) -1;
    1425:     n7 = size - (m*n-rank) - m * (n-1);
    1426:     n8 = size - (m*n-rank) - m * (n-1) +1;
    1427:   }

    1429:   if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
    1430:     n24 = rank - (size-m) -1;
    1431:     n25 = rank - (size-m);
    1432:     n26 = rank - (size-m) +1;
    1433:   }

    1435:   /* Check for Corners */
    1436:   if ((xs==0)   && (ys==0) && (zs==0)) { n0  = size -1;}
    1437:   if ((xs==0)   && (ys==0) && (ze==P)) { n18 = m*n-1;}
    1438:   if ((xs==0)   && (ye==N) && (zs==0)) { n6  = (size-1)-m*(n-1);}
    1439:   if ((xs==0)   && (ye==N) && (ze==P)) { n24 = m-1;}
    1440:   if ((xe==M*dof) && (ys==0) && (zs==0)) { n2  = size-m;}
    1441:   if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
    1442:   if ((xe==M*dof) && (ye==N) && (zs==0)) { n8  = size-m*n;}
    1443:   if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}

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

    1447:   /* If not X periodic */
    1448:   if (!DAXPeriodic(wrap)){
    1449:     if (xs==0)   {n0  = n3  = n6  = n9  = n12 = n15 = n18 = n21 = n24 = -2;}
    1450:     if (xe==M*dof) {n2  = n5  = n8  = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
    1451:   }

    1453:   /* If not Y periodic */
    1454:   if (!DAYPeriodic(wrap)){
    1455:     if (ys==0)   {n0  = n1  = n2  = n9  = n10 = n11 = n18 = n19 = n20 = -2;}
    1456:     if (ye==N)   {n6  = n7  = n8  = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
    1457:   }

    1459:   /* If not Z periodic */
    1460:   if (!DAZPeriodic(wrap)){
    1461:     if (zs==0)   {n0  = n1  = n2  = n3  = n4  = n5  = n6  = n7  = n8  = -2;}
    1462:     if (ze==P)   {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
    1463:   }

    1465:   nn = 0;

    1467:   /* Bottom Level */
    1468:   for (k=0; k<s_z; k++) {
    1469:     for (i=1; i<=s_y; i++) {
    1470:       if (n0 >= 0) { /* left below */
    1471:         x_t = lx[n0 % m]*dof;
    1472:         y_t = ly[(n0 % (m*n))/m];
    1473:         z_t = lz[n0 / (m*n)];
    1474:         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;
    1475:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1476:       }
    1477:       if (n1 >= 0) { /* directly below */
    1478:         x_t = x;
    1479:         y_t = ly[(n1 % (m*n))/m];
    1480:         z_t = lz[n1 / (m*n)];
    1481:         s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
    1482:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1483:       }
    1484:       if (n2 >= 0) { /* right below */
    1485:         x_t = lx[n2 % m]*dof;
    1486:         y_t = ly[(n2 % (m*n))/m];
    1487:         z_t = lz[n2 / (m*n)];
    1488:         s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
    1489:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1490:       }
    1491:     }

    1493:     for (i=0; i<y; i++) {
    1494:       if (n3 >= 0) { /* directly left */
    1495:         x_t = lx[n3 % m]*dof;
    1496:         y_t = y;
    1497:         z_t = lz[n3 / (m*n)];
    1498:         s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1499:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1500:       }

    1502:       if (n4 >= 0) { /* middle */
    1503:         x_t = x;
    1504:         y_t = y;
    1505:         z_t = lz[n4 / (m*n)];
    1506:         s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1507:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1508:       }

    1510:       if (n5 >= 0) { /* directly right */
    1511:         x_t = lx[n5 % m]*dof;
    1512:         y_t = y;
    1513:         z_t = lz[n5 / (m*n)];
    1514:         s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1515:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1516:       }
    1517:     }

    1519:     for (i=1; i<=s_y; i++) {
    1520:       if (n6 >= 0) { /* left above */
    1521:         x_t = lx[n6 % m]*dof;
    1522:         y_t = ly[(n6 % (m*n))/m];
    1523:         z_t = lz[n6 / (m*n)];
    1524:         s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1525:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1526:       }
    1527:       if (n7 >= 0) { /* directly above */
    1528:         x_t = x;
    1529:         y_t = ly[(n7 % (m*n))/m];
    1530:         z_t = lz[n7 / (m*n)];
    1531:         s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1532:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1533:       }
    1534:       if (n8 >= 0) { /* right above */
    1535:         x_t = lx[n8 % m]*dof;
    1536:         y_t = ly[(n8 % (m*n))/m];
    1537:         z_t = lz[n8 / (m*n)];
    1538:         s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
    1539:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1540:       }
    1541:     }
    1542:   }

    1544:   /* Middle Level */
    1545:   for (k=0; k<z; k++) {
    1546:     for (i=1; i<=s_y; i++) {
    1547:       if (n9 >= 0) { /* left below */
    1548:         x_t = lx[n9 % m]*dof;
    1549:         y_t = ly[(n9 % (m*n))/m];
    1550:         /* z_t = z; */
    1551:         s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    1552:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1553:       }
    1554:       if (n10 >= 0) { /* directly below */
    1555:         x_t = x;
    1556:         y_t = ly[(n10 % (m*n))/m];
    1557:         /* z_t = z; */
    1558:         s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1559:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1560:       }
    1561:       if (n11 >= 0) { /* right below */
    1562:         x_t = lx[n11 % m]*dof;
    1563:         y_t = ly[(n11 % (m*n))/m];
    1564:         /* z_t = z; */
    1565:         s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1566:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1567:       }
    1568:     }

    1570:     for (i=0; i<y; i++) {
    1571:       if (n12 >= 0) { /* directly left */
    1572:         x_t = lx[n12 % m]*dof;
    1573:         y_t = y;
    1574:         /* z_t = z; */
    1575:         s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
    1576:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1577:       }

    1579:       /* Interior */
    1580:       s_t = bases[rank] + i*x + k*x*y;
    1581:       for (j=0; j<x; j++) { idx[nn++] = s_t++;}

    1583:       if (n14 >= 0) { /* directly right */
    1584:         x_t = lx[n14 % m]*dof;
    1585:         y_t = y;
    1586:         /* z_t = z; */
    1587:         s_t = bases[n14] + i*x_t + k*x_t*y_t;
    1588:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1589:       }
    1590:     }

    1592:     for (i=1; i<=s_y; i++) {
    1593:       if (n15 >= 0) { /* left above */
    1594:         x_t = lx[n15 % m]*dof;
    1595:         y_t = ly[(n15 % (m*n))/m];
    1596:         /* z_t = z; */
    1597:         s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
    1598:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1599:       }
    1600:       if (n16 >= 0) { /* directly above */
    1601:         x_t = x;
    1602:         y_t = ly[(n16 % (m*n))/m];
    1603:         /* z_t = z; */
    1604:         s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
    1605:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1606:       }
    1607:       if (n17 >= 0) { /* right above */
    1608:         x_t = lx[n17 % m]*dof;
    1609:         y_t = ly[(n17 % (m*n))/m];
    1610:         /* z_t = z; */
    1611:         s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
    1612:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1613:       }
    1614:     }
    1615:   }
    1616:
    1617:   /* Upper Level */
    1618:   for (k=0; k<s_z; k++) {
    1619:     for (i=1; i<=s_y; i++) {
    1620:       if (n18 >= 0) { /* left below */
    1621:         x_t = lx[n18 % m]*dof;
    1622:         y_t = ly[(n18 % (m*n))/m];
    1623:         /* z_t = lz[n18 / (m*n)]; */
    1624:         s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
    1625:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1626:       }
    1627:       if (n19 >= 0) { /* directly below */
    1628:         x_t = x;
    1629:         y_t = ly[(n19 % (m*n))/m];
    1630:         /* z_t = lz[n19 / (m*n)]; */
    1631:         s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1632:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1633:       }
    1634:       if (n20 >= 0) { /* right belodof */
    1635:         x_t = lx[n20 % m]*dof;
    1636:         y_t = ly[(n20 % (m*n))/m];
    1637:         /* z_t = lz[n20 / (m*n)]; */
    1638:         s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
    1639:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1640:       }
    1641:     }

    1643:     for (i=0; i<y; i++) {
    1644:       if (n21 >= 0) { /* directly left */
    1645:         x_t = lx[n21 % m]*dof;
    1646:         y_t = y;
    1647:         /* z_t = lz[n21 / (m*n)]; */
    1648:         s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
    1649:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1650:       }

    1652:       if (n22 >= 0) { /* middle */
    1653:         x_t = x;
    1654:         y_t = y;
    1655:         /* z_t = lz[n22 / (m*n)]; */
    1656:         s_t = bases[n22] + i*x_t + k*x_t*y_t;
    1657:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1658:       }

    1660:       if (n23 >= 0) { /* directly right */
    1661:         x_t = lx[n23 % m]*dof;
    1662:         y_t = y;
    1663:         /* z_t = lz[n23 / (m*n)]; */
    1664:         s_t = bases[n23] + i*x_t + k*x_t*y_t;
    1665:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1666:       }
    1667:     }

    1669:     for (i=1; i<=s_y; i++) {
    1670:       if (n24 >= 0) { /* left above */
    1671:         x_t = lx[n24 % m]*dof;
    1672:         y_t = ly[(n24 % (m*n))/m];
    1673:         /* z_t = lz[n24 / (m*n)]; */
    1674:         s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
    1675:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1676:       }
    1677:       if (n25 >= 0) { /* directly above */
    1678:         x_t = x;
    1679:         y_t = ly[(n25 % (m*n))/m];
    1680:         /* z_t = lz[n25 / (m*n)]; */
    1681:         s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
    1682:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    1683:       }
    1684:       if (n26 >= 0) { /* right above */
    1685:         x_t = lx[n26 % m]*dof;
    1686:         y_t = ly[(n26 % (m*n))/m];
    1687:         /* z_t = lz[n26 / (m*n)]; */
    1688:         s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
    1689:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    1690:       }
    1691:     }
    1692:   }
    1693:   PetscFree(bases);
    1694:   da->gtol      = gtol;
    1695:   da->ltog      = ltog;
    1696:   da->idx       = idx;
    1697:   da->Nl        = nn;
    1698:   da->base      = base;
    1699:   da->ops->view = DAView_3d;
    1700:   da->wrap      = wrap;
    1701:   *inra = da;

    1703:   /*
    1704:      Set the local to global ordering in the global vector, this allows use
    1705:      of VecSetValuesLocal().
    1706:   */
    1707:   ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
    1708:   ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
    1709:   PetscLogObjectParent(da,da->ltogmap);

    1711:   da->ltol = PETSC_NULL;
    1712:   da->ao   = PETSC_NULL;

    1714:   if (!flx) {
    1715:     PetscMalloc(m*sizeof(PetscInt),&flx);
    1716:     PetscMemcpy(flx,lx,m*sizeof(PetscInt));
    1717:   }
    1718:   if (!fly) {
    1719:     PetscMalloc(n*sizeof(PetscInt),&fly);
    1720:     PetscMemcpy(fly,ly,n*sizeof(PetscInt));
    1721:   }
    1722:   if (!flz) {
    1723:     PetscMalloc(p*sizeof(PetscInt),&flz);
    1724:     PetscMemcpy(flz,lz,p*sizeof(PetscInt));
    1725:   }
    1726:   da->lx = flx;
    1727:   da->ly = fly;
    1728:   da->lz = flz;

    1730:   PetscOptionsHasName(PETSC_NULL,"-da_view",&flg1);
    1731:   if (flg1) {DAView(da,PETSC_VIEWER_STDOUT_(da->comm));}
    1732:   PetscOptionsHasName(PETSC_NULL,"-da_view_draw",&flg1);
    1733:   if (flg1) {DAView(da,PETSC_VIEWER_DRAW_(da->comm));}
    1734:   PetscOptionsHasName(PETSC_NULL,"-help",&flg1);
    1735:   if (flg1) {DAPrintHelp(da);}
    1736:   PetscPublishAll(da);

    1738:   return(0);
    1739: }

    1741: /*@C
    1742:    DACreate - Creates an object that will manage the communication of regular array data that is distributed across some processors
    1743:        in 1, 2 or 3 dimensions

    1745:    Collective on MPI_Comm

    1747:    See the manual pages for the routines for each dimension.

    1749:    Level: beginner

    1751:
    1752: .keywords: distributed array, create, three-dimensional

    1754: .seealso: DACreate1d(), DACreate2d(), DACreate3d()

    1756: @*/
    1757: PetscErrorCode PETSCDM_DLLEXPORT DACreate(MPI_Comm comm,PetscInt dim,DAPeriodicType wrap,DAStencilType stencil_type,PetscInt M,
    1758:          PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof,PetscInt s,PetscInt *lx,PetscInt *ly,PetscInt *lz,DA *inra)
    1759: {

    1762:   if (dim == 3) {
    1763:     DACreate3d(comm,wrap,stencil_type,M,N,P,m,n,p,dof,s,lx,ly,lz,inra);
    1764:   } else if (dim == 2) {
    1765:     DACreate2d(comm,wrap,stencil_type,M,N,m,n,dof,s,lx,ly,inra);
    1766:   } else if (dim == 1) {
    1767:     DACreate1d(comm,wrap,M,dof,s,lx,inra);
    1768:   }
    1769:   return(0);
    1770: }