Actual source code: stag3d.c

petsc-master 2020-08-15
Report Typos and Errors
  1: /* Functions specific to the 3-dimensional implementation of DMStag */
  2:  #include <petsc/private/dmstagimpl.h>

  4: /*@C
  5:   DMStagCreate3d - Create an object to manage data living on the faces, edges, and vertices of a parallelized regular 3D grid.

  7:   Collective

  9:   Input Parameters:
 10: + comm - MPI communicator
 11: . bndx,bndy,bndz - boundary type: DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, or DM_BOUNDARY_GHOSTED
 12: . M,N,P - global number of grid points in x,y directions
 13: . m,n,p - number of ranks in the x,y directions (may be PETSC_DECIDE)
 14: . dof0 - number of degrees of freedom per vertex/point/node/0-cell
 15: . dof1 - number of degrees of freedom per edge/1-cell
 16: . dof2 - number of degrees of freedom per face/2-cell
 17: . dof3 - number of degrees of freedom per element/3-cell
 18: . stencilType - ghost/halo region type: DMSTAG_STENCIL_NONE, DMSTAG_STENCIL_BOX, or DMSTAG_STENCIL_STAR
 19: . stencilWidth - width, in elements, of halo/ghost region
 20: - lx,ly,lz - array sof local x,y,z element counts, of length equal to m,n,p, summing to M,N,P

 22:   Output Parameter:
 23: . dm - the new DMStag object

 25:   Options Database Keys:
 26: + -dm_view - calls DMViewFromOptions() a the conclusion of DMSetUp()
 27: . -stag_grid_x <nx> - number of elements in the x direction
 28: . -stag_grid_y <ny> - number of elements in the y direction
 29: . -stag_grid_z <nz> - number of elements in the z direction
 30: . -stag_ranks_x <rx> - number of ranks in the x direction
 31: . -stag_ranks_y <ry> - number of ranks in the y direction
 32: . -stag_ranks_z <rz> - number of ranks in the z direction
 33: . -stag_ghost_stencil_width - width of ghost region, in elements
 34: . -stag_boundary_type x <none,ghosted,periodic> - DMBoundaryType value
 35: . -stag_boundary_type y <none,ghosted,periodic> - DMBoundaryType value
 36: - -stag_boundary_type z <none,ghosted,periodic> - DMBoundaryType value

 38:   Notes:
 39:   You must call DMSetUp() after this call before using the DM.
 40:   If you wish to use the options database (see the keys above) to change values in the DMStag, you must call
 41:   DMSetFromOptions() after this function but before DMSetUp().

 43:   Level: beginner

 45: .seealso: DMSTAG, DMStagCreate1d(), DMStagCreate2d(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateLocalVector(), DMLocalToGlobalBegin(), DMDACreate3d()
 46: @*/
 47: PETSC_EXTERN PetscErrorCode DMStagCreate3d(MPI_Comm comm,DMBoundaryType bndx,DMBoundaryType bndy,DMBoundaryType bndz,PetscInt M,PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DMStagStencilType stencilType,PetscInt stencilWidth,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[],DM* dm)
 48: {

 52:   DMCreate(comm,dm);
 53:   DMSetDimension(*dm,3);
 54:   DMStagInitialize(bndx,bndy,bndz,M,N,P,m,n,p,dof0,dof1,dof2,dof3,stencilType,stencilWidth,lx,ly,lz,*dm);
 55:   return(0);
 56: }

 58: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_3d(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
 59: {
 61:   DM_Stag        *stagCoord;
 62:   DM             dmCoord;
 63:   Vec            coordLocal,coord;
 64:   PetscReal      h[3],min[3];
 65:   PetscScalar    ****arr;
 66:   PetscInt       ind[3],start[3],n[3],nExtra[3],s,c;
 67:   PetscInt       ibackdownleft,ibackdown,ibackleft,iback,idownleft,idown,ileft,ielement;

 70:   DMGetCoordinateDM(dm,&dmCoord);
 71:   stagCoord = (DM_Stag*) dmCoord->data;
 72:   for (s=0; s<4; ++s) {
 73:     if (stagCoord->dof[s] !=0 && stagCoord->dof[s] != 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Coordinate DM in 3 dimensions must have 0 or 3 dof on each stratum, but stratum %d has %d dof",s,stagCoord->dof[s]);
 74:   }
 75:   DMGetLocalVector(dmCoord,&coordLocal);
 76:   DMStagVecGetArray(dmCoord,coordLocal,&arr);
 77:   if (stagCoord->dof[0]) {
 78:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN_LEFT,0,&ibackdownleft);
 79:   }
 80:   if (stagCoord->dof[1]) {
 81:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN     ,0,&ibackdown    );
 82:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_LEFT     ,0,&ibackleft    );
 83:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN_LEFT     ,0,&idownleft    );
 84:   }
 85:   if (stagCoord->dof[2]) {
 86:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK          ,0,&iback        );
 87:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN          ,0,&idown        );
 88:     DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT          ,0,&ileft        );
 89:   }
 90:   if (stagCoord->dof[3]) {
 91:     DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT       ,0,&ielement     );
 92:     DMStagGetCorners(dmCoord,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
 93:   }
 94:   min[0] = xmin; min[1]= ymin; min[2] = zmin;
 95:   h[0] = (xmax-xmin)/stagCoord->N[0];
 96:   h[1] = (ymax-ymin)/stagCoord->N[1];
 97:   h[2] = (zmax-zmin)/stagCoord->N[2];

 99:   for(ind[2]=start[2]; ind[2]<start[2] + n[2] + nExtra[2]; ++ind[2]) {
100:     for(ind[1]=start[1]; ind[1]<start[1] + n[1] + nExtra[1]; ++ind[1]) {
101:       for(ind[0]=start[0]; ind[0]<start[0] + n[0] + nExtra[0]; ++ind[0]) {
102:         if (stagCoord->dof[0]) {
103:           const PetscReal offs[3] = {0.0,0.0,0.0};
104:           for (c=0; c<3; ++c) {
105:             arr[ind[2]][ind[1]][ind[0]][ibackdownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
106:           }
107:         }
108:         if (stagCoord->dof[1]) {
109:           const PetscReal offs[3] = {0.5,0.0,0.0};
110:           for (c=0; c<3; ++c) {
111:             arr[ind[2]][ind[1]][ind[0]][ibackdown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
112:           }
113:         }
114:         if (stagCoord->dof[1]) {
115:           const PetscReal offs[3] = {0.0,0.5,0.0};
116:           for (c=0; c<3; ++c) {
117:             arr[ind[2]][ind[1]][ind[0]][ibackleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
118:           }
119:         }
120:         if (stagCoord->dof[2]) {
121:           const PetscReal offs[3] = {0.5,0.5,0.0};
122:           for (c=0; c<3; ++c) {
123:             arr[ind[2]][ind[1]][ind[0]][iback + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
124:           }
125:         }
126:         if (stagCoord->dof[1]) {
127:           const PetscReal offs[3] = {0.0,0.0,0.5};
128:           for (c=0; c<3; ++c) {
129:             arr[ind[2]][ind[1]][ind[0]][idownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
130:           }
131:         }
132:         if (stagCoord->dof[2]) {
133:           const PetscReal offs[3] = {0.5,0.0,0.5};
134:           for (c=0; c<3; ++c) {
135:             arr[ind[2]][ind[1]][ind[0]][idown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
136:           }
137:         }
138:         if (stagCoord->dof[2]) {
139:           const PetscReal offs[3] = {0.0,0.5,0.5};
140:           for (c=0; c<3; ++c) {
141:             arr[ind[2]][ind[1]][ind[0]][ileft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
142:           }
143:         }
144:         if (stagCoord->dof[3]) {
145:           const PetscReal offs[3] = {0.5,0.5,0.5};
146:           for (c=0; c<3; ++c) {
147:             arr[ind[2]][ind[1]][ind[0]][ielement + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
148:           }
149:         }
150:       }
151:     }
152:   }
153:   DMStagVecRestoreArray(dmCoord,coordLocal,&arr);
154:   DMCreateGlobalVector(dmCoord,&coord);
155:   DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);
156:   DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);
157:   DMSetCoordinates(dm,coord);
158:   PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);
159:   DMRestoreLocalVector(dmCoord,&coordLocal);
160:   VecDestroy(&coord);
161:   return(0);
162: }

164: /* Helper functions used in DMSetUp_Stag() */
165: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM);
166: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM);
167: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM,PetscInt**);
168: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM,const PetscInt*);
169: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM,const PetscInt*);
170: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM);

172: PETSC_INTERN PetscErrorCode DMSetUp_Stag_3d(DM dm)
173: {
174:   PetscErrorCode  ierr;
175:   DM_Stag * const stag = (DM_Stag*)dm->data;
176:   PetscMPIInt     rank;
177:   PetscInt        i,j,d;
178:   PetscInt        *globalOffsets;
179:   const PetscInt  dim = 3;

182:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

184:   /* Rank grid sizes (populates stag->nRanks) */
185:   DMStagSetUpBuildRankGrid_3d(dm);

187:   /* Determine location of rank in grid */
188:     stag->rank[0] = rank % stag->nRanks[0];
189:     stag->rank[1] = rank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0];
190:     stag->rank[2] = rank / (stag->nRanks[0] * stag->nRanks[1]);
191:     for(d=0; d<dim; ++d) {
192:       stag->firstRank[d] = PetscNot(stag->rank[d]);
193:       stag->lastRank[d]  = (PetscBool)(stag->rank[d] == stag->nRanks[d]-1);
194:     }

196:   /* Determine locally owned region (if not already prescribed).
197:    Divide equally, giving lower ranks in each dimension and extra element if needbe.
198:    Note that this uses O(P) storage. If this ever becomes an issue, this could
199:    be refactored to not keep this data around.  */
200:   for (i=0; i<dim; ++i) {
201:     if (!stag->l[i]) {
202:       const PetscInt Ni = stag->N[i], nRanksi = stag->nRanks[i];
203:       PetscMalloc1(stag->nRanks[i],&stag->l[i]);
204:       for (j=0; j<stag->nRanks[i]; ++j){
205:         stag->l[i][j] = Ni/nRanksi + ((Ni % nRanksi) > j);
206:       }
207:     }
208:   }

