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 elements, 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/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 - arrays of 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: {
 49:   DMCreate(comm,dm);
 50:   DMSetDimension(*dm,3);
 51:   DMStagInitialize(bndx,bndy,bndz,M,N,P,m,n,p,dof0,dof1,dof2,dof3,stencilType,stencilWidth,lx,ly,lz,*dm);
 52:   return 0;
 53: }

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

 65:   DMGetCoordinateDM(dm,&dmCoord);
 66:   stagCoord = (DM_Stag*) dmCoord->data;
 67:   for (s=0; s<4; ++s) {
 69:   }
 70:   DMCreateLocalVector(dmCoord,&coordLocal);
 71:   DMStagVecGetArray(dmCoord,coordLocal,&arr);
 72:   if (stagCoord->dof[0]) {
 73:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN_LEFT,0,&ibackdownleft);
 74:   }
 75:   if (stagCoord->dof[1]) {
 76:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN     ,0,&ibackdown);
 77:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_LEFT     ,0,&ibackleft);
 78:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN_LEFT     ,0,&idownleft);
 79:   }
 80:   if (stagCoord->dof[2]) {
 81:     DMStagGetLocationSlot(dmCoord,DMSTAG_BACK          ,0,&iback);
 82:     DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN          ,0,&idown);
 83:     DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT          ,0,&ileft);
 84:   }
 85:   if (stagCoord->dof[3]) {
 86:     DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT       ,0,&ielement);
 87:   }
 88:   DMStagGetGhostCorners(dmCoord,&start_ghost[0],&start_ghost[1],&start_ghost[2],&n_ghost[0],&n_ghost[1],&n_ghost[2]);
 89:   min[0] = xmin; min[1]= ymin; min[2] = zmin;
 90:   h[0] = (xmax-xmin)/stagCoord->N[0];
 91:   h[1] = (ymax-ymin)/stagCoord->N[1];
 92:   h[2] = (zmax-zmin)/stagCoord->N[2];

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

155: /* Helper functions used in DMSetUp_Stag() */
156: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM);
157: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM);
158: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM,PetscInt**);
159: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM,const PetscInt*);
160: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM,const PetscInt*);
161: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM);

163: PETSC_INTERN PetscErrorCode DMSetUp_Stag_3d(DM dm)
164: {
165:   DM_Stag * const stag = (DM_Stag*)dm->data;
166:   PetscMPIInt     rank;
167:   PetscInt        i,j,d;
168:   PetscInt        *globalOffsets;
169:   const PetscInt  dim = 3;

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

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

176:   /* Determine location of rank in grid */
177:     stag->rank[0] = rank % stag->nRanks[0];
178:     stag->rank[1] = rank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0];
179:     stag->rank[2] = rank / (stag->nRanks[0] * stag->nRanks[1]);
180:     for (d=0; d<dim; ++d) {
181:       stag->firstRank[d] = PetscNot(stag->rank[d]);
182:       stag->lastRank[d]  = (PetscBool)(stag->rank[d] == stag->nRanks[d]-1);
183:     }

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

199:   /* Retrieve local size in stag->n */
200:   for (i=0; i<dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
201:   if (PetscDefined(USE_DEBUG)) {
202:     for (i=0; i<dim; ++i) {
203:       PetscInt Ncheck,j;
204:       Ncheck = 0;
205:       for (j=0; j<stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
207:     }
208:   }

210:   /* Compute starting elements */
211:   for (i=0; i<dim; ++i) {
212:     stag->start[i] = 0;
213:     for (j=0;j<stag->rank[i];++j) {
214:       stag->start[i] += stag->l[i][j];
215:     }
216:   }

218:   /* Determine ranks of neighbors */
219:   DMStagSetUpBuildNeighbors_3d(dm);

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

242:   /* Check that we will not overflow 32 bit indices (slightly overconservative) */
243:   if (!PetscDefined(USE_64BIT_INDICES)) {
245:   }

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

249:     This again requires O(P) storage, which could be replaced with some global
250:     communication.
251:   */
252:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsets);


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

320:   /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping */
321:   DMStagSetUpBuildScatter_3d(dm,globalOffsets);
322:   DMStagSetUpBuildL2G_3d(dm,globalOffsets);

324:   /* In special cases, create a dedicated injective local-to-global map */
325:   if ((stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) ||
326:       (stag->boundaryType[1] == DM_BOUNDARY_PERIODIC && stag->nRanks[1] == 1) ||
327:       (stag->boundaryType[2] == DM_BOUNDARY_PERIODIC && stag->nRanks[2] == 1)) {
328:     DMStagPopulateLocalToGlobalInjective(dm);
329:   }

331:   /* Free global offsets */
332:   PetscFree(globalOffsets);

334:   /* Precompute location offsets */
335:   DMStagComputeLocationOffsets_3d(dm);

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

340:   return 0;
341: }

343: /* adapted from da3.c */
344: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM dm)
345: {
346:   PetscMPIInt     rank,size;
347:   PetscInt        m,n,p,pm;
348:   DM_Stag * const stag = (DM_Stag*)dm->data;
349:   const PetscInt M = stag->N[0];
350:   const PetscInt N = stag->N[1];
351:   const PetscInt P = stag->N[2];

353:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
354:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

356:   m = stag->nRanks[0];
357:   n = stag->nRanks[1];
358:   p = stag->nRanks[2];

360:   if (m != PETSC_DECIDE) {
363:   }
364:   if (n != PETSC_DECIDE) {
367:   }
368:   if (p != PETSC_DECIDE) {
371:   }

374:   /* Partition the array among the processors */
375:   if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
376:     m = size/(n*p);
377:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
378:     n = size/(m*p);
379:   } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
380:     p = size/(m*n);
381:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
382:     /* try for squarish distribution */
383:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
384:     if (!m) m = 1;
385:     while (m > 0) {
386:       n = size/(m*p);
387:       if (m*n*p == size) break;
388:       m--;
389:     }
391:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
392:   } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
393:     /* try for squarish distribution */
394:     m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
395:     if (!m) m = 1;
396:     while (m > 0) {
397:       p = size/(m*n);
398:       if (m*n*p == size) break;
399:       m--;
400:     }
402:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
403:   } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
404:     /* try for squarish distribution */
405:     n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
406:     if (!n) n = 1;
407:     while (n > 0) {
408:       p = size/(m*n);
409:       if (m*n*p == size) break;
410:       n--;
411:     }
413:     if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
414:   } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
415:     /* try for squarish distribution */
416:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
417:     if (!n) n = 1;
418:     while (n > 0) {
419:       pm = size/n;
420:       if (n*pm == size) break;
421:       n--;
422:     }
423:     if (!n) n = 1;
424:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
425:     if (!m) m = 1;
426:     while (m > 0) {
427:       p = size/(m*n);
428:       if (m*n*p == size) break;
429:       m--;
430:     }
431:     if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}


439:   stag->nRanks[0] = m;
440:   stag->nRanks[1] = n;
441:   stag->nRanks[2] = p;
442:   return 0;
443: }

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

447:         n24 n25 n26
448:         n21 n22 n23
449:         n18 n19 n20 (Front, bigger z)

