Actual source code: stag3d.c

petsc-3.11.3 2019-06-26
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()
 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_PERIODIC:
313:       case DM_BOUNDARY_GHOSTED:
314:         switch (stag->stencilType) {
315:           case DMSTAG_STENCIL_NONE : /* only the extra one on the right/top edges */
316:             stag->nGhost[d] = stag->n[d];
317:             stag->startGhost[d] = stag->start[d];
318:             break;
319:           case DMSTAG_STENCIL_STAR :
320:           case DMSTAG_STENCIL_BOX :
321:             stag->nGhost[d] = stag->n[d] + 2*stag->stencilWidth;
322:             stag->startGhost[d] = stag->start[d] - stag->stencilWidth;
323:             break;
324:           default :
325:             SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
326:         }
327:         break;
328:       default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported boundary type in dimension %D",d);
329:     }
330:   }
331:   stag->entriesGhost = stag->nGhost[0]*stag->nGhost[1]*stag->nGhost[2]*stag->entriesPerElement;

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

337:   /* Free global offsets */
338:   PetscFree(globalOffsets);

340:   /* Precompute location offsets */
341:   DMStagComputeLocationOffsets_3d(dm);

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

346:   return(0);
347: }

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

361:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
362:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

364:   m = stag->nRanks[0];
365:   n = stag->nRanks[1];
366:   p = stag->nRanks[2];

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

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

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

447:   stag->nRanks[0] = m;
448:   stag->nRanks[1] = n;
449:   stag->nRanks[2] = p;
450:   return(0);
451: }

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

455:         n24 n25 n26
456:         n21 n22 n23
457:         n18 n19 n20 (Front, bigger z)

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

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

478:   /* Assemble some convenience variables */
479:   for (d=0; d<dim; ++d) {
480:     per[d]   = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
481:     first[d] = stag->firstRank[d];
482:     last[d]  = stag->lastRank[d];
483:     r[d]     = stag->rank[d];
484:     n[d]     = stag->nRanks[d];
485:   }

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

