Actual source code: stag3d.c

petsc-master 2019-10-14
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: {
 50:   DM_Stag        *stag;

 53:   DMCreate(comm,dm);
 54:   DMSetDimension(*dm,3); /* Must precede DMSetType */
 55:   DMSetType(*dm,DMSTAG);
 56:   stag = (DM_Stag*)(*dm)->data;
 57:   DMSetDimension(*dm,3);
 58:   stag->boundaryType[0] = bndx;
 59:   stag->boundaryType[1] = bndy;
 60:   stag->boundaryType[2] = bndz;
 61:   stag->N[0]            = M;
 62:   stag->N[1]            = N;
 63:   stag->N[2]            = P;
 64:   stag->nRanks[0]       = m; /* Adjusted later in DMSetUp_Stag */
 65:   stag->nRanks[1]       = n; /* Adjusted later in DMSetUp_Stag */
 66:   stag->nRanks[2]       = p; /* Adjusted later in DMSetUp_Stag */
 67:   stag->stencilType     = stencilType;
 68:   stag->stencilWidth    = stencilWidth;
 69:   DMStagSetDOF(*dm,dof0,dof1,dof2,dof3);
 70:   DMStagSetOwnershipRanges(*dm,lx,ly,lz);
 71:   return(0);
 72: }

 74: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_3d(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
 75: {
 77:   DM_Stag        *stagCoord;
 78:   DM             dmCoord;
 79:   Vec            coordLocal,coord;
 80:   PetscReal      h[3],min[3];
 81:   PetscScalar    ****arr;
 82:   PetscInt       ind[3],start[3],n[3],nExtra[3],s,c;
 83:   PetscInt       ibackdownleft,ibackdown,ibackleft,iback,idownleft,idown,ileft,ielement;

 86:   DMGetCoordinateDM(dm,&dmCoord);
 87:   stagCoord = (DM_Stag*) dmCoord->data;
 88:   for (s=0; s<4; ++s) {
 89:     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]);
 90:   }
 91:   DMGetLocalVector(dmCoord,&coordLocal);
 92:   DMStagVecGetArrayDOF(dmCoord,coordLocal,&arr);
 93:   if (stagCoord->dof[0]) {
 94:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN_LEFT,0,&ibackdownleft);
 95:   }
 96:   if (stagCoord->dof[1]) {
 97:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN     ,0,&ibackdown    );
 98:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_LEFT     ,0,&ibackleft    );
 99:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN_LEFT     ,0,&idownleft    );
100:   }
101:   if (stagCoord->dof[2]) {
102:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK          ,0,&iback        );
103:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN          ,0,&idown        );
104:     DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT          ,0,&ileft        );
105:   }
106:   if (stagCoord->dof[3]) {
107:     DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT       ,0,&ielement     );
108:     DMStagGetCorners(dmCoord,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
109:   }
110:   min[0] = xmin; min[1]= ymin; min[2] = zmin;
111:   h[0] = (xmax-xmin)/stagCoord->N[0];
112:   h[1] = (ymax-ymin)/stagCoord->N[1];
113:   h[2] = (zmax-zmin)/stagCoord->N[2];

115:   for(ind[2]=start[2]; ind[2]<start[2] + n[2] + nExtra[2]; ++ind[2]) {
116:     for(ind[1]=start[1]; ind[1]<start[1] + n[1] + nExtra[1]; ++ind[1]) {
117:       for(ind[0]=start[0]; ind[0]<start[0] + n[0] + nExtra[0]; ++ind[0]) {
118:         if (stagCoord->dof[0]) {
119:           const PetscReal offs[3] = {0.0,0.0,0.0};
120:           for (c=0; c<3; ++c) {
121:             arr[ind[2]][ind[1]][ind[0]][ibackdownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
122:           }
123:         }
124:         if (stagCoord->dof[1]) {
125:           const PetscReal offs[3] = {0.5,0.0,0.0};
126:           for (c=0; c<3; ++c) {
127:             arr[ind[2]][ind[1]][ind[0]][ibackdown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
128:           }
129:         }
130:         if (stagCoord->dof[1]) {
131:           const PetscReal offs[3] = {0.0,0.5,0.0};
132:           for (c=0; c<3; ++c) {
133:             arr[ind[2]][ind[1]][ind[0]][ibackleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
134:           }
135:         }
136:         if (stagCoord->dof[2]) {
137:           const PetscReal offs[3] = {0.5,0.5,0.0};
138:           for (c=0; c<3; ++c) {
139:             arr[ind[2]][ind[1]][ind[0]][iback + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
140:           }
141:         }
142:         if (stagCoord->dof[1]) {
143:           const PetscReal offs[3] = {0.0,0.0,0.5};
144:           for (c=0; c<3; ++c) {
145:             arr[ind[2]][ind[1]][ind[0]][idownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
146:           }
147:         }
148:         if (stagCoord->dof[2]) {
149:           const PetscReal offs[3] = {0.5,0.0,0.5};
150:           for (c=0; c<3; ++c) {
151:             arr[ind[2]][ind[1]][ind[0]][idown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
152:           }
153:         }
154:         if (stagCoord->dof[2]) {
155:           const PetscReal offs[3] = {0.0,0.5,0.5};
156:           for (c=0; c<3; ++c) {
157:             arr[ind[2]][ind[1]][ind[0]][ileft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
158:           }
159:         }
160:         if (stagCoord->dof[3]) {
161:           const PetscReal offs[3] = {0.5,0.5,0.5};
162:           for (c=0; c<3; ++c) {
163:             arr[ind[2]][ind[1]][ind[0]][ielement + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
164:           }
165:         }
166:       }
167:     }
168:   }
169:   DMStagVecRestoreArrayDOF(dmCoord,coordLocal,&arr);
170:   DMCreateGlobalVector(dmCoord,&coord);
171:   DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);
172:   DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);
173:   DMSetCoordinates(dm,coord);
174:   PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);
175:   DMRestoreLocalVector(dmCoord,&coordLocal);
176:   VecDestroy(&coord);
177:   return(0);
178: }

180: /* Helper functions used in DMSetUp_Stag() */
181: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM);
182: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM);
183: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM,PetscInt**);
184: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM,const PetscInt*);
185: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM,const PetscInt*);
186: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM);

188: PETSC_INTERN PetscErrorCode DMSetUp_Stag_3d(DM dm)
189: {
190:   PetscErrorCode  ierr;
191:   DM_Stag * const stag = (DM_Stag*)dm->data;
192:   PetscMPIInt     rank;
193:   PetscInt        i,j,d;
194:   PetscInt        *globalOffsets;
195:   const PetscInt  dim = 3;

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

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

203:   /* Determine location of rank in grid */
204:     stag->rank[0] = rank % stag->nRanks[0];
205:     stag->rank[1] = rank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0];
206:     stag->rank[2] = rank / (stag->nRanks[0] * stag->nRanks[1]);
207:     for(d=0; d<dim; ++d) {
208:       stag->firstRank[d] = PetscNot(stag->rank[d]);
209:       stag->lastRank[d]  = (PetscBool)(stag->rank[d] == stag->nRanks[d]-1);
210:     }

212:   /* Determine locally owned region (if not already prescribed).
213:    Divide equally, giving lower ranks in each dimension and extra element if needbe.
214:    Note that this uses O(P) storage. If this ever becomes an issue, this could
215:    be refactored to not keep this data around.  */
216:   for (i=0; i<dim; ++i) {
217:     if (!stag->l[i]) {
218:       const PetscInt Ni = stag->N[i], nRanksi = stag->nRanks[i];
219:       PetscMalloc1(stag->nRanks[i],&stag->l[i]);
220:       for (j=0; j<stag->nRanks[i]; ++j){
221:         stag->l[i][j] = Ni/nRanksi + ((Ni % nRanksi) > j);
222:       }
223:     }
224:   }

226:   /* Retrieve local size in stag->n */
227:   for (i=0; i<dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
228: #if defined(PETSC_USE_DEBUG)
229:   for (i=0; i<dim; ++i) {
230:     PetscInt Ncheck,j;
231:     Ncheck = 0;
232:     for (j=0; j<stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
233:     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]);
234:   }
235: #endif

237:   /* Compute starting elements */
238:   for (i=0; i<dim; ++i) {
239:     stag->start[i] = 0;
240:     for(j=0;j<stag->rank[i];++j){
241:       stag->start[i] += stag->l[i][j];
242:     }
243:   }

245:   /* Determine ranks of neighbors */
246:   DMStagSetUpBuildNeighbors_3d(dm);