451:         n15 n16 n17
452:         n12     n14   ^ y
453:         n9  n10 n11   |
454:                       +--> x
455:         n6  n7  n8
456:         n3  n4  n5
457:         n0  n1  n2 (Back, smaller z) */
458: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM dm)
459: {
460:   DM_Stag * const stag = (DM_Stag*)dm->data;
461:   PetscInt        d,i;
462:   PetscBool       per[3],first[3],last[3];
463:   PetscInt        neighborRank[27][3],r[3],n[3];
464:   const PetscInt  dim = 3;


468:   /* Assemble some convenience variables */
469:   for (d=0; d<dim; ++d) {
470:     per[d]   = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
471:     first[d] = stag->firstRank[d];
472:     last[d]  = stag->lastRank[d];
473:     r[d]     = stag->rank[d];
474:     n[d]     = stag->nRanks[d];
475:   }

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

479:   neighborRank[0][0]  = first[0] ? (per[0] ?  n[0]-1 : -1) : r[0] - 1; /* left  down back  */
480:   neighborRank[0][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
481:   neighborRank[0][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

483:   neighborRank[1][0]  =                                      r[0]    ; /*       down back  */
484:   neighborRank[1][1]  = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
485:   neighborRank[1][2]  = first[2] ? (per[2] ?  n[2]-1 : -1) : r[2] - 1;

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

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

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

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

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

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

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

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

519:   neighborRank[10][0] =                                      r[0]    ; /*       down       */
520:   neighborRank[10][1] = first[1] ? (per[1] ?  n[1]-1 : -1) : r[1] - 1;
521:   neighborRank[10][2] =                                      r[2]    ;

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

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

531:   neighborRank[13][0] =                                      r[0]    ; /*                  */
532:   neighborRank[13][1] =                                      r[1]    ;
533:   neighborRank[13][2] =                                      r[2]    ;

535:   neighborRank[14][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right            */
536:   neighborRank[14][1] =                                      r[1]    ;
537:   neighborRank[14][2] =                                      r[2]    ;

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

543:   neighborRank[16][0] =                                      r[0]    ; /*       up         */
544:   neighborRank[16][1] = last[1]  ? (per[1] ?  0      : -1) : r[1] + 1;
545:   neighborRank[16][2] =                                      r[2]    ;

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

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

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

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

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

567:   neighborRank[22][0] =                                      r[0]    ; /*            front */
568:   neighborRank[22][1] =                                      r[1]    ;
569:   neighborRank[22][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

571:   neighborRank[23][0] = last[0]  ? (per[0] ?  0      : -1) : r[0] + 1; /* right      front */
572:   neighborRank[23][1] =                                      r[1]    ;
573:   neighborRank[23][2] = last[2]  ? (per[2] ?  0      : -1) : r[2] + 1;

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

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

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

587:   /* Then, compute the rank of each in the linear ordering */
588:   PetscMalloc1(27,&stag->neighbors);
589:   for (i=0; i<27; ++i) {
590:     if  (neighborRank[i][0] >= 0 && neighborRank[i][1] >=0 && neighborRank[i][2] >=0) {
591:       stag->neighbors[i] = neighborRank[i][0] + n[0]*neighborRank[i][1] + n[0]*n[1]*neighborRank[i][2];
592:     } else {
593:       stag->neighbors[i] = -1;
594:     }
595:   }

597:   return 0;
598: }

600: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM dm,PetscInt **pGlobalOffsets)
601: {
602:   const DM_Stag * const stag = (DM_Stag*)dm->data;
603:   PetscInt              *globalOffsets;
604:   PetscInt              i,j,k,d,entriesPerEdge,entriesPerFace,count;
605:   PetscMPIInt           size;
606:   PetscBool             extra[3];

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

696:   return 0;
697: }

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

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

836: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM dm,const PetscInt *globalOffsets)
837: {
838:   DM_Stag * const stag = (DM_Stag*)dm->data;
839:   PetscInt       d,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerCorner,entriesPerEdge,entriesPerFace,entriesToTransferTotal,count,eprGhost,eplGhost;
840:   PetscInt       *idxLocal,*idxGlobal;
841:   PetscMPIInt    rank;
842:   PetscInt       nNeighbors[27][3];
843:   PetscBool      star,nextToDummyEnd[3],dummyStart[3],dummyEnd[3];

845:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
846:   if (stag->stencilType != DMSTAG_STENCIL_NONE && (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth || stag->n[2] < stag->stencilWidth)) {
847:     SETERRQ(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);
848:   }

850:   /* Check stencil type */
852:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR || stag->stencilType == DMSTAG_STENCIL_NONE);

854:   /* Compute numbers of elements on each neighbor */
855:   {
856:     PetscInt i;
857:     for (i=0; i<27; ++i) {
858:       const PetscInt neighborRank = stag->neighbors[i];
859:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
860:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
861:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
862:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
863:       } /* else leave uninitialized - error if accessed */
864:     }
865:   }

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

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

886:   /* Convenience variables */
887:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
888:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
889:   entriesPerCorner                    = stag->dof[0];
890:   for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
891:   eprGhost                            = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row   */
892:   eplGhost                            = stag->nGhost[1]*eprGhost;                /* epl = entries per (element) layer */

894:   /* Compute the number of local entries which correspond to any global entry */
895:   {
896:     PetscInt nNonDummyGhost[3];
897:     for (d=0; d<3; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
898:     if (star) {
899:       entriesToTransferTotal = (
900:           nNonDummyGhost[0] * stag->n[1]        * stag->n[2]
901:         + stag->n[0]        * nNonDummyGhost[1] * stag->n[2]
902:         + stag->n[0]        * stag->n[1]        * nNonDummyGhost[2]
903:         - 2 * (stag->n[0] * stag->n[1] * stag->n[2])) * stag->entriesPerElement
904:         + (dummyEnd[0]                               ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace   : 0)
905:         + (dummyEnd[1]                               ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace   : 0)
906:         + (dummyEnd[2]                               ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace   : 0)
907:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
908:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
909:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
910:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
911:     } else {
912:       entriesToTransferTotal  =  nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
913:         + (dummyEnd[0]                               ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace   : 0)
914:         + (dummyEnd[1]                               ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace   : 0)
915:         + (dummyEnd[2]                               ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace   : 0)
916:         + (dummyEnd[0] && dummyEnd[1]                ? nNonDummyGhost[2]                     * entriesPerEdge   : 0)
917:         + (dummyEnd[2] && dummyEnd[0]                ? nNonDummyGhost[1]                     * entriesPerEdge   : 0)
918:         + (dummyEnd[1] && dummyEnd[2]                ? nNonDummyGhost[0]                     * entriesPerEdge   : 0)
919:         + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ?                                         entriesPerCorner : 0);
920:     }
921:   }

923:   /* Allocate arrays to populate */
924:   PetscMalloc1(entriesToTransferTotal,&idxLocal);
925:   PetscMalloc1(entriesToTransferTotal,&idxGlobal);

927:   /* Counts into idxLocal/idxGlobal */
928:   count = 0;

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

933:       1. The loop bounds (i/ighost,j/jghost,k/kghost)
934:       2. the strides used in the loop body, which depend on whether the current
935:       rank and/or its neighbor are a non-periodic right/top/front boundary, which has additional
936:       points in the global representation.
937:       - If the neighboring rank is a right/top/front boundary,
938:       then eprNeighbor (entries per element row on the neighbor) and/or eplNeighbor (entries per element layer on the neighbor)
939:       are different.
940:       - If this rank is a non-periodic right/top/front boundary (checking entries of stag->lastRank),
941:       there is an extra loop over 1-3 boundary faces)

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

948:   /* LEFT DOWN BACK */
949:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
950:     const PetscInt  neighbor     = 0;
951:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
952:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
953:     const PetscInt  epFaceRow    = -1;
954:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
955:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
956:     const PetscInt  start2       = nNeighbors[neighbor][2] - ghostOffsetStart[2];
957:     const PetscInt  startGhost0  = 0;
958:     const PetscInt  startGhost1  = 0;
959:     const PetscInt  startGhost2  = 0;
960:     const PetscInt  endGhost0    = ghostOffsetStart[0];
961:     const PetscInt  endGhost1    = ghostOffsetStart[1];
962:     const PetscInt  endGhost2    = ghostOffsetStart[2];
963:     const PetscBool extra0       = PETSC_FALSE;
964:     const PetscBool extra1       = PETSC_FALSE;
965:     const PetscBool extra2       = PETSC_FALSE;
966:     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);
967:   }

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

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

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

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

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

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

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

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

1137:   /* LEFT DOWN */
1138:   if (!star && !dummyStart[0] && !dummyStart[1]) {
1139:     const PetscInt  neighbor     = 9;
1140:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1141:     const PetscInt  eplNeighbor  = eprNeighbor  * nNeighbors[neighbor][1];
1142:     const PetscInt  epFaceRow    = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1143:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1144:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1145:     const PetscInt  start2       = 0;
1146:     const PetscInt  startGhost0  = 0;
1147:     const PetscInt  startGhost1  = 0;
1148:     const PetscInt  startGhost2  = ghostOffsetStart[2];
1149:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1150:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1151:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
1152:     const PetscBool extra0       = PETSC_FALSE;
1153:     const PetscBool extra1       = PETSC_FALSE;
1154:     const PetscBool extra2       = dummyEnd[2];
1155:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1156:   }

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

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

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

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

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

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

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

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

1326:   /* LEFT DOWN FRONT */
1327:   if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1328:     const PetscInt  neighbor     = 18;
1329:     const PetscInt  eprNeighbor  = stag->entriesPerElement * nNeighbors[neighbor][0];
1330:     const PetscInt  eplNeighbor  = eprNeighbor * nNeighbors[neighbor][1];
1331:     const PetscInt  epFaceRow    = -1;
1332:     const PetscInt  start0       = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1333:     const PetscInt  start1       = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1334:     const PetscInt  start2       = 0;
1335:     const PetscInt  startGhost0  = 0;
1336:     const PetscInt  startGhost1  = 0;
1337:     const PetscInt  startGhost2  = stag->nGhost[2] - ghostOffsetEnd[2];
1338:     const PetscInt  endGhost0    = ghostOffsetStart[0];
1339:     const PetscInt  endGhost1    = ghostOffsetStart[1];
1340:     const PetscInt  endGhost2    = stag->nGhost[2];
1341:     const PetscBool extra0       = PETSC_FALSE;
1342:     const PetscBool extra1       = PETSC_FALSE;
1343:     const PetscBool extra2       = PETSC_FALSE;
1344:     DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1345:   }

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

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

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

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

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

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

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

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


1517:   /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
1518:   {
1519:     Vec local,global;
1520:     IS isLocal,isGlobal;
1521:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
1522:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
1523:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
1524:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
1525:     VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
1526:     VecDestroy(&global);
1527:     VecDestroy(&local);
1528:     ISDestroy(&isLocal);  /* frees idxLocal */
1529:     ISDestroy(&isGlobal); /* free idxGlobal */
1530:   }

1532:   return 0;
1533: }

1535: /* Note: this function assumes that DMBoundary types of none, ghosted, and periodic are the only ones of interest.
1536: Adding support for others should be done very carefully.  */
1537: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM dm,const PetscInt *globalOffsets)
1538: {
1539:   const DM_Stag * const stag = (DM_Stag*)dm->data;
1540:   PetscInt              *idxGlobalAll;
1541:   PetscInt              d,count,ighost,jghost,kghost,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerFace,entriesPerEdge;
1542:   PetscInt              nNeighbors[27][3],eprNeighbor[27],eplNeighbor[27],globalOffsetNeighbor[27];
1543:   PetscBool             nextToDummyEnd[3],dummyStart[3],dummyEnd[3],star;


1546:   /* Check stencil type */
1548:   star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR || stag->stencilType == DMSTAG_STENCIL_NONE);

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

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

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

