Actual source code: stag3d.c

petsc-master 2020-10-22
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])) * stag->entriesPerElement
925:         + (dummyEnd[0]                               ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace   : 0)
926:         + (dummyEnd[1]                               ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace   : 0)
927:         + (dummyEnd[2]                               ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace   : 0)
928:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
929:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
930:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
931:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
932:     } else {
933:       entriesToTransferTotal  =  nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
934:         + (dummyEnd[0]                               ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace   : 0)
935:         + (dummyEnd[1]                               ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace   : 0)
936:         + (dummyEnd[2]                               ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace   : 0)
937:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
938:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
939:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
940:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
941:     }
942:   }

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

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

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

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

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

969:   /* LEFT DOWN BACK */
970:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
971:     const PetscInt  neighbor     = 0;
972:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
973:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
974:     const PetscInt  epFaceRow    = -1;
975:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
976:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
977:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
978:     const PetscInt  startGhost0  = 0;
979:     const PetscInt  startGhost1  = 0;
980:     const PetscInt  startGhost2  = 0;
981:     const PetscInt  endGhost0    = ghostOffsetStart[0];
982:     const PetscInt  endGhost1    = ghostOffsetStart[1];
983:     const PetscInt  endGhost2    = ghostOffsetStart[2];
984:     const PetscBool extra0       = PETSC_FALSE;
985:     const PetscBool extra1       = PETSC_FALSE;
986:     const PetscBool extra2       = PETSC_FALSE;
987:     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);
988:   }

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

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

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

1053:   /* BACK */
1054:   if (!dummyStart[2]) {
1055:     const PetscInt  neighbor     = 4;
1056:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may  be a right boundary */
1057:     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 */
1058:     const PetscInt  epFaceRow    = -1;
1059:     const PetscInt  start0       = 0;
1060:     const PetscInt  start1       = 0;
1061:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1062:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1063:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1064:     const PetscInt  startGhost2  = 0;
1065:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1066:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1067:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1068:     const PetscBool extra0       = dummyEnd[0];
1069:     const PetscBool extra1       = dummyEnd[1];
1070:     const PetscBool extra2       = PETSC_FALSE;
1071:     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);
1072:   }

1074:   /* RIGHT BACK */
1075:   if (!star && !dummyEnd[0] && !dummyStart[2]) {
1076:     const PetscInt  neighbor     = 5;
1077:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1078:     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 */
1079:     const PetscInt  epFaceRow    = -1;
1080:     const PetscInt  start0       = 0;
1081:     const PetscInt  start1       = 0;
1082:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1083:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1084:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1085:     const PetscInt  startGhost2  = 0;
1086:     const PetscInt  endGhost0    = stag->nGhost[0];
1087:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1088:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1089:     const PetscBool extra0       = PETSC_FALSE;
1090:     const PetscBool extra1       = dummyEnd[1];
1091:     const PetscBool extra2       = PETSC_FALSE;
1092:     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);
1093:   }

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

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

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

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

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

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

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

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

1263:   /* RIGHT */
1264:   if (!dummyEnd[0]) {
1265:     const PetscInt  neighbor     = 14;
1266:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1267:     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 */
1268:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1269:     const PetscInt  start0       = 0;
1270:     const PetscInt  start1       = 0;
1271:     const PetscInt  start2       = 0;
1272:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1273:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1274:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1275:     const PetscInt  endGhost0    = stag->nGhost[0];
1276:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1277:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1278:     const PetscBool extra0       = PETSC_FALSE;
1279:     const PetscBool extra1       = dummyEnd[1];
1280:     const PetscBool extra2       = dummyEnd[2];
1281:     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);
1282:   }

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

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

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

1347:   /* LEFT DOWN FRONT */
1348:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1349:     const PetscInt  neighbor     = 18;
1350:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1351:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1352:     const PetscInt  epFaceRow    = -1;
1353:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1354:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1355:     const PetscInt  start2       = 0;
1356:     const PetscInt  startGhost0  = 0;
1357:     const PetscInt  startGhost1  = 0;
1358:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1359:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1360:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1361:     const PetscInt  endGhost2    = stag->nGhost[2];
1362:     const PetscBool extra0       = PETSC_FALSE;
1363:     const PetscBool extra1       = PETSC_FALSE;
1364:     const PetscBool extra2       = PETSC_FALSE;
1365:     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);
1366:   }

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

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

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

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

1452:   /* RIGHT FRONT */
1453:   if (!star && !dummyEnd[0] && !dummyEnd[2]) {
1454:     const PetscInt  neighbor     = 23;
1455:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1456:     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 */
1457:     const PetscInt  epFaceRow    = -1;
1458:     const PetscInt  start0       = 0;
1459:     const PetscInt  start1       = 0;
1460:     const PetscInt  start2       = 0;
1461:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1462:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1463:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1464:     const PetscInt  endGhost0    = stag->nGhost[0];
1465:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1466:     const PetscInt  endGhost2    = stag->nGhost[2];
1467:     const PetscBool extra0       = PETSC_FALSE;
1468:     const PetscBool extra1       = dummyEnd[1];
1469:     const PetscBool extra2       = PETSC_FALSE;
1470:     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);
1471:   }

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

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

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

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

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

1553:   return(0);
1554: }

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


1569:   /* Check stencil type */
1570:   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]);
1571:   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");
1572:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

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

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

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

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

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

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

1717:   /* Loop over layers 1/3 : Back */
1718:   if (!dummyStart[2]) {
1719:     for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
1720:       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) */

1722:       /* Loop over rows 1/3: Down Back*/
1723:       if (!star && !dummyStart[1]) {
1724:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1725:           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)*/

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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