248:   /* Define useful sizes */
249:   {
250:     PetscInt entriesPerEdge,entriesPerFace,entriesPerCorner,entriesPerElementRow,entriesPerElementLayer;
251:     PetscBool dummyEnd[3];
252:     for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
253:     stag->entriesPerElement = stag->dof[0] + 3*stag->dof[1] + 3*stag->dof[2] + stag->dof[3];
254:     entriesPerFace          = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
255:     entriesPerEdge          = stag->dof[0] + stag->dof[1];
256:     entriesPerCorner        = stag->dof[0];
257:     entriesPerElementRow    = stag->n[0]*stag->entriesPerElement
258:                               + (dummyEnd[0]                               ? entriesPerFace                       : 0);
259:     entriesPerElementLayer  = stag->n[1]*entriesPerElementRow
260:                               + (dummyEnd[1]                               ? stag->n[0] * entriesPerFace          : 0)
261:                               + (dummyEnd[1] && dummyEnd[0]                ? entriesPerEdge                       : 0);
262:     stag->entries           = stag->n[2]*entriesPerElementLayer
263:                               + (dummyEnd[2]                               ? stag->n[0]*stag->n[1]*entriesPerFace : 0)
264:                               + (dummyEnd[2] && dummyEnd[0]                ? stag->n[1]*entriesPerEdge            : 0)
265:                               + (dummyEnd[2] && dummyEnd[1]                ? stag->n[0]*entriesPerEdge            : 0)
266:                               + (dummyEnd[2] && dummyEnd[1] && dummyEnd[0] ? entriesPerCorner                     : 0);
267:   }

269:   /* Check that we will not overflow 32 bit indices (slightly overconservative) */
270: #if !defined(PETSC_USE_64BIT_INDICES)
271:   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);
272: #endif

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

276:     This again requires O(P) storage, which could be replaced with some global
277:     communication.
278:   */
279:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsets);

281:   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");

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

347:   /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping */
348:   DMStagSetUpBuildScatter_3d(dm,globalOffsets);
349:   DMStagSetUpBuildL2G_3d(dm,globalOffsets);

351:   /* In special cases, create a dedicated injective local-to-global map */
352:   if ((stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) ||
353:       (stag->boundaryType[1] == DM_BOUNDARY_PERIODIC && stag->nRanks[1] == 1) ||
354:       (stag->boundaryType[2] == DM_BOUNDARY_PERIODIC && stag->nRanks[2] == 1)) {
355:     DMStagPopulateLocalToGlobalInjective(dm);
356:   }

358:   /* Free global offsets */
359:   PetscFree(globalOffsets);

361:   /* Precompute location offsets */
362:   DMStagComputeLocationOffsets_3d(dm);

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

367:   return(0);
368: }

370: /* adapted from da3.c */
371: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM dm)
372: {
373:   PetscErrorCode  ierr;
374:   PetscMPIInt     rank,size;
375:   PetscInt        m,n,p,pm;
376:   DM_Stag * const stag = (DM_Stag*)dm->data;
377:   const PetscInt M = stag->N[0];
378:   const PetscInt N = stag->N[1];
379:   const PetscInt P = stag->N[2];

382:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
383:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

385:   m = stag->nRanks[0];
386:   n = stag->nRanks[1];
387:   p = stag->nRanks[2];

389:   if (m != PETSC_DECIDE) {
390:     if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
391:     else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
392:   }
393:   if (n != PETSC_DECIDE) {
394:     if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
395:     else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
396:   }
397:   if (p != PETSC_DECIDE) {
398:     if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
399:     else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
400:   }
401:   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);

403:   /* Partition the array among the processors */
404:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
405:     m = size/(n*p);
406:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
407:     n = size/(m*p);
408:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
409:     p = size/(m*n);
410:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
411:     /* try for squarish distribution */
412:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
413:     if (!m) m = 1;
414:     while (m > 0) {
415:       n = size/(m*p);
416:       if (m*n*p == size) break;
417:       m--;
418:     }
419:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
420:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
421:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
422:     /* try for squarish distribution */
423:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
424:     if (!m) m = 1;
425:     while (m > 0) {
426:       p = size/(m*n);
427:       if (m*n*p == size) break;
428:       m--;
429:     }
430:     if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
431:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
432:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
433:     /* try for squarish distribution */
434:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
435:     if (!n) n = 1;
436:     while (n > 0) {
437:       p = size/(m*n);
438:       if (m*n*p == size) break;
439:       n--;
440:     }
441:     if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
442:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
443:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
444:     /* try for squarish distribution */
445:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
446:     if (!n) n = 1;
447:     while (n > 0) {
448:       pm = size/n;
449:       if (n*pm == size) break;
450:       n--;
451:     }
452:     if (!n) n = 1;
453:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
454:     if (!m) m = 1;
455:     while (m > 0) {
456:       p = size/(m*n);
457:       if (m*n*p == size) break;
458:       m--;
459:     }
460:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
461:   } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

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

468:   stag->nRanks[0] = m;
469:   stag->nRanks[1] = n;
470:   stag->nRanks[2] = p;
471:   return(0);
472: }

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

476:         n24 n25 n26
477:         n21 n22 n23
478:         n18 n19 n20 (Front, bigger z)

480:         n15 n16 n17
481:         n12     n14   ^ y
482:         n9  n10 n11   |
483:                       +--> x
484:         n6  n7  n8
485:         n3  n4  n5
486:         n0  n1  n2 (Back, smaller z) */
487: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM dm)
488: {
489:   PetscErrorCode  ierr;
490:   DM_Stag * const stag = (DM_Stag*)dm->data;
491:   PetscInt        d,i;
492:   PetscBool       per[3],first[3],last[3];
493:   PetscInt        neighborRank[27][3],r[3],n[3];
494:   const PetscInt  dim = 3;

497:   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]]);

499:   /* Assemble some convenience variables */
500:   for (d=0; d<dim; ++d) {
501:     per[d]   = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
502:     first[d] = stag->firstRank[d];
503:     last[d]  = stag->lastRank[d];
504:     r[d]     = stag->rank[d];
505:     n[d]     = stag->nRanks[d];
506:   }

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

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

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

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