1563:   /* Compute numbers of elements on each neighbor  and offset*/
1564:   {
1565:     PetscInt i;
1566:     for (i=0; i<27; ++i) {
1567:       const PetscInt neighborRank = stag->neighbors[i];
1568:       if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
1569:         nNeighbors[i][0] =  stag->l[0][neighborRank % stag->nRanks[0]];
1570:         nNeighbors[i][1] =  stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
1571:         nNeighbors[i][2] =  stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
1572:         globalOffsetNeighbor[i] = globalOffsets[stag->neighbors[i]];
1573:       } /* else leave uninitialized - error if accessed */
1574:     }
1575:   }

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

1689:   /* Populate idxGlobalAll */
1690:   PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
1691:   count = 0;

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

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

1703:           /* Loop over columns 1/3: Left Back Down*/
1704:           if (!dummyStart[0]) {
1705:             const PetscInt neighbor = 0;
1706:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1707:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1708:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1709:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1710:               }
1711:             }
1712:           } else {
1713:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1714:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1715:                 idxGlobalAll[count] = -1;
1716:               }
1717:             }
1718:           }

1720:           /* Loop over columns 2/3: (Middle) Down Back */
1721:           {
1722:             const PetscInt neighbor = 1;
1723:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1724:               const PetscInt i = ighost - ghostOffsetStart[0];
1725:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1726:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1727:               }
1728:             }
1729:           }

1731:           /* Loop over columns 3/3: Right Down Back */
1732:           if (!dummyEnd[0]) {
1733:             const PetscInt neighbor = 2;
1734:             PetscInt       i;
1735:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1736:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1737:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1738:               }
1739:             }
1740:           } else {
1741:             /* Partial dummy entries on (Middle) Down Back neighbor */
1742:             const PetscInt neighbor          = 1;
1743:             PetscInt       i,dLocal;
1744:             i = stag->n[0];
1745:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1746:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1747:             }
1748:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1749:               idxGlobalAll[count] = -1;
1750:             }
1751:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1752:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1753:             }
1754:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1755:               idxGlobalAll[count] = -1;
1756:             }
1757:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1758:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1759:             }
1760:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1761:               idxGlobalAll[count] = -1;
1762:             }
1763:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1764:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1765:             }
1766:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1767:               idxGlobalAll[count] = -1;
1768:             }
1769:             ++i;
1770:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1771:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1772:                 idxGlobalAll[count] = -1;
1773:               }
1774:             }
1775:           }
1776:         }
1777:       } else {
1778:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1779:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
1780:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
1781:               idxGlobalAll[count] = -1;
1782:             }
1783:           }
1784:         }
1785:       }

1787:       /* Loop over rows 2/3: (Middle) Back */
1788:       {
1789:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
1790:           const PetscInt j = jghost - ghostOffsetStart[1];

1792:           /* Loop over columns 1/3: Left (Middle) Back */
1793:           if (!star && !dummyStart[0]) {
1794:             const PetscInt neighbor = 3;
1795:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1796:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1797:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1798:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1799:               }
1800:             }
1801:           } else {
1802:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1803:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1804:                 idxGlobalAll[count] = -1;
1805:               }
1806:             }
1807:           }

1809:           /* Loop over columns 2/3: (Middle) (Middle) Back */
1810:           {
1811:             const PetscInt neighbor = 4;
1812:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1813:               const PetscInt i = ighost - ghostOffsetStart[0];
1814:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1815:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1816:               }
1817:             }
1818:           }

1820:           /* Loop over columns 3/3: Right (Middle) Back */
1821:           if (!star && !dummyEnd[0]) {
1822:             const PetscInt neighbor = 5;
1823:             PetscInt       i;
1824:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1825:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1826:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1827:               }
1828:             }
1829:           } else if (dummyEnd[0]) {
1830:             /* Partial dummy entries on (Middle) (Middle) Back rank */
1831:             const PetscInt neighbor = 4;
1832:             PetscInt       i,dLocal;
1833:             i = stag->n[0];
1834:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1835:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1836:             }
1837:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1838:               idxGlobalAll[count] = -1;
1839:             }
1840:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1841:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1842:             }
1843:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1844:               idxGlobalAll[count] = -1;
1845:             }
1846:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1847:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1848:             }
1849:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1850:               idxGlobalAll[count] = -1;
1851:             }
1852:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1853:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1854:             }
1855:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1856:               idxGlobalAll[count] = -1;
1857:             }
1858:             ++i;
1859:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1860:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1861:                 idxGlobalAll[count] = -1;
1862:               }
1863:             }
1864:           } else {
1865:             /* Right (Middle) Back dummies */
1866:             PetscInt       i;
1867:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1868:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1869:                 idxGlobalAll[count] = -1;
1870:               }
1871:             }
1872:           }
1873:         }
1874:       }

1876:       /* Loop over rows 3/3: Up Back */
1877:       if (!star && !dummyEnd[1]) {
1878:         PetscInt j;
1879:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
1880:           /* Loop over columns 1/3: Left Up Back*/
1881:           if (!dummyStart[0]) {
1882:             const PetscInt neighbor = 6;
1883:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1884:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1885:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1886:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1887:               }
1888:             }
1889:           } else {
1890:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1891:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1892:                 idxGlobalAll[count] = -1;
1893:               }
1894:             }
1895:           }

1897:           /* Loop over columns 2/3: (Middle) Up Back*/
1898:           {
1899:             const PetscInt neighbor = 7;
1900:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1901:               const PetscInt i = ighost - ghostOffsetStart[0];
1902:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1903:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1904:               }
1905:             }
1906:           }

1908:           /* Loop over columns 3/3: Right Up Back*/
1909:           if (!dummyEnd[0]) {
1910:             const PetscInt neighbor = 8;
1911:             PetscInt       i;
1912:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
1913:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1914:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1915:               }
1916:             }
1917:           } else {
1918:             /* Partial dummies on (Middle) Up Back */
1919:             const PetscInt neighbor = 7;
1920:             PetscInt       i,dLocal;
1921:             i = stag->n[0];
1922:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1923:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1924:             }
1925:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1926:               idxGlobalAll[count] = -1;
1927:             }
1928:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1929:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1930:             }
1931:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1932:               idxGlobalAll[count] = -1;
1933:             }
1934:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1935:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1936:             }
1937:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1938:               idxGlobalAll[count] = -1;
1939:             }
1940:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1941:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1942:             }
1943:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1944:               idxGlobalAll[count] = -1;
1945:             }
1946:             ++i;
1947:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1948:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
1949:                 idxGlobalAll[count] = -1;
1950:               }
1951:             }
1952:           }
1953:         }
1954:       } else if (dummyEnd[1]) {
1955:         /* Up Back partial dummy row */
1956:         PetscInt j = stag->n[1];

1958:         if (!star && !dummyStart[0]) {
1959:           /* Up Back partial dummy row: Loop over columns 1/3: Left Up Back, on Left (Middle) Back rank */
1960:           const PetscInt neighbor = 3;
1961:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1962:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1963:             PetscInt       dLocal;
1964:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
1965:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1966:             }

1968:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
1969:               idxGlobalAll[count] = -1;
1970:             }
1971:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
1972:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1973:             }
1974:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
1975:               idxGlobalAll[count] = -1;
1976:             }
1977:           }
1978:         } else {
1979:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1980:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
1981:               idxGlobalAll[count] = -1;
1982:             }
1983:           }
1984:         }

1986:         /* Up Back partial dummy row: Loop over columns 2/3: (Middle) Up Back, on (Middle) (Middle) Back rank */
1987:         {
1988:           const PetscInt neighbor = 4;
1989:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1990:             const PetscInt i = ighost - ghostOffsetStart[0];
1991:             PetscInt       dLocal;
1992:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
1993:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1994:             }
1995:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
1996:               idxGlobalAll[count] = -1;
1997:             }
1998:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
1999:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2000:             }
2001:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2002:               idxGlobalAll[count] = -1;
2003:             }
2004:           }
2005:         }

2007:         /* Up Back partial dummy row: Loop over columns 3/3: Right Up Back, on Right (Middle) Back rank */
2008:         if (!star && !dummyEnd[0]) {
2009:           const PetscInt neighbor = 5;
2010:           PetscInt       i;
2011:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2012:             PetscInt       dLocal;
2013:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2014:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2015:             }
2016:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2017:               idxGlobalAll[count] = -1;
2018:             }
2019:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2020:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2021:             }
2022:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2023:               idxGlobalAll[count] = -1;
2024:             }
2025:           }
2026:         } else if (dummyEnd[0]) {
2027:           /* Up Right Back partial dummy element, on (Middle) (Middle) Back rank */
2028:           const PetscInt neighbor = 4;
2029:           PetscInt       i,dLocal;
2030:           i = stag->n[0];
2031:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2032:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2033:           }
2034:           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 */
2035:             idxGlobalAll[count] = -1;
2036:           }
2037:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2038:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2039:           }
2040:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2041:             idxGlobalAll[count] = -1;
2042:           }
2043:           ++i;
2044:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2045:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2046:               idxGlobalAll[count] = -1;
2047:             }
2048:           }
2049:         } else {
2050:           /* Up Right Back dummies */
2051:           PetscInt i;
2052:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2053:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2054:               idxGlobalAll[count] = -1;
2055:             }
2056:           }
2057:         }
2058:         ++j;
2059:         /* Up Back additional dummy rows */
2060:         for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2061:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2062:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2063:               idxGlobalAll[count] = -1;
2064:             }
2065:           }
2066:         }
2067:       } else {
2068:         /* Up Back dummy rows */
2069:         PetscInt j;
2070:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2071:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2072:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2073:               idxGlobalAll[count] = -1;
2074:             }
2075:           }
2076:         }
2077:       }
2078:     }
2079:   } else {
2080:     for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
2081:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
2082:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2083:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2084:             idxGlobalAll[count] = -1;
2085:           }
2086:         }
2087:       }
2088:     }
2089:   }