210:   /* Retrieve local size in stag->n */
211:   for (i=0; i<dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
212:   if (PetscDefined(USE_DEBUG)) {
213:     for (i=0; i<dim; ++i) {
214:       PetscInt Ncheck,j;
215:       Ncheck = 0;
216:       for (j=0; j<stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
217:       if (Ncheck != stag->N[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local sizes in dimension %d don't add up. %d != %d\n",i,Ncheck,stag->N[i]);
218:     }
219:   }

221:   /* Compute starting elements */
222:   for (i=0; i<dim; ++i) {
223:     stag->start[i] = 0;
224:     for(j=0;j<stag->rank[i];++j){
225:       stag->start[i] += stag->l[i][j];
226:     }
227:   }

229:   /* Determine ranks of neighbors */
230:   DMStagSetUpBuildNeighbors_3d(dm);

232:   /* Define useful sizes */
233:   {
234:     PetscInt entriesPerEdge,entriesPerFace,entriesPerCorner,entriesPerElementRow,entriesPerElementLayer;
235:     PetscBool dummyEnd[3];
236:     for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
237:     stag->entriesPerElement = stag->dof[0] + 3*stag->dof[1] + 3*stag->dof[2] + stag->dof[3];
238:     entriesPerFace          = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
239:     entriesPerEdge          = stag->dof[0] + stag->dof[1];
240:     entriesPerCorner        = stag->dof[0];
241:     entriesPerElementRow    = stag->n[0]*stag->entriesPerElement
242:                               + (dummyEnd[0]                               ? entriesPerFace                       : 0);
243:     entriesPerElementLayer  = stag->n[1]*entriesPerElementRow
244:                               + (dummyEnd[1]                               ? stag->n[0] * entriesPerFace          : 0)
245:                               + (dummyEnd[1] && dummyEnd[0]                ? entriesPerEdge                       : 0);
246:     stag->entries           = stag->n[2]*entriesPerElementLayer
247:                               + (dummyEnd[2]                               ? stag->n[0]*stag->n[1]*entriesPerFace : 0)
248:                               + (dummyEnd[2] && dummyEnd[0]                ? stag->n[1]*entriesPerEdge            : 0)
249:                               + (dummyEnd[2] && dummyEnd[1]                ? stag->n[0]*entriesPerEdge            : 0)
250:                               + (dummyEnd[2] && dummyEnd[1] && dummyEnd[0] ? entriesPerCorner                     : 0);
251:   }

253:   /* Check that we will not overflow 32 bit indices (slightly overconservative) */
254:   if (!PetscDefined(USE_64BIT_INDICES)) {
255:     if (((PetscInt64) stag->n[0])*((PetscInt64) stag->n[1])*((PetscInt64) stag->n[2])*((PetscInt64) stag->entriesPerElement) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ4(PetscObjectComm((PetscObject)dm),PETSC_ERR_INT_OVERFLOW,"Mesh of %D x %D x %D with %D entries per (interior) element is likely too large for 32 bit indices",stag->n[0],stag->n[1],stag->n[2],stag->entriesPerElement);
256:   }

258:   /* Compute offsets for each rank into global vectors

260:     This again requires O(P) storage, which could be replaced with some global
261:     communication.
262:   */
263:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsets);

265:   for (d=0; d<dim; ++d) if (stag->boundaryType[d] != DM_BOUNDARY_NONE && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->boundaryType[d] != DM_BOUNDARY_GHOSTED) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported boundary type");

267:   /* Define ghosted/local sizes */
268:   for (d=0; d<dim; ++d) {
269:     switch (stag->boundaryType[d]) {
270:       case DM_BOUNDARY_NONE:
271:         /* Note: for a elements-only DMStag, the extra elements on the edges aren't necessary but we include them anyway */
272:         switch (stag->stencilType) {
273:           case DMSTAG_STENCIL_NONE : /* only the extra one on the right/top edges */
274:             stag->nGhost[d] = stag->n[d];
275:             stag->startGhost[d] = stag->start[d];
276:             if (stag->lastRank[d]) stag->nGhost[d] += 1;
277:             break;
278:           case DMSTAG_STENCIL_STAR : /* allocate the corners but don't use them */
279:           case DMSTAG_STENCIL_BOX :
280:             stag->nGhost[d] = stag->n[d];
281:             stag->startGhost[d] = stag->start[d];
282:             if (!stag->firstRank[d]) {
283:               stag->nGhost[d]     += stag->stencilWidth; /* add interior ghost elements */
284:               stag->startGhost[d] -= stag->stencilWidth;
285:             }
286:             if (!stag->lastRank[d]) {
287:               stag->nGhost[d] += stag->stencilWidth; /* add interior ghost elements */
288:             } else {
289:               stag->nGhost[d] += 1; /* one element on the boundary to complete blocking */
290:             }
291:             break;
292:           default :
293:             SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
294:         }
295:         break;
296:       case DM_BOUNDARY_GHOSTED:
297:         switch (stag->stencilType) {
298:           case DMSTAG_STENCIL_NONE :
299:             stag->startGhost[d] = stag->start[d];
300:             stag->nGhost[d]     = stag->n[d] + (stag->lastRank[d] ? 1 : 0);
301:             break;
302:           case DMSTAG_STENCIL_STAR :
303:           case DMSTAG_STENCIL_BOX :
304:             stag->startGhost[d] = stag->start[d] - stag->stencilWidth; /* This value may be negative */
305:             stag->nGhost[d]     = stag->n[d] + 2*stag->stencilWidth + (stag->lastRank[d] && stag->stencilWidth == 0 ? 1 : 0);
306:             break;
307:           default :
308:             SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
309:         }
310:         break;
311:       case DM_BOUNDARY_PERIODIC:
312:         switch (stag->stencilType) {
313:           case DMSTAG_STENCIL_NONE : /* only the extra one on the right/top edges */
314:             stag->nGhost[d] = stag->n[d];
315:             stag->startGhost[d] = stag->start[d];
316:             break;
317:           case DMSTAG_STENCIL_STAR :
318:           case DMSTAG_STENCIL_BOX :
319:             stag->nGhost[d] = stag->n[d] + 2*stag->stencilWidth;
320:             stag->startGhost[d] = stag->start[d] - stag->stencilWidth;
321:             break;
322:           default :
323:             SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
324:         }
325:         break;
326:       default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported boundary type in dimension %D",d);
327:     }
328:   }
329:   stag->entriesGhost = stag->nGhost[0]*stag->nGhost[1]*stag->nGhost[2]*stag->entriesPerElement;

331:   /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping */
332:   DMStagSetUpBuildScatter_3d(dm,globalOffsets);
333:   DMStagSetUpBuildL2G_3d(dm,globalOffsets);

335:   /* In special cases, create a dedicated injective local-to-global map */
336:   if ((stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) ||
337:       (stag->boundaryType[1] == DM_BOUNDARY_PERIODIC && stag->nRanks[1] == 1) ||
338:       (stag->boundaryType[2] == DM_BOUNDARY_PERIODIC && stag->nRanks[2] == 1)) {
339:     DMStagPopulateLocalToGlobalInjective(dm);
340:   }

342:   /* Free global offsets */
343:   PetscFree(globalOffsets);

345:   /* Precompute location offsets */
346:   DMStagComputeLocationOffsets_3d(dm);

348:   /* View from Options */
349:   DMViewFromOptions(dm,NULL,"-dm_view");

351:   return(0);
352: }

354: /* adapted from da3.c */
355: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM dm)
356: {
357:   PetscErrorCode  ierr;
358:   PetscMPIInt     rank,size;
359:   PetscInt        m,n,p,pm;
360:   DM_Stag * const stag = (DM_Stag*)dm->data;
361:   const PetscInt M = stag->N[0];
362:   const PetscInt N = stag->N[1];
363:   const PetscInt P = stag->N[2];

366:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
367:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

369:   m = stag->nRanks[0];
370:   n = stag->nRanks[1];
371:   p = stag->nRanks[2];

373:   if (m != PETSC_DECIDE) {
374:     if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
375:     else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
376:   }
377:   if (n != PETSC_DECIDE) {
378:     if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
379:     else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
380:   }
381:   if (p != PETSC_DECIDE) {
382:     if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
383:     else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
384:   }
385:   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);

387:   /* Partition the array among the processors */
388:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
389:     m = size/(n*p);
390:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
391:     n = size/(m*p);
392:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
393:     p = size/(m*n);
394:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
395:     /* try for squarish distribution */
396:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
397:     if (!m) m = 1;
398:     while (m > 0) {
399:       n = size/(m*p);
400:       if (m*n*p == size) break;
401:       m--;
402:     }
403:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
404:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
405:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
406:     /* try for squarish distribution */
407:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
408:     if (!m) m = 1;
409:     while (m > 0) {
410:       p = size/(m*n);
411:       if (m*n*p == size) break;
412:       m--;
413:     }
414:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
415:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
416:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
417:     /* try for squarish distribution */
418:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
419:     if (!n) n = 1;
420:     while (n > 0) {
421:       p = size/(m*n);
422:       if (m*n*p == size) break;
423:       n--;
424:     }
425:     if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
426:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
427:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
428:     /* try for squarish distribution */
429:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
430:     if (!n) n = 1;
431:     while (n > 0) {
432:       pm = size/n;
433:       if (n*pm == size) break;
434:       n--;
435:     }
436:     if (!n) n = 1;
437:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
438:     if (!m) m = 1;
439:     while (m > 0) {
440:       p = size/(m*n);
441:       if (m*n*p == size) break;
442:       m--;
443:     }
444:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
445:   } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

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

452:   stag->nRanks[0] = m;
453:   stag->nRanks[1] = n;
454:   stag->nRanks[2] = p;
455:   return(0);
456: }