522:   neighborRank[3][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left       back  */
523:   neighborRank[3][1]  =                                      r[1]    ;
524:   neighborRank[3][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

526:   neighborRank[4][0]  =                                      r[0]    ; /*            back  */
527:   neighborRank[4][1]  =                                      r[1]    ;
528:   neighborRank[4][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

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

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

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

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

546:   neighborRank[9][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down       */
547:   neighborRank[9][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
548:   neighborRank[9][2]  =                                      r[2]    ;

550:   neighborRank[10][0] =                                      r[0]    ; /*       down       */
551:   neighborRank[10][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
552:   neighborRank[10][2] =                                      r[2]    ;

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

558:   neighborRank[12][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left             */
559:   neighborRank[12][1] =                                      r[1]    ;
560:   neighborRank[12][2] =                                      r[2]    ;

562:   neighborRank[13][0] =                                      r[0]    ; /*                  */
563:   neighborRank[13][1] =                                      r[1]    ;
564:   neighborRank[13][2] =                                      r[2]    ;

566:   neighborRank[14][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right            */
567:   neighborRank[14][1] =                                      r[1]    ;
568:   neighborRank[14][2] =                                      r[2]    ;

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

574:   neighborRank[16][0] =                                      r[0]    ; /*       up         */
575:   neighborRank[16][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
576:   neighborRank[16][2] =                                      r[2]    ;

578:   neighborRank[17][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right up         */
579:   neighborRank[17][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
580:   neighborRank[17][2] =                                      r[2]    ;

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

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

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

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

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

602:   neighborRank[23][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right      front */
603:   neighborRank[23][1] =                                      r[1]    ;
604:   neighborRank[23][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

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

610:   neighborRank[25][0] =                                      r[0]    ; /*       up   front */
611:   neighborRank[25][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
612:   neighborRank[25][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

614:   neighborRank[26][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right up   front */
615:   neighborRank[26][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
616:   neighborRank[26][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

618:   /* Then, compute the rank of each in the linear ordering */
619:   PetscMalloc1(27,&stag->neighbors);
620:   for (i=0; i<27; ++i){
621:     if  (neighborRank[i][0] >= 0 && neighborRank[i][1] >=0 && neighborRank[i][2] >=0) {
622:       stag->neighbors[i] = neighborRank[i][0] + n[0]*neighborRank[i][1] + n[0]*n[1]*neighborRank[i][2];
623:     } else {
624:       stag->neighbors[i] = -1;
625:     }
626:   }

628:   return(0);
629: }

631: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM dm,PetscInt **pGlobalOffsets)
632: {
633:   PetscErrorCode        ierr;
634:   const DM_Stag * const stag = (DM_Stag*)dm->data;
635:   PetscInt              *globalOffsets;
636:   PetscInt              i,j,k,d,entriesPerEdge,entriesPerFace,count;
637:   PetscMPIInt           size;
638:   PetscBool             extra[3];

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

729:   return(0);
730: }

732: /* A helper function to reduce code duplication as we loop over various ranges.
733:    i,j,k refer to element numbers on the rank where an element lives in the global
734:    representation (without ghosts) to be offset by a global offset per rank.
735:    ig,jg,kg refer to element numbers in the local (this rank) ghosted numbering.
736:    Note that this function could easily be converted to a macro (it does not compute
737:    anything except loop indices and the values of idxGlobal and idxLocal).  */
738: 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)
739: {
740:   PetscInt ig,jg,kg,d,c;

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

870: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM dm,const PetscInt *globalOffsets)
871: {
873:   DM_Stag * const stag = (DM_Stag*)dm->data;
874:   PetscInt       d,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerCorner,entriesPerEdge,entriesPerFace,entriesToTransferTotal,count,eprGhost,eplGhost;
875:   PetscInt       *idxLocal,*idxGlobal;
876:   PetscMPIInt    rank;
877:   PetscInt       nNeighbors[27][3];
878:   PetscBool      star,nextToDummyEnd[3],dummyStart[3],dummyEnd[3];

881:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
882:   if (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth || stag->n[2] < stag->stencilWidth) {
883:     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);
884:   }

886:   /* Check stencil type */
887:   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]);
888:   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");
889:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

891:   /* Compute numbers of elements on each neighbor */
892:   {
893:     PetscInt i;
894:     for (i=0; i<27; ++i) {
895:       const PetscInt neighborRank = stag->neighbors[i];
896:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
897:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
898:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
899:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
900:       } /* else leave uninitialized - error if accessed */
901:     }
902:   }

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

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

923:   /* Convenience variables */
924:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
925:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
926:   entriesPerCorner                    = stag->dof[0];
927:   for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
928:   eprGhost                            = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row   */
929:   eplGhost                            = stag->nGhost[1]*eprGhost;                /* epl = entries per (element) layer */

931:   /* Compute the number of local entries which correspond to any global entry */
932:   {
933:     PetscInt nNonDummyGhost[3];
934:     for (d=0; d<3; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
935:     if (star) {
936:       entriesToTransferTotal = (
937:           nNonDummyGhost[0] * stag->n[1]        * stag->n[2]
938:         + stag->n[0]        * nNonDummyGhost[1] * stag->n[2]
939:         + stag->n[0]        * stag->n[1]        * nNonDummyGhost[2]
940:         - 2 * (stag->n[0] * stag->n[1] * stag->n[2])
941:         ) * stag->entriesPerElement
942:         + (dummyEnd[0]                               ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace   : 0)
943:         + (dummyEnd[1]                               ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace   : 0)
944:         + (dummyEnd[2]                               ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace   : 0)
945:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
946:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
947:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
948:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
949:     } else {
950:       entriesToTransferTotal  =  nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
951:         + (dummyEnd[0]                               ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace   : 0)
952:         + (dummyEnd[1]                               ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace   : 0)
953:         + (dummyEnd[2]                               ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace   : 0)
954:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
955:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
956:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
957:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
958:     }
959:   }

961:   /* Allocate arrays to populate */
962:   PetscMalloc1(entriesToTransferTotal,&idxLocal);
963:   PetscMalloc1(entriesToTransferTotal,&idxGlobal);

965:   /* Counts into idxLocal/idxGlobal */
966:   count = 0;

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

971:       1. The loop bounds (i/ighost,j/jghost,k/kghost)
972:       2. the strides used in the loop body, which depend on whether the current
973:       rank and/or its neighbor are a non-periodic right/top/front boundary, which has additional
974:       points in the global representation.
975:       - If the neighboring rank is a right/top/front boundary,
976:       then eprNeighbor (entries per element row on the neighbor) and/or eplNeighbor (entries per element layer on the neighbor)
977:       are different.
978:       - If this rank is a non-periodic right/top/front boundary (checking entries of stag->lastRank),
979:       there is an extra loop over 1-3 boundary faces)

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

986:   /* LEFT DOWN BACK */
987:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
988:     const PetscInt  neighbor     = 0;
989:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
990:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
991:     const PetscInt  epFaceRow    = -1;
992:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
993:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
994:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
995:     const PetscInt  startGhost0  = 0;
996:     const PetscInt  startGhost1  = 0;
997:     const PetscInt  startGhost2  = 0;
998:     const PetscInt  endGhost0    = ghostOffsetStart[0];
999:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1000:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1001:     const PetscBool extra0       = PETSC_FALSE;
1002:     const PetscBool extra1       = PETSC_FALSE;
1003:     const PetscBool extra2       = PETSC_FALSE;
1004:     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);
1005:   }

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

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

1049:   /* LEFT BACK */
1050:   if (!star && !dummyStart[0] && !dummyStart[2]) {
1051:     const PetscInt  neighbor     = 3;
1052:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1053:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* May be a top boundary */
1054:     const PetscInt  epFaceRow    = -1;
1055:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1056:     const PetscInt  start1       = 0;
1057:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1058:     const PetscInt  startGhost0  = 0;
1059:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1060:     const PetscInt  startGhost2  = 0;
1061:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1062:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1063:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1064:     const PetscBool extra0       = PETSC_FALSE;
1065:     const PetscBool extra1       = dummyEnd[1];
1066:     const PetscBool extra2       = PETSC_FALSE;
1067:     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);
1068:   }

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

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

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

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

1154:   /* RIGHT UP BACK */
1155:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyStart[2]) {
1156:     const PetscInt  neighbor     = 8;
1157:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1158:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1159:     const PetscInt  epFaceRow    = -1;
1160:     const PetscInt  start0       = 0;
1161:     const PetscInt  start1       = 0;
1162:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1163:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1164:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1165:     const PetscInt  startGhost2  = 0;
1166:     const PetscInt  endGhost0    = stag->nGhost[0];
1167:     const PetscInt  endGhost1    = stag->nGhost[1];
1168:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1169:     const PetscBool extra0       = PETSC_FALSE;
1170:     const PetscBool extra1       = PETSC_FALSE;
1171:     const PetscBool extra2       = PETSC_FALSE;
1172:     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);
1173:   }

1175:   /* LEFT DOWN */
1176:   if (!star && !dummyStart[0] && !dummyStart[1]) {
1177:     const PetscInt  neighbor     = 9;
1178:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1179:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1180:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1181:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1182:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1183:     const PetscInt  start2       = 0;
1184:     const PetscInt  startGhost0  = 0;
1185:     const PetscInt  startGhost1  = 0;
1186:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1187:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1188:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1189:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1190:     const PetscBool extra0       = PETSC_FALSE;
1191:     const PetscBool extra1       = PETSC_FALSE;
1192:     const PetscBool extra2       = dummyEnd[2];
1193:     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);
1194:   }

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

1217:   /* RIGHT DOWN */
1218:   if (!star && !dummyEnd[0] && !dummyStart[1]) {
1219:     const PetscInt  neighbor     = 11;
1220:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1221:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1222:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1223:     const PetscInt  start0       = 0;
1224:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1225:     const PetscInt  start2       = 0;
1226:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1227:     const PetscInt  startGhost1  = 0;
1228:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1229:     const PetscInt  endGhost0    = stag->nGhost[0];
1230:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1231:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1232:     const PetscBool extra0       = PETSC_FALSE;
1233:     const PetscBool extra1       = PETSC_FALSE;
1234:     const PetscBool extra2       = dummyEnd[2];
1235:     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);
1236:   }

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

1259:   /* (HERE) */
1260:   {
1261:     const PetscInt  neighbor     = 13;
1262:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1263:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1264:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
1265:     const PetscInt  start0       = 0;
1266:     const PetscInt  start1       = 0;
1267:     const PetscInt  start2       = 0;
1268:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1269:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1270:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1271:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1272:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1273:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1274:     const PetscBool extra0       = dummyEnd[0];
1275:     const PetscBool extra1       = dummyEnd[1];
1276:     const PetscBool extra2       = dummyEnd[2];
1277:     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);
1278:   }

1280:   /* RIGHT */
1281:   if (!dummyEnd[0]) {
1282:     const PetscInt  neighbor     = 14;
1283:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1284:     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 */
1285:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1286:     const PetscInt  start0       = 0;
1287:     const PetscInt  start1       = 0;
1288:     const PetscInt  start2       = 0;
1289:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1290:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1291:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1292:     const PetscInt  endGhost0    = stag->nGhost[0];
1293:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1294:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1295:     const PetscBool extra0       = PETSC_FALSE;
1296:     const PetscBool extra1       = dummyEnd[1];
1297:     const PetscBool extra2       = dummyEnd[2];
1298:     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);
1299:   }

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

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

1343:   /* RIGHT UP */
1344:   if (!star && !dummyEnd[0] && !dummyEnd[1]) {
1345:     const PetscInt  neighbor     = 17;
1346:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1347:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1348:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1349:     const PetscInt  start0       = 0;
1350:     const PetscInt  start1       = 0;
1351:     const PetscInt  start2       = 0;
1352:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1353:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1354:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1355:     const PetscInt  endGhost0    = stag->nGhost[0];
1356:     const PetscInt  endGhost1    = stag->nGhost[1];
1357:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1358:     const PetscBool extra0       = PETSC_FALSE;
1359:     const PetscBool extra1       = PETSC_FALSE;
1360:     const PetscBool extra2       = dummyEnd[2];
1361:     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);
1362:   }

1364:   /* LEFT DOWN FRONT */
1365:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1366:     const PetscInt  neighbor     = 18;
1367:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1368:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1369:     const PetscInt  epFaceRow    = -1;
1370:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1371:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1372:     const PetscInt  start2       = 0;
1373:     const PetscInt  startGhost0  = 0;
1374:     const PetscInt  startGhost1  = 0;
1375:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1376:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1377:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1378:     const PetscInt  endGhost2    = stag->nGhost[2];
1379:     const PetscBool extra0       = PETSC_FALSE;
1380:     const PetscBool extra1       = PETSC_FALSE;
1381:     const PetscBool extra2       = PETSC_FALSE;
1382:     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);
1383:   }

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

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