2091:   /* Loop over layers 2/3 : (Middle)  */
2092:   {
2093:     for (kghost = ghostOffsetStart[2]; kghost<stag->nGhost[2]-ghostOffsetEnd[2]; ++kghost) {
2094:       const PetscInt k = kghost - ghostOffsetStart[2];

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

2101:           /* Loop over columns 1/3: Left Down (Middle) */
2102:           if (!star && !dummyStart[0]) {
2103:             const PetscInt neighbor = 9;
2104:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2105:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2106:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2107:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2108:               }
2109:             }
2110:           } else {
2111:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2112:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2113:                 idxGlobalAll[count] = -1;
2114:               }
2115:             }
2116:           }

2118:           /* Loop over columns 2/3: (Middle) Down (Middle) */
2119:           {
2120:             const PetscInt neighbor = 10;
2121:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2122:               const PetscInt i = ighost - ghostOffsetStart[0];
2123:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2124:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2125:               }
2126:             }
2127:           }

2129:           /* Loop over columns 3/3: Right Down (Middle) */
2130:           if (!star && !dummyEnd[0]) {
2131:             const PetscInt neighbor = 11;
2132:             PetscInt       i;
2133:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2134:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2135:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2136:               }
2137:             }
2138:           } else if (dummyEnd[0]) {
2139:             /* Partial dummies on (Middle) Down (Middle) neighbor */
2140:             const PetscInt neighbor = 10;
2141:             PetscInt       i,dLocal;
2142:             i = stag->n[0];
2143:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2144:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2145:             }
2146:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2147:               idxGlobalAll[count] = -1;
2148:             }
2149:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2150:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2151:             }
2152:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2153:               idxGlobalAll[count] = -1;
2154:             }
2155:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2156:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2157:             }
2158:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2159:               idxGlobalAll[count] = -1;
2160:             }
2161:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2162:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2163:             }
2164:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2165:               idxGlobalAll[count] = -1;
2166:             }
2167:             ++i;
2168:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2169:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2170:                 idxGlobalAll[count] = -1;
2171:               }
2172:             }
2173:           } else {
2174:             /* Right Down (Middle) dummies */
2175:             PetscInt       i;
2176:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2177:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2178:                 idxGlobalAll[count] = -1;
2179:               }
2180:             }
2181:           }
2182:         }
2183:       } else {
2184:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2185:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2186:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2187:               idxGlobalAll[count] = -1;
2188:             }
2189:           }
2190:         }
2191:       }

2193:       /* Loop over rows 2/3: (Middle) (Middle) */
2194:       {
2195:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2196:           const PetscInt j = jghost - ghostOffsetStart[1];

2198:           /* Loop over columns 1/3: Left (Middle) (Middle) */
2199:           if (!dummyStart[0]) {
2200:             const PetscInt neighbor = 12;
2201:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2202:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2203:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2204:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2205:               }
2206:             }
2207:           } else {
2208:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2209:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2210:                 idxGlobalAll[count] = -1;
2211:               }
2212:             }
2213:           }

2215:           /* Loop over columns 2/3: (Middle) (Middle) (Middle) aka (Here) */
2216:           {
2217:             const PetscInt neighbor = 13;
2218:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2219:               const PetscInt i = ighost - ghostOffsetStart[0];
2220:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2221:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2222:               }
2223:             }
2224:           }

2226:           /* Loop over columns 3/3: Right (Middle) (Middle) */
2227:           if (!dummyEnd[0]) {
2228:             const PetscInt neighbor = 14;
2229:             PetscInt       i;
2230:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2231:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2232:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2233:               }
2234:             }
2235:           } else {
2236:             /* Partial dummies on this rank */
2237:             const PetscInt neighbor = 13;
2238:             PetscInt       i,dLocal;
2239:             i = stag->n[0];
2240:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2241:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2242:             }
2243:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2244:               idxGlobalAll[count] = -1;
2245:             }
2246:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2247:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2248:             }
2249:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2250:               idxGlobalAll[count] = -1;
2251:             }
2252:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2253:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2254:             }
2255:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2256:               idxGlobalAll[count] = -1;
2257:             }
2258:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2259:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2260:             }
2261:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2262:               idxGlobalAll[count] = -1;
2263:             }
2264:             ++i;
2265:             for (;i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2266:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2267:                 idxGlobalAll[count] = -1;
2268:               }
2269:             }
2270:           }
2271:         }
2272:       }

2274:       /* Loop over rows 3/3: Up (Middle) */
2275:       if (!dummyEnd[1]) {
2276:         PetscInt j;
2277:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2279:           /* Loop over columns 1/3: Left Up (Middle) */
2280:           if (!star && !dummyStart[0]) {
2281:             const PetscInt neighbor = 15;
2282:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2283:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2284:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2285:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2286:               }
2287:             }
2288:           } else {
2289:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2290:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2291:                 idxGlobalAll[count] = -1;
2292:               }
2293:             }
2294:           }

2296:           /* Loop over columns 2/3: (Middle) Up (Middle) */
2297:           {
2298:             const PetscInt neighbor = 16;
2299:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2300:               const PetscInt i = ighost - ghostOffsetStart[0];
2301:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2302:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2303:               }
2304:             }
2305:           }

2307:           /* Loop over columns 3/3: Right Up (Middle) */
2308:           if (!star && !dummyEnd[0]) {
2309:             const PetscInt neighbor = 17;
2310:             PetscInt       i;
2311:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2312:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2313:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2314:               }
2315:             }
2316:           } else if (dummyEnd[0]) {
2317:             /* Partial dummies on (Middle) Up (Middle) rank */
2318:             const PetscInt neighbor = 16;
2319:             PetscInt       i,dLocal;
2320:             i = stag->n[0];
2321:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2322:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2323:             }
2324:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2325:               idxGlobalAll[count] = -1;
2326:             }
2327:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2328:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2329:             }
2330:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2331:               idxGlobalAll[count] = -1;
2332:             }
2333:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2334:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2335:             }
2336:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2337:               idxGlobalAll[count] = -1;
2338:             }
2339:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2340:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2341:             }
2342:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2343:               idxGlobalAll[count] = -1;
2344:             }
2345:             ++i;
2346:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2347:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2348:                 idxGlobalAll[count] = -1;
2349:               }
2350:             }
2351:           } else {
2352:             /* Right Up (Middle) dumies */
2353:             PetscInt       i;
2354:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2355:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2356:                 idxGlobalAll[count] = -1;
2357:               }
2358:             }
2359:           }
2360:         }
2361:       } else {
2362:         /* Up (Middle) partial dummy row */
2363:         PetscInt j = stag->n[1];

2365:         /* Up (Middle) partial dummy row: columns 1/3: Left Up (Middle), on Left (Middle) (Middle) rank */
2366:         if (!dummyStart[0]) {
2367:           const PetscInt neighbor = 12;
2368:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2369:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2370:             PetscInt       dLocal;
2371:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2372:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2373:             }
2374:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2375:               idxGlobalAll[count] = -1;
2376:             }
2377:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2378:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2379:             }
2380:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2381:               idxGlobalAll[count] = -1;
2382:             }
2383:           }
2384:         } else {
2385:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2386:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2387:               idxGlobalAll[count] = -1;
2388:             }
2389:           }
2390:         }

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

2413:         if (!dummyEnd[0]) {
2414:           /* Up (Middle) partial dummy row: columns 3/3: Right Up (Middle), on Right (Middle) (Middle) rank */
2415:           const PetscInt neighbor = 14;
2416:           PetscInt       i;
2417:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2418:             PetscInt       dLocal;
2419:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2420:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2421:             }
2422:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2423:               idxGlobalAll[count] = -1;
2424:             }
2425:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2426:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2427:             }
2428:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2429:               idxGlobalAll[count] = -1;
2430:             }
2431:           }
2432:         } else {
2433:           /* Up (Middle) Right partial dummy element, on (Middle) (Middle) (Middle) = here rank */
2434:           const PetscInt neighbor = 13;
2435:           PetscInt       i,dLocal;
2436:           i = stag->n[0];
2437:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2438:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2439:           }
2440:           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 */
2441:             idxGlobalAll[count] = -1;
2442:           }
2443:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2444:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2445:           }
2446:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2447:             idxGlobalAll[count] = -1;
2448:           }
2449:           ++i;
2450:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2451:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2452:               idxGlobalAll[count] = -1;
2453:             }
2454:           }
2455:         }
2456:         ++j;
2457:         /* Up (Middle) additional dummy rows */
2458:         for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2459:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2460:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2461:               idxGlobalAll[count] = -1;
2462:             }
2463:           }
2464:         }
2465:       }
2466:     }
2467:   }

2469:   /* Loop over layers 3/3 : Front */
2470:   if (!dummyEnd[2]) {
2471:     PetscInt k;
2472:     for (k=0; k<ghostOffsetEnd[2]; ++k) {

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

2479:           /* Loop over columns 1/3: Left Down Front */
2480:           if (!dummyStart[0]) {
2481:             const PetscInt neighbor = 18;
2482:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2483:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2484:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2485:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2486:               }
2487:             }
2488:           } else {
2489:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2490:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2491:                 idxGlobalAll[count] = -1;
2492:               }
2493:             }
2494:           }

2496:           /* Loop over columns 2/3: (Middle) Down Front */
2497:           {
2498:             const PetscInt neighbor = 19;
2499:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2500:               const PetscInt i = ighost - ghostOffsetStart[0];
2501:               for (d=0; d<stag->entriesPerElement; ++d,++count) { /* vertex, 2 edges, and face associated with back face */
2502:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2503:               }
2504:             }
2505:           }