458: /* Determine ranks of neighbors, using DMDA's convention

460:         n24 n25 n26
461:         n21 n22 n23
462:         n18 n19 n20 (Front, bigger z)

464:         n15 n16 n17
465:         n12     n14   ^ y
466:         n9  n10 n11   |
467:                       +--> x
468:         n6  n7  n8
469:         n3  n4  n5
470:         n0  n1  n2 (Back, smaller z) */
471: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM dm)
472: {
473:   PetscErrorCode  ierr;
474:   DM_Stag * const stag = (DM_Stag*)dm->data;
475:   PetscInt        d,i;
476:   PetscBool       per[3],first[3],last[3];
477:   PetscInt        neighborRank[27][3],r[3],n[3];
478:   const PetscInt  dim = 3;

481:   for (d=0; d<dim; ++d) if (stag->boundaryType[d] != DM_BOUNDARY_NONE && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->boundaryType[d] != DM_BOUNDARY_GHOSTED) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Neighbor determination not implemented for %s",DMBoundaryTypes[stag->boundaryType[d]]);

483:   /* Assemble some convenience variables */
484:   for (d=0; d<dim; ++d) {
485:     per[d]   = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
486:     first[d] = stag->firstRank[d];
487:     last[d]  = stag->lastRank[d];
488:     r[d]     = stag->rank[d];
489:     n[d]     = stag->nRanks[d];
490:   }

492:   /* First, compute the position in the rank grid for all neighbors */

494:   neighborRank[0][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down back  */
495:   neighborRank[0][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
496:   neighborRank[0][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

498:   neighborRank[1][0]  =                                      r[0]    ; /*       down back  */
499:   neighborRank[1][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
500:   neighborRank[1][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

502:   neighborRank[2][0]  = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right down back  */
503:   neighborRank[2][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
504:   neighborRank[2][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

506:   neighborRank[3][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left       back  */
507:   neighborRank[3][1]  =                                      r[1]    ;
508:   neighborRank[3][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

510:   neighborRank[4][0]  =                                      r[0]    ; /*            back  */
511:   neighborRank[4][1]  =                                      r[1]    ;
512:   neighborRank[4][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

514:   neighborRank[5][0]  = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right      back  */
515:   neighborRank[5][1]  =                                      r[1]    ;
516:   neighborRank[5][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

518:   neighborRank[6][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  up   back  */
519:   neighborRank[6][1]  = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
520:   neighborRank[6][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

522:   neighborRank[7][0]  =                                      r[0]    ; /*       up   back  */
523:   neighborRank[7][1]  = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
524:   neighborRank[7][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

526:   neighborRank[8][0]  = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right up   back  */
527:   neighborRank[8][1]  = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
528:   neighborRank[8][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

530:   neighborRank[9][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down       */
531:   neighborRank[9][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
532:   neighborRank[9][2]  =                                      r[2]    ;

534:   neighborRank[10][0] =                                      r[0]    ; /*       down       */
535:   neighborRank[10][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
536:   neighborRank[10][2] =                                      r[2]    ;

538:   neighborRank[11][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right down       */
539:   neighborRank[11][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
540:   neighborRank[11][2] =                                      r[2]    ;

542:   neighborRank[12][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left             */
543:   neighborRank[12][1] =                                      r[1]    ;
544:   neighborRank[12][2] =                                      r[2]    ;

546:   neighborRank[13][0] =                                      r[0]    ; /*                  */
547:   neighborRank[13][1] =                                      r[1]    ;
548:   neighborRank[13][2] =                                      r[2]    ;

550:   neighborRank[14][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right            */
551:   neighborRank[14][1] =                                      r[1]    ;
552:   neighborRank[14][2] =                                      r[2]    ;

554:   neighborRank[15][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  up         */
555:   neighborRank[15][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
556:   neighborRank[15][2] =                                      r[2]    ;

558:   neighborRank[16][0] =                                      r[0]    ; /*       up         */
559:   neighborRank[16][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
560:   neighborRank[16][2] =                                      r[2]    ;

562:   neighborRank[17][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right up         */
563:   neighborRank[17][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
564:   neighborRank[17][2] =                                      r[2]    ;

566:   neighborRank[18][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down front */
567:   neighborRank[18][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
568:   neighborRank[18][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

570:   neighborRank[19][0] =                                      r[0]    ; /*       down front */
571:   neighborRank[19][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
572:   neighborRank[19][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

574:   neighborRank[20][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right down front */
575:   neighborRank[20][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
576:   neighborRank[20][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

578:   neighborRank[21][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left       front */
579:   neighborRank[21][1] =                                      r[1]    ;
580:   neighborRank[21][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

582:   neighborRank[22][0] =                                      r[0]    ; /*            front */
583:   neighborRank[22][1] =                                      r[1]    ;
584:   neighborRank[22][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

586:   neighborRank[23][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right      front */
587:   neighborRank[23][1] =                                      r[1]    ;
588:   neighborRank[23][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

590:   neighborRank[24][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  up   front */
591:   neighborRank[24][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
592:   neighborRank[24][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

594:   neighborRank[25][0] =                                      r[0]    ; /*       up   front */
595:   neighborRank[25][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
596:   neighborRank[25][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

598:   neighborRank[26][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right up   front */
599:   neighborRank[26][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
600:   neighborRank[26][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

602:   /* Then, compute the rank of each in the linear ordering */
603:   PetscMalloc1(27,&stag->neighbors);
604:   for (i=0; i<27; ++i){
605:     if  (neighborRank[i][0] >= 0 && neighborRank[i][1] >=0 && neighborRank[i][2] >=0) {
606:       stag->neighbors[i] = neighborRank[i][0] + n[0]*neighborRank[i][1] + n[0]*n[1]*neighborRank[i][2];
607:     } else {
608:       stag->neighbors[i] = -1;
609:     }
610:   }

612:   return(0);
613: }

615: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM dm,PetscInt **pGlobalOffsets)
616: {
617:   PetscErrorCode        ierr;
618:   const DM_Stag * const stag = (DM_Stag*)dm->data;
619:   PetscInt              *globalOffsets;
620:   PetscInt              i,j,k,d,entriesPerEdge,entriesPerFace,count;
621:   PetscMPIInt           size;
622:   PetscBool             extra[3];

625:   for (d=0; d<3; ++d) extra[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC); /* Extra points in global rep */
626:   entriesPerFace               = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
627:   entriesPerEdge               = stag->dof[0] + stag->dof[1];
628:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
629:   PetscMalloc1(size,pGlobalOffsets);
630:   globalOffsets = *pGlobalOffsets;
631:   globalOffsets[0] = 0;
632:   count = 1; /* note the count is offset by 1 here. We add the size of the previous rank */
633:   for (k=0; k<stag->nRanks[2]-1; ++k){
634:     const PetscInt nnk = stag->l[2][k];
635:     for (j=0; j<stag->nRanks[1]-1; ++j) {
636:       const PetscInt nnj = stag->l[1][j];
637:       for (i=0; i<stag->nRanks[0]-1; ++i) {
638:         const PetscInt nni = stag->l[0][i];
639:         /* Interior : No right/top/front boundaries */
640:         globalOffsets[count] = globalOffsets[count-1] + nni*nnj*nnk*stag->entriesPerElement;
641:         ++count;
642:       }
643:       {
644:         /* Right boundary - extra faces */
645:         /* i = stag->nRanks[0]-1; */
646:         const PetscInt nni = stag->l[0][i];
647:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
648:                                + (extra[0] ? nnj*nnk*entriesPerFace : 0);
649:         ++count;
650:       }
651:     }
652:     {
653:       /* j = stag->nRanks[1]-1; */
654:       const PetscInt nnj = stag->l[1][j];
655:       for (i=0; i<stag->nRanks[0]-1; ++i) {
656:         const PetscInt nni = stag->l[0][i];
657:         /* Up boundary - extra faces */
658:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
659:                                + (extra[1] ? nni*nnk*entriesPerFace : 0);
660:         ++count;
661:       }
662:       {
663:         /* i = stag->nRanks[0]-1; */
664:         const PetscInt nni = stag->l[0][i];
665:         /* Up right boundary - 2x extra faces and extra edges */
666:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
667:                                + (extra[0]             ? nnj*nnk*entriesPerFace : 0)
668:                                + (extra[1]             ? nni*nnk*entriesPerFace : 0)
669:                                + (extra[0] && extra[1] ? nnk*entriesPerEdge     : 0);
670:         ++count;
671:       }
672:     }
673:   }
674:   {
675:     /* k = stag->nRanks[2]-1; */
676:     const PetscInt nnk = stag->l[2][k];
677:     for (j=0; j<stag->nRanks[1]-1; ++j) {
678:       const PetscInt nnj = stag->l[1][j];
679:       for (i=0; i<stag->nRanks[0]-1; ++i) {
680:         const PetscInt nni = stag->l[0][i];
681:         /* Front boundary - extra faces */
682:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
683:                                + (extra[2] ? nni*nnj*entriesPerFace : 0);
684:         ++count;
685:       }
686:       {
687:         /* i = stag->nRanks[0]-1; */
688:         const PetscInt nni = stag->l[0][i];
689:         /* Front right boundary - 2x extra faces and extra edges */
690:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
691:                                + (extra[0]             ? nnk*nnj*entriesPerFace : 0)
692:                                + (extra[2]             ? nni*nnj*entriesPerFace : 0)
693:                                + (extra[0] && extra[2] ? nnj*entriesPerEdge     : 0);
694:         ++count;
695:       }
696:     }
697:     {
698:       /* j = stag->nRanks[1]-1; */
699:       const PetscInt nnj = stag->l[1][j];
700:       for (i=0; i<stag->nRanks[0]-1; ++i) {
701:         const PetscInt nni = stag->l[0][i];
702:         /* Front up boundary - 2x extra faces and extra edges */
703:         globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
704:                                + (extra[1]             ? nnk*nni*entriesPerFace : 0)
705:                                + (extra[2]             ? nnj*nni*entriesPerFace : 0)
706:                                + (extra[1] && extra[2] ? nni*entriesPerEdge     : 0);
707:         ++count;
708:       }
709:       /* Don't need to compute entries in last element */
710:     }
711:   }

713:   return(0);
714: }

716: /* A helper function to reduce code duplication as we loop over various ranges.
717:    i,j,k refer to element numbers on the rank where an element lives in the global
718:    representation (without ghosts) to be offset by a global offset per rank.
719:    ig,jg,kg refer to element numbers in the local (this rank) ghosted numbering.
720:    Note that this function could easily be converted to a macro (it does not compute
721:    anything except loop indices and the values of idxGlobal and idxLocal).  */
722: static PetscErrorCode DMStagSetUpBuildScatterPopulateIdx_3d(DM_Stag *stag,PetscInt *count,PetscInt *idxLocal,PetscInt *idxGlobal,PetscInt entriesPerEdge, PetscInt entriesPerFace,PetscInt eprNeighbor,PetscInt eplNeighbor,PetscInt eprGhost,PetscInt eplGhost,PetscInt epFaceRow,PetscInt globalOffset,PetscInt startx,PetscInt starty,PetscInt startz,PetscInt startGhostx,PetscInt startGhosty,PetscInt startGhostz,PetscInt endGhostx,PetscInt endGhosty,PetscInt endGhostz, PetscBool extrax,PetscBool extray,PetscBool extraz)
723: {
724:   PetscInt ig,jg,kg,d,c;

727:   c = *count;
728:   for(kg = startGhostz; kg < endGhostz; ++kg) {
729:     const PetscInt k = kg - startGhostz + startz;
730:     for (jg = startGhosty; jg < endGhosty; ++jg) {
731:       const PetscInt j = jg - startGhosty + starty;
732:       for (ig = startGhostx; ig<endGhostx; ++ig) {
733:         const PetscInt i = ig - startGhostx + startx;
734:         for (d=0; d<stag->entriesPerElement; ++d,++c) {
735:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
736:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + d;
737:         }
738:       }
739:       if (extrax) {
740:         PetscInt dLocal;
741:         const PetscInt i = endGhostx - startGhostx + startx;
742:         ig = endGhostx;
743:         for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
744:           idxGlobal[c] = globalOffset + k *eplNeighbor + j*eprNeighbor + i *stag->entriesPerElement + d;
745:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost   + ig*stag->entriesPerElement + dLocal;
746:         }
747:         dLocal += stag->dof[1]; /* Skip back down edge */
748:         for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++c) { /* Back left edge */
749:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
750:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
751:         }
752:         dLocal += stag->dof[2]; /* Skip back face */
753:         for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++c) { /* Down left edge */
754:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
755:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
756:         }
757:         dLocal += stag->dof[2]; /* Skip down face */
758:         for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++c) { /* Left face */
759:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
760:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
761:         }
762:         /* Skip element */
763:       }
764:     }
765:     if (extray) {
766:       const PetscInt j = endGhosty - startGhosty + starty;
767:       jg = endGhosty;
768:       for (ig = startGhostx; ig<endGhostx; ++ig) {
769:         const PetscInt i = ig - startGhostx + startx;
770:         /* Vertex and Back down edge */
771:         PetscInt dLocal;
772:         for (d=0,dLocal=0; d<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++c) { /* Vertex */
773:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace          + d; /* Note face increment in  x */
774:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
775:         }
776:         /* Skip back left edge and back face */
777:         dLocal += stag->dof[1] + stag->dof[2];
778:         /* Down face and down left edge */
779:         for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++c) { /* Back left edge */
780:           idxGlobal[c] = globalOffset + k *eplNeighbor + j*eprNeighbor + i *entriesPerFace          + d; /* Note face increment in x */
781:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost   + ig*stag->entriesPerElement + dLocal;
782:         }
783:         /* Skip left face and element */
784:       }
785:       if (extrax) {
786:         PetscInt dLocal;
787:         const PetscInt i = endGhostx - startGhostx + startx;
788:         ig = endGhostx;
789:         for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
790:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace          + d; /* Note face increment in x */
791:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
792:         }
793:         dLocal += 2*stag->dof[1]+stag->dof[2]; /* Skip back down edge, back face, and back left edge */
794:         for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++c) { /* Down left edge */
795:           idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace          + d; /* Note face increment in x */
796:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost    + ig*stag->entriesPerElement + dLocal;
797:         }
798:         /* Skip remaining entries */
799:       }
800:     }
801:   }
802:   if (extraz) {
803:     const PetscInt k = endGhostz - startGhostz + startz;
804:     kg = endGhostz;
805:     for (jg = startGhosty; jg < endGhosty; ++jg) {
806:       const PetscInt j = jg - startGhosty + starty;
807:       for (ig = startGhostx; ig<endGhostx; ++ig) {
808:         const PetscInt i = ig - startGhostx + startx;
809:         for (d=0; d<entriesPerFace; ++d, ++c) {
810:           idxGlobal[c] = globalOffset + k*eplNeighbor  + j *epFaceRow + i *entriesPerFace          + d; /* Note face-based x and y increments */
811:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost  + ig*stag->entriesPerElement + d;
812:         }
813:       }
814:       if (extrax) {
815:         PetscInt dLocal;
816:         const PetscInt i = endGhostx - startGhostx + startx;
817:         ig = endGhostx;
818:         for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
819:           idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerFace          + d; /* Note face-based x and y increments */
820:           idxLocal[c]  =                kg*eplGhost   + jg*eprGhost  + ig*stag->entriesPerElement + dLocal;
821:         }
822:         dLocal += stag->dof[1]; /* Skip back down edge */
823:         for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++c) { /* Back left edge */
824:           idxGlobal[c] = globalOffset + k*eplNeighbor  + j *epFaceRow + i*entriesPerFace           + d; /* Note face-based x and y increments */
825:           idxLocal[c]  =                kg*eplGhost    + jg*eprGhost  + ig*stag->entriesPerElement + dLocal;
826:         }
827:         /* Skip the rest */
828:       }
829:     }
830:     if (extray) {
831:       const PetscInt j = endGhosty - startGhosty + starty;
832:       jg = endGhosty;
833:       for (ig = startGhostx; ig<endGhostx; ++ig) {
834:         const PetscInt i = ig - startGhostx + startx;
835:         for (d=0; d<entriesPerEdge; ++d,++c) {
836:           idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerEdge          + d; /* Note face-based y increment and edge-based x increment */
837:           idxLocal[c]  =                kg*eplGhost   + jg*eprGhost  + ig*stag->entriesPerElement + d;
838:         }
839:       }
840:       if (extrax) {
841:         const PetscInt i = endGhostx - startGhostx + startx;
842:         ig = endGhostx;
843:         for (d=0; d<stag->dof[0]; ++d, ++c) { /* Vertex (only) */
844:           idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerEdge          + d; /* Note face-based y increment and edge-based x increment */
845:           idxLocal[c]  =                kg*eplGhost   + jg*eprGhost  + ig*stag->entriesPerElement + d;
846:         }
847:       }
848:     }
849:   }
850:   *count = c;
851:   return(0);
852: }

854: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM dm,const PetscInt *globalOffsets)
855: {
857:   DM_Stag * const stag = (DM_Stag*)dm->data;
858:   PetscInt       d,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerCorner,entriesPerEdge,entriesPerFace,entriesToTransferTotal,count,eprGhost,eplGhost;
859:   PetscInt       *idxLocal,*idxGlobal;
860:   PetscMPIInt    rank;
861:   PetscInt       nNeighbors[27][3];
862:   PetscBool      star,nextToDummyEnd[3],dummyStart[3],dummyEnd[3];

865:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
866:   if (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth || stag->n[2] < stag->stencilWidth) {
867:     SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMStag 3d setup does not support local sizes (%D x %D x %D) smaller than the elementwise stencil width (%D)",stag->n[0],stag->n[1],stag->n[2],stag->stencilWidth);
868:   }

870:   /* Check stencil type */
871:   if (stag->stencilType != DMSTAG_STENCIL_NONE && stag->stencilType != DMSTAG_STENCIL_BOX && stag->stencilType != DMSTAG_STENCIL_STAR) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported stencil type %s",DMStagStencilTypes[stag->stencilType]);
872:   if (stag->stencilType == DMSTAG_STENCIL_NONE && stag->stencilWidth != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMStag 3d setup requires stencil width 0 with stencil type none");
873:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

875:   /* Compute numbers of elements on each neighbor */
876:   {
877:     PetscInt i;
878:     for (i=0; i<27; ++i) {
879:       const PetscInt neighborRank = stag->neighbors[i];
880:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
881:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
882:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
883:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
884:       } /* else leave uninitialized - error if accessed */
885:     }
886:   }

888:   /* These offsets should always be non-negative, and describe how many
889:      ghost elements exist at each boundary. These are not always equal to the stencil width,
890:      because we may have different numbers of ghost elements at the boundaries. In particular,
891:      in the non-periodic casewe always have at least one ghost (dummy) element at the right/top/front. */
892:   for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
893:   for (d=0; d<3; ++d) ghostOffsetEnd[d]   = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);

895:   /* Determine whether the ghost region includes dummies or not. This is currently
896:      equivalent to having a non-periodic boundary. If not, then
897:      ghostOffset{Start,End}[d] elements correspond to elements on the neighbor.
898:      If true, then
899:      - at the start, there are ghostOffsetStart[d] ghost elements
900:      - at the end, there is a layer of extra "physical" points inside a layer of
901:        ghostOffsetEnd[d] ghost elements
902:      Note that this computation should be updated if any boundary types besides
903:      NONE, GHOSTED, and PERIODIC are supported.  */
904:   for (d=0; d<3; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
905:   for (d=0; d<3; ++d) dummyEnd[d]   = (PetscBool)(stag->lastRank[d]  && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);

907:   /* Convenience variables */
908:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
909:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
910:   entriesPerCorner                    = stag->dof[0];
911:   for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
912:   eprGhost                            = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row   */
913:   eplGhost                            = stag->nGhost[1]*eprGhost;                /* epl = entries per (element) layer */

915:   /* Compute the number of local entries which correspond to any global entry */
916:   {
917:     PetscInt nNonDummyGhost[3];
918:     for (d=0; d<3; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
919:     if (star) {
920:       entriesToTransferTotal = (
921:           nNonDummyGhost[0] * stag->n[1]        * stag->n[2]
922:         + stag->n[0]        * nNonDummyGhost[1] * stag->n[2]
923:         + stag->n[0]        * stag->n[1]        * nNonDummyGhost[2]
924:         - 2 * (stag->n[0] * stag->n[1] * stag->n[2])
925:         ) * stag->entriesPerElement
926:         + (dummyEnd[0]                               ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace   : 0)
927:         + (dummyEnd[1]                               ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace   : 0)
928:         + (dummyEnd[2]                               ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace   : 0)
929:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
930:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
931:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
932:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
933:     } else {
934:       entriesToTransferTotal  =  nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
935:         + (dummyEnd[0]                               ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace   : 0)
936:         + (dummyEnd[1]                               ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace   : 0)
937:         + (dummyEnd[2]                               ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace   : 0)
938:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
939:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
940:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
941:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
942:     }
943:   }

945:   /* Allocate arrays to populate */
946:   PetscMalloc1(entriesToTransferTotal,&idxLocal);
947:   PetscMalloc1(entriesToTransferTotal,&idxGlobal);

949:   /* Counts into idxLocal/idxGlobal */
950:   count = 0;

952:   /*  Loop over each of the 27 neighbor, populating idxLocal and idxGlobal. These
953:       cases are principally distinguished by

955:       1. The loop bounds (i/ighost,j/jghost,k/kghost)
956:       2. the strides used in the loop body, which depend on whether the current
957:       rank and/or its neighbor are a non-periodic right/top/front boundary, which has additional
958:       points in the global representation.
959:       - If the neighboring rank is a right/top/front boundary,
960:       then eprNeighbor (entries per element row on the neighbor) and/or eplNeighbor (entries per element layer on the neighbor)
961:       are different.
962:       - If this rank is a non-periodic right/top/front boundary (checking entries of stag->lastRank),
963:       there is an extra loop over 1-3 boundary faces)

965:       Here, we do not include "dummy" dof (points in the local representation which
966:       do not correspond to any global dof). This, and the fact that we iterate over points in terms of
967:       increasing global ordering, are the main two differences from the construction of
968:       the Local-to-global map, which iterates over points in local order, and does include dummy points. */

970:   /* LEFT DOWN BACK */
971:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
972:     const PetscInt  neighbor     = 0;
973:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
974:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
975:     const PetscInt  epFaceRow    = -1;
976:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
977:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
978:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
979:     const PetscInt  startGhost0  = 0;
980:     const PetscInt  startGhost1  = 0;
981:     const PetscInt  startGhost2  = 0;
982:     const PetscInt  endGhost0    = ghostOffsetStart[0];
983:     const PetscInt  endGhost1    = ghostOffsetStart[1];
984:     const PetscInt  endGhost2    = ghostOffsetStart[2];
985:     const PetscBool extra0       = PETSC_FALSE;
986:     const PetscBool extra1       = PETSC_FALSE;
987:     const PetscBool extra2       = PETSC_FALSE;
988:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
989:   }

991:   /* DOWN BACK */
992:   if (!star && !dummyStart[1] && !dummyStart[2]) {
993:     const PetscInt  neighbor     = 1;
994:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
995:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
996:     const PetscInt  epFaceRow    = -1;
997:     const PetscInt  start0       = 0;
998:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
999:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1000:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1001:     const PetscInt  startGhost1  = 0;
1002:     const PetscInt  startGhost2  = 0;
1003:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1004:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1005:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1006:     const PetscBool extra0       = dummyEnd[0];
1007:     const PetscBool extra1       = PETSC_FALSE;
1008:     const PetscBool extra2       = PETSC_FALSE;
1009:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1010:   }

1012:   /* RIGHT DOWN BACK */
1013:   if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyStart[2]) {
1014:     const PetscInt  neighbor     = 2;
1015:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1016:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1017:     const PetscInt  epFaceRow    = -1;
1018:     const PetscInt  start0       = 0;
1019:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1020:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1021:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1022:     const PetscInt  startGhost1  = 0;
1023:     const PetscInt  startGhost2  = 0;
1024:     const PetscInt  endGhost0    = stag->nGhost[0];
1025:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1026:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1027:     const PetscBool extra0       = PETSC_FALSE;
1028:     const PetscBool extra1       = PETSC_FALSE;
1029:     const PetscBool extra2       = PETSC_FALSE;
1030:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1031:   }

1033:   /* LEFT BACK */
1034:   if (!star && !dummyStart[0] && !dummyStart[2]) {
1035:     const PetscInt  neighbor     = 3;
1036:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1037:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* May be a top boundary */
1038:     const PetscInt  epFaceRow    = -1;
1039:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1040:     const PetscInt  start1       = 0;
1041:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1042:     const PetscInt  startGhost0  = 0;
1043:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1044:     const PetscInt  startGhost2  = 0;
1045:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1046:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1047:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1048:     const PetscBool extra0       = PETSC_FALSE;
1049:     const PetscBool extra1       = dummyEnd[1];
1050:     const PetscBool extra2       = PETSC_FALSE;
1051:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1052:   }

1054:   /* BACK */
1055:   if (!dummyStart[2]) {
1056:     const PetscInt  neighbor     = 4;
1057:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may  be a right boundary */
1058:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We+neighbor may be a top boundary */
1059:     const PetscInt  epFaceRow    = -1;
1060:     const PetscInt  start0       = 0;
1061:     const PetscInt  start1       = 0;
1062:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1063:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1064:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1065:     const PetscInt  startGhost2  = 0;
1066:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1067:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1068:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1069:     const PetscBool extra0       = dummyEnd[0];
1070:     const PetscBool extra1       = dummyEnd[1];
1071:     const PetscBool extra2       = PETSC_FALSE;
1072:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1073:   }

1075:   /* RIGHT BACK */
1076:   if (!star && !dummyEnd[0] && !dummyStart[2]) {
1077:     const PetscInt  neighbor     = 5;
1078:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1079:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1080:     const PetscInt  epFaceRow    = -1;
1081:     const PetscInt  start0       = 0;
1082:     const PetscInt  start1       = 0;
1083:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1084:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1085:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1086:     const PetscInt  startGhost2  = 0;
1087:     const PetscInt  endGhost0    = stag->nGhost[0];
1088:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1089:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1090:     const PetscBool extra0       = PETSC_FALSE;
1091:     const PetscBool extra1       = dummyEnd[1];
1092:     const PetscBool extra2       = PETSC_FALSE;
1093:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1094:   }

1096:   /* LEFT UP BACK */
1097:   if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyStart[2]) {
1098:     const PetscInt  neighbor     = 6;
1099:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1100:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1101:     const PetscInt  epFaceRow    = -1;
1102:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1103:     const PetscInt  start1       = 0;
1104:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1105:     const PetscInt  startGhost0  = 0;
1106:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1107:     const PetscInt  startGhost2  = 0;
1108:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1109:     const PetscInt  endGhost1    = stag->nGhost[1];
1110:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1111:     const PetscBool extra0       = PETSC_FALSE;
1112:     const PetscBool extra1       = PETSC_FALSE;
1113:     const PetscBool extra2       = PETSC_FALSE;
1114:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1115:   }

1117:   /* UP BACK */
1118:   if (!star && !dummyEnd[1] && !dummyStart[2]) {
1119:     const PetscInt  neighbor     = 7;
1120:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1121:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1122:     const PetscInt  epFaceRow    = -1;
1123:     const PetscInt  start0       = 0;
1124:     const PetscInt  start1       = 0;
1125:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1126:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1127:     const PetscInt  startGhost1  = stag->nGhost[1]-ghostOffsetEnd[1];
1128:     const PetscInt  startGhost2  = 0;
1129:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1130:     const PetscInt  endGhost1    = stag->nGhost[1];
1131:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1132:     const PetscBool extra0       = dummyEnd[0];
1133:     const PetscBool extra1       = PETSC_FALSE;
1134:     const PetscBool extra2       = PETSC_FALSE;
1135:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1136:   }

1138:   /* RIGHT UP BACK */
1139:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyStart[2]) {
1140:     const PetscInt  neighbor     = 8;
1141:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1142:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1143:     const PetscInt  epFaceRow    = -1;
1144:     const PetscInt  start0       = 0;
1145:     const PetscInt  start1       = 0;
1146:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1147:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1148:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1149:     const PetscInt  startGhost2  = 0;
1150:     const PetscInt  endGhost0    = stag->nGhost[0];
1151:     const PetscInt  endGhost1    = stag->nGhost[1];
1152:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1153:     const PetscBool extra0       = PETSC_FALSE;
1154:     const PetscBool extra1       = PETSC_FALSE;
1155:     const PetscBool extra2       = PETSC_FALSE;
1156:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1157:   }

1159:   /* LEFT DOWN */
1160:   if (!star && !dummyStart[0] && !dummyStart[1]) {
1161:     const PetscInt  neighbor     = 9;
1162:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1163:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1164:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1165:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1166:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1167:     const PetscInt  start2       = 0;
1168:     const PetscInt  startGhost0  = 0;
1169:     const PetscInt  startGhost1  = 0;
1170:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1171:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1172:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1173:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1174:     const PetscBool extra0       = PETSC_FALSE;
1175:     const PetscBool extra1       = PETSC_FALSE;
1176:     const PetscBool extra2       = dummyEnd[2];
1177:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1178:   }

1180:   /* DOWN */
1181:   if (!dummyStart[1]) {
1182:     const PetscInt  neighbor     = 10;
1183:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1184:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1185:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1186:     const PetscInt  start0       = 0;
1187:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1188:     const PetscInt  start2       = 0;
1189:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1190:     const PetscInt  startGhost1  = 0;
1191:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1192:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1193:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1194:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1195:     const PetscBool extra0       = dummyEnd[0];
1196:     const PetscBool extra1       = PETSC_FALSE;
1197:     const PetscBool extra2       = dummyEnd[2];
1198:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1199:   }

1201:   /* RIGHT DOWN */
1202:   if (!star && !dummyEnd[0] && !dummyStart[1]) {
1203:     const PetscInt  neighbor     = 11;
1204:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1205:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1206:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1207:     const PetscInt  start0       = 0;
1208:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1209:     const PetscInt  start2       = 0;
1210:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1211:     const PetscInt  startGhost1  = 0;
1212:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1213:     const PetscInt  endGhost0    = stag->nGhost[0];
1214:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1215:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1216:     const PetscBool extra0       = PETSC_FALSE;
1217:     const PetscBool extra1       = PETSC_FALSE;
1218:     const PetscBool extra2       = dummyEnd[2];
1219:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1220:   }

1222:   /* LEFT */
1223:   if (!dummyStart[0]) {
1224:     const PetscInt  neighbor     = 12;
1225:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1226:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1227:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0];
1228:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1229:     const PetscInt  start1       = 0;
1230:     const PetscInt  start2       = 0;
1231:     const PetscInt  startGhost0  = 0;
1232:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1233:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1234:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1235:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1236:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1237:     const PetscBool extra0       = PETSC_FALSE;
1238:     const PetscBool extra1       = dummyEnd[1];
1239:     const PetscBool extra2       = dummyEnd[2];
1240:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1241:   }

1243:   /* (HERE) */
1244:   {
1245:     const PetscInt  neighbor     = 13;
1246:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1247:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1248:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
1249:     const PetscInt  start0       = 0;
1250:     const PetscInt  start1       = 0;
1251:     const PetscInt  start2       = 0;
1252:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1253:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1254:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1255:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1256:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1257:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1258:     const PetscBool extra0       = dummyEnd[0];
1259:     const PetscBool extra1       = dummyEnd[1];
1260:     const PetscBool extra2       = dummyEnd[2];
1261:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1262:   }

1264:   /* RIGHT */
1265:   if (!dummyEnd[0]) {
1266:     const PetscInt  neighbor     = 14;
1267:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1268:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1269:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1270:     const PetscInt  start0       = 0;
1271:     const PetscInt  start1       = 0;
1272:     const PetscInt  start2       = 0;
1273:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1274:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1275:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1276:     const PetscInt  endGhost0    = stag->nGhost[0];
1277:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1278:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1279:     const PetscBool extra0       = PETSC_FALSE;
1280:     const PetscBool extra1       = dummyEnd[1];
1281:     const PetscBool extra2       = dummyEnd[2];
1282:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1283:   }

1285:   /* LEFT UP */
1286:   if (!star && !dummyStart[0] && !dummyEnd[1]) {
1287:     const PetscInt  neighbor     = 15;
1288:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1289:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1290:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0];
1291:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1292:     const PetscInt  start1       = 0;
1293:     const PetscInt  start2       = 0;
1294:     const PetscInt  startGhost0  = 0;
1295:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1296:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1297:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1298:     const PetscInt  endGhost1    = stag->nGhost[1];
1299:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1300:     const PetscBool extra0       = PETSC_FALSE;
1301:     const PetscBool extra1       = PETSC_FALSE;
1302:     const PetscBool extra2       = dummyEnd[2];
1303:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1304:   }

1306:   /* UP */
1307:   if (!dummyEnd[1]) {
1308:     const PetscInt  neighbor     = 16;
1309:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1310:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1311:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1312:     const PetscInt  start0       = 0;
1313:     const PetscInt  start1       = 0;
1314:     const PetscInt  start2       = 0;
1315:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1316:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1317:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1318:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1319:     const PetscInt  endGhost1    = stag->nGhost[1];
1320:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1321:     const PetscBool extra0       = dummyEnd[0];
1322:     const PetscBool extra1       = PETSC_FALSE;
1323:     const PetscBool extra2       = dummyEnd[2];
1324:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1325:   }

1327:   /* RIGHT UP */
1328:   if (!star && !dummyEnd[0] && !dummyEnd[1]) {
1329:     const PetscInt  neighbor     = 17;
1330:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1331:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1332:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1333:     const PetscInt  start0       = 0;
1334:     const PetscInt  start1       = 0;
1335:     const PetscInt  start2       = 0;
1336:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1337:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1338:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1339:     const PetscInt  endGhost0    = stag->nGhost[0];
1340:     const PetscInt  endGhost1    = stag->nGhost[1];
1341:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1342:     const PetscBool extra0       = PETSC_FALSE;
1343:     const PetscBool extra1       = PETSC_FALSE;
1344:     const PetscBool extra2       = dummyEnd[2];
1345:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1346:   }

1348:   /* LEFT DOWN FRONT */
1349:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1350:     const PetscInt  neighbor     = 18;
1351:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1352:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1353:     const PetscInt  epFaceRow    = -1;
1354:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1355:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1356:     const PetscInt  start2       = 0;
1357:     const PetscInt  startGhost0  = 0;
1358:     const PetscInt  startGhost1  = 0;
1359:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1360:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1361:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1362:     const PetscInt  endGhost2    = stag->nGhost[2];
1363:     const PetscBool extra0       = PETSC_FALSE;
1364:     const PetscBool extra1       = PETSC_FALSE;
1365:     const PetscBool extra2       = PETSC_FALSE;
1366:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1367:   }

1369:   /* DOWN FRONT */
1370:   if (!star && !dummyStart[1] && !dummyEnd[2]) {
1371:     const PetscInt  neighbor     = 19;
1372:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1373:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1374:     const PetscInt  epFaceRow    = -1;
1375:     const PetscInt  start0       = 0;
1376:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1377:     const PetscInt  start2       = 0;
1378:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1379:     const PetscInt  startGhost1  = 0;
1380:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1381:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1382:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1383:     const PetscInt  endGhost2    = stag->nGhost[2];
1384:     const PetscBool extra0       = dummyEnd[0];
1385:     const PetscBool extra1       = PETSC_FALSE;
1386:     const PetscBool extra2       = PETSC_FALSE;
1387:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1388:   }

1390:   /* RIGHT DOWN FRONT */
1391:   if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyEnd[2]) {
1392:     const PetscInt  neighbor     = 20;
1393:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1394:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1395:     const PetscInt  epFaceRow    = -1;
1396:     const PetscInt  start0       = 0;
1397:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1398:     const PetscInt  start2       = 0;
1399:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1400:     const PetscInt  startGhost1  = 0;
1401:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1402:     const PetscInt  endGhost0    = stag->nGhost[0];
1403:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1404:     const PetscInt  endGhost2    = stag->nGhost[2];
1405:     const PetscBool extra0       = PETSC_FALSE;
1406:     const PetscBool extra1       = PETSC_FALSE;
1407:     const PetscBool extra2       = PETSC_FALSE;
1408:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1409:   }

1411:   /* LEFT FRONT */
1412:   if (!star && !dummyStart[0] && !dummyEnd[2]) {
1413:     const PetscInt  neighbor     = 21;
1414:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1415:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1416:     const PetscInt  epFaceRow    = -1;
1417:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1418:     const PetscInt  start1       = 0;
1419:     const PetscInt  start2       = 0;
1420:     const PetscInt  startGhost0  = 0;
1421:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1422:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1423:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1424:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1425:     const PetscInt  endGhost2    = stag->nGhost[2];
1426:     const PetscBool extra0       = PETSC_FALSE;
1427:     const PetscBool extra1       = dummyEnd[1];
1428:     const PetscBool extra2       = PETSC_FALSE;
1429:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1430:   }

1432:   /* FRONT */
1433:   if (!dummyEnd[2]) {
1434:     const PetscInt  neighbor     = 22;
1435:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1436:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1437:     const PetscInt  epFaceRow    = -1;
1438:     const PetscInt  start0       = 0;
1439:     const PetscInt  start1       = 0;
1440:     const PetscInt  start2       = 0;
1441:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1442:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1443:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1444:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1445:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1446:     const PetscInt  endGhost2    = stag->nGhost[2];
1447:     const PetscBool extra0       = dummyEnd[0];
1448:     const PetscBool extra1       = dummyEnd[1];
1449:     const PetscBool extra2       = PETSC_FALSE;
1450:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1451:   }

1453:   /* RIGHT FRONT */
1454:   if (!star && !dummyEnd[0] && !dummyEnd[2]) {
1455:     const PetscInt  neighbor     = 23;
1456:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1457:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1458:     const PetscInt  epFaceRow    = -1;
1459:     const PetscInt  start0       = 0;
1460:     const PetscInt  start1       = 0;
1461:     const PetscInt  start2       = 0;
1462:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1463:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1464:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1465:     const PetscInt  endGhost0    = stag->nGhost[0];
1466:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1467:     const PetscInt  endGhost2    = stag->nGhost[2];
1468:     const PetscBool extra0       = PETSC_FALSE;
1469:     const PetscBool extra1       = dummyEnd[1];
1470:     const PetscBool extra2       = PETSC_FALSE;
1471:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1472:   }

1474:   /* LEFT UP FRONT */
1475:   if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyEnd[2]) {
1476:     const PetscInt  neighbor     = 24;
1477:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1478:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1479:     const PetscInt  epFaceRow    = -1;
1480:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1481:     const PetscInt  start1       = 0;
1482:     const PetscInt  start2       = 0;
1483:     const PetscInt  startGhost0  = 0;
1484:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1485:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1486:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1487:     const PetscInt  endGhost1    = stag->nGhost[1];
1488:     const PetscInt  endGhost2    = stag->nGhost[2];
1489:     const PetscBool extra0       = PETSC_FALSE;
1490:     const PetscBool extra1       = PETSC_FALSE;
1491:     const PetscBool extra2       = PETSC_FALSE;
1492:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1493:   }

1495:   /* UP FRONT */
1496:   if (!star && !dummyEnd[1] && !dummyEnd[2]) {
1497:     const PetscInt  neighbor     = 25;
1498:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1499:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1500:     const PetscInt  epFaceRow    = -1;
1501:     const PetscInt  start0       = 0;
1502:     const PetscInt  start1       = 0;
1503:     const PetscInt  start2       = 0;
1504:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1505:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1506:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1507:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1508:     const PetscInt  endGhost1    = stag->nGhost[1];
1509:     const PetscInt  endGhost2    = stag->nGhost[2];
1510:     const PetscBool extra0       = dummyEnd[0];
1511:     const PetscBool extra1       = PETSC_FALSE;
1512:     const PetscBool extra2       = PETSC_FALSE;
1513:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1514:   }

1516:   /* RIGHT UP FRONT */
1517:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyEnd[2]) {
1518:     const PetscInt  neighbor     = 26;
1519:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1520:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1521:     const PetscInt  epFaceRow    = -1;
1522:     const PetscInt  start0       = 0;
1523:     const PetscInt  start1       = 0;
1524:     const PetscInt  start2       = 0;
1525:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1526:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1527:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1528:     const PetscInt  endGhost0    = stag->nGhost[0];
1529:     const PetscInt  endGhost1    = stag->nGhost[1];
1530:     const PetscInt  endGhost2    = stag->nGhost[2];
1531:     const PetscBool extra0       = PETSC_FALSE;
1532:     const PetscBool extra1       = PETSC_FALSE;
1533:     const PetscBool extra2       = PETSC_FALSE;
1534:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1535:   }

1537:   if (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);

1539:   /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
1540:   {
1541:     Vec local,global;
1542:     IS isLocal,isGlobal;
1543:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
1544:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
1545:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
1546:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
1547:     VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
1548:     VecDestroy(&global);
1549:     VecDestroy(&local);
1550:     ISDestroy(&isLocal);  /* frees idxLocal */
1551:     ISDestroy(&isGlobal); /* free idxGlobal */
1552:   }

1554:   return(0);
1555: }

1557: /* Note: this function assumes that DMBoundary types of none, ghosted, and periodic are the only ones of interest.
1558: Adding support for others should be done very carefully.  */
1559: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM dm,const PetscInt *globalOffsets)
1560: {
1561:   PetscErrorCode        ierr;
1562:   const DM_Stag * const stag = (DM_Stag*)dm->data;
1563:   PetscInt              *idxGlobalAll;
1564:   PetscInt              d,count,ighost,jghost,kghost,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerFace,entriesPerEdge;
1565:   PetscInt              nNeighbors[27][3],eprNeighbor[27],eplNeighbor[27],globalOffsetNeighbor[27];
1566:   PetscBool             nextToDummyEnd[3],dummyStart[3],dummyEnd[3],star;


1570:   /* Check stencil type */
1571:   if (stag->stencilType != DMSTAG_STENCIL_NONE && stag->stencilType != DMSTAG_STENCIL_BOX && stag->stencilType != DMSTAG_STENCIL_STAR) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported stencil type %s",DMStagStencilTypes[stag->stencilType]);
1572:   if (stag->stencilType == DMSTAG_STENCIL_NONE && stag->stencilWidth != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMStag 3d setup requires stencil width 0 with stencil type none");
1573:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

1575:   /* Convenience variables */
1576:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
1577:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
1578:   for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);

1580:   /* Ghost offsets (may not be the ghost width, since we always have at least one ghost element on the right/top/front) */
1581:   for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
1582:   for (d=0; d<3; ++d) ghostOffsetEnd[d]   = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);

1584:   /* Whether the ghost region includes dummies. Currently equivalent to being a non-periodic boundary. */
1585:   for (d=0; d<3; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1586:   for (d=0; d<3; ++d) dummyEnd[d]   = (PetscBool)(stag->lastRank[d]  && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);

1588:   /* Compute numbers of elements on each neighbor  and offset*/
1589:   {
1590:     PetscInt i;
1591:     for (i=0; i<27; ++i) {
1592:       const PetscInt neighborRank = stag->neighbors[i];
1593:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
1594:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
1595:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
1596:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
1597:         globalOffsetNeighbor[i] = globalOffsets[stag->neighbors[i]];
1598:       } /* else leave uninitialized - error if accessed */
1599:     }
1600:   }

1602:   /* Precompute elements per row and layer on neighbor (zero unused) */
1603:   PetscMemzero(eprNeighbor,sizeof(eprNeighbor));
1604:   PetscMemzero(eplNeighbor,sizeof(eplNeighbor));
1605:   if (stag->neighbors[0] >= 0) {
1606:     eprNeighbor[0] = stag->entriesPerElement * nNeighbors[0][0];
1607:     eplNeighbor[0] = eprNeighbor[0] * nNeighbors[0][1];
1608:   }
1609:   if (stag->neighbors[1] >= 0) {
1610:     eprNeighbor[1] = stag->entriesPerElement * nNeighbors[1][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1611:     eplNeighbor[1] = eprNeighbor[1] * nNeighbors[1][1];
1612:   }
1613:   if (stag->neighbors[2] >= 0) {
1614:     eprNeighbor[2] = stag->entriesPerElement * nNeighbors[2][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1615:     eplNeighbor[2] = eprNeighbor[2] * nNeighbors[2][1];
1616:   }
1617:   if (stag->neighbors[3] >= 0) {
1618:     eprNeighbor[3] = stag->entriesPerElement * nNeighbors[3][0];
1619:     eplNeighbor[3] = eprNeighbor[3] * nNeighbors[3][1] + (dummyEnd[1] ? nNeighbors[3][0] * entriesPerFace : 0);  /* May be a top boundary */
1620:   }
1621:   if (stag->neighbors[4] >= 0) {
1622:     eprNeighbor[4] = stag->entriesPerElement * nNeighbors[4][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may  be a right boundary */
1623:     eplNeighbor[4] = eprNeighbor[4] * nNeighbors[4][1] + (dummyEnd[1] ? nNeighbors[4][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We+neighbor may be a top boundary */
1624:   }
1625:   if (stag->neighbors[5] >= 0) {
1626:     eprNeighbor[5] = stag->entriesPerElement * nNeighbors[5][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1627:     eplNeighbor[5] = eprNeighbor[5]  * nNeighbors[5][1] + (dummyEnd[1] ? nNeighbors[5][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1628:   }
1629:   if (stag->neighbors[6] >= 0) {
1630:     eprNeighbor[6] = stag->entriesPerElement * nNeighbors[6][0];
1631:     eplNeighbor[6] = eprNeighbor[6] * nNeighbors[6][1] + (nextToDummyEnd[1] ? nNeighbors[6][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1632:   }
1633:   if (stag->neighbors[7] >= 0) {
1634:     eprNeighbor[7] = stag->entriesPerElement * nNeighbors[7][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1635:     eplNeighbor[7] = eprNeighbor[7] * nNeighbors[7][1] + (nextToDummyEnd[1] ? nNeighbors[7][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1636:   }
1637:   if (stag->neighbors[8] >= 0) {
1638:     eprNeighbor[8] = stag->entriesPerElement * nNeighbors[8][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1639:     eplNeighbor[8] = eprNeighbor[8] * nNeighbors[8][1] + (nextToDummyEnd[1] ? nNeighbors[8][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1640:   }
1641:   if (stag->neighbors[9] >= 0) {
1642:     eprNeighbor[9] = stag->entriesPerElement * nNeighbors[9][0];
1643:     eplNeighbor[9] = eprNeighbor[9] * nNeighbors[9][1];
1644:   }
1645:   if (stag->neighbors[10] >= 0) {
1646:     eprNeighbor[10] = stag->entriesPerElement * nNeighbors[10][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1647:     eplNeighbor[10] = eprNeighbor[10] * nNeighbors[10][1];
1648:   }
1649:   if (stag->neighbors[11] >= 0) {
1650:     eprNeighbor[11] = stag->entriesPerElement * nNeighbors[11][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1651:     eplNeighbor[11] = eprNeighbor[11] * nNeighbors[11][1];
1652:   }
1653:   if (stag->neighbors[12] >= 0) {
1654:     eprNeighbor[12] = stag->entriesPerElement * nNeighbors[12][0];
1655:     eplNeighbor[12] = eprNeighbor[12] * nNeighbors[12][1] + (dummyEnd[1] ? nNeighbors[12][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1656:   }
1657:   if (stag->neighbors[13] >= 0) {
1658:     eprNeighbor[13] = stag->entriesPerElement * nNeighbors[13][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1659:     eplNeighbor[13] = eprNeighbor[13] * nNeighbors[13][1] + (dummyEnd[1] ? nNeighbors[13][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1660:   }
1661:   if (stag->neighbors[14] >= 0) {
1662:     eprNeighbor[14] = stag->entriesPerElement * nNeighbors[14][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1663:     eplNeighbor[14] = eprNeighbor[14] * nNeighbors[14][1] + (dummyEnd[1] ? nNeighbors[14][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1664:   }
1665:   if (stag->neighbors[15] >= 0) {
1666:     eprNeighbor[15] = stag->entriesPerElement * nNeighbors[15][0];
1667:     eplNeighbor[15] = eprNeighbor[15] * nNeighbors[15][1] + (nextToDummyEnd[1] ? nNeighbors[15][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1668:   }
1669:   if (stag->neighbors[16] >= 0) {
1670:     eprNeighbor[16] = stag->entriesPerElement * nNeighbors[16][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1671:     eplNeighbor[16] = eprNeighbor[16] * nNeighbors[16][1] + (nextToDummyEnd[1] ? nNeighbors[16][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1672:   }
1673:   if (stag->neighbors[17] >= 0) {
1674:     eprNeighbor[17] = stag->entriesPerElement * nNeighbors[17][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1675:     eplNeighbor[17] = eprNeighbor[17] * nNeighbors[17][1] + (nextToDummyEnd[1] ? nNeighbors[17][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1676:   }
1677:   if (stag->neighbors[18] >= 0) {
1678:     eprNeighbor[18] = stag->entriesPerElement * nNeighbors[18][0];
1679:     eplNeighbor[18] = eprNeighbor[18] * nNeighbors[18][1];
1680:   }
1681:   if (stag->neighbors[19] >= 0) {
1682:     eprNeighbor[19] = stag->entriesPerElement * nNeighbors[19][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1683:     eplNeighbor[19] = eprNeighbor[19] * nNeighbors[19][1];
1684:   }
1685:   if (stag->neighbors[20] >= 0) {
1686:     eprNeighbor[20] = stag->entriesPerElement * nNeighbors[20][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1687:     eplNeighbor[20] = eprNeighbor[20] * nNeighbors[20][1];
1688:   }
1689:   if (stag->neighbors[21] >= 0) {
1690:     eprNeighbor[21] = stag->entriesPerElement * nNeighbors[21][0];
1691:     eplNeighbor[21] = eprNeighbor[21] * nNeighbors[21][1] + (dummyEnd[1] ? nNeighbors[21][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1692:   }
1693:   if (stag->neighbors[22] >= 0) {
1694:     eprNeighbor[22] = stag->entriesPerElement * nNeighbors[22][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1695:     eplNeighbor[22] = eprNeighbor[22] * nNeighbors[22][1] + (dummyEnd[1] ? nNeighbors[22][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1696:   }
1697:   if (stag->neighbors[23] >= 0) {
1698:     eprNeighbor[23] = stag->entriesPerElement * nNeighbors[23][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1699:     eplNeighbor[23] = eprNeighbor[23] * nNeighbors[23][1] + (dummyEnd[1] ? nNeighbors[23][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* We and neighbor may be a top boundary */
1700:   }
1701:   if (stag->neighbors[24] >= 0) {
1702:     eprNeighbor[24] = stag->entriesPerElement * nNeighbors[24][0];
1703:     eplNeighbor[24] = eprNeighbor[24] * nNeighbors[24][1] + (nextToDummyEnd[1] ? nNeighbors[24][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1704:   }
1705:   if (stag->neighbors[25] >= 0) {
1706:     eprNeighbor[25] = stag->entriesPerElement * nNeighbors[25][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1707:     eplNeighbor[25] = eprNeighbor[25] * nNeighbors[25][1] + (nextToDummyEnd[1] ? nNeighbors[25][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1708:   }
1709:   if (stag->neighbors[26] >= 0) {
1710:     eprNeighbor[26] = stag->entriesPerElement * nNeighbors[26][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1711:     eplNeighbor[26] = eprNeighbor[26] * nNeighbors[26][1] + (nextToDummyEnd[1] ? nNeighbors[26][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1712:   }

1714:   /* Populate idxGlobalAll */
1715:   PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
1716:   count = 0;

1718:   /* Loop over layers 1/3 : Back */
1719:   if (!dummyStart[2]) {
1720:     for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
1721:       const PetscInt k = nNeighbors[4][2] - ghostOffsetStart[2] + kghost; /* Note: this is the same value for all neighbors in this layer (use neighbor 4 which will always exist if we're lookng at this layer) */

1723:       /* Loop over rows 1/3: Down Back*/
1724:       if (!star && !dummyStart[1]) {
1725:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1726:           const PetscInt j = nNeighbors[1][1] - ghostOffsetStart[1] + jghost; /* Note: this is the same value for all neighbors in this row (use neighbor 1, down back)*/

1728:           /* Loop over columns 1/3: Left Back Down*/
1729:           if (!dummyStart[0]) {
1730:             const PetscInt neighbor = 0;
1731:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1732:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1733:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1734:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1735:               }
1736:             }
1737:           } else {
1738:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1739:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1740:                 idxGlobalAll[count] = -1;
1741:               }
1742:             }
1743:           }

1745:           /* Loop over columns 2/3: (Middle) Down Back */
1746:           {
1747:             const PetscInt neighbor = 1;
1748:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1749:               const PetscInt i = ighost - ghostOffsetStart[0];
1750:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1751:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1752:               }
1753:             }
1754:           }

1756:           /* Loop over columns 3/3: Right Down Back */
1757:           if (!dummyEnd[0]) {
1758:             const PetscInt neighbor = 2;
1759:             PetscInt       i;
1760:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1761:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1762:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1763:               }
1764:             }
1765:           } else {
1766:             /* Partial dummy entries on (Middle) Down Back neighbor */
1767:             const PetscInt neighbor          = 1;
1768:             PetscInt       i,dLocal;
1769:             i = stag->n[0];
1770:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1771:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1772:             }
1773:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1774:               idxGlobalAll[count] = -1;
1775:             }
1776:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1777:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1778:             }
1779:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1780:               idxGlobalAll[count] = -1;
1781:             }
1782:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1783:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1784:             }
1785:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1786:               idxGlobalAll[count] = -1;
1787:             }
1788:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1789:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1790:             }
1791:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1792:               idxGlobalAll[count] = -1;
1793:             }
1794:             ++i;
1795:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1796:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1797:                 idxGlobalAll[count] = -1;
1798:               }
1799:             }
1800:           }
1801:         }
1802:       } else {
1803:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1804:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
1805:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
1806:               idxGlobalAll[count] = -1;
1807:             }
1808:           }
1809:         }
1810:       }

1812:       /* Loop over rows 2/3: (Middle) Back */
1813:       {
1814:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
1815:           const PetscInt j = jghost - ghostOffsetStart[1];

1817:           /* Loop over columns 1/3: Left (Middle) Back */
1818:           if (!star && !dummyStart[0]) {
1819:             const PetscInt neighbor = 3;
1820:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1821:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1822:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1823:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1824:               }
1825:             }
1826:           } else {
1827:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1828:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1829:                 idxGlobalAll[count] = -1;
1830:               }
1831:             }
1832:           }

1834:           /* Loop over columns 2/3: (Middle) (Middle) Back */
1835:           {
1836:             const PetscInt neighbor = 4;
1837:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1838:               const PetscInt i = ighost - ghostOffsetStart[0];
1839:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1840:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1841:               }
1842:             }
1843:           }

1845:           /* Loop over columns 3/3: Right (Middle) Back */
1846:           if (!star && !dummyEnd[0]) {
1847:             const PetscInt neighbor = 5;
1848:             PetscInt       i;
1849:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1850:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1851:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1852:               }
1853:             }
1854:           } else if (dummyEnd[0]) {
1855:             /* Partial dummy entries on (Middle) (Middle) Back rank */
1856:             const PetscInt neighbor = 4;
1857:             PetscInt       i,dLocal;
1858:             i = stag->n[0];
1859:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1860:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1861:             }
1862:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1863:               idxGlobalAll[count] = -1;
1864:             }
1865:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1866:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1867:             }
1868:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1869:               idxGlobalAll[count] = -1;
1870:             }
1871:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1872:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1873:             }
1874:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1875:               idxGlobalAll[count] = -1;
1876:             }
1877:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1878:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1879:             }
1880:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1881:               idxGlobalAll[count] = -1;
1882:             }
1883:             ++i;
1884:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1885:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1886:                 idxGlobalAll[count] = -1;
1887:               }
1888:             }
1889:           } else {
1890:             /* Right (Middle) Back dummies */
1891:             PetscInt       i;
1892:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1893:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1894:                 idxGlobalAll[count] = -1;
1895:               }
1896:             }
1897:           }
1898:         }
1899:       }

1901:       /* Loop over rows 3/3: Up Back */
1902:       if (!star && !dummyEnd[1]) {
1903:         PetscInt j;
1904:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
1905:           /* Loop over columns 1/3: Left Up Back*/
1906:           if (!dummyStart[0]) {
1907:             const PetscInt neighbor = 6;
1908:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1909:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1910:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1911:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1912:               }
1913:             }
1914:           } else {
1915:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1916:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1917:                 idxGlobalAll[count] = -1;
1918:               }
1919:             }
1920:           }

1922:           /* Loop over columns 2/3: (Middle) Up Back*/
1923:           {
1924:             const PetscInt neighbor = 7;
1925:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1926:               const PetscInt i = ighost - ghostOffsetStart[0];
1927:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1928:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1929:               }
1930:             }
1931:           }

1933:           /* Loop over columns 3/3: Right Up Back*/
1934:           if (!dummyEnd[0]) {
1935:             const PetscInt neighbor = 8;
1936:             PetscInt       i;
1937:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1938:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1939:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1940:               }
1941:             }
1942:           } else {
1943:             /* Partial dummies on (Middle) Up Back */
1944:             const PetscInt neighbor = 7;
1945:             PetscInt       i,dLocal;
1946:             i = stag->n[0];
1947:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1948:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1949:             }
1950:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1951:               idxGlobalAll[count] = -1;
1952:             }
1953:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1954:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1955:             }
1956:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1957:               idxGlobalAll[count] = -1;
1958:             }
1959:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1960:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1961:             }
1962:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1963:               idxGlobalAll[count] = -1;
1964:             }
1965:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1966:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1967:             }
1968:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1969:               idxGlobalAll[count] = -1;
1970:             }
1971:             ++i;
1972:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1973:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1974:                 idxGlobalAll[count] = -1;
1975:               }
1976:             }
1977:           }
1978:         }
1979:       } else if (dummyEnd[1]) {
1980:         /* Up Back partial dummy row */
1981:         PetscInt j = stag->n[1];

1983:         if (!star && !dummyStart[0]) {
1984:           /* Up Back partial dummy row: Loop over columns 1/3: Left Up Back, on Left (Middle) Back rank */
1985:           const PetscInt neighbor = 3;
1986:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1987:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1988:             PetscInt       dLocal;
1989:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
1990:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1991:             }

1993:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
1994:               idxGlobalAll[count] = -1;
1995:             }
1996:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
1997:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1998:             }
1999:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2000:               idxGlobalAll[count] = -1;
2001:             }
2002:           }
2003:         } else {
2004:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2005:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2006:               idxGlobalAll[count] = -1;
2007:             }
2008:           }
2009:         }

2011:         /* Up Back partial dummy row: Loop over columns 2/3: (Middle) Up Back, on (Middle) (Middle) Back rank */
2012:         {
2013:           const PetscInt neighbor = 4;
2014:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2015:             const PetscInt i = ighost - ghostOffsetStart[0];
2016:             PetscInt       dLocal;
2017:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2018:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2019:             }
2020:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2021:               idxGlobalAll[count] = -1;
2022:             }
2023:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2024:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2025:             }
2026:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2027:               idxGlobalAll[count] = -1;
2028:             }
2029:           }
2030:         }

2032:         /* Up Back partial dummy row: Loop over columns 3/3: Right Up Back, on Right (Middle) Back rank */
2033:         if (!star && !dummyEnd[0]) {
2034:           const PetscInt neighbor = 5;
2035:           PetscInt       i;
2036:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2037:             PetscInt       dLocal;
2038:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2039:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2040:             }
2041:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2042:               idxGlobalAll[count] = -1;
2043:             }
2044:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2045:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2046:             }
2047:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2048:               idxGlobalAll[count] = -1;
2049:             }
2050:           }
2051:         } else if (dummyEnd[0]) {
2052:           /* Up Right Back partial dummy element, on (Middle) (Middle) Back rank */
2053:           const PetscInt neighbor = 4;
2054:           PetscInt       i,dLocal;
2055:           i = stag->n[0];
2056:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2057:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2058:           }
2059:           for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back down edge, back face and back left edge */
2060:             idxGlobalAll[count] = -1;
2061:           }
2062:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2063:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2064:           }
2065:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2066:             idxGlobalAll[count] = -1;
2067:           }
2068:           ++i;
2069:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2070:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2071:               idxGlobalAll[count] = -1;
2072:             }
2073:           }
2074:         } else {
2075:           /* Up Right Back dummies */
2076:           PetscInt i;
2077:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2078:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2079:               idxGlobalAll[count] = -1;
2080:             }
2081:           }
2082:         }
2083:         ++j;
2084:         /* Up Back additional dummy rows */
2085:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2086:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2087:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2088:               idxGlobalAll[count] = -1;
2089:             }
2090:           }
2091:         }
2092:       } else {
2093:         /* Up Back dummy rows */
2094:         PetscInt j;
2095:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2096:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2097:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2098:               idxGlobalAll[count] = -1;
2099:             }
2100:           }
2101:         }
2102:       }
2103:     }
2104:   } else {
2105:     for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
2106:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
2107:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2108:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2109:             idxGlobalAll[count] = -1;
2110:           }
2111:         }
2112:       }
2113:     }
2114:   }