1427:   /* LEFT FRONT */
1428:   if (!star && !dummyStart[0] && !dummyEnd[2]) {
1429:     const PetscInt  neighbor     = 21;
1430:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1431:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1432:     const PetscInt  epFaceRow    = -1;
1433:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1434:     const PetscInt  start1       = 0;
1435:     const PetscInt  start2       = 0;
1436:     const PetscInt  startGhost0  = 0;
1437:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1438:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1439:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1440:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1441:     const PetscInt  endGhost2    = stag->nGhost[2];
1442:     const PetscBool extra0       = PETSC_FALSE;
1443:     const PetscBool extra1       = dummyEnd[1];
1444:     const PetscBool extra2       = PETSC_FALSE;
1445:     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);
1446:   }

1448:   /* FRONT */
1449:   if (!dummyEnd[2]) {
1450:     const PetscInt  neighbor     = 22;
1451:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1452:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1453:     const PetscInt  epFaceRow    = -1;
1454:     const PetscInt  start0       = 0;
1455:     const PetscInt  start1       = 0;
1456:     const PetscInt  start2       = 0;
1457:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1458:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1459:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1460:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1461:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1462:     const PetscInt  endGhost2    = stag->nGhost[2];
1463:     const PetscBool extra0       = dummyEnd[0];
1464:     const PetscBool extra1       = dummyEnd[1];
1465:     const PetscBool extra2       = PETSC_FALSE;
1466:     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);
1467:   }

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

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

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

1532:   /* RIGHT UP FRONT */
1533:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyEnd[2]) {
1534:     const PetscInt  neighbor     = 26;
1535:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1536:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1537:     const PetscInt  epFaceRow    = -1;
1538:     const PetscInt  start0       = 0;
1539:     const PetscInt  start1       = 0;
1540:     const PetscInt  start2       = 0;
1541:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1542:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1543:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1544:     const PetscInt  endGhost0    = stag->nGhost[0];
1545:     const PetscInt  endGhost1    = stag->nGhost[1];
1546:     const PetscInt  endGhost2    = stag->nGhost[2];
1547:     const PetscBool extra0       = PETSC_FALSE;
1548:     const PetscBool extra1       = PETSC_FALSE;
1549:     const PetscBool extra2       = PETSC_FALSE;
1550:     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);
1551:   }

1553: #if defined(PETSC_USE_DEBUG)
1554:   if (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);
1555: #endif

1557:   /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
1558:   {
1559:     Vec local,global;
1560:     IS isLocal,isGlobal;
1561:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
1562:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
1563:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
1564:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
1565:     VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
1566:     VecDestroy(&global);
1567:     VecDestroy(&local);
1568:     ISDestroy(&isLocal);  /* frees idxLocal */
1569:     ISDestroy(&isGlobal); /* free idxGlobal */
1570:   }

1572:   return(0);
1573: }

1575: /* Note: this function assumes that DMBoundary types of none, ghosted, and periodic are the only ones of interest.
1576: Adding support for others should be done very carefully.  */
1577: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM dm,const PetscInt *globalOffsets)
1578: {
1579:   PetscErrorCode        ierr;
1580:   const DM_Stag * const stag = (DM_Stag*)dm->data;
1581:   PetscInt              *idxGlobalAll;
1582:   PetscInt              d,count,ighost,jghost,kghost,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerFace,entriesPerEdge;
1583:   PetscInt              nNeighbors[27][3],eprNeighbor[27],eplNeighbor[27],globalOffsetNeighbor[27];
1584:   PetscBool             nextToDummyEnd[3],dummyStart[3],dummyEnd[3],star;


1588:   /* Check stencil type */
1589:   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]);
1590:   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");
1591:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

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

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

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

1606:   /* Compute numbers of elements on each neighbor  and offset*/
1607:   {
1608:     PetscInt i;
1609:     for (i=0; i<27; ++i) {
1610:       const PetscInt neighborRank = stag->neighbors[i];
1611:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
1612:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
1613:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
1614:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
1615:         globalOffsetNeighbor[i] = globalOffsets[stag->neighbors[i]];
1616:       } /* else leave uninitialized - error if accessed */
1617:     }
1618:   }

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

1732:   /* Populate idxGlobalAll */
1733:   PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
1734:   count = 0;

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

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

1746:           /* Loop over columns 1/3: Left Back Down*/
1747:           if (!dummyStart[0]) {
1748:             const PetscInt neighbor = 0;
1749:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1750:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1751:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1752:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1753:               }
1754:             }
1755:           } else {
1756:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1757:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1758:                 idxGlobalAll[count] = -1;
1759:               }
1760:             }
1761:           }

1763:           /* Loop over columns 2/3: (Middle) Down Back */
1764:           {
1765:             const PetscInt neighbor = 1;
1766:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1767:               const PetscInt i = ighost - ghostOffsetStart[0];
1768:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1769:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1770:               }
1771:             }
1772:           }

1774:           /* Loop over columns 3/3: Right Down Back */
1775:           if (!dummyEnd[0]) {
1776:             const PetscInt neighbor = 2;
1777:             PetscInt       i;
1778:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1779:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1780:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1781:               }
1782:             }
1783:           } else {
1784:             /* Partial dummy entries on (Middle) Down Back neighbor */
1785:             const PetscInt neighbor          = 1;
1786:             PetscInt       i,dLocal;
1787:             i = stag->n[0];
1788:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1789:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1790:             }
1791:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1792:               idxGlobalAll[count] = -1;
1793:             }
1794:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1795:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1796:             }
1797:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1798:               idxGlobalAll[count] = -1;
1799:             }
1800:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1801:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1802:             }
1803:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1804:               idxGlobalAll[count] = -1;
1805:             }
1806:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1807:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1808:             }
1809:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1810:               idxGlobalAll[count] = -1;
1811:             }
1812:             ++i;
1813:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1814:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1815:                 idxGlobalAll[count] = -1;
1816:               }
1817:             }
1818:           }
1819:         }
1820:       } else {
1821:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1822:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
1823:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
1824:               idxGlobalAll[count] = -1;
1825:             }
1826:           }
1827:         }
1828:       }

1830:       /* Loop over rows 2/3: (Middle) Back */
1831:       {
1832:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
1833:           const PetscInt j = jghost - ghostOffsetStart[1];

1835:           /* Loop over columns 1/3: Left (Middle) Back */
1836:           if (!star && !dummyStart[0]) {
1837:             const PetscInt neighbor = 3;
1838:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1839:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1840:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1841:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1842:               }
1843:             }
1844:           } else {
1845:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1846:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1847:                 idxGlobalAll[count] = -1;
1848:               }
1849:             }
1850:           }

1852:           /* Loop over columns 2/3: (Middle) (Middle) Back */
1853:           {
1854:             const PetscInt neighbor = 4;
1855:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1856:               const PetscInt i = ighost - ghostOffsetStart[0];
1857:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1858:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1859:               }
1860:             }
1861:           }

1863:           /* Loop over columns 3/3: Right (Middle) Back */
1864:           if (!star && !dummyEnd[0]) {
1865:             const PetscInt neighbor = 5;
1866:             PetscInt       i;
1867:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1868:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1869:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1870:               }
1871:             }
1872:           } else if (dummyEnd[0]) {
1873:             /* Partial dummy entries on (Middle) (Middle) Back rank */
1874:             const PetscInt neighbor = 4;
1875:             PetscInt       i,dLocal;
1876:             i = stag->n[0];
1877:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1878:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1879:             }
1880:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1881:               idxGlobalAll[count] = -1;
1882:             }
1883:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1884:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1885:             }
1886:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1887:               idxGlobalAll[count] = -1;
1888:             }
1889:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1890:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1891:             }
1892:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1893:               idxGlobalAll[count] = -1;
1894:             }
1895:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1896:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1897:             }
1898:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1899:               idxGlobalAll[count] = -1;
1900:             }
1901:             ++i;
1902:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1903:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1904:                 idxGlobalAll[count] = -1;
1905:               }
1906:             }
1907:           } else {
1908:             /* Right (Middle) Back dummies */
1909:             PetscInt       i;
1910:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1911:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1912:                 idxGlobalAll[count] = -1;
1913:               }
1914:             }
1915:           }
1916:         }
1917:       }

1919:       /* Loop over rows 3/3: Up Back */
1920:       if (!star && !dummyEnd[1]) {
1921:         PetscInt j;
1922:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
1923:           /* Loop over columns 1/3: Left Up Back*/
1924:           if (!dummyStart[0]) {
1925:             const PetscInt neighbor = 6;
1926:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1927:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1928:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1929:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1930:               }
1931:             }
1932:           } else {
1933:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1934:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1935:                 idxGlobalAll[count] = -1;
1936:               }
1937:             }
1938:           }