2507:           /* Loop over columns 3/3: Right Down Front */
2508:           if (!dummyEnd[0]) {
2509:             const PetscInt neighbor = 20;
2510:             PetscInt       i;
2511:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2512:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2513:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2514:               }
2515:             }
2516:           } else {
2517:             /* Partial dummy element on (Middle) Down Front rank */
2518:             const PetscInt neighbor = 19;
2519:             PetscInt       i,dLocal;
2520:             i = stag->n[0];
2521:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2522:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2523:             }
2524:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2525:               idxGlobalAll[count] = -1;
2526:             }
2527:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2528:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2529:             }
2530:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2531:               idxGlobalAll[count] = -1;
2532:             }
2533:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2534:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2535:             }
2536:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2537:               idxGlobalAll[count] = -1;
2538:             }
2539:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2540:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2541:             }
2542:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2543:               idxGlobalAll[count] = -1;
2544:             }
2545:             ++i;
2546:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2547:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2548:                 idxGlobalAll[count] = -1;
2549:               }
2550:             }
2551:           }
2552:         }
2553:       } else {
2554:         for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2555:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2556:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2557:               idxGlobalAll[count] = -1;
2558:             }
2559:           }
2560:         }
2561:       }

2563:       /* Loop over rows 2/3: (Middle) Front */
2564:       {
2565:         for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2566:           const PetscInt j = jghost - ghostOffsetStart[1];

2568:           /* Loop over columns 1/3: Left (Middle) Front*/
2569:           if (!star && !dummyStart[0]) {
2570:             const PetscInt neighbor = 21;
2571:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2572:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2573:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2574:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2575:               }
2576:             }
2577:           } else {
2578:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2579:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2580:                 idxGlobalAll[count] = -1;
2581:               }
2582:             }
2583:           }

2585:           /* Loop over columns 2/3: (Middle) (Middle) Front*/
2586:           {
2587:             const PetscInt neighbor = 22;
2588:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2589:               const PetscInt i = ighost - ghostOffsetStart[0];
2590:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2591:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2592:               }
2593:             }
2594:           }

2596:           /* Loop over columns 3/3: Right (Middle) Front*/
2597:           if (!star && !dummyEnd[0]) {
2598:             const PetscInt neighbor = 23;
2599:             PetscInt       i;
2600:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2601:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2602:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2603:               }
2604:             }
2605:           } else if (dummyEnd[0]) {
2606:             /* Partial dummy element on (Middle) (Middle) Front element */
2607:             const PetscInt neighbor = 22;
2608:             PetscInt       i,dLocal;
2609:             i = stag->n[0];
2610:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2611:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2612:             }
2613:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2614:               idxGlobalAll[count] = -1;
2615:             }
2616:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2617:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2618:             }
2619:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2620:               idxGlobalAll[count] = -1;
2621:             }
2622:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2623:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2624:             }
2625:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2626:               idxGlobalAll[count] = -1;
2627:             }
2628:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2629:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2630:             }
2631:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2632:               idxGlobalAll[count] = -1;
2633:             }
2634:             ++i;
2635:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2636:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2637:                 idxGlobalAll[count] = -1;
2638:               }
2639:             }
2640:           } else {
2641:             /* Right (Middle) Front dummies */
2642:             PetscInt       i;
2643:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2644:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2645:                 idxGlobalAll[count] = -1;
2646:               }
2647:             }
2648:           }
2649:         }
2650:       }

2652:       /* Loop over rows 3/3: Up Front */
2653:       if (!star && !dummyEnd[1]) {
2654:         PetscInt j;
2655:         for (j=0; j<ghostOffsetEnd[1]; ++j) {

2657:           /* Loop over columns 1/3: Left Up Front */
2658:           if (!dummyStart[0]) {
2659:             const PetscInt neighbor = 24;
2660:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2661:               const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2662:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2663:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2664:               }
2665:             }
2666:           } else {
2667:             for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2668:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2669:                 idxGlobalAll[count] = -1;
2670:               }
2671:             }
2672:           }

2674:           /* Loop over columns 2/3: (Middle) Up Front */
2675:           {
2676:             const PetscInt neighbor = 25;
2677:             for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2678:               const PetscInt i = ighost - ghostOffsetStart[0];
2679:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2680:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2681:               }
2682:             }
2683:           }

2685:           /* Loop over columns 3/3: Right Up Front */
2686:           if (!dummyEnd[0]) {
2687:             const PetscInt neighbor = 26;
2688:             PetscInt     i;
2689:             for (i=0; i<ghostOffsetEnd[0]; ++i) {
2690:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2691:                 idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2692:               }
2693:             }
2694:           } else {
2695:             /* Partial dummy element on (Middle) Up Front rank */
2696:             const PetscInt neighbor = 25;
2697:             PetscInt       i,dLocal;
2698:             i = stag->n[0];
2699:             for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2700:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2701:             }
2702:             for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2703:               idxGlobalAll[count] = -1;
2704:             }
2705:             for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2706:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2707:             }
2708:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2709:               idxGlobalAll[count] = -1;
2710:             }
2711:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2712:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2713:             }
2714:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2715:               idxGlobalAll[count] = -1;
2716:             }
2717:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2718:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2719:             }
2720:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2721:               idxGlobalAll[count] = -1;
2722:             }
2723:             ++i;
2724:             for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2725:               for (d=0; d<stag->entriesPerElement; ++d,++count) {
2726:                 idxGlobalAll[count] = -1;
2727:               }
2728:             }
2729:           }
2730:         }
2731:       } else if (dummyEnd[1]) {
2732:         /* Up Front partial dummy row */
2733:         PetscInt j = stag->n[1];

2735:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) Front rank */
2736:         if (!star && !dummyStart[0]) {
2737:           const PetscInt neighbor = 21;
2738:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2739:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2740:             PetscInt       dLocal;
2741:             for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2742:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2743:             }
2744:             for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2745:               idxGlobalAll[count] = -1;
2746:             }
2747:             for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2748:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2749:             }
2750:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2751:               idxGlobalAll[count] = -1;
2752:             }
2753:           }
2754:         } else {
2755:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2756:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2757:               idxGlobalAll[count] = -1;
2758:             }
2759:           }
2760:         }

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

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

2803:         } else if (dummyEnd[0]) {
2804:           /* Partial Right Up Front dummy element, on (Middle) (Middle) Front rank */
2805:           const PetscInt neighbor = 22;
2806:           PetscInt       i,dLocal;
2807:           i = stag->n[0];
2808:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2809:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2810:           }
2811:           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 */
2812:             idxGlobalAll[count] = -1;
2813:           }
2814:           for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2815:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2816:           }
2817:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2818:             idxGlobalAll[count] = -1;
2819:           }
2820:           ++i;
2821:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2822:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2823:               idxGlobalAll[count] = -1;
2824:             }
2825:           }
2826:         } else {
2827:           /* Right Up Front dummies */
2828:           PetscInt i;
2829:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2830:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2831:               idxGlobalAll[count] = -1;
2832:             }
2833:           }
2834:         }
2835:         ++j;
2836:         /* Up Front additional dummy rows */
2837:         for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2838:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2839:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2840:               idxGlobalAll[count] = -1;
2841:             }
2842:           }
2843:         }
2844:       } else {
2845:         /* Up Front dummies */
2846:         PetscInt j;
2847:         for (j=0; j<ghostOffsetEnd[1]; ++j) {
2848:           for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2849:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2850:               idxGlobalAll[count] = -1;
2851:             }
2852:           }
2853:         }
2854:       }
2855:     }
2856:   } else {
2857:     PetscInt k = stag->n[2];

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

2864:         /* Down Front partial ghost row: columns 1/3: Left Down Front, on  Left Down (Middle) */
2865:         if (!star && !dummyStart[0]) {
2866:           const PetscInt neighbor = 9;
2867:           const PetscInt epFaceRow         = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
2868:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2869:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2870:             PetscInt  dLocal;
2871:             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 */
2872:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2873:             }
2874:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2875:               idxGlobalAll[count] = -1;
2876:             }
2877:           }
2878:         } else {
2879:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2880:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2881:               idxGlobalAll[count] = -1;
2882:             }
2883:           }
2884:         }

2886:         /* Down Front partial ghost row: columns 2/3: (Middle) Down Front, on (Middle) Down (Middle) */
2887:         {
2888:           const PetscInt neighbor = 10;
2889:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2890:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2891:             const PetscInt i = ighost - ghostOffsetStart[0];
2892:             PetscInt  dLocal;
2893:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2894:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2895:             }
2896:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2897:               idxGlobalAll[count] = -1;
2898:             }
2899:           }
2900:         }

2902:         if (!star && !dummyEnd[0]) {
2903:           /* Down Front partial dummy row: columns 3/3: Right Down Front, on Right Down (Middle) */
2904:           const PetscInt neighbor = 11;
2905:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
2906:           PetscInt       i;
2907:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2908:             PetscInt  dLocal;
2909:             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 */
2910:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2911:             }
2912:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2913:               idxGlobalAll[count] = -1;
2914:             }
2915:           }
2916:         } else if (dummyEnd[0]) {
2917:           /* Right Down Front partial dummy element, living on (Middle) Down (Middle) rank */
2918:           const PetscInt neighbor = 10;
2919:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2920:           PetscInt       i,dLocal;
2921:           i = stag->n[0];
2922:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2923:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2924:           }
2925:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2926:             idxGlobalAll[count] = -1;
2927:           }
2928:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
2929:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2930:           }
2931:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
2932:             idxGlobalAll[count] = -1;
2933:           }
2934:           ++i;
2935:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2936:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2937:               idxGlobalAll[count] = -1;
2938:             }
2939:           }
2940:         } else {
2941:           PetscInt i;
2942:           /* Right Down Front dummies */
2943:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
2944:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2945:               idxGlobalAll[count] = -1;
2946:             }
2947:           }
2948:         }
2949:       }
2950:     } else {
2951:       for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2952:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2953:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
2954:             idxGlobalAll[count] = -1;
2955:           }
2956:         }
2957:       }
2958:     }

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