2116:   /* Loop over layers 2/3 : (Middle)  */
2117:   {
2118:     for (kghost = ghostOffsetStart[2]; kghost<stag->nGhost[2]-ghostOffsetEnd[2]; ++kghost) {
2119:       const PetscInt k = kghost - ghostOffsetStart[2];

2121:       /* Loop over rows 1/3: Down (Middle) */
2122:       if (!dummyStart[1]) {
2123:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2124:           const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* down neighbor (10) always exists */

2126:           /* Loop over columns 1/3: Left Down (Middle) */
2127:           if (!star && !dummyStart[0]) {
2128:             const PetscInt neighbor = 9;
2129:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2130:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2131:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2132:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2133:               }
2134:             }
2135:           } else {
2136:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2137:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2138:                 idxGlobalAll[count] = -1;
2139:               }
2140:             }
2141:           }

2143:           /* Loop over columns 2/3: (Middle) Down (Middle) */
2144:           {
2145:             const PetscInt neighbor = 10;
2146:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2147:               const PetscInt i = ighost - ghostOffsetStart[0];
2148:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2149:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2150:               }
2151:             }
2152:           }

2154:           /* Loop over columns 3/3: Right Down (Middle) */
2155:           if (!star && !dummyEnd[0]) {
2156:             const PetscInt neighbor = 11;
2157:             PetscInt       i;
2158:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2159:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2160:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2161:               }
2162:             }
2163:           } else if (dummyEnd[0]) {
2164:             /* Partial dummies on (Middle) Down (Middle) neighbor */
2165:             const PetscInt neighbor = 10;
2166:             PetscInt       i,dLocal;
2167:             i = stag->n[0];
2168:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2169:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2170:             }
2171:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2172:               idxGlobalAll[count] = -1;
2173:             }
2174:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2175:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2176:             }
2177:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2178:               idxGlobalAll[count] = -1;
2179:             }
2180:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2181:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2182:             }
2183:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2184:               idxGlobalAll[count] = -1;
2185:             }
2186:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2187:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2188:             }
2189:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2190:               idxGlobalAll[count] = -1;
2191:             }
2192:             ++i;
2193:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2194:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2195:                 idxGlobalAll[count] = -1;
2196:               }
2197:             }
2198:           } else {
2199:             /* Right Down (Middle) dummies */
2200:             PetscInt       i;
2201:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2202:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2203:                 idxGlobalAll[count] = -1;
2204:               }
2205:             }
2206:           }
2207:         }
2208:       } else {
2209:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2210:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2211:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2212:               idxGlobalAll[count] = -1;
2213:             }
2214:           }
2215:         }
2216:       }