1940:           /* Loop over columns 2/3: (Middle) Up Back*/
1941:           {
1942:             const PetscInt neighbor = 7;
1943:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1944:               const PetscInt i = ighost - ghostOffsetStart[0];
1945:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1946:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1947:               }
1948:             }
1949:           }

1951:           /* Loop over columns 3/3: Right Up Back*/
1952:           if (!dummyEnd[0]) {
1953:             const PetscInt neighbor = 8;
1954:             PetscInt       i;
1955:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1956:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1957:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1958:               }
1959:             }
1960:           } else {
1961:             /* Partial dummies on (Middle) Up Back */
1962:             const PetscInt neighbor = 7;
1963:             PetscInt       i,dLocal;
1964:             i = stag->n[0];
1965:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1966:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1967:             }
1968:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1969:               idxGlobalAll[count] = -1;
1970:             }
1971:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1972:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1973:             }
1974:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1975:               idxGlobalAll[count] = -1;
1976:             }
1977:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1978:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1979:             }
1980:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1981:               idxGlobalAll[count] = -1;
1982:             }
1983:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1984:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1985:             }
1986:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1987:               idxGlobalAll[count] = -1;
1988:             }
1989:             ++i;
1990:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1991:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1992:                 idxGlobalAll[count] = -1;
1993:               }
1994:             }
1995:           }
1996:         }
1997:       } else if (dummyEnd[1]) {
1998:         /* Up Back partial dummy row */
1999:         PetscInt j = stag->n[1];

2001:         if (!star && !dummyStart[0]) {
2002:           /* Up Back partial dummy row: Loop over columns 1/3: Left Up Back, on Left (Middle) Back rank */
2003:           const PetscInt neighbor = 3;
2004:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2005:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2006:             PetscInt       dLocal;
2007:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2008:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2009:             }

2011:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2012:               idxGlobalAll[count] = -1;
2013:             }
2014:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2015:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2016:             }
2017:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2018:               idxGlobalAll[count] = -1;
2019:             }
2020:           }
2021:         } else {
2022:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2023:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2024:               idxGlobalAll[count] = -1;
2025:             }
2026:           }
2027:         }

2029:         /* Up Back partial dummy row: Loop over columns 2/3: (Middle) Up Back, on (Middle) (Middle) Back rank */
2030:         {
2031:           const PetscInt neighbor = 4;
2032:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2033:             const PetscInt i = ighost - ghostOffsetStart[0];
2034:             PetscInt       dLocal;
2035:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2036:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2037:             }
2038:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2039:               idxGlobalAll[count] = -1;
2040:             }
2041:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2042:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2043:             }
2044:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2045:               idxGlobalAll[count] = -1;
2046:             }
2047:           }
2048:         }

2050:         /* Up Back partial dummy row: Loop over columns 3/3: Right Up Back, on Right (Middle) Back rank */
2051:         if (!star && !dummyEnd[0]) {
2052:           const PetscInt neighbor = 5;
2053:           PetscInt       i;
2054:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2055:             PetscInt       dLocal;
2056:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2057:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2058:             }
2059:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2060:               idxGlobalAll[count] = -1;
2061:             }
2062:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2063:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2064:             }
2065:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2066:               idxGlobalAll[count] = -1;
2067:             }
2068:           }
2069:         } else if (dummyEnd[0]) {
2070:           /* Up Right Back partial dummy element, on (Middle) (Middle) Back rank */
2071:           const PetscInt neighbor = 4;
2072:           PetscInt       i,dLocal;
2073:           i = stag->n[0];
2074:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2075:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2076:           }
2077:           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 */
2078:             idxGlobalAll[count] = -1;
2079:           }
2080:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2081:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2082:           }
2083:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2084:             idxGlobalAll[count] = -1;
2085:           }
2086:           ++i;
2087:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2088:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2089:               idxGlobalAll[count] = -1;
2090:             }
2091:           }
2092:         } else {
2093:           /* Up Right Back dummies */
2094:           PetscInt i;
2095:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2096:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2097:               idxGlobalAll[count] = -1;
2098:             }
2099:           }
2100:         }
2101:         ++j;
2102:         /* Up Back additional dummy rows */
2103:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2104:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2105:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2106:               idxGlobalAll[count] = -1;
2107:             }
2108:           }
2109:         }
2110:       } else {
2111:         /* Up Back dummy rows */
2112:         PetscInt j;
2113:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2114:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2115:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2116:               idxGlobalAll[count] = -1;
2117:             }
2118:           }
2119:         }
2120:       }
2121:     }
2122:   } else {
2123:     for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
2124:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
2125:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2126:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2127:             idxGlobalAll[count] = -1;
2128:           }
2129:         }
2130:       }
2131:     }
2132:   }

2134:   /* Loop over layers 2/3 : (Middle)  */
2135:   {
2136:     for (kghost = ghostOffsetStart[2]; kghost<stag->nGhost[2]-ghostOffsetEnd[2]; ++kghost) {
2137:       const PetscInt k = kghost - ghostOffsetStart[2];

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

2144:           /* Loop over columns 1/3: Left Down (Middle) */
2145:           if (!star && !dummyStart[0]) {
2146:             const PetscInt neighbor = 9;
2147:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2148:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2149:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2150:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2151:               }
2152:             }
2153:           } else {
2154:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2155:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2156:                 idxGlobalAll[count] = -1;
2157:               }
2158:             }
2159:           }

2161:           /* Loop over columns 2/3: (Middle) Down (Middle) */
2162:           {
2163:             const PetscInt neighbor = 10;
2164:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2165:               const PetscInt i = ighost - ghostOffsetStart[0];
2166:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2167:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2168:               }
2169:             }
2170:           }

2172:           /* Loop over columns 3/3: Right Down (Middle) */
2173:           if (!star && !dummyEnd[0]) {
2174:             const PetscInt neighbor = 11;
2175:             PetscInt       i;
2176:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2177:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2178:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2179:               }
2180:             }
2181:           } else if (dummyEnd[0]) {
2182:             /* Partial dummies on (Middle) Down (Middle) neighbor */
2183:             const PetscInt neighbor = 10;
2184:             PetscInt       i,dLocal;
2185:             i = stag->n[0];
2186:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2187:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2188:             }
2189:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2190:               idxGlobalAll[count] = -1;
2191:             }
2192:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2193:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2194:             }
2195:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2196:               idxGlobalAll[count] = -1;
2197:             }
2198:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2199:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2200:             }
2201:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2202:               idxGlobalAll[count] = -1;
2203:             }
2204:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2205:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2206:             }
2207:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2208:               idxGlobalAll[count] = -1;
2209:             }
2210:             ++i;
2211:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2212:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2213:                 idxGlobalAll[count] = -1;
2214:               }
2215:             }
2216:           } else {
2217:             /* Right Down (Middle) dummies */
2218:             PetscInt       i;
2219:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2220:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2221:                 idxGlobalAll[count] = -1;
2222:               }
2223:             }
2224:           }
2225:         }
2226:       } else {
2227:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2228:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2229:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2230:               idxGlobalAll[count] = -1;
2231:             }
2232:           }
2233:         }
2234:       }

2236:       /* Loop over rows 2/3: (Middle) (Middle) */
2237:       {
2238:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2239:           const PetscInt j = jghost - ghostOffsetStart[1];

2241:           /* Loop over columns 1/3: Left (Middle) (Middle) */
2242:           if (!dummyStart[0]) {
2243:             const PetscInt neighbor = 12;
2244:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2245:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2246:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2247:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2248:               }
2249:             }
2250:           } else {
2251:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2252:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2253:                 idxGlobalAll[count] = -1;
2254:               }
2255:             }
2256:           }

2258:           /* Loop over columns 2/3: (Middle) (Middle) (Middle) aka (Here) */
2259:           {
2260:             const PetscInt neighbor = 13;
2261:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2262:               const PetscInt i = ighost - ghostOffsetStart[0];
2263:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2264:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2265:               }
2266:             }
2267:           }

2269:           /* Loop over columns 3/3: Right (Middle) (Middle) */
2270:           if (!dummyEnd[0]) {
2271:             const PetscInt neighbor = 14;
2272:             PetscInt       i;
2273:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2274:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2275:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2276:               }
2277:             }
2278:           } else {
2279:             /* Partial dummies on this rank */
2280:             const PetscInt neighbor = 13;
2281:             PetscInt       i,dLocal;
2282:             i = stag->n[0];
2283:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2284:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2285:             }
2286:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2287:               idxGlobalAll[count] = -1;
2288:             }
2289:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2290:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2291:             }
2292:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2293:               idxGlobalAll[count] = -1;
2294:             }
2295:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2296:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2297:             }
2298:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2299:               idxGlobalAll[count] = -1;
2300:             }
2301:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2302:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2303:             }
2304:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2305:               idxGlobalAll[count] = -1;
2306:             }
2307:             ++i;
2308:             for(;i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2309:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2310:                 idxGlobalAll[count] = -1;
2311:               }
2312:             }
2313:           }
2314:         }
2315:       }