2965:         /* (Middle) Front partial dummy row: columns 1/3: Left (Middle) Front, on Left (Middle) (Middle) */
2966:         if (!dummyStart[0]) {
2967:           const PetscInt neighbor = 12;
2968:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
2969:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2970:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2971:             PetscInt  dLocal;
2972:             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 */
2973:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2974:             }
2975:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2976:               idxGlobalAll[count] = -1;
2977:             }
2978:           }
2979:         } else {
2980:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2981:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
2982:               idxGlobalAll[count] = -1;
2983:             }
2984:           }
2985:         }

2987:         /* (Middle) Front partial dummy row: columns 2/3: (Middle) (Middle) Front, on (Middle) (Middle) (Middle) */
2988:         {
2989:           const PetscInt neighbor = 13;
2990:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
2991:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2992:             const PetscInt i = ighost - ghostOffsetStart[0];
2993:             PetscInt  dLocal;
2994:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2995:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2996:             }
2997:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2998:               idxGlobalAll[count] = -1;
2999:             }
3000:           }
3001:         }

3003:         if (!dummyEnd[0]) {
3004:           /* (Middle) Front partial dummy row: columns 3/3: Right (Middle) Front, on Right (Middle) (Middle) */
3005:           const PetscInt neighbor = 14;
3006:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3007:           PetscInt i;
3008:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3009:             PetscInt  dLocal;
3010:             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 */
3011:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3012:             }
3013:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3014:               idxGlobalAll[count] = -1;
3015:             }
3016:           }
3017:         } else {
3018:           /* Right (Middle) Front partial dummy element, on (Middle) (Middle) (Middle) rank */
3019:           const PetscInt neighbor = 13;
3020:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3021:           PetscInt       i,dLocal;
3022:           i = stag->n[0];
3023:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3024:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3025:           }
3026:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3027:             idxGlobalAll[count] = -1;
3028:           }
3029:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3030:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3031:           }
3032:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3033:             idxGlobalAll[count] = -1;
3034:           }
3035:           ++i;
3036:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3037:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3038:               idxGlobalAll[count] = -1;
3039:             }
3040:           }
3041:         }
3042:       }
3043:     }

3045:     /* Front partial ghost layer: rows 3/3: Up Front, on Up (Middle) */
3046:     if (!dummyEnd[1]) {
3047:       PetscInt j;
3048:       for (j=0; j<ghostOffsetEnd[1]; ++j) {

3050:         /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left Up (Middle) */
3051:         if (!star && !dummyStart[0]) {
3052:           const PetscInt neighbor = 15;
3053:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3054:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3055:             const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3056:             PetscInt  dLocal;
3057:             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 */
3058:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3059:             }
3060:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3061:               idxGlobalAll[count] = -1;
3062:             }
3063:           }
3064:         } else {
3065:           for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3066:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3067:               idxGlobalAll[count] = -1;
3068:             }
3069:           }
3070:         }

3072:         /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) Up (Middle) */
3073:         {
3074:           const PetscInt neighbor = 16;
3075:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3076:           for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3077:             const PetscInt i = ighost - ghostOffsetStart[0];
3078:             PetscInt  dLocal;
3079:             for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3080:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3081:             }
3082:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3083:               idxGlobalAll[count] = -1;
3084:             }
3085:           }
3086:         }

3088:         if (!star && !dummyEnd[0]) {
3089:           /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right Up (Middle) */
3090:           const PetscInt neighbor = 17;
3091:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3092:           PetscInt       i;
3093:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3094:             PetscInt  dLocal;
3095:             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 */
3096:               idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3097:             }
3098:             for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3099:               idxGlobalAll[count] = -1;
3100:             }
3101:           }
3102:         } else if (dummyEnd[0]) {
3103:           /* Right Up Front partial dummy element, on (Middle) Up (Middle) */
3104:           const PetscInt neighbor = 16;
3105:           const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3106:           PetscInt       i,dLocal;
3107:           i = stag->n[0];
3108:           for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3109:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3110:           }
3111:           for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3112:             idxGlobalAll[count] = -1;
3113:           }
3114:           for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3115:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3116:           }
3117:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3118:             idxGlobalAll[count] = -1;
3119:           }
3120:           ++i;
3121:           for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3122:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3123:               idxGlobalAll[count] = -1;
3124:             }
3125:           }
3126:         } else {
3127:           /* Right Up Front dummies */
3128:           PetscInt i;
3129:           /* Right Down Front dummies */
3130:           for (i=0; i<ghostOffsetEnd[0]; ++i) {
3131:             for (d=0; d<stag->entriesPerElement; ++d,++count) {
3132:               idxGlobalAll[count] = -1;
3133:             }
3134:           }
3135:         }
3136:       }
3137:     } else {
3138:       /* Up Front partial dummy row */
3139:       PetscInt j = stag->n[1];

3141:       /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) (Middle) */
3142:       if (!dummyStart[0]) {
3143:         const PetscInt neighbor = 12;
3144:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3145:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3146:           const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3147:           PetscInt       dLocal;
3148:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3149:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3150:           }
3151:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3152:             idxGlobalAll[count] = -1;
3153:           }
3154:         }
3155:       } else {
3156:         for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3157:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3158:             idxGlobalAll[count] = -1;
3159:           }
3160:         }
3161:       }

3163:       /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) (Middle) */
3164:       {
3165:         const PetscInt neighbor = 13;
3166:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3167:         for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3168:           const PetscInt i = ighost - ghostOffsetStart[0];
3169:           PetscInt       dLocal;
3170:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3171:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3172:           }
3173:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3174:             idxGlobalAll[count] = -1;
3175:           }
3176:         }
3177:       }

3179:       if (!dummyEnd[0]) {
3180:         /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) (Middle) */
3181:         const PetscInt neighbor = 14;
3182:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3183:         PetscInt       i;
3184:         for (i=0; i<ghostOffsetEnd[0]; ++i) {
3185:           PetscInt       dLocal;
3186:           for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex  and down back edge */
3187:             idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3188:           }
3189:           for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3190:             idxGlobalAll[count] = -1;
3191:           }
3192:         }
3193:       } else {
3194:         /* Right Up Front partial dummy element, on this rank */
3195:         const PetscInt neighbor = 13;
3196:         const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3197:         PetscInt       i,dLocal;
3198:         i = stag->n[0];
3199:         for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3200:           idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3201:         }
3202:         for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3203:           idxGlobalAll[count] = -1;
3204:         }
3205:         ++i;
3206:         for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3207:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3208:             idxGlobalAll[count] = -1;
3209:           }
3210:         }
3211:       }
3212:       ++j;
3213:       /* Up Front additional dummy rows */
3214:       for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
3215:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3216:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3217:             idxGlobalAll[count] = -1;
3218:           }
3219:         }
3220:       }
3221:     }
3222:     /* Front additional dummy layers */
3223:     ++k;
3224:     for (; k<stag->n[2] + ghostOffsetEnd[2]; ++k) {
3225:       for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
3226:         for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3227:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
3228:             idxGlobalAll[count] = -1;
3229:           }
3230:         }
3231:       }
3232:     }
3233:   }

3235:   /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
3236:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
3237:   PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
3238:   return 0;
3239: }

3241: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM dm)
3242: {
3243:   DM_Stag * const stag = (DM_Stag*)dm->data;
3244:   const PetscInt epe = stag->entriesPerElement;
3245:   const PetscInt epr = stag->nGhost[0]*epe;
3246:   const PetscInt epl = stag->nGhost[1]*epr;

3248:   PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
3249:   stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]   = 0;
3250:   stag->locationOffsets[DMSTAG_BACK_DOWN]        = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + stag->dof[0];
3251:   stag->locationOffsets[DMSTAG_BACK_DOWN_RIGHT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epe;
3252:   stag->locationOffsets[DMSTAG_BACK_LEFT]        = stag->locationOffsets[DMSTAG_BACK_DOWN]       + stag->dof[1];
3253:   stag->locationOffsets[DMSTAG_BACK]             = stag->locationOffsets[DMSTAG_BACK_LEFT]       + stag->dof[1];
3254:   stag->locationOffsets[DMSTAG_BACK_RIGHT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epe;
3255:   stag->locationOffsets[DMSTAG_BACK_UP_LEFT]     = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epr;
3256:   stag->locationOffsets[DMSTAG_BACK_UP]          = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epr;
3257:   stag->locationOffsets[DMSTAG_BACK_UP_RIGHT]    = stag->locationOffsets[DMSTAG_BACK_UP_LEFT]    + epe;
3258:   stag->locationOffsets[DMSTAG_DOWN_LEFT]        = stag->locationOffsets[DMSTAG_BACK]            + stag->dof[2];
3259:   stag->locationOffsets[DMSTAG_DOWN]             = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + stag->dof[1];
3260:   stag->locationOffsets[DMSTAG_DOWN_RIGHT]       = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + epe;
3261:   stag->locationOffsets[DMSTAG_LEFT]             = stag->locationOffsets[DMSTAG_DOWN]            + stag->dof[2];
3262:   stag->locationOffsets[DMSTAG_ELEMENT]          = stag->locationOffsets[DMSTAG_LEFT]            + stag->dof[2];
3263:   stag->locationOffsets[DMSTAG_RIGHT]            = stag->locationOffsets[DMSTAG_LEFT]            + epe;
3264:   stag->locationOffsets[DMSTAG_UP_LEFT]          = stag->locationOffsets[DMSTAG_DOWN_LEFT]       + epr;
3265:   stag->locationOffsets[DMSTAG_UP]               = stag->locationOffsets[DMSTAG_DOWN]            + epr;
3266:   stag->locationOffsets[DMSTAG_UP_RIGHT]         = stag->locationOffsets[DMSTAG_UP_LEFT]         + epe;
3267:   stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT]  = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT]  + epl;
3268:   stag->locationOffsets[DMSTAG_FRONT_DOWN]       = stag->locationOffsets[DMSTAG_BACK_DOWN]       + epl;
3269:   stag->locationOffsets[DMSTAG_FRONT_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epe;
3270:   stag->locationOffsets[DMSTAG_FRONT_LEFT]       = stag->locationOffsets[DMSTAG_BACK_LEFT]       + epl;
3271:   stag->locationOffsets[DMSTAG_FRONT]            = stag->locationOffsets[DMSTAG_BACK]            + epl;
3272:   stag->locationOffsets[DMSTAG_FRONT_RIGHT]      = stag->locationOffsets[DMSTAG_FRONT_LEFT]      + epe;
3273:   stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]    = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epr;
3274:   stag->locationOffsets[DMSTAG_FRONT_UP]         = stag->locationOffsets[DMSTAG_FRONT_DOWN]      + epr;
3275:   stag->locationOffsets[DMSTAG_FRONT_UP_RIGHT]   = stag->locationOffsets[DMSTAG_FRONT_UP_LEFT]   + epe;
3276:   return 0;
3277: }