2218:       /* Loop over rows 2/3: (Middle) (Middle) */
2219:       {
2220:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2221:           const PetscInt j = jghost - ghostOffsetStart[1];

2223:           /* Loop over columns 1/3: Left (Middle) (Middle) */
2224:           if (!dummyStart[0]) {
2225:             const PetscInt neighbor = 12;
2226:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2227:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2228:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2229:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2230:               }
2231:             }
2232:           } else {
2233:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2234:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2235:                 idxGlobalAll[count] = -1;
2236:               }
2237:             }
2238:           }

2240:           /* Loop over columns 2/3: (Middle) (Middle) (Middle) aka (Here) */
2241:           {
2242:             const PetscInt neighbor = 13;
2243:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2244:               const PetscInt i = ighost - ghostOffsetStart[0];
2245:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2246:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2247:               }
2248:             }
2249:           }

2251:           /* Loop over columns 3/3: Right (Middle) (Middle) */
2252:           if (!dummyEnd[0]) {
2253:             const PetscInt neighbor = 14;
2254:             PetscInt       i;
2255:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2256:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2257:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2258:               }
2259:             }
2260:           } else {
2261:             /* Partial dummies on this rank */
2262:             const PetscInt neighbor = 13;
2263:             PetscInt       i,dLocal;
2264:             i = stag->n[0];
2265:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2266:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2267:             }
2268:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2269:               idxGlobalAll[count] = -1;
2270:             }
2271:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2272:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2273:             }
2274:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2275:               idxGlobalAll[count] = -1;
2276:             }
2277:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2278:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2279:             }
2280:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2281:               idxGlobalAll[count] = -1;
2282:             }
2283:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2284:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2285:             }
2286:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2287:               idxGlobalAll[count] = -1;
2288:             }
2289:             ++i;
2290:             for(;i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2291:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2292:                 idxGlobalAll[count] = -1;
2293:               }
2294:             }
2295:           }
2296:         }
2297:       }