2317:       /* Loop over rows 3/3: Up (Middle) */
2318:       if (!dummyEnd[1]) {
2319:         PetscInt j;
2320:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2322:           /* Loop over columns 1/3: Left Up (Middle) */
2323:           if (!star && !dummyStart[0]) {
2324:             const PetscInt neighbor = 15;
2325:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2326:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2327:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2328:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2329:               }
2330:             }
2331:           } else {
2332:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2333:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2334:                 idxGlobalAll[count] = -1;
2335:               }
2336:             }
2337:           }

2339:           /* Loop over columns 2/3: (Middle) Up (Middle) */
2340:           {
2341:             const PetscInt neighbor = 16;
2342:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2343:               const PetscInt i = ighost - ghostOffsetStart[0];
2344:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2345:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2346:               }
2347:             }
2348:           }

2350:           /* Loop over columns 3/3: Right Up (Middle) */
2351:           if (!star && !dummyEnd[0]) {
2352:             const PetscInt neighbor = 17;
2353:             PetscInt       i;
2354:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2355:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2356:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2357:               }
2358:             }
2359:           } else if (dummyEnd[0]) {
2360:             /* Partial dummies on (Middle) Up (Middle) rank */
2361:             const PetscInt neighbor = 16;
2362:             PetscInt       i,dLocal;
2363:             i = stag->n[0];
2364:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2365:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2366:             }
2367:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2368:               idxGlobalAll[count] = -1;
2369:             }
2370:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2371:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2372:             }
2373:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2374:               idxGlobalAll[count] = -1;
2375:             }
2376:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2377:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2378:             }
2379:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2380:               idxGlobalAll[count] = -1;
2381:             }
2382:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2383:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2384:             }
2385:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2386:               idxGlobalAll[count] = -1;
2387:             }
2388:             ++i;
2389:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2390:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2391:                 idxGlobalAll[count] = -1;
2392:               }
2393:             }
2394:           } else {
2395:             /* Right Up (Middle) dumies */
2396:             PetscInt       i;
2397:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2398:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2399:                 idxGlobalAll[count] = -1;
2400:               }
2401:             }
2402:           }
2403:         }
2404:       } else {
2405:         /* Up (Middle) partial dummy row */
2406:         PetscInt j = stag->n[1];

2408:         /* Up (Middle) partial dummy row: colums 1/3: Left Up (Middle), on Left (Middle) (Middle) rank */
2409:         if (!dummyStart[0]) {
2410:           const PetscInt neighbor = 12;
2411:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2412:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2413:             PetscInt       dLocal;
2414:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2415:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2416:             }
2417:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2418:               idxGlobalAll[count] = -1;
2419:             }
2420:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2421:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2422:             }
2423:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2424:               idxGlobalAll[count] = -1;
2425:             }
2426:           }
2427:         } else {
2428:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2429:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2430:               idxGlobalAll[count] = -1;
2431:             }
2432:           }
2433:         }

2435:         /* Up (Middle) partial dummy row: columns 2/3: (Middle) Up (Middle), on (Middle) (Middle) (Middle) = here rank */
2436:         {
2437:           const PetscInt neighbor = 13;
2438:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2439:             const PetscInt i = ighost - ghostOffsetStart[0];
2440:             PetscInt       dLocal;
2441:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2442:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2443:             }
2444:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2445:               idxGlobalAll[count] = -1;
2446:             }
2447:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2448:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2449:             }
2450:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2451:               idxGlobalAll[count] = -1;
2452:             }
2453:           }
2454:         }

2456:         if (!dummyEnd[0]) {
2457:           /* Up (Middle) partial dummy row: columns 3/3: Right Up (Middle), on Right (Middle) (Middle) rank */
2458:           const PetscInt neighbor = 14;
2459:           PetscInt       i;
2460:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2461:             PetscInt       dLocal;
2462:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2463:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2464:             }
2465:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2466:               idxGlobalAll[count] = -1;
2467:             }
2468:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2469:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2470:             }
2471:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2472:               idxGlobalAll[count] = -1;
2473:             }
2474:           }
2475:         } else {
2476:           /* Up (Middle) Right partial dummy element, on (Middle) (Middle) (Middle) = here rank */
2477:           const PetscInt neighbor = 13;
2478:           PetscInt       i,dLocal;
2479:           i = stag->n[0];
2480:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2481:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2482:           }
2483:           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 */
2484:             idxGlobalAll[count] = -1;
2485:           }
2486:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2487:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2488:           }
2489:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2490:             idxGlobalAll[count] = -1;
2491:           }
2492:           ++i;
2493:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2494:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2495:               idxGlobalAll[count] = -1;
2496:             }
2497:           }
2498:         }
2499:         ++j;
2500:         /* Up (Middle) additional dummy rows */
2501:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2502:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2503:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2504:               idxGlobalAll[count] = -1;
2505:             }
2506:           }
2507:         }
2508:       }
2509:     }
2510:   }

2512:   /* Loop over layers 3/3 : Front */
2513:   if (!dummyEnd[2]) {
2514:     PetscInt k;
2515:     for (k=0; k<ghostOffsetEnd[2]; ++k) {

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

2522:           /* Loop over columns 1/3: Left Down Front */
2523:           if (!dummyStart[0]) {
2524:             const PetscInt neighbor = 18;
2525:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2526:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2527:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2528:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2529:               }
2530:             }
2531:           } else {
2532:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2533:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2534:                 idxGlobalAll[count] = -1;
2535:               }
2536:             }
2537:           }

2539:           /* Loop over columns 2/3: (Middle) Down Front */
2540:           {
2541:             const PetscInt neighbor = 19;
2542:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2543:               const PetscInt i = ighost - ghostOffsetStart[0];
2544:               for (d=0; d<stag->entriesPerElement; ++d,++count) { /* vertex, 2 edges, and face associated with back face */
2545:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2546:               }
2547:             }
2548:           }

2550:           /* Loop over columns 3/3: Right Down Front */
2551:           if (!dummyEnd[0]) {
2552:             const PetscInt neighbor = 20;
2553:             PetscInt       i;
2554:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2555:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2556:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2557:               }
2558:             }
2559:           } else {
2560:             /* Partial dummy element on (Middle) Down Front rank */
2561:             const PetscInt neighbor = 19;
2562:             PetscInt       i,dLocal;
2563:             i = stag->n[0];
2564:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2565:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2566:             }
2567:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2568:               idxGlobalAll[count] = -1;
2569:             }
2570:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2571:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2572:             }
2573:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2574:               idxGlobalAll[count] = -1;
2575:             }
2576:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2577:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2578:             }
2579:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2580:               idxGlobalAll[count] = -1;
2581:             }
2582:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2583:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2584:             }
2585:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2586:               idxGlobalAll[count] = -1;
2587:             }
2588:             ++i;
2589:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2590:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2591:                 idxGlobalAll[count] = -1;
2592:               }
2593:             }
2594:           }
2595:         }
2596:       } else {
2597:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2598:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2599:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2600:               idxGlobalAll[count] = -1;
2601:             }
2602:           }
2603:         }
2604:       }

2606:       /* Loop over rows 2/3: (Middle) Front */
2607:       {
2608:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2609:           const PetscInt j = jghost - ghostOffsetStart[1];

2611:           /* Loop over columns 1/3: Left (Middle) Front*/
2612:           if (!star && !dummyStart[0]) {
2613:             const PetscInt neighbor = 21;
2614:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2615:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2616:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2617:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2618:               }
2619:             }
2620:           } else {
2621:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2622:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2623:                 idxGlobalAll[count] = -1;
2624:               }
2625:             }
2626:           }

2628:           /* Loop over columns 2/3: (Middle) (Middle) Front*/
2629:           {
2630:             const PetscInt neighbor = 22;
2631:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2632:               const PetscInt i = ighost - ghostOffsetStart[0];
2633:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2634:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2635:               }
2636:             }
2637:           }

2639:           /* Loop over columns 3/3: Right (Middle) Front*/
2640:           if (!star && !dummyEnd[0]) {
2641:             const PetscInt neighbor = 23;
2642:             PetscInt       i;
2643:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2644:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2645:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2646:               }
2647:             }
2648:           } else if (dummyEnd[0]) {
2649:             /* Partial dummy element on (Middle) (Middle) Front element */
2650:             const PetscInt neighbor = 22;
2651:             PetscInt       i,dLocal;
2652:             i = stag->n[0];
2653:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2654:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2655:             }
2656:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2657:               idxGlobalAll[count] = -1;
2658:             }
2659:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2660:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2661:             }
2662:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2663:               idxGlobalAll[count] = -1;
2664:             }
2665:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2666:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2667:             }
2668:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2669:               idxGlobalAll[count] = -1;
2670:             }
2671:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2672:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2673:             }
2674:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2675:               idxGlobalAll[count] = -1;
2676:             }
2677:             ++i;
2678:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2679:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2680:                 idxGlobalAll[count] = -1;
2681:               }
2682:             }
2683:           } else {
2684:             /* Right (Middle) Front dummies */
2685:             PetscInt       i;
2686:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2687:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2688:                 idxGlobalAll[count] = -1;
2689:               }
2690:             }
2691:           }
2692:         }
2693:       }