3279: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_3d(DM dm)
3280: {
3281:   DM_Stag * const stag = (DM_Stag*)dm->data;
3282:   PetscInt        *idxLocal,*idxGlobal,*globalOffsetsRecomputed;
3283:   const PetscInt  *globalOffsets;
3284:   PetscInt        count,d,entriesPerEdge,entriesPerFace,eprGhost,eplGhost,ghostOffsetStart[3],ghostOffsetEnd[3];
3285:   IS              isLocal,isGlobal;
3286:   PetscBool       dummyEnd[3];

3288:   DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsetsRecomputed); /* note that we don't actually use all of these. An available optimization is to pass them, when available */
3289:   globalOffsets = globalOffsetsRecomputed;
3290:   PetscMalloc1(stag->entries,&idxLocal);
3291:   PetscMalloc1(stag->entries,&idxGlobal);
3292:   for (d=0; d<3; ++d) dummyEnd[d]   = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
3293:   entriesPerFace                      = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
3294:   entriesPerEdge                      = stag->dof[0] + stag->dof[1];
3295:   eprGhost                            = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row   */
3296:   eplGhost                            = stag->nGhost[1]*eprGhost;                /* epl = entries per (element) layer */
3297:   count = 0;
3298:   for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
3299:   for (d=0; d<3; ++d) ghostOffsetEnd[d]   = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
3300:   {
3301:     const PetscInt  neighbor     = 13;
3302:     const PetscInt  epr          = stag->entriesPerElement * stag->n[0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
3303:     const PetscInt  epl          = epr                     * stag->n[1] + (dummyEnd[1] ? stag->n[0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
3304:     const PetscInt  epFaceRow    = entriesPerFace          * stag->n[0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3305:     const PetscInt  start0       = 0;
3306:     const PetscInt  start1       = 0;
3307:     const PetscInt  start2       = 0;
3308:     const PetscInt  startGhost0  = ghostOffsetStart[0];
3309:     const PetscInt  startGhost1  = ghostOffsetStart[1];
3310:     const PetscInt  startGhost2  = ghostOffsetStart[2];
3311:     const PetscInt  endGhost0    = stag->nGhost[0] - ghostOffsetEnd[0];
3312:     const PetscInt  endGhost1    = stag->nGhost[1] - ghostOffsetEnd[1];
3313:     const PetscInt  endGhost2    = stag->nGhost[2] - ghostOffsetEnd[2];
3314:     const PetscBool extra0       = dummyEnd[0];
3315:     const PetscBool extra1       = dummyEnd[1];
3316:     const PetscBool extra2       = dummyEnd[2];
3317:     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);
3318:   }
3319:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
3320:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
3321:   {
3322:     Vec local,global;
3323:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
3324:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
3325:     VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
3326:     VecDestroy(&global);
3327:     VecDestroy(&local);
3328:   }
3329:   ISDestroy(&isLocal);
3330:   ISDestroy(&isGlobal);
3331:   if (globalOffsetsRecomputed) {
3332:     PetscFree(globalOffsetsRecomputed);
3333:   }
3334:   return 0;
3335: }

3337: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_3D_AIJ_Assemble(DM dm,Mat A)
3338: {
3339:   PetscInt          dof[DMSTAG_MAX_STRATA],epe,stencil_width,N[3],start[3],n[3],n_extra[3];
3340:   DMStagStencilType stencil_type;
3341:   DMBoundaryType    boundary_type[3];

3343:   /* This implementation gives a very dense stencil, which is likely unsuitable for
3344:      (typical) applications which have fewer couplings */
3345:   DMStagGetDOF(dm,&dof[0],&dof[1],&dof[2],&dof[3]);
3346:   DMStagGetStencilType(dm,&stencil_type);
3347:   DMStagGetStencilWidth(dm,&stencil_width);
3348:   DMStagGetEntriesPerElement(dm,&epe);
3349:   DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&n_extra[0],&n_extra[1],&n_extra[2]);
3350:   DMStagGetGlobalSizes(dm,&N[0],&N[1],&N[2]);
3351:   DMStagGetBoundaryTypes(dm,&boundary_type[0],&boundary_type[1],&boundary_type[2]);

3353:   if (stencil_type == DMSTAG_STENCIL_NONE) {
3354:     /* Couple all DOF at each location to each other */
3355:     DMStagStencil *row_vertex,*row_edge_down_left,*row_edge_back_down,*row_edge_back_left,*row_face_down,*row_face_left,*row_face_back,*row_element;

3357:     PetscMalloc1(dof[0],&row_vertex);
3358:     for (PetscInt c=0; c<dof[0]; ++c) {
3359:       row_vertex[c].loc = DMSTAG_BACK_DOWN_LEFT;
3360:       row_vertex[c].c = c;
3361:     }

3363:     PetscMalloc1(dof[1],&row_edge_down_left);
3364:     for (PetscInt c=0; c<dof[1]; ++c) {
3365:       row_edge_down_left[c].loc = DMSTAG_DOWN_LEFT;
3366:       row_edge_down_left[c].c = c;
3367:     }

3369:     PetscMalloc1(dof[1],&row_edge_back_left);
3370:     for (PetscInt c=0; c<dof[1]; ++c) {
3371:       row_edge_back_left[c].loc = DMSTAG_BACK_LEFT;
3372:       row_edge_back_left[c].c = c;
3373:     }

3375:     PetscMalloc1(dof[1],&row_edge_back_down);
3376:     for (PetscInt c=0; c<dof[1]; ++c) {
3377:       row_edge_back_down[c].loc = DMSTAG_BACK_DOWN;
3378:       row_edge_back_down[c].c = c;
3379:     }

3381:     PetscMalloc1(dof[2],&row_face_left);
3382:     for (PetscInt c=0; c<dof[2]; ++c) {
3383:       row_face_left[c].loc = DMSTAG_LEFT;
3384:       row_face_left[c].c = c;
3385:     }

3387:     PetscMalloc1(dof[2],&row_face_down);
3388:     for (PetscInt c=0; c<dof[2]; ++c) {
3389:       row_face_down[c].loc = DMSTAG_DOWN;
3390:       row_face_down[c].c = c;
3391:     }

3393:     PetscMalloc1(dof[2],&row_face_back);
3394:     for (PetscInt c=0; c<dof[2]; ++c) {
3395:       row_face_back[c].loc = DMSTAG_BACK;
3396:       row_face_back[c].c = c;
3397:     }

3399:     PetscMalloc1(dof[3],&row_element);
3400:     for (PetscInt c=0; c<dof[3]; ++c) {
3401:       row_element[c].loc = DMSTAG_ELEMENT;
3402:       row_element[c].c = c;
3403:     }

3405:     for (PetscInt ez=start[2]; ez<start[2]+n[2]+n_extra[2]; ++ez) {
3406:       for (PetscInt ey=start[1]; ey<start[1]+n[1]+n_extra[1]; ++ey) {
3407:         for (PetscInt ex=start[0]; ex<start[0]+n[0]+n_extra[0]; ++ex) {
3408:           for (PetscInt c=0; c<dof[0]; ++c) {
3409:             row_vertex[c].i = ex;
3410:             row_vertex[c].j = ey;
3411:             row_vertex[c].k = ez;
3412:           }
3413:           DMStagMatSetValuesStencil(dm,A,dof[0],row_vertex,dof[0],row_vertex,NULL,INSERT_VALUES);

3415:           if (ez < N[2]) {
3416:             for (PetscInt c=0; c<dof[1]; ++c) {
3417:               row_edge_down_left[c].i = ex;
3418:               row_edge_down_left[c].j = ey;
3419:               row_edge_down_left[c].k = ez;
3420:             }
3421:             DMStagMatSetValuesStencil(dm,A,dof[1],row_edge_down_left,dof[1],row_edge_down_left,NULL,INSERT_VALUES);
3422:           }

3424:           if (ey < N[1]) {
3425:             for (PetscInt c=0; c<dof[1]; ++c) {
3426:               row_edge_back_left[c].i = ex;
3427:               row_edge_back_left[c].j = ey;
3428:               row_edge_back_left[c].k = ez;
3429:             }
3430:             DMStagMatSetValuesStencil(dm,A,dof[1],row_edge_back_left,dof[1],row_edge_back_left,NULL,INSERT_VALUES);
3431:           }

3433:           if (ey < N[0]) {
3434:             for (PetscInt c=0; c<dof[1]; ++c) {
3435:               row_edge_back_down[c].i = ex;
3436:               row_edge_back_down[c].j = ey;
3437:               row_edge_back_down[c].k = ez;
3438:             }
3439:             DMStagMatSetValuesStencil(dm,A,dof[1],row_edge_back_down,dof[1],row_edge_back_down,NULL,INSERT_VALUES);
3440:           }

3442:           if (ey < N[1] && ez < N[2]) {
3443:             for (PetscInt c=0; c<dof[2]; ++c) {
3444:               row_face_left[c].i = ex;
3445:               row_face_left[c].j = ey;
3446:               row_face_left[c].k = ez;
3447:             }
3448:             DMStagMatSetValuesStencil(dm,A,dof[2],row_face_left,dof[2],row_face_left,NULL,INSERT_VALUES);
3449:           }

3451:           if (ex < N[0] && ez < N[2]) {
3452:             for (PetscInt c=0; c<dof[2]; ++c) {
3453:               row_face_down[c].i = ex;
3454:               row_face_down[c].j = ey;
3455:               row_face_down[c].k = ez;
3456:             }
3457:             DMStagMatSetValuesStencil(dm,A,dof[2],row_face_down,dof[2],row_face_down,NULL,INSERT_VALUES);
3458:           }

3460:           if (ex < N[0] && ey < N[1]) {
3461:             for (PetscInt c=0; c<dof[2]; ++c) {
3462:               row_face_back[c].i = ex;
3463:               row_face_back[c].j = ey;
3464:               row_face_back[c].k = ez;
3465:             }
3466:             DMStagMatSetValuesStencil(dm,A,dof[2],row_face_back,dof[2],row_face_back,NULL,INSERT_VALUES);
3467:           }

3469:           if (ex < N[0] && ey < N[1] && ez < N[2]) {
3470:             for (PetscInt c=0; c<dof[3]; ++c){
3471:               row_element[c].i = ex;
3472:               row_element[c].j = ey;
3473:               row_element[c].k = ez;
3474:             }
3475:             DMStagMatSetValuesStencil(dm,A,dof[3],row_element,dof[3],row_element,NULL,INSERT_VALUES);
3476:           }
3477:         }
3478:       }
3479:     }
3480:     PetscFree(row_vertex);
3481:     PetscFree(row_edge_back_left);
3482:     PetscFree(row_edge_back_down);
3483:     PetscFree(row_edge_down_left);
3484:     PetscFree(row_face_left);
3485:     PetscFree(row_face_back);
3486:     PetscFree(row_face_down);
3487:     PetscFree(row_element);
3488:   } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
3489:     DMStagStencil *col,*row;

3491:     PetscMalloc1(epe,&row);
3492:     {
3493:       PetscInt nrows = 0;

3495:       for (PetscInt c=0; c<dof[0]; ++c) {
3496:         row[nrows].c = c;
3497:         row[nrows].loc = DMSTAG_BACK_DOWN_LEFT;
3498:         ++nrows;
3499:       }
3500:       for (PetscInt c=0; c<dof[1]; ++c) {
3501:         row[nrows].c = c;
3502:         row[nrows].loc = DMSTAG_DOWN_LEFT;
3503:         ++nrows;
3504:       }
3505:       for (PetscInt c=0; c<dof[1]; ++c) {
3506:         row[nrows].c = c;
3507:         row[nrows].loc = DMSTAG_BACK_LEFT;
3508:         ++nrows;
3509:       }
3510:       for (PetscInt c=0; c<dof[1]; ++c) {
3511:         row[nrows].c = c;
3512:         row[nrows].loc = DMSTAG_BACK_DOWN;
3513:         ++nrows;
3514:       }
3515:       for (PetscInt c=0; c<dof[2]; ++c) {
3516:         row[nrows].c = c;
3517:         row[nrows].loc = DMSTAG_LEFT;
3518:         ++nrows;
3519:       }
3520:       for (PetscInt c=0; c<dof[2]; ++c) {
3521:         row[nrows].c = c;
3522:         row[nrows].loc = DMSTAG_DOWN;
3523:         ++nrows;
3524:       }
3525:       for (PetscInt c=0; c<dof[2]; ++c) {
3526:         row[nrows].c = c;
3527:         row[nrows].loc = DMSTAG_BACK;
3528:         ++nrows;
3529:       }
3530:       for (PetscInt c=0; c<dof[3]; ++c) {
3531:         row[nrows].c = c;
3532:         row[nrows].loc = DMSTAG_ELEMENT;
3533:         ++nrows;
3534:       }
3535:     }

3537:     PetscMalloc1(epe,&col);
3538:     {
3539:       PetscInt ncols = 0;

3541:       for (PetscInt c=0; c<dof[0]; ++c) {
3542:         col[ncols].c = c;
3543:         col[ncols].loc = DMSTAG_BACK_DOWN_LEFT;
3544:         ++ncols;
3545:       }
3546:       for (PetscInt c=0; c<dof[1]; ++c) {
3547:         col[ncols].c = c;
3548:         col[ncols].loc = DMSTAG_DOWN_LEFT;
3549:         ++ncols;
3550:       }
3551:       for (PetscInt c=0; c<dof[1]; ++c) {
3552:         col[ncols].c = c;
3553:         col[ncols].loc = DMSTAG_BACK_LEFT;
3554:         ++ncols;
3555:       }
3556:       for (PetscInt c=0; c<dof[1]; ++c) {
3557:         col[ncols].c = c;
3558:         col[ncols].loc = DMSTAG_BACK_DOWN;
3559:         ++ncols;
3560:       }
3561:       for (PetscInt c=0; c<dof[2]; ++c) {
3562:         col[ncols].c = c;
3563:         col[ncols].loc = DMSTAG_LEFT;
3564:         ++ncols;
3565:       }
3566:       for (PetscInt c=0; c<dof[2]; ++c) {
3567:         col[ncols].c = c;
3568:         col[ncols].loc = DMSTAG_DOWN;
3569:         ++ncols;
3570:       }
3571:       for (PetscInt c=0; c<dof[2]; ++c) {
3572:         col[ncols].c = c;
3573:         col[ncols].loc = DMSTAG_BACK;
3574:         ++ncols;
3575:       }
3576:       for (PetscInt c=0; c<dof[3]; ++c) {
3577:         col[ncols].c = c;
3578:         col[ncols].loc = DMSTAG_ELEMENT;
3579:         ++ncols;
3580:       }
3581:     }

3583:     for (PetscInt ez=start[2]; ez<start[2]+n[2]+n_extra[2]; ++ez) {
3584:       for (PetscInt ey=start[1]; ey<start[1]+n[1]+n_extra[1]; ++ey) {
3585:         for (PetscInt ex=start[0]; ex<start[0]+n[0]+n_extra[0]; ++ex) {
3586:           for (PetscInt i=0; i<epe; ++i) {
3587:             row[i].i = ex;
3588:             row[i].j = ey;
3589:             row[i].k = ez;
3590:           }
3591:           for (PetscInt offset_z = -stencil_width; offset_z<=stencil_width; ++offset_z) {
3592:             const PetscInt ez_offset = ez + offset_z;
3593:             for (PetscInt offset_y = -stencil_width; offset_y<=stencil_width; ++offset_y) {
3594:               const PetscInt ey_offset = ey + offset_y;
3595:               for (PetscInt offset_x = -stencil_width; offset_x<=stencil_width; ++offset_x) {
3596:                 const PetscInt ex_offset = ex + offset_x;
3597:                 const PetscBool is_star_point = (PetscBool) (((offset_x == 0) && (offset_y == 0 || offset_z == 0)) || (offset_y == 0 && offset_z == 0));
3598:                 /* Only set values corresponding to elements which can have non-dummy entries,
3599:                    meaning those that map to unknowns in the global representation. In the periodic
3600:                    case, this is the entire stencil, but in all other cases, only includes a single
3601:                    "extra" element which is partially outside the physical domain (those points in the
3602:                    global representation */
3603:                 if ((stencil_type == DMSTAG_STENCIL_BOX || is_star_point) &&
3604:                     (boundary_type[0] == DM_BOUNDARY_PERIODIC || (ex_offset < N[0]+1 && ex_offset >= 0)) &&
3605:                     (boundary_type[1] == DM_BOUNDARY_PERIODIC || (ey_offset < N[1]+1 && ey_offset >= 0)) &&
3606:                     (boundary_type[2] == DM_BOUNDARY_PERIODIC || (ez_offset < N[2]+1 && ez_offset >= 0)))
3607:                 {
3608:                   for (PetscInt i=0; i<epe; ++i) {
3609:                     col[i].i = ex_offset;
3610:                     col[i].j = ey_offset;
3611:                     col[i].k = ez_offset;
3612:                   }
3613:                   DMStagMatSetValuesStencil(dm,A,epe,row,epe,col,NULL,INSERT_VALUES);
3614:                 }
3615:               }
3616:             }
3617:           }
3618:         }
3619:       }
3620:     }
3621:     PetscFree(row);
3622:     PetscFree(col);
3623:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported stencil type %s",DMStagStencilTypes[stencil_type]);
3624:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
3625:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
3626:   return 0;
3627: }