2299:       /* Loop over rows 3/3: Up (Middle) */
2300:       if (!dummyEnd[1]) {
2301:         PetscInt j;
2302:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2304:           /* Loop over columns 1/3: Left Up (Middle) */
2305:           if (!star && !dummyStart[0]) {
2306:             const PetscInt neighbor = 15;
2307:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2308:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2309:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2310:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2311:               }
2312:             }
2313:           } else {
2314:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2315:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2316:                 idxGlobalAll[count] = -1;
2317:               }
2318:             }
2319:           }

2321:           /* Loop over columns 2/3: (Middle) Up (Middle) */
2322:           {
2323:             const PetscInt neighbor = 16;
2324:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2325:               const PetscInt i = ighost - ghostOffsetStart[0];
2326:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2327:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2328:               }
2329:             }
2330:           }

2332:           /* Loop over columns 3/3: Right Up (Middle) */
2333:           if (!star && !dummyEnd[0]) {
2334:             const PetscInt neighbor = 17;
2335:             PetscInt       i;
2336:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2337:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2338:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2339:               }
2340:             }
2341:           } else if (dummyEnd[0]) {
2342:             /* Partial dummies on (Middle) Up (Middle) rank */
2343:             const PetscInt neighbor = 16;
2344:             PetscInt       i,dLocal;
2345:             i = stag->n[0];
2346:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2347:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2348:             }
2349:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2350:               idxGlobalAll[count] = -1;
2351:             }
2352:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2353:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2354:             }
2355:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2356:               idxGlobalAll[count] = -1;
2357:             }
2358:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2359:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2360:             }
2361:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2362:               idxGlobalAll[count] = -1;
2363:             }
2364:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2365:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2366:             }
2367:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2368:               idxGlobalAll[count] = -1;
2369:             }
2370:             ++i;
2371:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2372:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2373:                 idxGlobalAll[count] = -1;
2374:               }
2375:             }
2376:           } else {
2377:             /* Right Up (Middle) dumies */
2378:             PetscInt       i;
2379:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2380:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2381:                 idxGlobalAll[count] = -1;
2382:               }
2383:             }
2384:           }
2385:         }
2386:       } else {
2387:         /* Up (Middle) partial dummy row */
2388:         PetscInt j = stag->n[1];

2390:         /* Up (Middle) partial dummy row: colums 1/3: Left Up (Middle), on Left (Middle) (Middle) rank */
2391:         if (!dummyStart[0]) {
2392:           const PetscInt neighbor = 12;
2393:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2394:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2395:             PetscInt       dLocal;
2396:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2397:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2398:             }
2399:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2400:               idxGlobalAll[count] = -1;
2401:             }
2402:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2403:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2404:             }
2405:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2406:               idxGlobalAll[count] = -1;
2407:             }
2408:           }
2409:         } else {
2410:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2411:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2412:               idxGlobalAll[count] = -1;
2413:             }
2414:           }
2415:         }