2695:       /* Loop over rows 3/3: Up Front */
2696:       if (!star && !dummyEnd[1]) {
2697:         PetscInt j;
2698:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2700:           /* Loop over columns 1/3: Left Up Front */
2701:           if (!dummyStart[0]) {
2702:             const PetscInt neighbor = 24;
2703:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2704:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2705:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2706:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2707:               }
2708:             }
2709:           } else {
2710:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2711:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2712:                 idxGlobalAll[count] = -1;
2713:               }
2714:             }
2715:           }

2717:           /* Loop over columns 2/3: (Middle) Up Front */
2718:           {
2719:             const PetscInt neighbor = 25;
2720:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2721:               const PetscInt i = ighost - ghostOffsetStart[0];
2722:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2723:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2724:               }
2725:             }
2726:           }

2728:           /* Loop over columns 3/3: Right Up Front */
2729:           if (!dummyEnd[0]) {
2730:             const PetscInt neighbor = 26;
2731:             PetscInt     i;
2732:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2733:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2734:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2735:               }
2736:             }
2737:           } else {
2738:             /* Partial dummy element on (Middle) Up Front rank */
2739:             const PetscInt neighbor = 25;
2740:             PetscInt       i,dLocal;
2741:             i = stag->n[0];
2742:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2743:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2744:             }
2745:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2746:               idxGlobalAll[count] = -1;
2747:             }
2748:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2749:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2750:             }
2751:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2752:               idxGlobalAll[count] = -1;
2753:             }
2754:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2755:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2756:             }
2757:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2758:               idxGlobalAll[count] = -1;
2759:             }
2760:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2761:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2762:             }
2763:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2764:               idxGlobalAll[count] = -1;
2765:             }
2766:             ++i;
2767:             for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2768:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2769:                 idxGlobalAll[count] = -1;
2770:               }
2771:             }
2772:           }
2773:         }
2774:       } else if (dummyEnd[1]) {
2775:         /* Up Front partial dummy row */
2776:         PetscInt j = stag->n[1];

2778:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) Front rank */
2779:         if (!star && !dummyStart[0]) {
2780:           const PetscInt neighbor = 21;
2781:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2782:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2783:             PetscInt       dLocal;
2784:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2785:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2786:             }
2787:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2788:               idxGlobalAll[count] = -1;
2789:             }
2790:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2791:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2792:             }
2793:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2794:               idxGlobalAll[count] = -1;
2795:             }
2796:           }
2797:         } else {
2798:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2799:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2800:               idxGlobalAll[count] = -1;
2801:             }
2802:           }
2803:         }

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

2826:         if (!star && !dummyEnd[0]) {
2827:           /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) Front rank */
2828:           const PetscInt neighbor = 23;
2829:           PetscInt       i;
2830:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2831:             PetscInt       dLocal;
2832:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
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 face and back left edge */
2836:               idxGlobalAll[count] = -1;
2837:             }
2838:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
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 left face and element */
2842:               idxGlobalAll[count] = -1;
2843:             }
2844:           }

2846:         } else if (dummyEnd[0]) {
2847:           /* Partial Right Up Front dummy element, on (Middle) (Middle) Front rank */
2848:           const PetscInt neighbor = 22;
2849:           PetscInt       i,dLocal;
2850:           i = stag->n[0];
2851:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2852:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2853:           }
2854:           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 */
2855:             idxGlobalAll[count] = -1;
2856:           }
2857:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2858:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2859:           }
2860:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2861:             idxGlobalAll[count] = -1;
2862:           }
2863:           ++i;
2864:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2865:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2866:               idxGlobalAll[count] = -1;
2867:             }
2868:           }
2869:         } else {
2870:           /* Right Up Front dummies */
2871:           PetscInt i;
2872:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2873:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2874:               idxGlobalAll[count] = -1;
2875:             }
2876:           }
2877:         }
2878:         ++j;
2879:         /* Up Front additional dummy rows */
2880:         for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2881:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2882:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2883:               idxGlobalAll[count] = -1;
2884:             }
2885:           }
2886:         }
2887:       } else {
2888:         /* Up Front dummies */
2889:         PetscInt j;
2890:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2891:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2892:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2893:               idxGlobalAll[count] = -1;
2894:             }
2895:           }
2896:         }
2897:       }
2898:     }
2899:   } else {
2900:     PetscInt k = stag->n[2];

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

2907:         /* Down Front partial ghost row: colums 1/3: Left Down Front, on  Left Down (Middle) */
2908:         if (!star && !dummyStart[0]) {
2909:           const PetscInt neighbor = 9;
2910:           const PetscInt epFaceRow         = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
2911:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2912:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2913:             PetscInt  dLocal;
2914:             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 */
2915:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2916:             }
2917:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2918:               idxGlobalAll[count] = -1;
2919:             }
2920:           }
2921:         } else {
2922:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2923:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2924:               idxGlobalAll[count] = -1;
2925:             }
2926:           }
2927:         }

2929:         /* Down Front partial ghost row: columns 2/3: (Middle) Down Front, on (Middle) Down (Middle) */
2930:         {
2931:           const PetscInt neighbor = 10;
2932:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2933:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2934:             const PetscInt i = ighost - ghostOffsetStart[0];
2935:             PetscInt  dLocal;
2936:             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 */
2937:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2938:             }
2939:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2940:               idxGlobalAll[count] = -1;
2941:             }
2942:           }
2943:         }

2945:         if (!star && !dummyEnd[0]) {
2946:           /* Down Front partial dummy row: columns 3/3: Right Down Front, on Right Down (Middle) */
2947:           const PetscInt neighbor = 11;
2948:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
2949:           PetscInt       i;
2950:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2951:             PetscInt  dLocal;
2952:             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 */
2953:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2954:             }
2955:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2956:               idxGlobalAll[count] = -1;
2957:             }
2958:           }
2959:         } else if (dummyEnd[0]) {
2960:           /* Right Down Front partial dummy element, living on (Middle) Down (Middle) rank */
2961:           const PetscInt neighbor = 10;
2962:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2963:           PetscInt       i,dLocal;
2964:           i = stag->n[0];
2965:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2966:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2967:           }
2968:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2969:             idxGlobalAll[count] = -1;
2970:           }
2971:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
2972:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2973:           }
2974:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
2975:             idxGlobalAll[count] = -1;
2976:           }
2977:           ++i;
2978:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2979:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2980:               idxGlobalAll[count] = -1;
2981:             }
2982:           }
2983:         } else {
2984:           PetscInt i;
2985:           /* Right Down Front dummies */
2986:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
2987:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2988:               idxGlobalAll[count] = -1;
2989:             }
2990:           }
2991:         }
2992:       }
2993:     } else {
2994:       for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2995:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2996:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2997:             idxGlobalAll[count] = -1;
2998:           }
2999:         }
3000:       }
3001:     }

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

3008:         /* (Middle) Front partial dummy row: columns 1/3: Left (Middle) Front, on Left (Middle) (Middle) */
3009:         if (!dummyStart[0]) {
3010:           const PetscInt neighbor = 12;
3011:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3012:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3013:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3014:             PetscInt  dLocal;
3015:             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 */
3016:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3017:             }
3018:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3019:               idxGlobalAll[count] = -1;
3020:             }
3021:           }
3022:         } else {
3023:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3024:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3025:               idxGlobalAll[count] = -1;
3026:             }
3027:           }
3028:         }

3030:         /* (Middle) Front partial dummy row: columns 2/3: (Middle) (Middle) Front, on (Middle) (Middle) (Middle) */
3031:         {
3032:           const PetscInt neighbor = 13;
3033:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3034:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3035:             const PetscInt i = ighost - ghostOffsetStart[0];
3036:             PetscInt  dLocal;
3037:             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 */
3038:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3039:             }
3040:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3041:               idxGlobalAll[count] = -1;
3042:             }
3043:           }
3044:         }

3046:         if (!dummyEnd[0]) {
3047:           /* (Middle) Front partial dummy row: columns 3/3: Right (Middle) Front, on Right (Middle) (Middle) */
3048:           const PetscInt neighbor = 14;
3049:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3050:           PetscInt i;
3051:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3052:             PetscInt  dLocal;
3053:             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 */
3054:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3055:             }
3056:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3057:               idxGlobalAll[count] = -1;
3058:             }
3059:           }
3060:         } else {
3061:           /* Right (Middle) Front partial dummy element, on (Middle) (Middle) (Middle) rank */
3062:           const PetscInt neighbor = 13;
3063:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3064:           PetscInt       i,dLocal;
3065:           i = stag->n[0];
3066:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3067:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3068:           }
3069:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3070:             idxGlobalAll[count] = -1;
3071:           }
3072:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3073:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3074:           }
3075:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3076:             idxGlobalAll[count] = -1;
3077:           }
3078:           ++i;
3079:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3080:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3081:               idxGlobalAll[count] = -1;
3082:             }
3083:           }
3084:         }
3085:       }
3086:     }

