Actual source code: stag3d.c

  1: /* Functions specific to the 3-dimensional implementation of DMStag */
  2: #include <petsc/private/dmstagimpl.h>

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

  7:   Collective

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

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

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

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

 43:   Level: beginner

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

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

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

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

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

160: /* Helper functions used in DMSetUp_Stag() */
161: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM);
162: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM);
163: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM,PetscInt**);
164: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM,const PetscInt*);
165: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM,const PetscInt*);
166: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM);

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

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

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

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

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

206:   /* Retrieve local size in stag->n */
207:   for (i=0; i<dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
208:   if (PetscDefined(USE_DEBUG)) {
209:     for (i=0; i<dim; ++i) {
210:       PetscInt Ncheck,j;
211:       Ncheck = 0;
212:       for (j=0; j<stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
213:       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]);
214:     }
215:   }

217:   /* Compute starting elements */
218:   for (i=0; i<dim; ++i) {
219:     stag->start[i] = 0;
220:     for (j=0;j<stag->rank[i];++j) {
221:       stag->start[i] += stag->l[i][j];
222:     }
223:   }

225:   /* Determine ranks of neighbors */
226:   DMStagSetUpBuildNeighbors_3d(dm);

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

249:   /* Check that we will not overflow 32 bit indices (slightly overconservative) */
250:   if (!PetscDefined(USE_64BIT_INDICES)) {
251:     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);
252:   }

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

256:     This again requires O(P) storage, which could be replaced with some global
257:     communication.
258:   */
259:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsets);

261:   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(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported boundary type");

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

327:   /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping */
328:   DMStagSetUpBuildScatter_3d(dm,globalOffsets);
329:   DMStagSetUpBuildL2G_3d(dm,globalOffsets);

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

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

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

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

347:   return(0);
348: }

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

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

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

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

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

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

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

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

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

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

477:   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(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Neighbor determination not implemented for %s",DMBoundaryTypes[stag->boundaryType[d]]);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

608:   return(0);
609: }

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

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

709:   return(0);
710: }

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

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

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

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

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

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

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

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

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

911:   /* Compute the number of local entries which correspond to any global entry */
912:   {
913:     PetscInt nNonDummyGhost[3];
914:     for (d=0; d<3; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
915:     if (star) {
916:       entriesToTransferTotal = (
917:           nNonDummyGhost[0] * stag->n[1]        * stag->n[2]
918:         + stag->n[0]        * nNonDummyGhost[1] * stag->n[2]
919:         + stag->n[0]        * stag->n[1]        * nNonDummyGhost[2]
920:         - 2 * (stag->n[0] * stag->n[1] * stag->n[2])) * 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 (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);

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

1549:   return(0);
1550: }

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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