2417:         /* Up (Middle) partial dummy row: columns 2/3: (Middle) Up (Middle), on (Middle) (Middle) (Middle) = here rank */
2418:         {
2419:           const PetscInt neighbor = 13;
2420:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2421:             const PetscInt i = ighost - ghostOffsetStart[0];
2422:             PetscInt       dLocal;
2423:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2424:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2425:             }
2426:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2427:               idxGlobalAll[count] = -1;
2428:             }
2429:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2430:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2431:             }
2432:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2433:               idxGlobalAll[count] = -1;
2434:             }
2435:           }
2436:         }

2438:         if (!dummyEnd[0]) {
2439:           /* Up (Middle) partial dummy row: columns 3/3: Right Up (Middle), on Right (Middle) (Middle) rank */
2440:           const PetscInt neighbor = 14;
2441:           PetscInt       i;
2442:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2443:             PetscInt       dLocal;
2444:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2445:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2446:             }
2447:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2448:               idxGlobalAll[count] = -1;
2449:             }
2450:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2451:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2452:             }
2453:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2454:               idxGlobalAll[count] = -1;
2455:             }
2456:           }
2457:         } else {
2458:           /* Up (Middle) Right partial dummy element, on (Middle) (Middle) (Middle) = here rank */
2459:           const PetscInt neighbor = 13;
2460:           PetscInt       i,dLocal;
2461:           i = stag->n[0];
2462:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2463:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2464:           }
2465:           for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back down edge, back face and back left edge */
2466:             idxGlobalAll[count] = -1;
2467:           }
2468:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2469:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2470:           }
2471:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2472:             idxGlobalAll[count] = -1;
2473:           }
2474:           ++i;
2475:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2476:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2477:               idxGlobalAll[count] = -1;
2478:             }
2479:           }
2480:         }
2481:         ++j;
2482:         /* Up (Middle) additional dummy rows */
2483:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2484:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2485:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2486:               idxGlobalAll[count] = -1;
2487:             }
2488:           }
2489:         }
2490:       }
2491:     }
2492:   }

2494:   /* Loop over layers 3/3 : Front */
2495:   if (!dummyEnd[2]) {
2496:     PetscInt k;
2497:     for (k=0; k<ghostOffsetEnd[2]; ++k) {

2499:       /* Loop over rows 1/3: Down Front */
2500:       if (!star && !dummyStart[1]) {
2501:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2502:           const PetscInt j = nNeighbors[19][1] - ghostOffsetStart[1] + jghost; /* Constant for whole row (use down front neighbor) */

2504:           /* Loop over columns 1/3: Left Down Front */
2505:           if (!dummyStart[0]) {
2506:             const PetscInt neighbor = 18;
2507:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2508:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2509:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2510:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2511:               }
2512:             }
2513:           } else {
2514:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2515:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2516:                 idxGlobalAll[count] = -1;
2517:               }
2518:             }
2519:           }

2521:           /* Loop over columns 2/3: (Middle) Down Front */
2522:           {
2523:             const PetscInt neighbor = 19;
2524:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2525:               const PetscInt i = ighost - ghostOffsetStart[0];
2526:               for (d=0; d<stag->entriesPerElement; ++d,++count) { /* vertex, 2 edges, and face associated with back face */
2527:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2528:               }
2529:             }
2530:           }

2532:           /* Loop over columns 3/3: Right Down Front */
2533:           if (!dummyEnd[0]) {
2534:             const PetscInt neighbor = 20;
2535:             PetscInt       i;
2536:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2537:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2538:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2539:               }
2540:             }
2541:           } else {
2542:             /* Partial dummy element on (Middle) Down Front rank */
2543:             const PetscInt neighbor = 19;
2544:             PetscInt       i,dLocal;
2545:             i = stag->n[0];
2546:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2547:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2548:             }
2549:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2550:               idxGlobalAll[count] = -1;
2551:             }
2552:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2553:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2554:             }
2555:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2556:               idxGlobalAll[count] = -1;
2557:             }
2558:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2559:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2560:             }
2561:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2562:               idxGlobalAll[count] = -1;
2563:             }
2564:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2565:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2566:             }
2567:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2568:               idxGlobalAll[count] = -1;
2569:             }
2570:             ++i;
2571:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2572:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2573:                 idxGlobalAll[count] = -1;
2574:               }
2575:             }
2576:           }
2577:         }
2578:       } else {
2579:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2580:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2581:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2582:               idxGlobalAll[count] = -1;
2583:             }
2584:           }
2585:         }
2586:       }

2588:       /* Loop over rows 2/3: (Middle) Front */
2589:       {
2590:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2591:           const PetscInt j = jghost - ghostOffsetStart[1];

2593:           /* Loop over columns 1/3: Left (Middle) Front*/
2594:           if (!star && !dummyStart[0]) {
2595:             const PetscInt neighbor = 21;
2596:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2597:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2598:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2599:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2600:               }
2601:             }
2602:           } else {
2603:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2604:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2605:                 idxGlobalAll[count] = -1;
2606:               }
2607:             }
2608:           }

2610:           /* Loop over columns 2/3: (Middle) (Middle) Front*/
2611:           {
2612:             const PetscInt neighbor = 22;
2613:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2614:               const PetscInt i = ighost - ghostOffsetStart[0];
2615:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2616:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2617:               }
2618:             }
2619:           }

2621:           /* Loop over columns 3/3: Right (Middle) Front*/
2622:           if (!star && !dummyEnd[0]) {
2623:             const PetscInt neighbor = 23;
2624:             PetscInt       i;
2625:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2626:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2627:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2628:               }
2629:             }
2630:           } else if (dummyEnd[0]) {
2631:             /* Partial dummy element on (Middle) (Middle) Front element */
2632:             const PetscInt neighbor = 22;
2633:             PetscInt       i,dLocal;
2634:             i = stag->n[0];
2635:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2636:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2637:             }
2638:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2639:               idxGlobalAll[count] = -1;
2640:             }
2641:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2642:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2643:             }
2644:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2645:               idxGlobalAll[count] = -1;
2646:             }
2647:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2648:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2649:             }
2650:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2651:               idxGlobalAll[count] = -1;
2652:             }
2653:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2654:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2655:             }
2656:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2657:               idxGlobalAll[count] = -1;
2658:             }
2659:             ++i;
2660:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2661:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2662:                 idxGlobalAll[count] = -1;
2663:               }
2664:             }
2665:           } else {
2666:             /* Right (Middle) Front dummies */
2667:             PetscInt       i;
2668:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2669:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2670:                 idxGlobalAll[count] = -1;
2671:               }
2672:             }
2673:           }
2674:         }
2675:       }

2677:       /* Loop over rows 3/3: Up Front */
2678:       if (!star && !dummyEnd[1]) {
2679:         PetscInt j;
2680:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2682:           /* Loop over columns 1/3: Left Up Front */
2683:           if (!dummyStart[0]) {
2684:             const PetscInt neighbor = 24;
2685:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2686:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2687:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2688:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2689:               }
2690:             }
2691:           } else {
2692:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2693:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2694:                 idxGlobalAll[count] = -1;
2695:               }
2696:             }
2697:           }

2699:           /* Loop over columns 2/3: (Middle) Up Front */
2700:           {
2701:             const PetscInt neighbor = 25;
2702:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2703:               const PetscInt i = ighost - ghostOffsetStart[0];
2704:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2705:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2706:               }
2707:             }
2708:           }

2710:           /* Loop over columns 3/3: Right Up Front */
2711:           if (!dummyEnd[0]) {
2712:             const PetscInt neighbor = 26;
2713:             PetscInt     i;
2714:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2715:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2716:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2717:               }
2718:             }
2719:           } else {
2720:             /* Partial dummy element on (Middle) Up Front rank */
2721:             const PetscInt neighbor = 25;
2722:             PetscInt       i,dLocal;
2723:             i = stag->n[0];
2724:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2725:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2726:             }
2727:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2728:               idxGlobalAll[count] = -1;
2729:             }
2730:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2731:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2732:             }
2733:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2734:               idxGlobalAll[count] = -1;
2735:             }
2736:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2737:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2738:             }
2739:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2740:               idxGlobalAll[count] = -1;
2741:             }
2742:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2743:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2744:             }
2745:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2746:               idxGlobalAll[count] = -1;
2747:             }
2748:             ++i;
2749:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2750:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2751:                 idxGlobalAll[count] = -1;
2752:               }
2753:             }
2754:           }
2755:         }
2756:       } else if (dummyEnd[1]) {
2757:         /* Up Front partial dummy row */
2758:         PetscInt j = stag->n[1];

2760:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) Front rank */
2761:         if (!star && !dummyStart[0]) {
2762:           const PetscInt neighbor = 21;
2763:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2764:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2765:             PetscInt       dLocal;
2766:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2767:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2768:             }
2769:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2770:               idxGlobalAll[count] = -1;
2771:             }
2772:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2773:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2774:             }
2775:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2776:               idxGlobalAll[count] = -1;
2777:             }
2778:           }
2779:         } else {
2780:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2781:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2782:               idxGlobalAll[count] = -1;
2783:             }
2784:           }
2785:         }

2787:         /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) Front rank */
2788:         {
2789:           const PetscInt neighbor = 22;
2790:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2791:             const PetscInt i = ighost - ghostOffsetStart[0];
2792:             PetscInt       dLocal;
2793:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2794:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2795:             }
2796:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2797:               idxGlobalAll[count] = -1;
2798:             }
2799:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2800:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2801:             }
2802:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2803:               idxGlobalAll[count] = -1;
2804:             }
2805:           }
2806:         }

2808:         if (!star && !dummyEnd[0]) {
2809:           /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) Front rank */
2810:           const PetscInt neighbor = 23;
2811:           PetscInt       i;
2812:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2813:             PetscInt       dLocal;
2814:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2815:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2816:             }
2817:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2818:               idxGlobalAll[count] = -1;
2819:             }
2820:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2821:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2822:             }
2823:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2824:               idxGlobalAll[count] = -1;
2825:             }
2826:           }

2828:         } else if (dummyEnd[0]) {
2829:           /* Partial Right Up Front dummy element, on (Middle) (Middle) Front rank */
2830:           const PetscInt neighbor = 22;
2831:           PetscInt       i,dLocal;
2832:           i = stag->n[0];
2833:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2834:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2835:           }
2836:           for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back down edge, back face and back left edge */
2837:             idxGlobalAll[count] = -1;
2838:           }
2839:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2840:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2841:           }
2842:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2843:             idxGlobalAll[count] = -1;
2844:           }
2845:           ++i;
2846:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2847:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2848:               idxGlobalAll[count] = -1;
2849:             }
2850:           }
2851:         } else {
2852:           /* Right Up Front dummies */
2853:           PetscInt i;
2854:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2855:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2856:               idxGlobalAll[count] = -1;
2857:             }
2858:           }
2859:         }
2860:         ++j;
2861:         /* Up Front additional dummy rows */
2862:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2863:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2864:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2865:               idxGlobalAll[count] = -1;
2866:             }
2867:           }
2868:         }
2869:       } else {
2870:         /* Up Front dummies */
2871:         PetscInt j;
2872:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2873:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2874:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2875:               idxGlobalAll[count] = -1;
2876:             }
2877:           }
2878:         }
2879:       }
2880:     }
2881:   } else {
2882:     PetscInt k = stag->n[2];

2884:     /* Front partial ghost layer: rows 1/3: Down Front, on Down (Middle) */
2885:     if (!dummyStart[1]) {
2886:       for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2887:         const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* wrt down neighbor (10) */

2889:         /* Down Front partial ghost row: colums 1/3: Left Down Front, on  Left Down (Middle) */
2890:         if (!star && !dummyStart[0]) {
2891:           const PetscInt neighbor = 9;
2892:           const PetscInt epFaceRow         = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
2893:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2894:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2895:             PetscInt  dLocal;
2896:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2897:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2898:             }
2899:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2900:               idxGlobalAll[count] = -1;
2901:             }
2902:           }
2903:         } else {
2904:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2905:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2906:               idxGlobalAll[count] = -1;
2907:             }
2908:           }
2909:         }

2911:         /* Down Front partial ghost row: columns 2/3: (Middle) Down Front, on (Middle) Down (Middle) */
2912:         {
2913:           const PetscInt neighbor = 10;
2914:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2915:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2916:             const PetscInt i = ighost - ghostOffsetStart[0];
2917:             PetscInt  dLocal;
2918:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2919:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2920:             }
2921:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2922:               idxGlobalAll[count] = -1;
2923:             }
2924:           }
2925:         }