3088:     /* Front partial ghost layer: rows 3/3: Up Front, on Up (Middle) */
3089:     if (!dummyEnd[1]) {
3090:       PetscInt j;
3091:       for (j=0; j<ghostOffsetEnd[1]; ++j) {

3093:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left Up (Middle) */
3094:         if (!star && !dummyStart[0]) {
3095:           const PetscInt neighbor = 15;
3096:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3097:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3098:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3099:             PetscInt  dLocal;
3100:             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 */
3101:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3102:             }
3103:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3104:               idxGlobalAll[count] = -1;
3105:             }
3106:           }
3107:         } else {
3108:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3109:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3110:               idxGlobalAll[count] = -1;
3111:             }
3112:           }
3113:         }

3115:         /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) Up (Middle) */
3116:         {
3117:           const PetscInt neighbor = 16;
3118:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3119:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3120:             const PetscInt i = ighost - ghostOffsetStart[0];
3121:             PetscInt  dLocal;
3122:             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 */
3123:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3124:             }
3125:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3126:               idxGlobalAll[count] = -1;
3127:             }
3128:           }
3129:         }

3131:         if (!star && !dummyEnd[0]) {
3132:           /* Up Front partial dummy row: columsn 3/3: Right Up Front, on Right Up (Middle) */
3133:           const PetscInt neighbor = 17;
3134:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3135:           PetscInt       i;
3136:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3137:             PetscInt  dLocal;
3138:             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 */
3139:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3140:             }
3141:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3142:               idxGlobalAll[count] = -1;
3143:             }
3144:           }
3145:         } else if (dummyEnd[0]) {
3146:           /* Right Up Front partial dummy element, on (Middle) Up (Middle) */
3147:           const PetscInt neighbor = 16;
3148:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3149:           PetscInt       i,dLocal;
3150:           i = stag->n[0];
3151:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3152:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3153:           }
3154:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3155:             idxGlobalAll[count] = -1;
3156:           }
3157:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3158:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3159:           }
3160:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3161:             idxGlobalAll[count] = -1;
3162:           }
3163:           ++i;
3164:           for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3165:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3166:               idxGlobalAll[count] = -1;
3167:             }
3168:           }
3169:         } else {
3170:           /* Right Up Front dummies */
3171:           PetscInt i;
3172:           /* Right Down Front dummies */
3173:           for(i=0; i<ghostOffsetEnd[0]; ++i) {
3174:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3175:               idxGlobalAll[count] = -1;
3176:             }
3177:           }
3178:         }
3179:       }
3180:     } else {
3181:       /* Up Front partial dummy row */
3182:       PetscInt j = stag->n[1];

3184:       /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) (Middle) */
3185:       if (!dummyStart[0]) {
3186:         const PetscInt neighbor = 12;
3187:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3188:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3189:           const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3190:           PetscInt       dLocal;
3191:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3192:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3193:           }
3194:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3195:             idxGlobalAll[count] = -1;
3196:           }
3197:         }
3198:       } else {
3199:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3200:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3201:             idxGlobalAll[count] = -1;
3202:           }
3203:         }
3204:       }

3206:       /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) (Middle) */
3207:       {
3208:         const PetscInt neighbor = 13;
3209:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3210:         for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3211:           const PetscInt i = ighost - ghostOffsetStart[0];
3212:           PetscInt       dLocal;
3213:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3214:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3215:           }
3216:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3217:             idxGlobalAll[count] = -1;
3218:           }
3219:         }
3220:       }

3222:       if (!dummyEnd[0]) {
3223:         /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) (Middle) */
3224:         const PetscInt neighbor = 14;
3225:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3226:         PetscInt       i;
3227:         for (i=0; i<ghostOffsetEnd[0]; ++i) {
3228:           PetscInt       dLocal;
3229:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3230:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3231:           }
3232:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3233:             idxGlobalAll[count] = -1;
3234:           }
3235:         }
3236:       } else {
3237:         /* Right Up Front partial dummy element, on this rank */
3238:         const PetscInt neighbor = 13;
3239:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3240:         PetscInt       i,dLocal;
3241:         i = stag->n[0];
3242:         for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3243:           idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3244:         }
3245:         for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3246:           idxGlobalAll[count] = -1;
3247:         }
3248:         ++i;
3249:         for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3250:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3251:             idxGlobalAll[count] = -1;
3252:           }
3253:         }
3254:       }
3255:       ++j;
3256:       /* Up Front additional dummy rows */
3257:       for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
3258:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3259:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3260:             idxGlobalAll[count] = -1;
3261:           }
3262:         }
3263:       }
3264:     }
3265:     /* Front additional dummy layers */
3266:     ++k;
3267:     for (; k<stag->n[2] + ghostOffsetEnd[2]; ++k) {
3268:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
3269:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3270:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3271:             idxGlobalAll[count] = -1;
3272:           }
3273:         }
3274:       }
3275:     }
3276:   }

3278:   /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
3279:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
3280:   PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
3281:   return(0);
3282: }

3284: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM dm)
3285: {
3286:   PetscErrorCode  ierr;
3287:   DM_Stag * const stag = (DM_Stag*)dm->data;
3288:   const PetscInt epe = stag->entriesPerElement;
3289:   const PetscInt epr = stag->nGhost[0]*epe;
3290:   const PetscInt epl = stag->nGhost[1]*epr;

3293:   PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
3294:   stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]   = 0;
3295:   stag->locationOffsets[DMSTAG_BACK_DOWN]        = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + stag->dof[0];
3296:   stag->locationOffsets[DMSTAG_BACK_DOWN_RIGHT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epe;
3297:   stag->locationOffsets[DMSTAG_BACK_LEFT]        = stag->locationOffsets[DMSTAG_BACK_DOWN]       + stag->dof[1];
3298:   stag->locationOffsets[DMSTAG_BACK]             = stag->locationOffsets[DMSTAG_BACK_LEFT]       + stag->dof[1];
3299:   stag->locationOffsets[DMSTAG_BACK_RIGHT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epe;
3300:   stag->locationOffsets[DMSTAG_BACK_UP_LEFT]     = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epr;
3301:   stag->locationOffsets[DMSTAG_BACK_UP]          = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epr;
3302:   stag->locationOffsets[DMSTAG_BACK_UP_RIGHT]    = stag->locationOffsets[DMSTAG_BACK_UP_LEFT]    + epe;
3303:   stag->locationOffsets[DMSTAG_DOWN_LEFT]        = stag->locationOffsets[DMSTAG_BACK]            + stag->dof[2];
3304:   stag->locationOffsets[DMSTAG_DOWN]             = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + stag->dof[1];
3305:   stag->locationOffsets[DMSTAG_DOWN_RIGHT]       = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + epe;
3306:   stag->locationOffsets[DMSTAG_LEFT]             = stag->locationOffsets[DMSTAG_DOWN]            + stag->dof[2];
3307:   stag->locationOffsets[DMSTAG_ELEMENT]          = stag->locationOffsets[DMSTAG_LEFT]            + stag->dof[2];
3308:   stag->locationOffsets[DMSTAG_RIGHT]            = stag->locationOffsets[DMSTAG_LEFT]            + epe;
3309:   stag->locationOffsets[DMSTAG_UP_LEFT]          = stag->locationOffsets[DMSTAG_DOWN]            + epr;
3310:   stag->locationOffsets[DMSTAG_UP]               = stag->locationOffsets[DMSTAG_DOWN]            + epr;
3311:   stag->locationOffsets[DMSTAG_UP_RIGHT]         = stag->locationOffsets[DMSTAG_UP_LEFT]         + epe;
3312:   stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epl;
3313:   stag->locationOffsets[DMSTAG_FRONT_DOWN]       = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epl;
3314:   stag->locationOffsets[DMSTAG_FRONT_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epe;
3315:   stag->locationOffsets[DMSTAG_FRONT_LEFT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epl;
3316:   stag->locationOffsets[DMSTAG_FRONT]            = stag->locationOffsets[DMSTAG_BACK]            + epl;
3317:   stag->locationOffsets[DMSTAG_FRONT_RIGHT]      = stag->locationOffsets[DMSTAG_FRONT_LEFT]      + epe;
3318:   stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]    = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epr;
3319:   stag->locationOffsets[DMSTAG_FRONT_UP]         = stag->locationOffsets[DMSTAG_FRONT_DOWN]      + epr;
3320:   stag->locationOffsets[DMSTAG_FRONT_UP_RIGHT]   = stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]   + epe;
3321:   return(0);
3322: }

3324: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_3d(DM dm)
3325: {
3326:   PetscErrorCode  ierr;
3327:   DM_Stag * const stag = (DM_Stag*)dm->data;
3328:   PetscInt        *idxLocal,*idxGlobal,*globalOffsetsRecomputed;
3329:   const PetscInt  *globalOffsets;
3330:   PetscInt        count,d,entriesPerEdge,entriesPerFace,eprGhost,eplGhost,ghostOffsetStart[3],ghostOffsetEnd[3];
3331:   IS              isLocal,isGlobal;
3332:   PetscBool       dummyEnd[3];

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