489:   neighborRank[0][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down back  */
490:   neighborRank[0][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
491:   neighborRank[0][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

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

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

501:   neighborRank[3][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left       back  */
502:   neighborRank[3][1]  =                                      r[1]    ;
503:   neighborRank[3][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

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

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

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

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

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

525:   neighborRank[9][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down       */
526:   neighborRank[9][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
527:   neighborRank[9][2]  =                                      r[2]    ;

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

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

537:   neighborRank[12][0] = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left             */
538:   neighborRank[12][1] =                                      r[1]    ;
539:   neighborRank[12][2] =                                      r[2]    ;

541:   neighborRank[13][0] =                                      r[0]    ; /*                  */
542:   neighborRank[13][1] =                                      r[1]    ;
543:   neighborRank[13][2] =                                      r[2]    ;

545:   neighborRank[14][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right            */
546:   neighborRank[14][1] =                                      r[1]    ;
547:   neighborRank[14][2] =                                      r[2]    ;

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

553:   neighborRank[16][0] =                                      r[0]    ; /*       up         */
554:   neighborRank[16][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
555:   neighborRank[16][2] =                                      r[2]    ;

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

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

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

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

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

577:   neighborRank[22][0] =                                      r[0]    ; /*            front */
578:   neighborRank[22][1] =                                      r[1]    ;
579:   neighborRank[22][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

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

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

589:   neighborRank[25][0] =                                      r[0]    ; /*       up   front */
590:   neighborRank[25][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
591:   neighborRank[25][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

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

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

607:   return(0);
608: }

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

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

708:   return(0);
709: }

711: /* A helper function to reduce code duplicaion as we loop over various ranges.
712:    i,j,k refer to element numbers on the rank where an element lives in the global
713:    representation (without ghosts) to be offset by a global offset per rank.
714:    ig,jg,kg refer to element numbers in the local (this rank) ghosted numbering.
715:    Note that this function could easily be converted to a macro (it does not compute
716:    anything except loop indices and the values of idxGlobal and idxLocal).  */
717: 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)
718: {
719:   PetscInt ig,jg,kg,d,c;

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

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

860:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
861:   if (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth || stag->n[2] < stag->stencilWidth) {
862:     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);
863:   }

865:   /* Check stencil type */
866:   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]);
867:   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");
868:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);

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

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

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

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

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

940:   /* Allocate arrays to populate */
941:   PetscMalloc1(entriesToTransferTotal,&idxLocal);
942:   PetscMalloc1(entriesToTransferTotal,&idxGlobal);

944:   /* Counts into idxLocal/idxGlobal */
945:   count = 0;

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

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

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

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

986:   /* DOWN BACK */
987:   if (!star && !dummyStart[1] && !dummyStart[2]) {
988:     const PetscInt  neighbor     = 1;
989:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
990:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
991:     const PetscInt  epFaceRow    = -1;
992:     const PetscInt  start0       = 0;
993:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
994:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
995:     const PetscInt  startGhost0  = ghostOffsetStart[0];
996:     const PetscInt  startGhost1  = 0;
997:     const PetscInt  startGhost2  = 0;
998:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
999:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1000:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1001:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT DOWN BACK */
1008:   if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyStart[2]) {
1009:     const PetscInt  neighbor     = 2;
1010:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* 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  = stag->nGhost[0] - ghostOffsetEnd[0];
1017:     const PetscInt  startGhost1  = 0;
1018:     const PetscInt  startGhost2  = 0;
1019:     const PetscInt  endGhost0    = stag->nGhost[0];
1020:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1021:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1022:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT BACK */
1029:   if (!star && !dummyStart[0] && !dummyStart[2]) {
1030:     const PetscInt  neighbor     = 3;
1031:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1032:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* May be a top boundary */
1033:     const PetscInt  epFaceRow    = -1;
1034:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1035:     const PetscInt  start1       = 0;
1036:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1037:     const PetscInt  startGhost0  = 0;
1038:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1039:     const PetscInt  startGhost2  = 0;
1040:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1041:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1042:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1043:     const PetscBool extra0       = PETSC_FALSE;
1044:     const PetscBool extra1       = dummyEnd[1];
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:   /* BACK */
1050:   if (!dummyStart[2]) {
1051:     const PetscInt  neighbor     = 4;
1052:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may  be a right boundary */
1053:     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 */
1054:     const PetscInt  epFaceRow    = -1;
1055:     const PetscInt  start0       = 0;
1056:     const PetscInt  start1       = 0;
1057:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1058:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1059:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1060:     const PetscInt  startGhost2  = 0;
1061:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1062:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1063:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1064:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT BACK */
1071:   if (!star && !dummyEnd[0] && !dummyStart[2]) {
1072:     const PetscInt  neighbor     = 5;
1073:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1074:     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 */
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  = stag->nGhost[0] - ghostOffsetEnd[0];
1080:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1081:     const PetscInt  startGhost2  = 0;
1082:     const PetscInt  endGhost0    = stag->nGhost[0];
1083:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1084:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1085:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT UP BACK */
1092:   if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyStart[2]) {
1093:     const PetscInt  neighbor     = 6;
1094:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1095:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1096:     const PetscInt  epFaceRow    = -1;
1097:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1098:     const PetscInt  start1       = 0;
1099:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1100:     const PetscInt  startGhost0  = 0;
1101:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1102:     const PetscInt  startGhost2  = 0;
1103:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1104:     const PetscInt  endGhost1    = stag->nGhost[1];
1105:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1106:     const PetscBool extra0       = PETSC_FALSE;
1107:     const PetscBool extra1       = PETSC_FALSE;
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:   /* UP BACK */
1113:   if (!star && !dummyEnd[1] && !dummyStart[2]) {
1114:     const PetscInt  neighbor     = 7;
1115:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1116:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1117:     const PetscInt  epFaceRow    = -1;
1118:     const PetscInt  start0       = 0;
1119:     const PetscInt  start1       = 0;
1120:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1121:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1122:     const PetscInt  startGhost1  = stag->nGhost[1]-ghostOffsetEnd[1];
1123:     const PetscInt  startGhost2  = 0;
1124:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1125:     const PetscInt  endGhost1    = stag->nGhost[1];
1126:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1127:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT UP BACK */
1134:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyStart[2]) {
1135:     const PetscInt  neighbor     = 8;
1136:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1137:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[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  = stag->nGhost[0] - ghostOffsetEnd[0];
1143:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1144:     const PetscInt  startGhost2  = 0;
1145:     const PetscInt  endGhost0    = stag->nGhost[0];
1146:     const PetscInt  endGhost1    = stag->nGhost[1];
1147:     const PetscInt  endGhost2    = ghostOffsetStart[2];
1148:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT DOWN */
1155:   if (!star && !dummyStart[0] && !dummyStart[1]) {
1156:     const PetscInt  neighbor     = 9;
1157:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1158:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1159:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1160:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1161:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1162:     const PetscInt  start2       = 0;
1163:     const PetscInt  startGhost0  = 0;
1164:     const PetscInt  startGhost1  = 0;
1165:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1166:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1167:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1168:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1169:     const PetscBool extra0       = PETSC_FALSE;
1170:     const PetscBool extra1       = PETSC_FALSE;
1171:     const PetscBool extra2       = dummyEnd[2];
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:   /* DOWN */
1176:   if (!dummyStart[1]) {
1177:     const PetscInt  neighbor     = 10;
1178:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1179:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1180:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1181:     const PetscInt  start0       = 0;
1182:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1183:     const PetscInt  start2       = 0;
1184:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1185:     const PetscInt  startGhost1  = 0;
1186:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1187:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1188:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1189:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1190:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT DOWN */
1197:   if (!star && !dummyEnd[0] && !dummyStart[1]) {
1198:     const PetscInt  neighbor     = 11;
1199:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1200:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1201:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* 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  = stag->nGhost[0] - ghostOffsetEnd[0];
1206:     const PetscInt  startGhost1  = 0;
1207:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1208:     const PetscInt  endGhost0    = stag->nGhost[0];
1209:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1210:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1211:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT */
1218:   if (!dummyStart[0]) {
1219:     const PetscInt  neighbor     = 12;
1220:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1221:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1222:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0];
1223:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1224:     const PetscInt  start1       = 0;
1225:     const PetscInt  start2       = 0;
1226:     const PetscInt  startGhost0  = 0;
1227:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1228:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1229:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1230:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1231:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1232:     const PetscBool extra0       = PETSC_FALSE;
1233:     const PetscBool extra1       = dummyEnd[1];
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:   /* (HERE) */
1239:   {
1240:     const PetscInt  neighbor     = 13;
1241:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1242:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1243:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
1244:     const PetscInt  start0       = 0;
1245:     const PetscInt  start1       = 0;
1246:     const PetscInt  start2       = 0;
1247:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1248:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1249:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1250:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1251:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1252:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1253:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT */
1260:   if (!dummyEnd[0]) {
1261:     const PetscInt  neighbor     = 14;
1262:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1263:     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 */
1264:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1265:     const PetscInt  start0       = 0;
1266:     const PetscInt  start1       = 0;
1267:     const PetscInt  start2       = 0;
1268:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1269:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1270:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1271:     const PetscInt  endGhost0    = stag->nGhost[0];
1272:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1273:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1274:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT UP */
1281:   if (!star && !dummyStart[0] && !dummyEnd[1]) {
1282:     const PetscInt  neighbor     = 15;
1283:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1284:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1285:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0];
1286:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1287:     const PetscInt  start1       = 0;
1288:     const PetscInt  start2       = 0;
1289:     const PetscInt  startGhost0  = 0;
1290:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1291:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1292:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1293:     const PetscInt  endGhost1    = stag->nGhost[1];
1294:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1295:     const PetscBool extra0       = PETSC_FALSE;
1296:     const PetscBool extra1       = PETSC_FALSE;
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:   /* UP */
1302:   if (!dummyEnd[1]) {
1303:     const PetscInt  neighbor     = 16;
1304:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1305:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1306:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1307:     const PetscInt  start0       = 0;
1308:     const PetscInt  start1       = 0;
1309:     const PetscInt  start2       = 0;
1310:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1311:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1312:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1313:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1314:     const PetscInt  endGhost1    = stag->nGhost[1];
1315:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1316:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT UP */
1323:   if (!star && !dummyEnd[0] && !dummyEnd[1]) {
1324:     const PetscInt  neighbor     = 17;
1325:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1326:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1327:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* 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  = stag->nGhost[0] - ghostOffsetEnd[0];
1332:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1333:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1334:     const PetscInt  endGhost0    = stag->nGhost[0];
1335:     const PetscInt  endGhost1    = stag->nGhost[1];
1336:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1337:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT DOWN FRONT */
1344:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1345:     const PetscInt  neighbor     = 18;
1346:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1347:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1348:     const PetscInt  epFaceRow    = -1;
1349:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1350:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1351:     const PetscInt  start2       = 0;
1352:     const PetscInt  startGhost0  = 0;
1353:     const PetscInt  startGhost1  = 0;
1354:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1355:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1356:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1357:     const PetscInt  endGhost2    = stag->nGhost[2];
1358:     const PetscBool extra0       = PETSC_FALSE;
1359:     const PetscBool extra1       = PETSC_FALSE;
1360:     const PetscBool extra2       = PETSC_FALSE;
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:   /* DOWN FRONT */
1365:   if (!star && !dummyStart[1] && !dummyEnd[2]) {
1366:     const PetscInt  neighbor     = 19;
1367:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1368:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1369:     const PetscInt  epFaceRow    = -1;
1370:     const PetscInt  start0       = 0;
1371:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1372:     const PetscInt  start2       = 0;
1373:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1374:     const PetscInt  startGhost1  = 0;
1375:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1376:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1377:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1378:     const PetscInt  endGhost2    = stag->nGhost[2];
1379:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT DOWN FRONT */
1386:   if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyEnd[2]) {
1387:     const PetscInt  neighbor     = 20;
1388:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* 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  = stag->nGhost[0] - ghostOffsetEnd[0];
1395:     const PetscInt  startGhost1  = 0;
1396:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1397:     const PetscInt  endGhost0    = stag->nGhost[0];
1398:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1399:     const PetscInt  endGhost2    = stag->nGhost[2];
1400:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT FRONT */
1407:   if (!star && !dummyStart[0] && !dummyEnd[2]) {
1408:     const PetscInt  neighbor     = 21;
1409:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1410:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0);  /* We+neighbor may be a top boundary */
1411:     const PetscInt  epFaceRow    = -1;
1412:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1413:     const PetscInt  start1       = 0;
1414:     const PetscInt  start2       = 0;
1415:     const PetscInt  startGhost0  = 0;
1416:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1417:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1418:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1419:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1420:     const PetscInt  endGhost2    = stag->nGhost[2];
1421:     const PetscBool extra0       = PETSC_FALSE;
1422:     const PetscBool extra1       = dummyEnd[1];
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:   /* FRONT */
1428:   if (!dummyEnd[2]) {
1429:     const PetscInt  neighbor     = 22;
1430:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1431:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1432:     const PetscInt  epFaceRow    = -1;
1433:     const PetscInt  start0       = 0;
1434:     const PetscInt  start1       = 0;
1435:     const PetscInt  start2       = 0;
1436:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1437:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1438:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1439:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1440:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1441:     const PetscInt  endGhost2    = stag->nGhost[2];
1442:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT FRONT */
1449:   if (!star && !dummyEnd[0] && !dummyEnd[2]) {
1450:     const PetscInt  neighbor     = 23;
1451:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1452:     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 */
1453:     const PetscInt  epFaceRow    = -1;
1454:     const PetscInt  start0       = 0;
1455:     const PetscInt  start1       = 0;
1456:     const PetscInt  start2       = 0;
1457:     const PetscInt  startGhost0  = stag->nGhost[0] - ghostOffsetEnd[0];
1458:     const PetscInt  startGhost1  = ghostOffsetStart[1];
1459:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1460:     const PetscInt  endGhost0    = stag->nGhost[0];
1461:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
1462:     const PetscInt  endGhost2    = stag->nGhost[2];
1463:     const PetscBool extra0       = PETSC_FALSE;
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:   /* LEFT UP FRONT */
1470:   if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyEnd[2]) {
1471:     const PetscInt  neighbor     = 24;
1472:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1473:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1474:     const PetscInt  epFaceRow    = -1;
1475:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1476:     const PetscInt  start1       = 0;
1477:     const PetscInt  start2       = 0;
1478:     const PetscInt  startGhost0  = 0;
1479:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1480:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1481:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1482:     const PetscInt  endGhost1    = stag->nGhost[1];
1483:     const PetscInt  endGhost2    = stag->nGhost[2];
1484:     const PetscBool extra0       = PETSC_FALSE;
1485:     const PetscBool extra1       = PETSC_FALSE;
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:   /* UP FRONT */
1491:   if (!star && !dummyEnd[1] && !dummyEnd[2]) {
1492:     const PetscInt  neighbor     = 25;
1493:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1494:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1495:     const PetscInt  epFaceRow    = -1;
1496:     const PetscInt  start0       = 0;
1497:     const PetscInt  start1       = 0;
1498:     const PetscInt  start2       = 0;
1499:     const PetscInt  startGhost0  = ghostOffsetStart[0];
1500:     const PetscInt  startGhost1  = stag->nGhost[1] - ghostOffsetEnd[1];
1501:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1502:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
1503:     const PetscInt  endGhost1    = stag->nGhost[1];
1504:     const PetscInt  endGhost2    = stag->nGhost[2];
1505:     const PetscBool extra0       = dummyEnd[0];
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:   /* RIGHT UP FRONT */
1512:   if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyEnd[2]) {
1513:     const PetscInt  neighbor     = 26;
1514:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1515:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[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  = stag->nGhost[0] - ghostOffsetEnd[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];
1524:     const PetscInt  endGhost1    = stag->nGhost[1];
1525:     const PetscInt  endGhost2    = stag->nGhost[2];
1526:     const PetscBool extra0       = PETSC_FALSE;
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: #if defined(PETSC_USE_DEBUG)
1533:   if (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);
1534: #endif

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

1551:   return(0);
1552: }

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2805:         if (!star && !dummyEnd[0]) {
2806:           /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) Front rank */
2807:           const PetscInt neighbor = 23;
2808:           PetscInt       i;
2809:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
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:           }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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