2927:         if (!star && !dummyEnd[0]) {
2928:           /* Down Front partial dummy row: columns 3/3: Right Down Front, on Right Down (Middle) */
2929:           const PetscInt neighbor = 11;
2930:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
2931:           PetscInt       i;
2932:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2933:             PetscInt  dLocal;
2934:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2935:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2936:             }
2937:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2938:               idxGlobalAll[count] = -1;
2939:             }
2940:           }
2941:         } else if (dummyEnd[0]) {
2942:           /* Right Down Front partial dummy element, living on (Middle) Down (Middle) rank */
2943:           const PetscInt neighbor = 10;
2944:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2945:           PetscInt       i,dLocal;
2946:           i = stag->n[0];
2947:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2948:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2949:           }
2950:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2951:             idxGlobalAll[count] = -1;
2952:           }
2953:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
2954:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2955:           }
2956:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
2957:             idxGlobalAll[count] = -1;
2958:           }
2959:           ++i;
2960:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2961:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2962:               idxGlobalAll[count] = -1;
2963:             }
2964:           }
2965:         } else {
2966:           PetscInt i;
2967:           /* Right Down Front dummies */
2968:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2969:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2970:               idxGlobalAll[count] = -1;
2971:             }
2972:           }
2973:         }
2974:       }
2975:     } else {
2976:       for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2977:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2978:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2979:             idxGlobalAll[count] = -1;
2980:           }
2981:         }
2982:       }
2983:     }

2985:     /* Front partial ghost layer: rows 2/3: (Middle) Front, on (Middle) (Middle) */
2986:     {
2987:       for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2988:         const PetscInt j = jghost - ghostOffsetStart[1];

2990:         /* (Middle) Front partial dummy row: columns 1/3: Left (Middle) Front, on Left (Middle) (Middle) */
2991:         if (!dummyStart[0]) {
2992:           const PetscInt neighbor = 12;
2993:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
2994:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2995:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2996:             PetscInt  dLocal;
2997:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2998:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2999:             }
3000:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3001:               idxGlobalAll[count] = -1;
3002:             }
3003:           }
3004:         } else {
3005:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3006:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3007:               idxGlobalAll[count] = -1;
3008:             }
3009:           }
3010:         }

3012:         /* (Middle) Front partial dummy row: columns 2/3: (Middle) (Middle) Front, on (Middle) (Middle) (Middle) */
3013:         {
3014:           const PetscInt neighbor = 13;
3015:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3016:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3017:             const PetscInt i = ighost - ghostOffsetStart[0];
3018:             PetscInt  dLocal;
3019:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3020:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3021:             }
3022:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3023:               idxGlobalAll[count] = -1;
3024:             }
3025:           }
3026:         }

3028:         if (!dummyEnd[0]) {
3029:           /* (Middle) Front partial dummy row: columns 3/3: Right (Middle) Front, on Right (Middle) (Middle) */
3030:           const PetscInt neighbor = 14;
3031:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3032:           PetscInt i;
3033:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3034:             PetscInt  dLocal;
3035:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3036:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3037:             }
3038:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3039:               idxGlobalAll[count] = -1;
3040:             }
3041:           }
3042:         } else {
3043:           /* Right (Middle) Front partial dummy element, on (Middle) (Middle) (Middle) rank */
3044:           const PetscInt neighbor = 13;
3045:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3046:           PetscInt       i,dLocal;
3047:           i = stag->n[0];
3048:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3049:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3050:           }
3051:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3052:             idxGlobalAll[count] = -1;
3053:           }
3054:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3055:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3056:           }
3057:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3058:             idxGlobalAll[count] = -1;
3059:           }
3060:           ++i;
3061:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3062:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3063:               idxGlobalAll[count] = -1;
3064:             }
3065:           }
3066:         }
3067:       }
3068:     }

3070:     /* Front partial ghost layer: rows 3/3: Up Front, on Up (Middle) */
3071:     if (!dummyEnd[1]) {
3072:       PetscInt j;
3073:       for (j=0; j<ghostOffsetEnd[1]; ++j) {

3075:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left Up (Middle) */
3076:         if (!star && !dummyStart[0]) {
3077:           const PetscInt neighbor = 15;
3078:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3079:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3080:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3081:             PetscInt  dLocal;
3082:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3083:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3084:             }
3085:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3086:               idxGlobalAll[count] = -1;
3087:             }
3088:           }
3089:         } else {
3090:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3091:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3092:               idxGlobalAll[count] = -1;
3093:             }
3094:           }
3095:         }

3097:         /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) Up (Middle) */
3098:         {
3099:           const PetscInt neighbor = 16;
3100:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3101:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3102:             const PetscInt i = ighost - ghostOffsetStart[0];
3103:             PetscInt  dLocal;
3104:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3105:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3106:             }
3107:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3108:               idxGlobalAll[count] = -1;
3109:             }
3110:           }
3111:         }

3113:         if (!star && !dummyEnd[0]) {
3114:           /* Up Front partial dummy row: columsn 3/3: Right Up Front, on Right Up (Middle) */
3115:           const PetscInt neighbor = 17;
3116:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3117:           PetscInt       i;
3118:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3119:             PetscInt  dLocal;
3120:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3121:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3122:             }
3123:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3124:               idxGlobalAll[count] = -1;
3125:             }
3126:           }
3127:         } else if (dummyEnd[0]) {
3128:           /* Right Up Front partial dummy element, on (Middle) Up (Middle) */
3129:           const PetscInt neighbor = 16;
3130:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3131:           PetscInt       i,dLocal;
3132:           i = stag->n[0];
3133:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3134:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3135:           }
3136:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3137:             idxGlobalAll[count] = -1;
3138:           }
3139:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3140:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3141:           }
3142:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3143:             idxGlobalAll[count] = -1;
3144:           }
3145:           ++i;
3146:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3147:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3148:               idxGlobalAll[count] = -1;
3149:             }
3150:           }
3151:         } else {
3152:           /* Right Up Front dummies */
3153:           PetscInt i;
3154:           /* Right Down Front dummies */
3155:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
3156:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3157:               idxGlobalAll[count] = -1;
3158:             }
3159:           }
3160:         }
3161:       }
3162:     } else {
3163:       /* Up Front partial dummy row */
3164:       PetscInt j = stag->n[1];

3166:       /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) (Middle) */
3167:       if (!dummyStart[0]) {
3168:         const PetscInt neighbor = 12;
3169:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3170:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3171:           const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3172:           PetscInt       dLocal;
3173:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3174:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3175:           }
3176:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3177:             idxGlobalAll[count] = -1;
3178:           }
3179:         }
3180:       } else {
3181:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3182:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3183:             idxGlobalAll[count] = -1;
3184:           }
3185:         }
3186:       }

3188:       /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) (Middle) */
3189:       {
3190:         const PetscInt neighbor = 13;
3191:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3192:         for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3193:           const PetscInt i = ighost - ghostOffsetStart[0];
3194:           PetscInt       dLocal;
3195:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3196:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3197:           }
3198:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3199:             idxGlobalAll[count] = -1;
3200:           }
3201:         }
3202:       }

3204:       if (!dummyEnd[0]) {
3205:         /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) (Middle) */
3206:         const PetscInt neighbor = 14;
3207:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3208:         PetscInt       i;
3209:         for (i=0; i<ghostOffsetEnd[0]; ++i) {
3210:           PetscInt       dLocal;
3211:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3212:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3213:           }
3214:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3215:             idxGlobalAll[count] = -1;
3216:           }
3217:         }
3218:       } else {
3219:         /* Right Up Front partial dummy element, on this rank */
3220:         const PetscInt neighbor = 13;
3221:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3222:         PetscInt       i,dLocal;
3223:         i = stag->n[0];
3224:         for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3225:           idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3226:         }
3227:         for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3228:           idxGlobalAll[count] = -1;
3229:         }
3230:         ++i;
3231:         for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3232:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3233:             idxGlobalAll[count] = -1;
3234:           }
3235:         }
3236:       }
3237:       ++j;
3238:       /* Up Front additional dummy rows */
3239:       for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
3240:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3241:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3242:             idxGlobalAll[count] = -1;
3243:           }
3244:         }
3245:       }
3246:     }
3247:     /* Front additional dummy layers */
3248:     ++k;
3249:     for (; k<stag->n[2] + ghostOffsetEnd[2]; ++k) {
3250:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
3251:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3252:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3253:             idxGlobalAll[count] = -1;
3254:           }
3255:         }
3256:       }
3257:     }
3258:   }

3260:   /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
3261:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
3262:   PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
3263:   return(0);
3264: }

3266: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM dm)
3267: {
3268:   PetscErrorCode  ierr;
3269:   DM_Stag * const stag = (DM_Stag*)dm->data;
3270:   const PetscInt epe = stag->entriesPerElement;
3271:   const PetscInt epr = stag->nGhost[0]*epe;
3272:   const PetscInt epl = stag->nGhost[1]*epr;

3275:   PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
3276:   stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]   = 0;
3277:   stag->locationOffsets[DMSTAG_BACK_DOWN]        = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + stag->dof[0];
3278:   stag->locationOffsets[DMSTAG_BACK_DOWN_RIGHT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epe;
3279:   stag->locationOffsets[DMSTAG_BACK_LEFT]        = stag->locationOffsets[DMSTAG_BACK_DOWN]       + stag->dof[1];
3280:   stag->locationOffsets[DMSTAG_BACK]             = stag->locationOffsets[DMSTAG_BACK_LEFT]       + stag->dof[1];
3281:   stag->locationOffsets[DMSTAG_BACK_RIGHT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epe;
3282:   stag->locationOffsets[DMSTAG_BACK_UP_LEFT]     = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epr;
3283:   stag->locationOffsets[DMSTAG_BACK_UP]          = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epr;
3284:   stag->locationOffsets[DMSTAG_BACK_UP_RIGHT]    = stag->locationOffsets[DMSTAG_BACK_UP_LEFT]    + epe;
3285:   stag->locationOffsets[DMSTAG_DOWN_LEFT]        = stag->locationOffsets[DMSTAG_BACK]            + stag->dof[2];
3286:   stag->locationOffsets[DMSTAG_DOWN]             = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + stag->dof[1];
3287:   stag->locationOffsets[DMSTAG_DOWN_RIGHT]       = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + epe;
3288:   stag->locationOffsets[DMSTAG_LEFT]             = stag->locationOffsets[DMSTAG_DOWN]            + stag->dof[2];
3289:   stag->locationOffsets[DMSTAG_ELEMENT]          = stag->locationOffsets[DMSTAG_LEFT]            + stag->dof[2];
3290:   stag->locationOffsets[DMSTAG_RIGHT]            = stag->locationOffsets[DMSTAG_LEFT]            + epe;
3291:   stag->locationOffsets[DMSTAG_UP_LEFT]          = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + epr;
3292:   stag->locationOffsets[DMSTAG_UP]               = stag->locationOffsets[DMSTAG_DOWN]            + epr;
3293:   stag->locationOffsets[DMSTAG_UP_RIGHT]         = stag->locationOffsets[DMSTAG_UP_LEFT]         + epe;
3294:   stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epl;
3295:   stag->locationOffsets[DMSTAG_FRONT_DOWN]       = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epl;
3296:   stag->locationOffsets[DMSTAG_FRONT_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epe;
3297:   stag->locationOffsets[DMSTAG_FRONT_LEFT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epl;
3298:   stag->locationOffsets[DMSTAG_FRONT]            = stag->locationOffsets[DMSTAG_BACK]            + epl;
3299:   stag->locationOffsets[DMSTAG_FRONT_RIGHT]      = stag->locationOffsets[DMSTAG_FRONT_LEFT]      + epe;
3300:   stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]    = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epr;
3301:   stag->locationOffsets[DMSTAG_FRONT_UP]         = stag->locationOffsets[DMSTAG_FRONT_DOWN]      + epr;
3302:   stag->locationOffsets[DMSTAG_FRONT_UP_RIGHT]   = stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]   + epe;
3303:   return(0);
3304: }

3306: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_3d(DM dm)
3307: {
3308:   PetscErrorCode  ierr;
3309:   DM_Stag * const stag = (DM_Stag*)dm->data;
3310:   PetscInt        *idxLocal,*idxGlobal,*globalOffsetsRecomputed;
3311:   const PetscInt  *globalOffsets;
3312:   PetscInt        count,d,entriesPerEdge,entriesPerFace,eprGhost,eplGhost,ghostOffsetStart[3],ghostOffsetEnd[3];
3313:   IS              isLocal,isGlobal;
3314:   PetscBool       dummyEnd[3];

3317:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsetsRecomputed); /* note that we don't actually use all of these. An available optimization is to pass them, when available */
3318:   globalOffsets = globalOffsetsRecomputed;
3319:   PetscMalloc1(stag->entries,&idxLocal);
3320:   PetscMalloc1(stag->entries,&idxGlobal);
3321:   for (d=0; d<3; ++d) dummyEnd[d]   = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
3322:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
3323:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
3324:   eprGhost                            = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row   */
3325:   eplGhost                            = stag->nGhost[1]*eprGhost;                /* epl = entries per (element) layer */
3326:   count = 0;
3327:   for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
3328:   for (d=0; d<3; ++d) ghostOffsetEnd[d]   = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
3329:   {
3330:     const PetscInt  neighbor     = 13;
3331:     const PetscInt  epr          = stag->entriesPerElement * stag->n[0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
3332:     const PetscInt  epl          = epr                     * stag->n[1] + (dummyEnd[1] ? stag->n[0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
3333:     const PetscInt  epFaceRow    = entriesPerFace          * stag->n[0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3334:     const PetscInt  start0       = 0;
3335:     const PetscInt  start1       = 0;
3336:     const PetscInt  start2       = 0;
3337:     const PetscInt  startGhost0  = ghostOffsetStart[0];
3338:     const PetscInt  startGhost1  = ghostOffsetStart[1];
3339:     const PetscInt  startGhost2  = ghostOffsetStart[2];
3340:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
3341:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
3342:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
3343:     const PetscBool extra0       = dummyEnd[0];
3344:     const PetscBool extra1       = dummyEnd[1];
3345:     const PetscBool extra2       = dummyEnd[2];
3346:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,epr,epl,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
3347:   }
3348:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
3349:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
3350:   {
3351:     Vec local,global;
3352:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
3353:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
3354:     VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
3355:     VecDestroy(&global);
3356:     VecDestroy(&local);
3357:   }
3358:   ISDestroy(&isLocal);
3359:   ISDestroy(&isGlobal);
3360:   if (globalOffsetsRecomputed) {
3361:     PetscFree(globalOffsetsRecomputed);
3362:   }
3363:   return(0);
3364: }