Actual source code: stag3d.c
petsc-3.13.6 2020-09-29
1: /* Functions specific to the 3-dimensional implementation of DMStag */
2: #include <petsc/private/dmstagimpl.h>
4: /*@C
5: DMStagCreate3d - Create an object to manage data living on the faces, edges, and vertices of a parallelized regular 3D grid.
7: Collective
9: Input Parameters:
10: + comm - MPI communicator
11: . bndx,bndy,bndz - boundary type: DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, or DM_BOUNDARY_GHOSTED
12: . M,N,P - global number of grid points in x,y directions
13: . m,n,p - number of ranks in the x,y directions (may be PETSC_DECIDE)
14: . dof0 - number of degrees of freedom per vertex/point/node/0-cell
15: . dof1 - number of degrees of freedom per edge/1-cell
16: . dof2 - number of degrees of freedom per face/2-cell
17: . dof3 - number of degrees of freedom per element/3-cell
18: . stencilType - ghost/halo region type: DMSTAG_STENCIL_NONE, DMSTAG_STENCIL_BOX, or DMSTAG_STENCIL_STAR
19: . stencilWidth - width, in elements, of halo/ghost region
20: - lx,ly,lz - array sof local x,y,z element counts, of length equal to m,n,p, summing to M,N,P
22: Output Parameter:
23: . dm - the new DMStag object
25: Options Database Keys:
26: + -dm_view - calls DMViewFromOptions() a the conclusion of DMSetUp()
27: . -stag_grid_x <nx> - number of elements in the x direction
28: . -stag_grid_y <ny> - number of elements in the y direction
29: . -stag_grid_z <nz> - number of elements in the z direction
30: . -stag_ranks_x <rx> - number of ranks in the x direction
31: . -stag_ranks_y <ry> - number of ranks in the y direction
32: . -stag_ranks_z <rz> - number of ranks in the z direction
33: . -stag_ghost_stencil_width - width of ghost region, in elements
34: . -stag_boundary_type x <none,ghosted,periodic> - DMBoundaryType value
35: . -stag_boundary_type y <none,ghosted,periodic> - DMBoundaryType value
36: - -stag_boundary_type z <none,ghosted,periodic> - DMBoundaryType value
38: Notes:
39: You must call DMSetUp() after this call before using the DM.
40: If you wish to use the options database (see the keys above) to change values in the DMStag, you must call
41: DMSetFromOptions() after this function but before DMSetUp().
43: Level: beginner
45: .seealso: DMSTAG, DMStagCreate1d(), DMStagCreate2d(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateLocalVector(), DMLocalToGlobalBegin(), DMDACreate3d()
46: @*/
47: PETSC_EXTERN PetscErrorCode DMStagCreate3d(MPI_Comm comm,DMBoundaryType bndx,DMBoundaryType bndy,DMBoundaryType bndz,PetscInt M,PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DMStagStencilType stencilType,PetscInt stencilWidth,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[],DM* dm)
48: {
52: DMCreate(comm,dm);
53: DMSetDimension(*dm,3);
54: DMStagInitialize(bndx,bndy,bndz,M,N,P,m,n,p,dof0,dof1,dof2,dof3,stencilType,stencilWidth,lx,ly,lz,*dm);
55: return(0);
56: }
58: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_3d(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
59: {
61: DM_Stag *stagCoord;
62: DM dmCoord;
63: Vec coordLocal,coord;
64: PetscReal h[3],min[3];
65: PetscScalar ****arr;
66: PetscInt ind[3],start[3],n[3],nExtra[3],s,c;
67: PetscInt ibackdownleft,ibackdown,ibackleft,iback,idownleft,idown,ileft,ielement;
70: DMGetCoordinateDM(dm,&dmCoord);
71: stagCoord = (DM_Stag*) dmCoord->data;
72: for (s=0; s<4; ++s) {
73: if (stagCoord->dof[s] !=0 && stagCoord->dof[s] != 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Coordinate DM in 3 dimensions must have 0 or 3 dof on each stratum, but stratum %d has %d dof",s,stagCoord->dof[s]);
74: }
75: DMGetLocalVector(dmCoord,&coordLocal);
76: DMStagVecGetArray(dmCoord,coordLocal,&arr);
77: if (stagCoord->dof[0]) {
78: DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN_LEFT,0,&ibackdownleft);
79: }
80: if (stagCoord->dof[1]) {
81: DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_DOWN ,0,&ibackdown );
82: DMStagGetLocationSlot(dmCoord,DMSTAG_BACK_LEFT ,0,&ibackleft );
83: DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN_LEFT ,0,&idownleft );
84: }
85: if (stagCoord->dof[2]) {
86: DMStagGetLocationSlot(dmCoord,DMSTAG_BACK ,0,&iback );
87: DMStagGetLocationSlot(dmCoord,DMSTAG_DOWN ,0,&idown );
88: DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT ,0,&ileft );
89: }
90: if (stagCoord->dof[3]) {
91: DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT ,0,&ielement );
92: DMStagGetCorners(dmCoord,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
93: }
94: min[0] = xmin; min[1]= ymin; min[2] = zmin;
95: h[0] = (xmax-xmin)/stagCoord->N[0];
96: h[1] = (ymax-ymin)/stagCoord->N[1];
97: h[2] = (zmax-zmin)/stagCoord->N[2];
99: for(ind[2]=start[2]; ind[2]<start[2] + n[2] + nExtra[2]; ++ind[2]) {
100: for(ind[1]=start[1]; ind[1]<start[1] + n[1] + nExtra[1]; ++ind[1]) {
101: for(ind[0]=start[0]; ind[0]<start[0] + n[0] + nExtra[0]; ++ind[0]) {
102: if (stagCoord->dof[0]) {
103: const PetscReal offs[3] = {0.0,0.0,0.0};
104: for (c=0; c<3; ++c) {
105: arr[ind[2]][ind[1]][ind[0]][ibackdownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
106: }
107: }
108: if (stagCoord->dof[1]) {
109: const PetscReal offs[3] = {0.5,0.0,0.0};
110: for (c=0; c<3; ++c) {
111: arr[ind[2]][ind[1]][ind[0]][ibackdown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
112: }
113: }
114: if (stagCoord->dof[1]) {
115: const PetscReal offs[3] = {0.0,0.5,0.0};
116: for (c=0; c<3; ++c) {
117: arr[ind[2]][ind[1]][ind[0]][ibackleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
118: }
119: }
120: if (stagCoord->dof[2]) {
121: const PetscReal offs[3] = {0.5,0.5,0.0};
122: for (c=0; c<3; ++c) {
123: arr[ind[2]][ind[1]][ind[0]][iback + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
124: }
125: }
126: if (stagCoord->dof[1]) {
127: const PetscReal offs[3] = {0.0,0.0,0.5};
128: for (c=0; c<3; ++c) {
129: arr[ind[2]][ind[1]][ind[0]][idownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
130: }
131: }
132: if (stagCoord->dof[2]) {
133: const PetscReal offs[3] = {0.5,0.0,0.5};
134: for (c=0; c<3; ++c) {
135: arr[ind[2]][ind[1]][ind[0]][idown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
136: }
137: }
138: if (stagCoord->dof[2]) {
139: const PetscReal offs[3] = {0.0,0.5,0.5};
140: for (c=0; c<3; ++c) {
141: arr[ind[2]][ind[1]][ind[0]][ileft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
142: }
143: }
144: if (stagCoord->dof[3]) {
145: const PetscReal offs[3] = {0.5,0.5,0.5};
146: for (c=0; c<3; ++c) {
147: arr[ind[2]][ind[1]][ind[0]][ielement + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
148: }
149: }
150: }
151: }
152: }
153: DMStagVecRestoreArray(dmCoord,coordLocal,&arr);
154: DMCreateGlobalVector(dmCoord,&coord);
155: DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);
156: DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);
157: DMSetCoordinates(dm,coord);
158: PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);
159: DMRestoreLocalVector(dmCoord,&coordLocal);
160: VecDestroy(&coord);
161: return(0);
162: }
164: /* Helper functions used in DMSetUp_Stag() */
165: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM);
166: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM);
167: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM,PetscInt**);
168: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM,const PetscInt*);
169: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM,const PetscInt*);
170: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM);
172: PETSC_INTERN PetscErrorCode DMSetUp_Stag_3d(DM dm)
173: {
174: PetscErrorCode ierr;
175: DM_Stag * const stag = (DM_Stag*)dm->data;
176: PetscMPIInt rank;
177: PetscInt i,j,d;
178: PetscInt *globalOffsets;
179: const PetscInt dim = 3;
182: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
184: /* Rank grid sizes (populates stag->nRanks) */
185: DMStagSetUpBuildRankGrid_3d(dm);
187: /* Determine location of rank in grid */
188: stag->rank[0] = rank % stag->nRanks[0];
189: stag->rank[1] = rank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0];
190: stag->rank[2] = rank / (stag->nRanks[0] * stag->nRanks[1]);
191: for(d=0; d<dim; ++d) {
192: stag->firstRank[d] = PetscNot(stag->rank[d]);
193: stag->lastRank[d] = (PetscBool)(stag->rank[d] == stag->nRanks[d]-1);
194: }
196: /* Determine locally owned region (if not already prescribed).
197: Divide equally, giving lower ranks in each dimension and extra element if needbe.
198: Note that this uses O(P) storage. If this ever becomes an issue, this could
199: be refactored to not keep this data around. */
200: for (i=0; i<dim; ++i) {
201: if (!stag->l[i]) {
202: const PetscInt Ni = stag->N[i], nRanksi = stag->nRanks[i];
203: PetscMalloc1(stag->nRanks[i],&stag->l[i]);
204: for (j=0; j<stag->nRanks[i]; ++j){
205: stag->l[i][j] = Ni/nRanksi + ((Ni % nRanksi) > j);
206: }
207: }
208: }
210: /* Retrieve local size in stag->n */
211: for (i=0; i<dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
212: #if defined(PETSC_USE_DEBUG)
213: for (i=0; i<dim; ++i) {
214: PetscInt Ncheck,j;
215: Ncheck = 0;
216: for (j=0; j<stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
217: if (Ncheck != stag->N[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local sizes in dimension %d don't add up. %d != %d\n",i,Ncheck,stag->N[i]);
218: }
219: #endif
221: /* Compute starting elements */
222: for (i=0; i<dim; ++i) {
223: stag->start[i] = 0;
224: for(j=0;j<stag->rank[i];++j){
225: stag->start[i] += stag->l[i][j];
226: }
227: }
229: /* Determine ranks of neighbors */
230: DMStagSetUpBuildNeighbors_3d(dm);
232: /* Define useful sizes */
233: {
234: PetscInt entriesPerEdge,entriesPerFace,entriesPerCorner,entriesPerElementRow,entriesPerElementLayer;
235: PetscBool dummyEnd[3];
236: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
237: stag->entriesPerElement = stag->dof[0] + 3*stag->dof[1] + 3*stag->dof[2] + stag->dof[3];
238: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
239: entriesPerEdge = stag->dof[0] + stag->dof[1];
240: entriesPerCorner = stag->dof[0];
241: entriesPerElementRow = stag->n[0]*stag->entriesPerElement
242: + (dummyEnd[0] ? entriesPerFace : 0);
243: entriesPerElementLayer = stag->n[1]*entriesPerElementRow
244: + (dummyEnd[1] ? stag->n[0] * entriesPerFace : 0)
245: + (dummyEnd[1] && dummyEnd[0] ? entriesPerEdge : 0);
246: stag->entries = stag->n[2]*entriesPerElementLayer
247: + (dummyEnd[2] ? stag->n[0]*stag->n[1]*entriesPerFace : 0)
248: + (dummyEnd[2] && dummyEnd[0] ? stag->n[1]*entriesPerEdge : 0)
249: + (dummyEnd[2] && dummyEnd[1] ? stag->n[0]*entriesPerEdge : 0)
250: + (dummyEnd[2] && dummyEnd[1] && dummyEnd[0] ? entriesPerCorner : 0);
251: }
253: /* Check that we will not overflow 32 bit indices (slightly overconservative) */
254: #if !defined(PETSC_USE_64BIT_INDICES)
255: if (((PetscInt64) stag->n[0])*((PetscInt64) stag->n[1])*((PetscInt64) stag->n[2])*((PetscInt64) stag->entriesPerElement) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ4(PetscObjectComm((PetscObject)dm),PETSC_ERR_INT_OVERFLOW,"Mesh of %D x %D x %D with %D entries per (interior) element is likely too large for 32 bit indices",stag->n[0],stag->n[1],stag->n[2],stag->entriesPerElement);
256: #endif
258: /* Compute offsets for each rank into global vectors
260: This again requires O(P) storage, which could be replaced with some global
261: communication.
262: */
263: DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsets);
265: for (d=0; d<dim; ++d) if (stag->boundaryType[d] != DM_BOUNDARY_NONE && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->boundaryType[d] != DM_BOUNDARY_GHOSTED) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported boundary type");
267: /* Define ghosted/local sizes */
268: for (d=0; d<dim; ++d) {
269: switch (stag->boundaryType[d]) {
270: case DM_BOUNDARY_NONE:
271: /* Note: for a elements-only DMStag, the extra elements on the edges aren't necessary but we include them anyway */
272: switch (stag->stencilType) {
273: case DMSTAG_STENCIL_NONE : /* only the extra one on the right/top edges */
274: stag->nGhost[d] = stag->n[d];
275: stag->startGhost[d] = stag->start[d];
276: if (stag->lastRank[d]) stag->nGhost[d] += 1;
277: break;
278: case DMSTAG_STENCIL_STAR : /* allocate the corners but don't use them */
279: case DMSTAG_STENCIL_BOX :
280: stag->nGhost[d] = stag->n[d];
281: stag->startGhost[d] = stag->start[d];
282: if (!stag->firstRank[d]) {
283: stag->nGhost[d] += stag->stencilWidth; /* add interior ghost elements */
284: stag->startGhost[d] -= stag->stencilWidth;
285: }
286: if (!stag->lastRank[d]) {
287: stag->nGhost[d] += stag->stencilWidth; /* add interior ghost elements */
288: } else {
289: stag->nGhost[d] += 1; /* one element on the boundary to complete blocking */
290: }
291: break;
292: default :
293: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
294: }
295: break;
296: case DM_BOUNDARY_GHOSTED:
297: switch (stag->stencilType) {
298: case DMSTAG_STENCIL_NONE :
299: stag->startGhost[d] = stag->start[d];
300: stag->nGhost[d] = stag->n[d] + (stag->lastRank[d] ? 1 : 0);
301: break;
302: case DMSTAG_STENCIL_STAR :
303: case DMSTAG_STENCIL_BOX :
304: stag->startGhost[d] = stag->start[d] - stag->stencilWidth; /* This value may be negative */
305: stag->nGhost[d] = stag->n[d] + 2*stag->stencilWidth + (stag->lastRank[d] && stag->stencilWidth == 0 ? 1 : 0);
306: break;
307: default :
308: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
309: }
310: break;
311: case DM_BOUNDARY_PERIODIC:
312: switch (stag->stencilType) {
313: case DMSTAG_STENCIL_NONE : /* only the extra one on the right/top edges */
314: stag->nGhost[d] = stag->n[d];
315: stag->startGhost[d] = stag->start[d];
316: break;
317: case DMSTAG_STENCIL_STAR :
318: case DMSTAG_STENCIL_BOX :
319: stag->nGhost[d] = stag->n[d] + 2*stag->stencilWidth;
320: stag->startGhost[d] = stag->start[d] - stag->stencilWidth;
321: break;
322: default :
323: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
324: }
325: break;
326: default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported boundary type in dimension %D",d);
327: }
328: }
329: stag->entriesGhost = stag->nGhost[0]*stag->nGhost[1]*stag->nGhost[2]*stag->entriesPerElement;
331: /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping */
332: DMStagSetUpBuildScatter_3d(dm,globalOffsets);
333: DMStagSetUpBuildL2G_3d(dm,globalOffsets);
335: /* In special cases, create a dedicated injective local-to-global map */
336: if ((stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) ||
337: (stag->boundaryType[1] == DM_BOUNDARY_PERIODIC && stag->nRanks[1] == 1) ||
338: (stag->boundaryType[2] == DM_BOUNDARY_PERIODIC && stag->nRanks[2] == 1)) {
339: DMStagPopulateLocalToGlobalInjective(dm);
340: }
342: /* Free global offsets */
343: PetscFree(globalOffsets);
345: /* Precompute location offsets */
346: DMStagComputeLocationOffsets_3d(dm);
348: /* View from Options */
349: DMViewFromOptions(dm,NULL,"-dm_view");
351: return(0);
352: }
354: /* adapted from da3.c */
355: static PetscErrorCode DMStagSetUpBuildRankGrid_3d(DM dm)
356: {
357: PetscErrorCode ierr;
358: PetscMPIInt rank,size;
359: PetscInt m,n,p,pm;
360: DM_Stag * const stag = (DM_Stag*)dm->data;
361: const PetscInt M = stag->N[0];
362: const PetscInt N = stag->N[1];
363: const PetscInt P = stag->N[2];
366: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
367: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
369: m = stag->nRanks[0];
370: n = stag->nRanks[1];
371: p = stag->nRanks[2];
373: if (m != PETSC_DECIDE) {
374: if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
375: else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
376: }
377: if (n != PETSC_DECIDE) {
378: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
379: else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
380: }
381: if (p != PETSC_DECIDE) {
382: if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p);
383: else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size);
384: }
385: if ((m > 0) && (n > 0) && (p > 0) && (m*n*p != size)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m %D * n %D * p %D != size %d",m,n,p,size);
387: /* Partition the array among the processors */
388: if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
389: m = size/(n*p);
390: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
391: n = size/(m*p);
392: } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
393: p = size/(m*n);
394: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
395: /* try for squarish distribution */
396: m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
397: if (!m) m = 1;
398: while (m > 0) {
399: n = size/(m*p);
400: if (m*n*p == size) break;
401: m--;
402: }
403: if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p);
404: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
405: } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
406: /* try for squarish distribution */
407: m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
408: if (!m) m = 1;
409: while (m > 0) {
410: p = size/(m*n);
411: if (m*n*p == size) break;
412: m--;
413: }
414: if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n);
415: if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
416: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
417: /* try for squarish distribution */
418: n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
419: if (!n) n = 1;
420: while (n > 0) {
421: p = size/(m*n);
422: if (m*n*p == size) break;
423: n--;
424: }
425: if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n);
426: if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;}
427: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
428: /* try for squarish distribution */
429: n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));
430: if (!n) n = 1;
431: while (n > 0) {
432: pm = size/n;
433: if (n*pm == size) break;
434: n--;
435: }
436: if (!n) n = 1;
437: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
438: if (!m) m = 1;
439: while (m > 0) {
440: p = size/(m*n);
441: if (m*n*p == size) break;
442: m--;
443: }
444: if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;}
445: } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
447: if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Could not find good partition");
448: if (M < m) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
449: if (N < n) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
450: if (P < p) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Partition in z direction is too fine! %D %D",P,p);
452: stag->nRanks[0] = m;
453: stag->nRanks[1] = n;
454: stag->nRanks[2] = p;
455: return(0);
456: }
458: /* Determine ranks of neighbors, using DMDA's convention
460: n24 n25 n26
461: n21 n22 n23
462: n18 n19 n20 (Front, bigger z)
464: n15 n16 n17
465: n12 n14 ^ y
466: n9 n10 n11 |
467: +--> x
468: n6 n7 n8
469: n3 n4 n5
470: n0 n1 n2 (Back, smaller z) */
471: static PetscErrorCode DMStagSetUpBuildNeighbors_3d(DM dm)
472: {
473: PetscErrorCode ierr;
474: DM_Stag * const stag = (DM_Stag*)dm->data;
475: PetscInt d,i;
476: PetscBool per[3],first[3],last[3];
477: PetscInt neighborRank[27][3],r[3],n[3];
478: const PetscInt dim = 3;
481: for (d=0; d<dim; ++d) if (stag->boundaryType[d] != DM_BOUNDARY_NONE && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->boundaryType[d] != DM_BOUNDARY_GHOSTED) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Neighbor determination not implemented for %s",DMBoundaryTypes[stag->boundaryType[d]]);
483: /* Assemble some convenience variables */
484: for (d=0; d<dim; ++d) {
485: per[d] = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
486: first[d] = stag->firstRank[d];
487: last[d] = stag->lastRank[d];
488: r[d] = stag->rank[d];
489: n[d] = stag->nRanks[d];
490: }
492: /* First, compute the position in the rank grid for all neighbors */
494: neighborRank[0][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left down back */
495: neighborRank[0][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
496: neighborRank[0][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
498: neighborRank[1][0] = r[0] ; /* down back */
499: neighborRank[1][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
500: neighborRank[1][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
502: neighborRank[2][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right down back */
503: neighborRank[2][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
504: neighborRank[2][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
506: neighborRank[3][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left back */
507: neighborRank[3][1] = r[1] ;
508: neighborRank[3][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
510: neighborRank[4][0] = r[0] ; /* back */
511: neighborRank[4][1] = r[1] ;
512: neighborRank[4][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
514: neighborRank[5][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right back */
515: neighborRank[5][1] = r[1] ;
516: neighborRank[5][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
518: neighborRank[6][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left up back */
519: neighborRank[6][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
520: neighborRank[6][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
522: neighborRank[7][0] = r[0] ; /* up back */
523: neighborRank[7][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
524: neighborRank[7][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
526: neighborRank[8][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right up back */
527: neighborRank[8][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
528: neighborRank[8][2] = first[2] ? (per[2] ? n[2]-1 : -1) : r[2] - 1;
530: neighborRank[9][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left down */
531: neighborRank[9][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
532: neighborRank[9][2] = r[2] ;
534: neighborRank[10][0] = r[0] ; /* down */
535: neighborRank[10][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
536: neighborRank[10][2] = r[2] ;
538: neighborRank[11][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right down */
539: neighborRank[11][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
540: neighborRank[11][2] = r[2] ;
542: neighborRank[12][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left */
543: neighborRank[12][1] = r[1] ;
544: neighborRank[12][2] = r[2] ;
546: neighborRank[13][0] = r[0] ; /* */
547: neighborRank[13][1] = r[1] ;
548: neighborRank[13][2] = r[2] ;
550: neighborRank[14][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right */
551: neighborRank[14][1] = r[1] ;
552: neighborRank[14][2] = r[2] ;
554: neighborRank[15][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left up */
555: neighborRank[15][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
556: neighborRank[15][2] = r[2] ;
558: neighborRank[16][0] = r[0] ; /* up */
559: neighborRank[16][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
560: neighborRank[16][2] = r[2] ;
562: neighborRank[17][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right up */
563: neighborRank[17][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
564: neighborRank[17][2] = r[2] ;
566: neighborRank[18][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left down front */
567: neighborRank[18][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
568: neighborRank[18][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
570: neighborRank[19][0] = r[0] ; /* down front */
571: neighborRank[19][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
572: neighborRank[19][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
574: neighborRank[20][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right down front */
575: neighborRank[20][1] = first[1] ? (per[1] ? n[1]-1 : -1) : r[1] - 1;
576: neighborRank[20][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
578: neighborRank[21][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left front */
579: neighborRank[21][1] = r[1] ;
580: neighborRank[21][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
582: neighborRank[22][0] = r[0] ; /* front */
583: neighborRank[22][1] = r[1] ;
584: neighborRank[22][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
586: neighborRank[23][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right front */
587: neighborRank[23][1] = r[1] ;
588: neighborRank[23][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
590: neighborRank[24][0] = first[0] ? (per[0] ? n[0]-1 : -1) : r[0] - 1; /* left up front */
591: neighborRank[24][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
592: neighborRank[24][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
594: neighborRank[25][0] = r[0] ; /* up front */
595: neighborRank[25][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
596: neighborRank[25][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
598: neighborRank[26][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right up front */
599: neighborRank[26][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
600: neighborRank[26][2] = last[2] ? (per[2] ? 0 : -1) : r[2] + 1;
602: /* Then, compute the rank of each in the linear ordering */
603: PetscMalloc1(27,&stag->neighbors);
604: for (i=0; i<27; ++i){
605: if (neighborRank[i][0] >= 0 && neighborRank[i][1] >=0 && neighborRank[i][2] >=0) {
606: stag->neighbors[i] = neighborRank[i][0] + n[0]*neighborRank[i][1] + n[0]*n[1]*neighborRank[i][2];
607: } else {
608: stag->neighbors[i] = -1;
609: }
610: }
612: return(0);
613: }
615: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_3d(DM dm,PetscInt **pGlobalOffsets)
616: {
617: PetscErrorCode ierr;
618: const DM_Stag * const stag = (DM_Stag*)dm->data;
619: PetscInt *globalOffsets;
620: PetscInt i,j,k,d,entriesPerEdge,entriesPerFace,count;
621: PetscMPIInt size;
622: PetscBool extra[3];
625: for (d=0; d<3; ++d) extra[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC); /* Extra points in global rep */
626: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
627: entriesPerEdge = stag->dof[0] + stag->dof[1];
628: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
629: PetscMalloc1(size,pGlobalOffsets);
630: globalOffsets = *pGlobalOffsets;
631: globalOffsets[0] = 0;
632: count = 1; /* note the count is offset by 1 here. We add the size of the previous rank */
633: for (k=0; k<stag->nRanks[2]-1; ++k){
634: const PetscInt nnk = stag->l[2][k];
635: for (j=0; j<stag->nRanks[1]-1; ++j) {
636: const PetscInt nnj = stag->l[1][j];
637: for (i=0; i<stag->nRanks[0]-1; ++i) {
638: const PetscInt nni = stag->l[0][i];
639: /* Interior : No right/top/front boundaries */
640: globalOffsets[count] = globalOffsets[count-1] + nni*nnj*nnk*stag->entriesPerElement;
641: ++count;
642: }
643: {
644: /* Right boundary - extra faces */
645: /* i = stag->nRanks[0]-1; */
646: const PetscInt nni = stag->l[0][i];
647: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
648: + (extra[0] ? nnj*nnk*entriesPerFace : 0);
649: ++count;
650: }
651: }
652: {
653: /* j = stag->nRanks[1]-1; */
654: const PetscInt nnj = stag->l[1][j];
655: for (i=0; i<stag->nRanks[0]-1; ++i) {
656: const PetscInt nni = stag->l[0][i];
657: /* Up boundary - extra faces */
658: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
659: + (extra[1] ? nni*nnk*entriesPerFace : 0);
660: ++count;
661: }
662: {
663: /* i = stag->nRanks[0]-1; */
664: const PetscInt nni = stag->l[0][i];
665: /* Up right boundary - 2x extra faces and extra edges */
666: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
667: + (extra[0] ? nnj*nnk*entriesPerFace : 0)
668: + (extra[1] ? nni*nnk*entriesPerFace : 0)
669: + (extra[0] && extra[1] ? nnk*entriesPerEdge : 0);
670: ++count;
671: }
672: }
673: }
674: {
675: /* k = stag->nRanks[2]-1; */
676: const PetscInt nnk = stag->l[2][k];
677: for (j=0; j<stag->nRanks[1]-1; ++j) {
678: const PetscInt nnj = stag->l[1][j];
679: for (i=0; i<stag->nRanks[0]-1; ++i) {
680: const PetscInt nni = stag->l[0][i];
681: /* Front boundary - extra faces */
682: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
683: + (extra[2] ? nni*nnj*entriesPerFace : 0);
684: ++count;
685: }
686: {
687: /* i = stag->nRanks[0]-1; */
688: const PetscInt nni = stag->l[0][i];
689: /* Front right boundary - 2x extra faces and extra edges */
690: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
691: + (extra[0] ? nnk*nnj*entriesPerFace : 0)
692: + (extra[2] ? nni*nnj*entriesPerFace : 0)
693: + (extra[0] && extra[2] ? nnj*entriesPerEdge : 0);
694: ++count;
695: }
696: }
697: {
698: /* j = stag->nRanks[1]-1; */
699: const PetscInt nnj = stag->l[1][j];
700: for (i=0; i<stag->nRanks[0]-1; ++i) {
701: const PetscInt nni = stag->l[0][i];
702: /* Front up boundary - 2x extra faces and extra edges */
703: globalOffsets[count] = globalOffsets[count-1] + nnj*nni*nnk*stag->entriesPerElement
704: + (extra[1] ? nnk*nni*entriesPerFace : 0)
705: + (extra[2] ? nnj*nni*entriesPerFace : 0)
706: + (extra[1] && extra[2] ? nni*entriesPerEdge : 0);
707: ++count;
708: }
709: /* Don't need to compute entries in last element */
710: }
711: }
713: return(0);
714: }
716: /* A helper function to reduce code duplication as we loop over various ranges.
717: i,j,k refer to element numbers on the rank where an element lives in the global
718: representation (without ghosts) to be offset by a global offset per rank.
719: ig,jg,kg refer to element numbers in the local (this rank) ghosted numbering.
720: Note that this function could easily be converted to a macro (it does not compute
721: anything except loop indices and the values of idxGlobal and idxLocal). */
722: static PetscErrorCode DMStagSetUpBuildScatterPopulateIdx_3d(DM_Stag *stag,PetscInt *count,PetscInt *idxLocal,PetscInt *idxGlobal,PetscInt entriesPerEdge, PetscInt entriesPerFace,PetscInt eprNeighbor,PetscInt eplNeighbor,PetscInt eprGhost,PetscInt eplGhost,PetscInt epFaceRow,PetscInt globalOffset,PetscInt startx,PetscInt starty,PetscInt startz,PetscInt startGhostx,PetscInt startGhosty,PetscInt startGhostz,PetscInt endGhostx,PetscInt endGhosty,PetscInt endGhostz, PetscBool extrax,PetscBool extray,PetscBool extraz)
723: {
724: PetscInt ig,jg,kg,d,c;
727: c = *count;
728: for(kg = startGhostz; kg < endGhostz; ++kg) {
729: const PetscInt k = kg - startGhostz + startz;
730: for (jg = startGhosty; jg < endGhosty; ++jg) {
731: const PetscInt j = jg - startGhosty + starty;
732: for (ig = startGhostx; ig<endGhostx; ++ig) {
733: const PetscInt i = ig - startGhostx + startx;
734: for (d=0; d<stag->entriesPerElement; ++d,++c) {
735: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
736: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + d;
737: }
738: }
739: if (extrax) {
740: PetscInt dLocal;
741: const PetscInt i = endGhostx - startGhostx + startx;
742: ig = endGhostx;
743: for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
744: idxGlobal[c] = globalOffset + k *eplNeighbor + j*eprNeighbor + i *stag->entriesPerElement + d;
745: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
746: }
747: dLocal += stag->dof[1]; /* Skip back down edge */
748: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++c) { /* Back left edge */
749: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
750: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
751: }
752: dLocal += stag->dof[2]; /* Skip back face */
753: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++c) { /* Down left edge */
754: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
755: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
756: }
757: dLocal += stag->dof[2]; /* Skip down face */
758: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++c) { /* Left face */
759: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *stag->entriesPerElement + d;
760: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
761: }
762: /* Skip element */
763: }
764: }
765: if (extray) {
766: const PetscInt j = endGhosty - startGhosty + starty;
767: jg = endGhosty;
768: for (ig = startGhostx; ig<endGhostx; ++ig) {
769: const PetscInt i = ig - startGhostx + startx;
770: /* Vertex and Back down edge */
771: PetscInt dLocal;
772: for (d=0,dLocal=0; d<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++c) { /* Vertex */
773: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace + d; /* Note face increment in x */
774: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
775: }
776: /* Skip back left edge and back face */
777: dLocal += stag->dof[1] + stag->dof[2];
778: /* Down face and down left edge */
779: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++c) { /* Back left edge */
780: idxGlobal[c] = globalOffset + k *eplNeighbor + j*eprNeighbor + i *entriesPerFace + d; /* Note face increment in x */
781: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
782: }
783: /* Skip left face and element */
784: }
785: if (extrax) {
786: PetscInt dLocal;
787: const PetscInt i = endGhostx - startGhostx + startx;
788: ig = endGhostx;
789: for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
790: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace + d; /* Note face increment in x */
791: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
792: }
793: dLocal += 2*stag->dof[1]+stag->dof[2]; /* Skip back down edge, back face, and back left edge */
794: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++c) { /* Down left edge */
795: idxGlobal[c] = globalOffset + k *eplNeighbor + j *eprNeighbor + i *entriesPerFace + d; /* Note face increment in x */
796: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
797: }
798: /* Skip remaining entries */
799: }
800: }
801: }
802: if (extraz) {
803: const PetscInt k = endGhostz - startGhostz + startz;
804: kg = endGhostz;
805: for (jg = startGhosty; jg < endGhosty; ++jg) {
806: const PetscInt j = jg - startGhosty + starty;
807: for (ig = startGhostx; ig<endGhostx; ++ig) {
808: const PetscInt i = ig - startGhostx + startx;
809: for (d=0; d<entriesPerFace; ++d, ++c) {
810: idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerFace + d; /* Note face-based x and y increments */
811: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + d;
812: }
813: }
814: if (extrax) {
815: PetscInt dLocal;
816: const PetscInt i = endGhostx - startGhostx + startx;
817: ig = endGhostx;
818: for (d=0,dLocal=0; d<stag->dof[0]; ++d, ++dLocal, ++c) { /* Vertex */
819: idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerFace + d; /* Note face-based x and y increments */
820: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
821: }
822: dLocal += stag->dof[1]; /* Skip back down edge */
823: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++c) { /* Back left edge */
824: idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i*entriesPerFace + d; /* Note face-based x and y increments */
825: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + dLocal;
826: }
827: /* Skip the rest */
828: }
829: }
830: if (extray) {
831: const PetscInt j = endGhosty - startGhosty + starty;
832: jg = endGhosty;
833: for (ig = startGhostx; ig<endGhostx; ++ig) {
834: const PetscInt i = ig - startGhostx + startx;
835: for (d=0; d<entriesPerEdge; ++d,++c) {
836: idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerEdge + d; /* Note face-based y increment and edge-based x increment */
837: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + d;
838: }
839: }
840: if (extrax) {
841: const PetscInt i = endGhostx - startGhostx + startx;
842: ig = endGhostx;
843: for (d=0; d<stag->dof[0]; ++d, ++c) { /* Vertex (only) */
844: idxGlobal[c] = globalOffset + k*eplNeighbor + j *epFaceRow + i *entriesPerEdge + d; /* Note face-based y increment and edge-based x increment */
845: idxLocal[c] = kg*eplGhost + jg*eprGhost + ig*stag->entriesPerElement + d;
846: }
847: }
848: }
849: }
850: *count = c;
851: return(0);
852: }
854: static PetscErrorCode DMStagSetUpBuildScatter_3d(DM dm,const PetscInt *globalOffsets)
855: {
857: DM_Stag * const stag = (DM_Stag*)dm->data;
858: PetscInt d,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerCorner,entriesPerEdge,entriesPerFace,entriesToTransferTotal,count,eprGhost,eplGhost;
859: PetscInt *idxLocal,*idxGlobal;
860: PetscMPIInt rank;
861: PetscInt nNeighbors[27][3];
862: PetscBool star,nextToDummyEnd[3],dummyStart[3],dummyEnd[3];
865: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
866: if (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth || stag->n[2] < stag->stencilWidth) {
867: SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMStag 3d setup does not support local sizes (%D x %D x %D) smaller than the elementwise stencil width (%D)",stag->n[0],stag->n[1],stag->n[2],stag->stencilWidth);
868: }
870: /* Check stencil type */
871: if (stag->stencilType != DMSTAG_STENCIL_NONE && stag->stencilType != DMSTAG_STENCIL_BOX && stag->stencilType != DMSTAG_STENCIL_STAR) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported stencil type %s",DMStagStencilTypes[stag->stencilType]);
872: if (stag->stencilType == DMSTAG_STENCIL_NONE && stag->stencilWidth != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMStag 3d setup requires stencil width 0 with stencil type none");
873: star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);
875: /* Compute numbers of elements on each neighbor */
876: {
877: PetscInt i;
878: for (i=0; i<27; ++i) {
879: const PetscInt neighborRank = stag->neighbors[i];
880: if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
881: nNeighbors[i][0] = stag->l[0][neighborRank % stag->nRanks[0]];
882: nNeighbors[i][1] = stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
883: nNeighbors[i][2] = stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
884: } /* else leave uninitialized - error if accessed */
885: }
886: }
888: /* These offsets should always be non-negative, and describe how many
889: ghost elements exist at each boundary. These are not always equal to the stencil width,
890: because we may have different numbers of ghost elements at the boundaries. In particular,
891: in the non-periodic casewe always have at least one ghost (dummy) element at the right/top/front. */
892: for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
893: for (d=0; d<3; ++d) ghostOffsetEnd[d] = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
895: /* Determine whether the ghost region includes dummies or not. This is currently
896: equivalent to having a non-periodic boundary. If not, then
897: ghostOffset{Start,End}[d] elements correspond to elements on the neighbor.
898: If true, then
899: - at the start, there are ghostOffsetStart[d] ghost elements
900: - at the end, there is a layer of extra "physical" points inside a layer of
901: ghostOffsetEnd[d] ghost elements
902: Note that this computation should be updated if any boundary types besides
903: NONE, GHOSTED, and PERIODIC are supported. */
904: for (d=0; d<3; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
905: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
907: /* Convenience variables */
908: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
909: entriesPerEdge = stag->dof[0] + stag->dof[1];
910: entriesPerCorner = stag->dof[0];
911: for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
912: eprGhost = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row */
913: eplGhost = stag->nGhost[1]*eprGhost; /* epl = entries per (element) layer */
915: /* Compute the number of local entries which correspond to any global entry */
916: {
917: PetscInt nNonDummyGhost[3];
918: for (d=0; d<3; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
919: if (star) {
920: entriesToTransferTotal = (
921: nNonDummyGhost[0] * stag->n[1] * stag->n[2]
922: + stag->n[0] * nNonDummyGhost[1] * stag->n[2]
923: + stag->n[0] * stag->n[1] * nNonDummyGhost[2]
924: - 2 * (stag->n[0] * stag->n[1] * stag->n[2])
925: ) * stag->entriesPerElement
926: + (dummyEnd[0] ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace : 0)
927: + (dummyEnd[1] ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace : 0)
928: + (dummyEnd[2] ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace : 0)
929: + (dummyEnd[0] && dummyEnd[1] ? nNonDummyGhost[2] * entriesPerEdge : 0)
930: + (dummyEnd[2] && dummyEnd[0] ? nNonDummyGhost[1] * entriesPerEdge : 0)
931: + (dummyEnd[1] && dummyEnd[2] ? nNonDummyGhost[0] * entriesPerEdge : 0)
932: + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ? entriesPerCorner : 0);
933: } else {
934: entriesToTransferTotal = nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
935: + (dummyEnd[0] ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace : 0)
936: + (dummyEnd[1] ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace : 0)
937: + (dummyEnd[2] ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace : 0)
938: + (dummyEnd[0] && dummyEnd[1] ? nNonDummyGhost[2] * entriesPerEdge : 0)
939: + (dummyEnd[2] && dummyEnd[0] ? nNonDummyGhost[1] * entriesPerEdge : 0)
940: + (dummyEnd[1] && dummyEnd[2] ? nNonDummyGhost[0] * entriesPerEdge : 0)
941: + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ? entriesPerCorner : 0);
942: }
943: }
945: /* Allocate arrays to populate */
946: PetscMalloc1(entriesToTransferTotal,&idxLocal);
947: PetscMalloc1(entriesToTransferTotal,&idxGlobal);
949: /* Counts into idxLocal/idxGlobal */
950: count = 0;
952: /* Loop over each of the 27 neighbor, populating idxLocal and idxGlobal. These
953: cases are principally distinguished by
955: 1. The loop bounds (i/ighost,j/jghost,k/kghost)
956: 2. the strides used in the loop body, which depend on whether the current
957: rank and/or its neighbor are a non-periodic right/top/front boundary, which has additional
958: points in the global representation.
959: - If the neighboring rank is a right/top/front boundary,
960: then eprNeighbor (entries per element row on the neighbor) and/or eplNeighbor (entries per element layer on the neighbor)
961: are different.
962: - If this rank is a non-periodic right/top/front boundary (checking entries of stag->lastRank),
963: there is an extra loop over 1-3 boundary faces)
965: Here, we do not include "dummy" dof (points in the local representation which
966: do not correspond to any global dof). This, and the fact that we iterate over points in terms of
967: increasing global ordering, are the main two differences from the construction of
968: the Local-to-global map, which iterates over points in local order, and does include dummy points. */
970: /* LEFT DOWN BACK */
971: if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
972: const PetscInt neighbor = 0;
973: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
974: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
975: const PetscInt epFaceRow = -1;
976: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
977: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
978: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
979: const PetscInt startGhost0 = 0;
980: const PetscInt startGhost1 = 0;
981: const PetscInt startGhost2 = 0;
982: const PetscInt endGhost0 = ghostOffsetStart[0];
983: const PetscInt endGhost1 = ghostOffsetStart[1];
984: const PetscInt endGhost2 = ghostOffsetStart[2];
985: const PetscBool extra0 = PETSC_FALSE;
986: const PetscBool extra1 = PETSC_FALSE;
987: const PetscBool extra2 = PETSC_FALSE;
988: 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);
989: }
991: /* DOWN BACK */
992: if (!star && !dummyStart[1] && !dummyStart[2]) {
993: const PetscInt neighbor = 1;
994: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
995: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
996: const PetscInt epFaceRow = -1;
997: const PetscInt start0 = 0;
998: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
999: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1000: const PetscInt startGhost0 = ghostOffsetStart[0];
1001: const PetscInt startGhost1 = 0;
1002: const PetscInt startGhost2 = 0;
1003: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1004: const PetscInt endGhost1 = ghostOffsetStart[1];
1005: const PetscInt endGhost2 = ghostOffsetStart[2];
1006: const PetscBool extra0 = dummyEnd[0];
1007: const PetscBool extra1 = PETSC_FALSE;
1008: const PetscBool extra2 = PETSC_FALSE;
1009: 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);
1010: }
1012: /* RIGHT DOWN BACK */
1013: if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyStart[2]) {
1014: const PetscInt neighbor = 2;
1015: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1016: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1017: const PetscInt epFaceRow = -1;
1018: const PetscInt start0 = 0;
1019: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1020: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1021: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1022: const PetscInt startGhost1 = 0;
1023: const PetscInt startGhost2 = 0;
1024: const PetscInt endGhost0 = stag->nGhost[0];
1025: const PetscInt endGhost1 = ghostOffsetStart[1];
1026: const PetscInt endGhost2 = ghostOffsetStart[2];
1027: const PetscBool extra0 = PETSC_FALSE;
1028: const PetscBool extra1 = PETSC_FALSE;
1029: const PetscBool extra2 = PETSC_FALSE;
1030: 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);
1031: }
1033: /* LEFT BACK */
1034: if (!star && !dummyStart[0] && !dummyStart[2]) {
1035: const PetscInt neighbor = 3;
1036: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1037: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* May be a top boundary */
1038: const PetscInt epFaceRow = -1;
1039: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1040: const PetscInt start1 = 0;
1041: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1042: const PetscInt startGhost0 = 0;
1043: const PetscInt startGhost1 = ghostOffsetStart[1];
1044: const PetscInt startGhost2 = 0;
1045: const PetscInt endGhost0 = ghostOffsetStart[0];
1046: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1047: const PetscInt endGhost2 = ghostOffsetStart[2];
1048: const PetscBool extra0 = PETSC_FALSE;
1049: const PetscBool extra1 = dummyEnd[1];
1050: const PetscBool extra2 = PETSC_FALSE;
1051: 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);
1052: }
1054: /* BACK */
1055: if (!dummyStart[2]) {
1056: const PetscInt neighbor = 4;
1057: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1058: 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 */
1059: const PetscInt epFaceRow = -1;
1060: const PetscInt start0 = 0;
1061: const PetscInt start1 = 0;
1062: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1063: const PetscInt startGhost0 = ghostOffsetStart[0];
1064: const PetscInt startGhost1 = ghostOffsetStart[1];
1065: const PetscInt startGhost2 = 0;
1066: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1067: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1068: const PetscInt endGhost2 = ghostOffsetStart[2];
1069: const PetscBool extra0 = dummyEnd[0];
1070: const PetscBool extra1 = dummyEnd[1];
1071: const PetscBool extra2 = PETSC_FALSE;
1072: 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);
1073: }
1075: /* RIGHT BACK */
1076: if (!star && !dummyEnd[0] && !dummyStart[2]) {
1077: const PetscInt neighbor = 5;
1078: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1079: 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 */
1080: const PetscInt epFaceRow = -1;
1081: const PetscInt start0 = 0;
1082: const PetscInt start1 = 0;
1083: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1084: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1085: const PetscInt startGhost1 = ghostOffsetStart[1];
1086: const PetscInt startGhost2 = 0;
1087: const PetscInt endGhost0 = stag->nGhost[0];
1088: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1089: const PetscInt endGhost2 = ghostOffsetStart[2];
1090: const PetscBool extra0 = PETSC_FALSE;
1091: const PetscBool extra1 = dummyEnd[1];
1092: const PetscBool extra2 = PETSC_FALSE;
1093: 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);
1094: }
1096: /* LEFT UP BACK */
1097: if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyStart[2]) {
1098: const PetscInt neighbor = 6;
1099: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1100: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1101: const PetscInt epFaceRow = -1;
1102: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1103: const PetscInt start1 = 0;
1104: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1105: const PetscInt startGhost0 = 0;
1106: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1107: const PetscInt startGhost2 = 0;
1108: const PetscInt endGhost0 = ghostOffsetStart[0];
1109: const PetscInt endGhost1 = stag->nGhost[1];
1110: const PetscInt endGhost2 = ghostOffsetStart[2];
1111: const PetscBool extra0 = PETSC_FALSE;
1112: const PetscBool extra1 = PETSC_FALSE;
1113: const PetscBool extra2 = PETSC_FALSE;
1114: 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);
1115: }
1117: /* UP BACK */
1118: if (!star && !dummyEnd[1] && !dummyStart[2]) {
1119: const PetscInt neighbor = 7;
1120: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1121: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1122: const PetscInt epFaceRow = -1;
1123: const PetscInt start0 = 0;
1124: const PetscInt start1 = 0;
1125: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1126: const PetscInt startGhost0 = ghostOffsetStart[0];
1127: const PetscInt startGhost1 = stag->nGhost[1]-ghostOffsetEnd[1];
1128: const PetscInt startGhost2 = 0;
1129: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1130: const PetscInt endGhost1 = stag->nGhost[1];
1131: const PetscInt endGhost2 = ghostOffsetStart[2];
1132: const PetscBool extra0 = dummyEnd[0];
1133: const PetscBool extra1 = PETSC_FALSE;
1134: const PetscBool extra2 = PETSC_FALSE;
1135: 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);
1136: }
1138: /* RIGHT UP BACK */
1139: if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyStart[2]) {
1140: const PetscInt neighbor = 8;
1141: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1142: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1143: const PetscInt epFaceRow = -1;
1144: const PetscInt start0 = 0;
1145: const PetscInt start1 = 0;
1146: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1147: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1148: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1149: const PetscInt startGhost2 = 0;
1150: const PetscInt endGhost0 = stag->nGhost[0];
1151: const PetscInt endGhost1 = stag->nGhost[1];
1152: const PetscInt endGhost2 = ghostOffsetStart[2];
1153: const PetscBool extra0 = PETSC_FALSE;
1154: const PetscBool extra1 = PETSC_FALSE;
1155: const PetscBool extra2 = PETSC_FALSE;
1156: 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);
1157: }
1159: /* LEFT DOWN */
1160: if (!star && !dummyStart[0] && !dummyStart[1]) {
1161: const PetscInt neighbor = 9;
1162: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1163: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1164: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1165: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1166: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1167: const PetscInt start2 = 0;
1168: const PetscInt startGhost0 = 0;
1169: const PetscInt startGhost1 = 0;
1170: const PetscInt startGhost2 = ghostOffsetStart[2];
1171: const PetscInt endGhost0 = ghostOffsetStart[0];
1172: const PetscInt endGhost1 = ghostOffsetStart[1];
1173: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1174: const PetscBool extra0 = PETSC_FALSE;
1175: const PetscBool extra1 = PETSC_FALSE;
1176: const PetscBool extra2 = dummyEnd[2];
1177: 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);
1178: }
1180: /* DOWN */
1181: if (!dummyStart[1]) {
1182: const PetscInt neighbor = 10;
1183: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1184: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1185: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1186: const PetscInt start0 = 0;
1187: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1188: const PetscInt start2 = 0;
1189: const PetscInt startGhost0 = ghostOffsetStart[0];
1190: const PetscInt startGhost1 = 0;
1191: const PetscInt startGhost2 = ghostOffsetStart[2];
1192: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1193: const PetscInt endGhost1 = ghostOffsetStart[1];
1194: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1195: const PetscBool extra0 = dummyEnd[0];
1196: const PetscBool extra1 = PETSC_FALSE;
1197: const PetscBool extra2 = dummyEnd[2];
1198: 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);
1199: }
1201: /* RIGHT DOWN */
1202: if (!star && !dummyEnd[0] && !dummyStart[1]) {
1203: const PetscInt neighbor = 11;
1204: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1205: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1206: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1207: const PetscInt start0 = 0;
1208: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1209: const PetscInt start2 = 0;
1210: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1211: const PetscInt startGhost1 = 0;
1212: const PetscInt startGhost2 = ghostOffsetStart[2];
1213: const PetscInt endGhost0 = stag->nGhost[0];
1214: const PetscInt endGhost1 = ghostOffsetStart[1];
1215: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1216: const PetscBool extra0 = PETSC_FALSE;
1217: const PetscBool extra1 = PETSC_FALSE;
1218: const PetscBool extra2 = dummyEnd[2];
1219: 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);
1220: }
1222: /* LEFT */
1223: if (!dummyStart[0]) {
1224: const PetscInt neighbor = 12;
1225: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1226: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1227: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
1228: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1229: const PetscInt start1 = 0;
1230: const PetscInt start2 = 0;
1231: const PetscInt startGhost0 = 0;
1232: const PetscInt startGhost1 = ghostOffsetStart[1];
1233: const PetscInt startGhost2 = ghostOffsetStart[2];
1234: const PetscInt endGhost0 = ghostOffsetStart[0];
1235: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1236: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1237: const PetscBool extra0 = PETSC_FALSE;
1238: const PetscBool extra1 = dummyEnd[1];
1239: const PetscBool extra2 = dummyEnd[2];
1240: 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);
1241: }
1243: /* (HERE) */
1244: {
1245: const PetscInt neighbor = 13;
1246: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1247: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1248: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
1249: const PetscInt start0 = 0;
1250: const PetscInt start1 = 0;
1251: const PetscInt start2 = 0;
1252: const PetscInt startGhost0 = ghostOffsetStart[0];
1253: const PetscInt startGhost1 = ghostOffsetStart[1];
1254: const PetscInt startGhost2 = ghostOffsetStart[2];
1255: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1256: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1257: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1258: const PetscBool extra0 = dummyEnd[0];
1259: const PetscBool extra1 = dummyEnd[1];
1260: const PetscBool extra2 = dummyEnd[2];
1261: 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);
1262: }
1264: /* RIGHT */
1265: if (!dummyEnd[0]) {
1266: const PetscInt neighbor = 14;
1267: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1268: 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 */
1269: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1270: const PetscInt start0 = 0;
1271: const PetscInt start1 = 0;
1272: const PetscInt start2 = 0;
1273: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1274: const PetscInt startGhost1 = ghostOffsetStart[1];
1275: const PetscInt startGhost2 = ghostOffsetStart[2];
1276: const PetscInt endGhost0 = stag->nGhost[0];
1277: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1278: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1279: const PetscBool extra0 = PETSC_FALSE;
1280: const PetscBool extra1 = dummyEnd[1];
1281: const PetscBool extra2 = dummyEnd[2];
1282: 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);
1283: }
1285: /* LEFT UP */
1286: if (!star && !dummyStart[0] && !dummyEnd[1]) {
1287: const PetscInt neighbor = 15;
1288: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1289: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1290: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
1291: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1292: const PetscInt start1 = 0;
1293: const PetscInt start2 = 0;
1294: const PetscInt startGhost0 = 0;
1295: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1296: const PetscInt startGhost2 = ghostOffsetStart[2];
1297: const PetscInt endGhost0 = ghostOffsetStart[0];
1298: const PetscInt endGhost1 = stag->nGhost[1];
1299: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1300: const PetscBool extra0 = PETSC_FALSE;
1301: const PetscBool extra1 = PETSC_FALSE;
1302: const PetscBool extra2 = dummyEnd[2];
1303: 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);
1304: }
1306: /* UP */
1307: if (!dummyEnd[1]) {
1308: const PetscInt neighbor = 16;
1309: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1310: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1311: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1312: const PetscInt start0 = 0;
1313: const PetscInt start1 = 0;
1314: const PetscInt start2 = 0;
1315: const PetscInt startGhost0 = ghostOffsetStart[0];
1316: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1317: const PetscInt startGhost2 = ghostOffsetStart[2];
1318: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1319: const PetscInt endGhost1 = stag->nGhost[1];
1320: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1321: const PetscBool extra0 = dummyEnd[0];
1322: const PetscBool extra1 = PETSC_FALSE;
1323: const PetscBool extra2 = dummyEnd[2];
1324: 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);
1325: }
1327: /* RIGHT UP */
1328: if (!star && !dummyEnd[0] && !dummyEnd[1]) {
1329: const PetscInt neighbor = 17;
1330: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1331: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1332: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1333: const PetscInt start0 = 0;
1334: const PetscInt start1 = 0;
1335: const PetscInt start2 = 0;
1336: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1337: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1338: const PetscInt startGhost2 = ghostOffsetStart[2];
1339: const PetscInt endGhost0 = stag->nGhost[0];
1340: const PetscInt endGhost1 = stag->nGhost[1];
1341: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1342: const PetscBool extra0 = PETSC_FALSE;
1343: const PetscBool extra1 = PETSC_FALSE;
1344: const PetscBool extra2 = dummyEnd[2];
1345: 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);
1346: }
1348: /* LEFT DOWN FRONT */
1349: if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1350: const PetscInt neighbor = 18;
1351: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1352: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1353: const PetscInt epFaceRow = -1;
1354: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1355: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1356: const PetscInt start2 = 0;
1357: const PetscInt startGhost0 = 0;
1358: const PetscInt startGhost1 = 0;
1359: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1360: const PetscInt endGhost0 = ghostOffsetStart[0];
1361: const PetscInt endGhost1 = ghostOffsetStart[1];
1362: const PetscInt endGhost2 = stag->nGhost[2];
1363: const PetscBool extra0 = PETSC_FALSE;
1364: const PetscBool extra1 = PETSC_FALSE;
1365: const PetscBool extra2 = PETSC_FALSE;
1366: 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);
1367: }
1369: /* DOWN FRONT */
1370: if (!star && !dummyStart[1] && !dummyEnd[2]) {
1371: const PetscInt neighbor = 19;
1372: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1373: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1374: const PetscInt epFaceRow = -1;
1375: const PetscInt start0 = 0;
1376: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1377: const PetscInt start2 = 0;
1378: const PetscInt startGhost0 = ghostOffsetStart[0];
1379: const PetscInt startGhost1 = 0;
1380: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1381: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1382: const PetscInt endGhost1 = ghostOffsetStart[1];
1383: const PetscInt endGhost2 = stag->nGhost[2];
1384: const PetscBool extra0 = dummyEnd[0];
1385: const PetscBool extra1 = PETSC_FALSE;
1386: const PetscBool extra2 = PETSC_FALSE;
1387: 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);
1388: }
1390: /* RIGHT DOWN FRONT */
1391: if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyEnd[2]) {
1392: const PetscInt neighbor = 20;
1393: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1394: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1395: const PetscInt epFaceRow = -1;
1396: const PetscInt start0 = 0;
1397: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1398: const PetscInt start2 = 0;
1399: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1400: const PetscInt startGhost1 = 0;
1401: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1402: const PetscInt endGhost0 = stag->nGhost[0];
1403: const PetscInt endGhost1 = ghostOffsetStart[1];
1404: const PetscInt endGhost2 = stag->nGhost[2];
1405: const PetscBool extra0 = PETSC_FALSE;
1406: const PetscBool extra1 = PETSC_FALSE;
1407: const PetscBool extra2 = PETSC_FALSE;
1408: 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);
1409: }
1411: /* LEFT FRONT */
1412: if (!star && !dummyStart[0] && !dummyEnd[2]) {
1413: const PetscInt neighbor = 21;
1414: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1415: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1416: const PetscInt epFaceRow = -1;
1417: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1418: const PetscInt start1 = 0;
1419: const PetscInt start2 = 0;
1420: const PetscInt startGhost0 = 0;
1421: const PetscInt startGhost1 = ghostOffsetStart[1];
1422: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1423: const PetscInt endGhost0 = ghostOffsetStart[0];
1424: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1425: const PetscInt endGhost2 = stag->nGhost[2];
1426: const PetscBool extra0 = PETSC_FALSE;
1427: const PetscBool extra1 = dummyEnd[1];
1428: const PetscBool extra2 = PETSC_FALSE;
1429: 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);
1430: }
1432: /* FRONT */
1433: if (!dummyEnd[2]) {
1434: const PetscInt neighbor = 22;
1435: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1436: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1437: const PetscInt epFaceRow = -1;
1438: const PetscInt start0 = 0;
1439: const PetscInt start1 = 0;
1440: const PetscInt start2 = 0;
1441: const PetscInt startGhost0 = ghostOffsetStart[0];
1442: const PetscInt startGhost1 = ghostOffsetStart[1];
1443: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1444: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1445: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1446: const PetscInt endGhost2 = stag->nGhost[2];
1447: const PetscBool extra0 = dummyEnd[0];
1448: const PetscBool extra1 = dummyEnd[1];
1449: const PetscBool extra2 = PETSC_FALSE;
1450: 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);
1451: }
1453: /* RIGHT FRONT */
1454: if (!star && !dummyEnd[0] && !dummyEnd[2]) {
1455: const PetscInt neighbor = 23;
1456: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1457: 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 */
1458: const PetscInt epFaceRow = -1;
1459: const PetscInt start0 = 0;
1460: const PetscInt start1 = 0;
1461: const PetscInt start2 = 0;
1462: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1463: const PetscInt startGhost1 = ghostOffsetStart[1];
1464: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1465: const PetscInt endGhost0 = stag->nGhost[0];
1466: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1467: const PetscInt endGhost2 = stag->nGhost[2];
1468: const PetscBool extra0 = PETSC_FALSE;
1469: const PetscBool extra1 = dummyEnd[1];
1470: const PetscBool extra2 = PETSC_FALSE;
1471: 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);
1472: }
1474: /* LEFT UP FRONT */
1475: if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyEnd[2]) {
1476: const PetscInt neighbor = 24;
1477: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1478: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1479: const PetscInt epFaceRow = -1;
1480: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1481: const PetscInt start1 = 0;
1482: const PetscInt start2 = 0;
1483: const PetscInt startGhost0 = 0;
1484: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1485: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1486: const PetscInt endGhost0 = ghostOffsetStart[0];
1487: const PetscInt endGhost1 = stag->nGhost[1];
1488: const PetscInt endGhost2 = stag->nGhost[2];
1489: const PetscBool extra0 = PETSC_FALSE;
1490: const PetscBool extra1 = PETSC_FALSE;
1491: const PetscBool extra2 = PETSC_FALSE;
1492: 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);
1493: }
1495: /* UP FRONT */
1496: if (!star && !dummyEnd[1] && !dummyEnd[2]) {
1497: const PetscInt neighbor = 25;
1498: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1499: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1500: const PetscInt epFaceRow = -1;
1501: const PetscInt start0 = 0;
1502: const PetscInt start1 = 0;
1503: const PetscInt start2 = 0;
1504: const PetscInt startGhost0 = ghostOffsetStart[0];
1505: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1506: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1507: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1508: const PetscInt endGhost1 = stag->nGhost[1];
1509: const PetscInt endGhost2 = stag->nGhost[2];
1510: const PetscBool extra0 = dummyEnd[0];
1511: const PetscBool extra1 = PETSC_FALSE;
1512: const PetscBool extra2 = PETSC_FALSE;
1513: 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);
1514: }
1516: /* RIGHT UP FRONT */
1517: if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyEnd[2]) {
1518: const PetscInt neighbor = 26;
1519: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1520: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1521: const PetscInt epFaceRow = -1;
1522: const PetscInt start0 = 0;
1523: const PetscInt start1 = 0;
1524: const PetscInt start2 = 0;
1525: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1526: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1527: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1528: const PetscInt endGhost0 = stag->nGhost[0];
1529: const PetscInt endGhost1 = stag->nGhost[1];
1530: const PetscInt endGhost2 = stag->nGhost[2];
1531: const PetscBool extra0 = PETSC_FALSE;
1532: const PetscBool extra1 = PETSC_FALSE;
1533: const PetscBool extra2 = PETSC_FALSE;
1534: 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);
1535: }
1537: #if defined(PETSC_USE_DEBUG)
1538: if (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);
1539: #endif
1541: /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
1542: {
1543: Vec local,global;
1544: IS isLocal,isGlobal;
1545: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
1546: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
1547: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
1548: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
1549: VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
1550: VecDestroy(&global);
1551: VecDestroy(&local);
1552: ISDestroy(&isLocal); /* frees idxLocal */
1553: ISDestroy(&isGlobal); /* free idxGlobal */
1554: }
1556: return(0);
1557: }
1559: /* Note: this function assumes that DMBoundary types of none, ghosted, and periodic are the only ones of interest.
1560: Adding support for others should be done very carefully. */
1561: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM dm,const PetscInt *globalOffsets)
1562: {
1563: PetscErrorCode ierr;
1564: const DM_Stag * const stag = (DM_Stag*)dm->data;
1565: PetscInt *idxGlobalAll;
1566: PetscInt d,count,ighost,jghost,kghost,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerFace,entriesPerEdge;
1567: PetscInt nNeighbors[27][3],eprNeighbor[27],eplNeighbor[27],globalOffsetNeighbor[27];
1568: PetscBool nextToDummyEnd[3],dummyStart[3],dummyEnd[3],star;
1572: /* Check stencil type */
1573: if (stag->stencilType != DMSTAG_STENCIL_NONE && stag->stencilType != DMSTAG_STENCIL_BOX && stag->stencilType != DMSTAG_STENCIL_STAR) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported stencil type %s",DMStagStencilTypes[stag->stencilType]);
1574: if (stag->stencilType == DMSTAG_STENCIL_NONE && stag->stencilWidth != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMStag 3d setup requires stencil width 0 with stencil type none");
1575: star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);
1577: /* Convenience variables */
1578: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
1579: entriesPerEdge = stag->dof[0] + stag->dof[1];
1580: for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
1582: /* Ghost offsets (may not be the ghost width, since we always have at least one ghost element on the right/top/front) */
1583: for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
1584: for (d=0; d<3; ++d) ghostOffsetEnd[d] = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
1586: /* Whether the ghost region includes dummies. Currently equivalent to being a non-periodic boundary. */
1587: for (d=0; d<3; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1588: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1590: /* Compute numbers of elements on each neighbor and offset*/
1591: {
1592: PetscInt i;
1593: for (i=0; i<27; ++i) {
1594: const PetscInt neighborRank = stag->neighbors[i];
1595: if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
1596: nNeighbors[i][0] = stag->l[0][neighborRank % stag->nRanks[0]];
1597: nNeighbors[i][1] = stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
1598: nNeighbors[i][2] = stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
1599: globalOffsetNeighbor[i] = globalOffsets[stag->neighbors[i]];
1600: } /* else leave uninitialized - error if accessed */
1601: }
1602: }
1604: /* Precompute elements per row and layer on neighbor (zero unused) */
1605: PetscMemzero(eprNeighbor,sizeof(eprNeighbor));
1606: PetscMemzero(eplNeighbor,sizeof(eplNeighbor));
1607: if (stag->neighbors[0] >= 0) {
1608: eprNeighbor[0] = stag->entriesPerElement * nNeighbors[0][0];
1609: eplNeighbor[0] = eprNeighbor[0] * nNeighbors[0][1];
1610: }
1611: if (stag->neighbors[1] >= 0) {
1612: eprNeighbor[1] = stag->entriesPerElement * nNeighbors[1][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1613: eplNeighbor[1] = eprNeighbor[1] * nNeighbors[1][1];
1614: }
1615: if (stag->neighbors[2] >= 0) {
1616: eprNeighbor[2] = stag->entriesPerElement * nNeighbors[2][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1617: eplNeighbor[2] = eprNeighbor[2] * nNeighbors[2][1];
1618: }
1619: if (stag->neighbors[3] >= 0) {
1620: eprNeighbor[3] = stag->entriesPerElement * nNeighbors[3][0];
1621: eplNeighbor[3] = eprNeighbor[3] * nNeighbors[3][1] + (dummyEnd[1] ? nNeighbors[3][0] * entriesPerFace : 0); /* May be a top boundary */
1622: }
1623: if (stag->neighbors[4] >= 0) {
1624: eprNeighbor[4] = stag->entriesPerElement * nNeighbors[4][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1625: 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 */
1626: }
1627: if (stag->neighbors[5] >= 0) {
1628: eprNeighbor[5] = stag->entriesPerElement * nNeighbors[5][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1629: 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 */
1630: }
1631: if (stag->neighbors[6] >= 0) {
1632: eprNeighbor[6] = stag->entriesPerElement * nNeighbors[6][0];
1633: eplNeighbor[6] = eprNeighbor[6] * nNeighbors[6][1] + (nextToDummyEnd[1] ? nNeighbors[6][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1634: }
1635: if (stag->neighbors[7] >= 0) {
1636: eprNeighbor[7] = stag->entriesPerElement * nNeighbors[7][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1637: eplNeighbor[7] = eprNeighbor[7] * nNeighbors[7][1] + (nextToDummyEnd[1] ? nNeighbors[7][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1638: }
1639: if (stag->neighbors[8] >= 0) {
1640: eprNeighbor[8] = stag->entriesPerElement * nNeighbors[8][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1641: eplNeighbor[8] = eprNeighbor[8] * nNeighbors[8][1] + (nextToDummyEnd[1] ? nNeighbors[8][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1642: }
1643: if (stag->neighbors[9] >= 0) {
1644: eprNeighbor[9] = stag->entriesPerElement * nNeighbors[9][0];
1645: eplNeighbor[9] = eprNeighbor[9] * nNeighbors[9][1];
1646: }
1647: if (stag->neighbors[10] >= 0) {
1648: eprNeighbor[10] = stag->entriesPerElement * nNeighbors[10][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1649: eplNeighbor[10] = eprNeighbor[10] * nNeighbors[10][1];
1650: }
1651: if (stag->neighbors[11] >= 0) {
1652: eprNeighbor[11] = stag->entriesPerElement * nNeighbors[11][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1653: eplNeighbor[11] = eprNeighbor[11] * nNeighbors[11][1];
1654: }
1655: if (stag->neighbors[12] >= 0) {
1656: eprNeighbor[12] = stag->entriesPerElement * nNeighbors[12][0];
1657: eplNeighbor[12] = eprNeighbor[12] * nNeighbors[12][1] + (dummyEnd[1] ? nNeighbors[12][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1658: }
1659: if (stag->neighbors[13] >= 0) {
1660: eprNeighbor[13] = stag->entriesPerElement * nNeighbors[13][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1661: eplNeighbor[13] = eprNeighbor[13] * nNeighbors[13][1] + (dummyEnd[1] ? nNeighbors[13][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1662: }
1663: if (stag->neighbors[14] >= 0) {
1664: eprNeighbor[14] = stag->entriesPerElement * nNeighbors[14][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1665: 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 */
1666: }
1667: if (stag->neighbors[15] >= 0) {
1668: eprNeighbor[15] = stag->entriesPerElement * nNeighbors[15][0];
1669: eplNeighbor[15] = eprNeighbor[15] * nNeighbors[15][1] + (nextToDummyEnd[1] ? nNeighbors[15][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1670: }
1671: if (stag->neighbors[16] >= 0) {
1672: eprNeighbor[16] = stag->entriesPerElement * nNeighbors[16][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1673: eplNeighbor[16] = eprNeighbor[16] * nNeighbors[16][1] + (nextToDummyEnd[1] ? nNeighbors[16][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1674: }
1675: if (stag->neighbors[17] >= 0) {
1676: eprNeighbor[17] = stag->entriesPerElement * nNeighbors[17][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1677: eplNeighbor[17] = eprNeighbor[17] * nNeighbors[17][1] + (nextToDummyEnd[1] ? nNeighbors[17][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1678: }
1679: if (stag->neighbors[18] >= 0) {
1680: eprNeighbor[18] = stag->entriesPerElement * nNeighbors[18][0];
1681: eplNeighbor[18] = eprNeighbor[18] * nNeighbors[18][1];
1682: }
1683: if (stag->neighbors[19] >= 0) {
1684: eprNeighbor[19] = stag->entriesPerElement * nNeighbors[19][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1685: eplNeighbor[19] = eprNeighbor[19] * nNeighbors[19][1];
1686: }
1687: if (stag->neighbors[20] >= 0) {
1688: eprNeighbor[20] = stag->entriesPerElement * nNeighbors[20][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1689: eplNeighbor[20] = eprNeighbor[20] * nNeighbors[20][1];
1690: }
1691: if (stag->neighbors[21] >= 0) {
1692: eprNeighbor[21] = stag->entriesPerElement * nNeighbors[21][0];
1693: eplNeighbor[21] = eprNeighbor[21] * nNeighbors[21][1] + (dummyEnd[1] ? nNeighbors[21][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1694: }
1695: if (stag->neighbors[22] >= 0) {
1696: eprNeighbor[22] = stag->entriesPerElement * nNeighbors[22][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1697: eplNeighbor[22] = eprNeighbor[22] * nNeighbors[22][1] + (dummyEnd[1] ? nNeighbors[22][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1698: }
1699: if (stag->neighbors[23] >= 0) {
1700: eprNeighbor[23] = stag->entriesPerElement * nNeighbors[23][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1701: 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 */
1702: }
1703: if (stag->neighbors[24] >= 0) {
1704: eprNeighbor[24] = stag->entriesPerElement * nNeighbors[24][0];
1705: eplNeighbor[24] = eprNeighbor[24] * nNeighbors[24][1] + (nextToDummyEnd[1] ? nNeighbors[24][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1706: }
1707: if (stag->neighbors[25] >= 0) {
1708: eprNeighbor[25] = stag->entriesPerElement * nNeighbors[25][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1709: eplNeighbor[25] = eprNeighbor[25] * nNeighbors[25][1] + (nextToDummyEnd[1] ? nNeighbors[25][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1710: }
1711: if (stag->neighbors[26] >= 0) {
1712: eprNeighbor[26] = stag->entriesPerElement * nNeighbors[26][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1713: eplNeighbor[26] = eprNeighbor[26] * nNeighbors[26][1] + (nextToDummyEnd[1] ? nNeighbors[26][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1714: }
1716: /* Populate idxGlobalAll */
1717: PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
1718: count = 0;
1720: /* Loop over layers 1/3 : Back */
1721: if (!dummyStart[2]) {
1722: for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
1723: 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) */
1725: /* Loop over rows 1/3: Down Back*/
1726: if (!star && !dummyStart[1]) {
1727: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1728: 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)*/
1730: /* Loop over columns 1/3: Left Back Down*/
1731: if (!dummyStart[0]) {
1732: const PetscInt neighbor = 0;
1733: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1734: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1735: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1736: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1737: }
1738: }
1739: } else {
1740: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1741: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1742: idxGlobalAll[count] = -1;
1743: }
1744: }
1745: }
1747: /* Loop over columns 2/3: (Middle) Down Back */
1748: {
1749: const PetscInt neighbor = 1;
1750: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1751: const PetscInt i = ighost - ghostOffsetStart[0];
1752: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1753: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1754: }
1755: }
1756: }
1758: /* Loop over columns 3/3: Right Down Back */
1759: if (!dummyEnd[0]) {
1760: const PetscInt neighbor = 2;
1761: PetscInt i;
1762: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1763: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1764: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1765: }
1766: }
1767: } else {
1768: /* Partial dummy entries on (Middle) Down Back neighbor */
1769: const PetscInt neighbor = 1;
1770: PetscInt i,dLocal;
1771: i = stag->n[0];
1772: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1773: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1774: }
1775: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1776: idxGlobalAll[count] = -1;
1777: }
1778: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1779: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1780: }
1781: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1782: idxGlobalAll[count] = -1;
1783: }
1784: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1785: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1786: }
1787: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1788: idxGlobalAll[count] = -1;
1789: }
1790: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1791: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1792: }
1793: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1794: idxGlobalAll[count] = -1;
1795: }
1796: ++i;
1797: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1798: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1799: idxGlobalAll[count] = -1;
1800: }
1801: }
1802: }
1803: }
1804: } else {
1805: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1806: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
1807: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1808: idxGlobalAll[count] = -1;
1809: }
1810: }
1811: }
1812: }
1814: /* Loop over rows 2/3: (Middle) Back */
1815: {
1816: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
1817: const PetscInt j = jghost - ghostOffsetStart[1];
1819: /* Loop over columns 1/3: Left (Middle) Back */
1820: if (!star && !dummyStart[0]) {
1821: const PetscInt neighbor = 3;
1822: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1823: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1824: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1825: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1826: }
1827: }
1828: } else {
1829: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1830: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1831: idxGlobalAll[count] = -1;
1832: }
1833: }
1834: }
1836: /* Loop over columns 2/3: (Middle) (Middle) Back */
1837: {
1838: const PetscInt neighbor = 4;
1839: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1840: const PetscInt i = ighost - ghostOffsetStart[0];
1841: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1842: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1843: }
1844: }
1845: }
1847: /* Loop over columns 3/3: Right (Middle) Back */
1848: if (!star && !dummyEnd[0]) {
1849: const PetscInt neighbor = 5;
1850: PetscInt i;
1851: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1852: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1853: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1854: }
1855: }
1856: } else if (dummyEnd[0]) {
1857: /* Partial dummy entries on (Middle) (Middle) Back rank */
1858: const PetscInt neighbor = 4;
1859: PetscInt i,dLocal;
1860: i = stag->n[0];
1861: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1862: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1863: }
1864: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1865: idxGlobalAll[count] = -1;
1866: }
1867: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1868: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1869: }
1870: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1871: idxGlobalAll[count] = -1;
1872: }
1873: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1874: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1875: }
1876: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1877: idxGlobalAll[count] = -1;
1878: }
1879: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1880: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1881: }
1882: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1883: idxGlobalAll[count] = -1;
1884: }
1885: ++i;
1886: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1887: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1888: idxGlobalAll[count] = -1;
1889: }
1890: }
1891: } else {
1892: /* Right (Middle) Back dummies */
1893: PetscInt i;
1894: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1895: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1896: idxGlobalAll[count] = -1;
1897: }
1898: }
1899: }
1900: }
1901: }
1903: /* Loop over rows 3/3: Up Back */
1904: if (!star && !dummyEnd[1]) {
1905: PetscInt j;
1906: for (j=0; j<ghostOffsetEnd[1]; ++j) {
1907: /* Loop over columns 1/3: Left Up Back*/
1908: if (!dummyStart[0]) {
1909: const PetscInt neighbor = 6;
1910: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1911: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1912: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1913: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1914: }
1915: }
1916: } else {
1917: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1918: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1919: idxGlobalAll[count] = -1;
1920: }
1921: }
1922: }
1924: /* Loop over columns 2/3: (Middle) Up Back*/
1925: {
1926: const PetscInt neighbor = 7;
1927: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1928: const PetscInt i = ighost - ghostOffsetStart[0];
1929: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1930: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1931: }
1932: }
1933: }
1935: /* Loop over columns 3/3: Right Up Back*/
1936: if (!dummyEnd[0]) {
1937: const PetscInt neighbor = 8;
1938: PetscInt i;
1939: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1940: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1941: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1942: }
1943: }
1944: } else {
1945: /* Partial dummies on (Middle) Up Back */
1946: const PetscInt neighbor = 7;
1947: PetscInt i,dLocal;
1948: i = stag->n[0];
1949: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1950: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1951: }
1952: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1953: idxGlobalAll[count] = -1;
1954: }
1955: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1956: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1957: }
1958: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1959: idxGlobalAll[count] = -1;
1960: }
1961: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1962: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1963: }
1964: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1965: idxGlobalAll[count] = -1;
1966: }
1967: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1968: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1969: }
1970: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1971: idxGlobalAll[count] = -1;
1972: }
1973: ++i;
1974: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1975: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1976: idxGlobalAll[count] = -1;
1977: }
1978: }
1979: }
1980: }
1981: } else if (dummyEnd[1]) {
1982: /* Up Back partial dummy row */
1983: PetscInt j = stag->n[1];
1985: if (!star && !dummyStart[0]) {
1986: /* Up Back partial dummy row: Loop over columns 1/3: Left Up Back, on Left (Middle) Back rank */
1987: const PetscInt neighbor = 3;
1988: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1989: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1990: PetscInt dLocal;
1991: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
1992: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1993: }
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: } else {
2006: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2007: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2008: idxGlobalAll[count] = -1;
2009: }
2010: }
2011: }
2013: /* Up Back partial dummy row: Loop over columns 2/3: (Middle) Up Back, on (Middle) (Middle) Back rank */
2014: {
2015: const PetscInt neighbor = 4;
2016: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2017: const PetscInt i = ighost - ghostOffsetStart[0];
2018: PetscInt dLocal;
2019: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2020: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2021: }
2022: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2023: idxGlobalAll[count] = -1;
2024: }
2025: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2026: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2027: }
2028: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2029: idxGlobalAll[count] = -1;
2030: }
2031: }
2032: }
2034: /* Up Back partial dummy row: Loop over columns 3/3: Right Up Back, on Right (Middle) Back rank */
2035: if (!star && !dummyEnd[0]) {
2036: const PetscInt neighbor = 5;
2037: PetscInt i;
2038: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2039: PetscInt dLocal;
2040: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2041: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2042: }
2043: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2044: idxGlobalAll[count] = -1;
2045: }
2046: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2047: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2048: }
2049: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2050: idxGlobalAll[count] = -1;
2051: }
2052: }
2053: } else if (dummyEnd[0]) {
2054: /* Up Right Back partial dummy element, on (Middle) (Middle) Back rank */
2055: const PetscInt neighbor = 4;
2056: PetscInt i,dLocal;
2057: i = stag->n[0];
2058: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2059: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2060: }
2061: 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 */
2062: idxGlobalAll[count] = -1;
2063: }
2064: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2065: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2066: }
2067: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2068: idxGlobalAll[count] = -1;
2069: }
2070: ++i;
2071: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2072: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2073: idxGlobalAll[count] = -1;
2074: }
2075: }
2076: } else {
2077: /* Up Right Back dummies */
2078: PetscInt i;
2079: for(i=0; i<ghostOffsetEnd[0]; ++i) {
2080: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2081: idxGlobalAll[count] = -1;
2082: }
2083: }
2084: }
2085: ++j;
2086: /* Up Back additional dummy rows */
2087: for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2088: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2089: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2090: idxGlobalAll[count] = -1;
2091: }
2092: }
2093: }
2094: } else {
2095: /* Up Back dummy rows */
2096: PetscInt j;
2097: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2098: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2099: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2100: idxGlobalAll[count] = -1;
2101: }
2102: }
2103: }
2104: }
2105: }
2106: } else {
2107: for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
2108: for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
2109: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2110: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2111: idxGlobalAll[count] = -1;
2112: }
2113: }
2114: }
2115: }
2116: }
2118: /* Loop over layers 2/3 : (Middle) */
2119: {
2120: for (kghost = ghostOffsetStart[2]; kghost<stag->nGhost[2]-ghostOffsetEnd[2]; ++kghost) {
2121: const PetscInt k = kghost - ghostOffsetStart[2];
2123: /* Loop over rows 1/3: Down (Middle) */
2124: if (!dummyStart[1]) {
2125: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2126: const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* down neighbor (10) always exists */
2128: /* Loop over columns 1/3: Left Down (Middle) */
2129: if (!star && !dummyStart[0]) {
2130: const PetscInt neighbor = 9;
2131: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2132: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2133: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2134: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2135: }
2136: }
2137: } else {
2138: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2139: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2140: idxGlobalAll[count] = -1;
2141: }
2142: }
2143: }
2145: /* Loop over columns 2/3: (Middle) Down (Middle) */
2146: {
2147: const PetscInt neighbor = 10;
2148: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2149: const PetscInt i = ighost - ghostOffsetStart[0];
2150: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2151: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2152: }
2153: }
2154: }
2156: /* Loop over columns 3/3: Right Down (Middle) */
2157: if (!star && !dummyEnd[0]) {
2158: const PetscInt neighbor = 11;
2159: PetscInt i;
2160: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2161: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2162: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2163: }
2164: }
2165: } else if (dummyEnd[0]) {
2166: /* Partial dummies on (Middle) Down (Middle) neighbor */
2167: const PetscInt neighbor = 10;
2168: PetscInt i,dLocal;
2169: i = stag->n[0];
2170: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2171: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2172: }
2173: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2174: idxGlobalAll[count] = -1;
2175: }
2176: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2177: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2178: }
2179: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2180: idxGlobalAll[count] = -1;
2181: }
2182: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2183: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2184: }
2185: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2186: idxGlobalAll[count] = -1;
2187: }
2188: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2189: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2190: }
2191: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2192: idxGlobalAll[count] = -1;
2193: }
2194: ++i;
2195: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2196: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2197: idxGlobalAll[count] = -1;
2198: }
2199: }
2200: } else {
2201: /* Right Down (Middle) dummies */
2202: PetscInt i;
2203: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2204: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2205: idxGlobalAll[count] = -1;
2206: }
2207: }
2208: }
2209: }
2210: } else {
2211: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2212: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2213: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2214: idxGlobalAll[count] = -1;
2215: }
2216: }
2217: }
2218: }
2220: /* Loop over rows 2/3: (Middle) (Middle) */
2221: {
2222: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2223: const PetscInt j = jghost - ghostOffsetStart[1];
2225: /* Loop over columns 1/3: Left (Middle) (Middle) */
2226: if (!dummyStart[0]) {
2227: const PetscInt neighbor = 12;
2228: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2229: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2230: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2231: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2232: }
2233: }
2234: } else {
2235: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2236: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2237: idxGlobalAll[count] = -1;
2238: }
2239: }
2240: }
2242: /* Loop over columns 2/3: (Middle) (Middle) (Middle) aka (Here) */
2243: {
2244: const PetscInt neighbor = 13;
2245: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2246: const PetscInt i = ighost - ghostOffsetStart[0];
2247: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2248: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2249: }
2250: }
2251: }
2253: /* Loop over columns 3/3: Right (Middle) (Middle) */
2254: if (!dummyEnd[0]) {
2255: const PetscInt neighbor = 14;
2256: PetscInt i;
2257: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2258: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2259: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2260: }
2261: }
2262: } else {
2263: /* Partial dummies on this rank */
2264: const PetscInt neighbor = 13;
2265: PetscInt i,dLocal;
2266: i = stag->n[0];
2267: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2268: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2269: }
2270: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2271: idxGlobalAll[count] = -1;
2272: }
2273: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2274: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2275: }
2276: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2277: idxGlobalAll[count] = -1;
2278: }
2279: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2280: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2281: }
2282: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2283: idxGlobalAll[count] = -1;
2284: }
2285: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2286: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2287: }
2288: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2289: idxGlobalAll[count] = -1;
2290: }
2291: ++i;
2292: for(;i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2293: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2294: idxGlobalAll[count] = -1;
2295: }
2296: }
2297: }
2298: }
2299: }
2301: /* Loop over rows 3/3: Up (Middle) */
2302: if (!dummyEnd[1]) {
2303: PetscInt j;
2304: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2306: /* Loop over columns 1/3: Left Up (Middle) */
2307: if (!star && !dummyStart[0]) {
2308: const PetscInt neighbor = 15;
2309: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2310: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2311: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2312: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2313: }
2314: }
2315: } else {
2316: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2317: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2318: idxGlobalAll[count] = -1;
2319: }
2320: }
2321: }
2323: /* Loop over columns 2/3: (Middle) Up (Middle) */
2324: {
2325: const PetscInt neighbor = 16;
2326: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2327: const PetscInt i = ighost - ghostOffsetStart[0];
2328: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2329: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2330: }
2331: }
2332: }
2334: /* Loop over columns 3/3: Right Up (Middle) */
2335: if (!star && !dummyEnd[0]) {
2336: const PetscInt neighbor = 17;
2337: PetscInt i;
2338: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2339: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2340: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2341: }
2342: }
2343: } else if (dummyEnd[0]) {
2344: /* Partial dummies on (Middle) Up (Middle) rank */
2345: const PetscInt neighbor = 16;
2346: PetscInt i,dLocal;
2347: i = stag->n[0];
2348: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2349: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2350: }
2351: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2352: idxGlobalAll[count] = -1;
2353: }
2354: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2355: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2356: }
2357: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2358: idxGlobalAll[count] = -1;
2359: }
2360: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2361: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2362: }
2363: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2364: idxGlobalAll[count] = -1;
2365: }
2366: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2367: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2368: }
2369: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2370: idxGlobalAll[count] = -1;
2371: }
2372: ++i;
2373: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2374: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2375: idxGlobalAll[count] = -1;
2376: }
2377: }
2378: } else {
2379: /* Right Up (Middle) dumies */
2380: PetscInt i;
2381: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2382: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2383: idxGlobalAll[count] = -1;
2384: }
2385: }
2386: }
2387: }
2388: } else {
2389: /* Up (Middle) partial dummy row */
2390: PetscInt j = stag->n[1];
2392: /* Up (Middle) partial dummy row: colums 1/3: Left Up (Middle), on Left (Middle) (Middle) rank */
2393: if (!dummyStart[0]) {
2394: const PetscInt neighbor = 12;
2395: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2396: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
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: } else {
2412: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2413: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2414: idxGlobalAll[count] = -1;
2415: }
2416: }
2417: }
2419: /* Up (Middle) partial dummy row: columns 2/3: (Middle) Up (Middle), on (Middle) (Middle) (Middle) = here rank */
2420: {
2421: const PetscInt neighbor = 13;
2422: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2423: const PetscInt i = ighost - ghostOffsetStart[0];
2424: PetscInt dLocal;
2425: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2426: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2427: }
2428: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2429: idxGlobalAll[count] = -1;
2430: }
2431: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2432: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2433: }
2434: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2435: idxGlobalAll[count] = -1;
2436: }
2437: }
2438: }
2440: if (!dummyEnd[0]) {
2441: /* Up (Middle) partial dummy row: columns 3/3: Right Up (Middle), on Right (Middle) (Middle) rank */
2442: const PetscInt neighbor = 14;
2443: PetscInt i;
2444: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2445: PetscInt dLocal;
2446: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2447: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2448: }
2449: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2450: idxGlobalAll[count] = -1;
2451: }
2452: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2453: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2454: }
2455: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2456: idxGlobalAll[count] = -1;
2457: }
2458: }
2459: } else {
2460: /* Up (Middle) Right partial dummy element, on (Middle) (Middle) (Middle) = here rank */
2461: const PetscInt neighbor = 13;
2462: PetscInt i,dLocal;
2463: i = stag->n[0];
2464: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2465: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2466: }
2467: 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 */
2468: idxGlobalAll[count] = -1;
2469: }
2470: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2471: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2472: }
2473: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2474: idxGlobalAll[count] = -1;
2475: }
2476: ++i;
2477: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2478: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2479: idxGlobalAll[count] = -1;
2480: }
2481: }
2482: }
2483: ++j;
2484: /* Up (Middle) additional dummy rows */
2485: for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2486: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2487: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2488: idxGlobalAll[count] = -1;
2489: }
2490: }
2491: }
2492: }
2493: }
2494: }
2496: /* Loop over layers 3/3 : Front */
2497: if (!dummyEnd[2]) {
2498: PetscInt k;
2499: for (k=0; k<ghostOffsetEnd[2]; ++k) {
2501: /* Loop over rows 1/3: Down Front */
2502: if (!star && !dummyStart[1]) {
2503: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2504: const PetscInt j = nNeighbors[19][1] - ghostOffsetStart[1] + jghost; /* Constant for whole row (use down front neighbor) */
2506: /* Loop over columns 1/3: Left Down Front */
2507: if (!dummyStart[0]) {
2508: const PetscInt neighbor = 18;
2509: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2510: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2511: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2512: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2513: }
2514: }
2515: } else {
2516: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2517: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2518: idxGlobalAll[count] = -1;
2519: }
2520: }
2521: }
2523: /* Loop over columns 2/3: (Middle) Down Front */
2524: {
2525: const PetscInt neighbor = 19;
2526: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2527: const PetscInt i = ighost - ghostOffsetStart[0];
2528: for (d=0; d<stag->entriesPerElement; ++d,++count) { /* vertex, 2 edges, and face associated with back face */
2529: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2530: }
2531: }
2532: }
2534: /* Loop over columns 3/3: Right Down Front */
2535: if (!dummyEnd[0]) {
2536: const PetscInt neighbor = 20;
2537: PetscInt i;
2538: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2539: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2540: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2541: }
2542: }
2543: } else {
2544: /* Partial dummy element on (Middle) Down Front rank */
2545: const PetscInt neighbor = 19;
2546: PetscInt i,dLocal;
2547: i = stag->n[0];
2548: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2549: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2550: }
2551: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2552: idxGlobalAll[count] = -1;
2553: }
2554: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2555: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2556: }
2557: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2558: idxGlobalAll[count] = -1;
2559: }
2560: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2561: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2562: }
2563: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2564: idxGlobalAll[count] = -1;
2565: }
2566: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2567: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2568: }
2569: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2570: idxGlobalAll[count] = -1;
2571: }
2572: ++i;
2573: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2574: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2575: idxGlobalAll[count] = -1;
2576: }
2577: }
2578: }
2579: }
2580: } else {
2581: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2582: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2583: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2584: idxGlobalAll[count] = -1;
2585: }
2586: }
2587: }
2588: }
2590: /* Loop over rows 2/3: (Middle) Front */
2591: {
2592: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2593: const PetscInt j = jghost - ghostOffsetStart[1];
2595: /* Loop over columns 1/3: Left (Middle) Front*/
2596: if (!star && !dummyStart[0]) {
2597: const PetscInt neighbor = 21;
2598: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2599: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2600: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2601: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2602: }
2603: }
2604: } else {
2605: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2606: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2607: idxGlobalAll[count] = -1;
2608: }
2609: }
2610: }
2612: /* Loop over columns 2/3: (Middle) (Middle) Front*/
2613: {
2614: const PetscInt neighbor = 22;
2615: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2616: const PetscInt i = ighost - ghostOffsetStart[0];
2617: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2618: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2619: }
2620: }
2621: }
2623: /* Loop over columns 3/3: Right (Middle) Front*/
2624: if (!star && !dummyEnd[0]) {
2625: const PetscInt neighbor = 23;
2626: PetscInt i;
2627: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2628: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2629: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2630: }
2631: }
2632: } else if (dummyEnd[0]) {
2633: /* Partial dummy element on (Middle) (Middle) Front element */
2634: const PetscInt neighbor = 22;
2635: PetscInt i,dLocal;
2636: i = stag->n[0];
2637: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2638: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2639: }
2640: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2641: idxGlobalAll[count] = -1;
2642: }
2643: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2644: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2645: }
2646: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2647: idxGlobalAll[count] = -1;
2648: }
2649: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2650: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2651: }
2652: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2653: idxGlobalAll[count] = -1;
2654: }
2655: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2656: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2657: }
2658: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2659: idxGlobalAll[count] = -1;
2660: }
2661: ++i;
2662: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2663: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2664: idxGlobalAll[count] = -1;
2665: }
2666: }
2667: } else {
2668: /* Right (Middle) Front dummies */
2669: PetscInt i;
2670: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2671: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2672: idxGlobalAll[count] = -1;
2673: }
2674: }
2675: }
2676: }
2677: }
2679: /* Loop over rows 3/3: Up Front */
2680: if (!star && !dummyEnd[1]) {
2681: PetscInt j;
2682: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2684: /* Loop over columns 1/3: Left Up Front */
2685: if (!dummyStart[0]) {
2686: const PetscInt neighbor = 24;
2687: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2688: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2689: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2690: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2691: }
2692: }
2693: } else {
2694: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2695: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2696: idxGlobalAll[count] = -1;
2697: }
2698: }
2699: }
2701: /* Loop over columns 2/3: (Middle) Up Front */
2702: {
2703: const PetscInt neighbor = 25;
2704: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2705: const PetscInt i = ighost - ghostOffsetStart[0];
2706: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2707: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2708: }
2709: }
2710: }
2712: /* Loop over columns 3/3: Right Up Front */
2713: if (!dummyEnd[0]) {
2714: const PetscInt neighbor = 26;
2715: PetscInt i;
2716: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2717: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2718: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2719: }
2720: }
2721: } else {
2722: /* Partial dummy element on (Middle) Up Front rank */
2723: const PetscInt neighbor = 25;
2724: PetscInt i,dLocal;
2725: i = stag->n[0];
2726: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2727: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2728: }
2729: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2730: idxGlobalAll[count] = -1;
2731: }
2732: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2733: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2734: }
2735: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2736: idxGlobalAll[count] = -1;
2737: }
2738: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2739: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2740: }
2741: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2742: idxGlobalAll[count] = -1;
2743: }
2744: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2745: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2746: }
2747: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2748: idxGlobalAll[count] = -1;
2749: }
2750: ++i;
2751: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2752: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2753: idxGlobalAll[count] = -1;
2754: }
2755: }
2756: }
2757: }
2758: } else if (dummyEnd[1]) {
2759: /* Up Front partial dummy row */
2760: PetscInt j = stag->n[1];
2762: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) Front rank */
2763: if (!star && !dummyStart[0]) {
2764: const PetscInt neighbor = 21;
2765: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2766: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
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: } else {
2782: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2783: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2784: idxGlobalAll[count] = -1;
2785: }
2786: }
2787: }
2789: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) Front rank */
2790: {
2791: const PetscInt neighbor = 22;
2792: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2793: const PetscInt i = ighost - ghostOffsetStart[0];
2794: PetscInt dLocal;
2795: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2796: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2797: }
2798: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2799: idxGlobalAll[count] = -1;
2800: }
2801: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2802: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2803: }
2804: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2805: idxGlobalAll[count] = -1;
2806: }
2807: }
2808: }
2810: if (!star && !dummyEnd[0]) {
2811: /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) Front rank */
2812: const PetscInt neighbor = 23;
2813: PetscInt i;
2814: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2815: PetscInt dLocal;
2816: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2817: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2818: }
2819: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2820: idxGlobalAll[count] = -1;
2821: }
2822: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2823: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2824: }
2825: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2826: idxGlobalAll[count] = -1;
2827: }
2828: }
2830: } else if (dummyEnd[0]) {
2831: /* Partial Right Up Front dummy element, on (Middle) (Middle) Front rank */
2832: const PetscInt neighbor = 22;
2833: PetscInt i,dLocal;
2834: i = stag->n[0];
2835: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2836: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2837: }
2838: 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 */
2839: idxGlobalAll[count] = -1;
2840: }
2841: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2842: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2843: }
2844: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2845: idxGlobalAll[count] = -1;
2846: }
2847: ++i;
2848: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2849: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2850: idxGlobalAll[count] = -1;
2851: }
2852: }
2853: } else {
2854: /* Right Up Front dummies */
2855: PetscInt i;
2856: for(i=0; i<ghostOffsetEnd[0]; ++i) {
2857: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2858: idxGlobalAll[count] = -1;
2859: }
2860: }
2861: }
2862: ++j;
2863: /* Up Front additional dummy rows */
2864: for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2865: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2866: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2867: idxGlobalAll[count] = -1;
2868: }
2869: }
2870: }
2871: } else {
2872: /* Up Front dummies */
2873: PetscInt j;
2874: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2875: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2876: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2877: idxGlobalAll[count] = -1;
2878: }
2879: }
2880: }
2881: }
2882: }
2883: } else {
2884: PetscInt k = stag->n[2];
2886: /* Front partial ghost layer: rows 1/3: Down Front, on Down (Middle) */
2887: if (!dummyStart[1]) {
2888: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2889: const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* wrt down neighbor (10) */
2891: /* Down Front partial ghost row: colums 1/3: Left Down Front, on Left Down (Middle) */
2892: if (!star && !dummyStart[0]) {
2893: const PetscInt neighbor = 9;
2894: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
2895: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2896: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2897: PetscInt dLocal;
2898: 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 */
2899: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2900: }
2901: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2902: idxGlobalAll[count] = -1;
2903: }
2904: }
2905: } else {
2906: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2907: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2908: idxGlobalAll[count] = -1;
2909: }
2910: }
2911: }
2913: /* Down Front partial ghost row: columns 2/3: (Middle) Down Front, on (Middle) Down (Middle) */
2914: {
2915: const PetscInt neighbor = 10;
2916: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2917: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2918: const PetscInt i = ighost - ghostOffsetStart[0];
2919: PetscInt dLocal;
2920: 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 */
2921: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2922: }
2923: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2924: idxGlobalAll[count] = -1;
2925: }
2926: }
2927: }
2929: if (!star && !dummyEnd[0]) {
2930: /* Down Front partial dummy row: columns 3/3: Right Down Front, on Right Down (Middle) */
2931: const PetscInt neighbor = 11;
2932: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
2933: PetscInt i;
2934: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2935: PetscInt dLocal;
2936: for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
2937: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2938: }
2939: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2940: idxGlobalAll[count] = -1;
2941: }
2942: }
2943: } else if (dummyEnd[0]) {
2944: /* Right Down Front partial dummy element, living on (Middle) Down (Middle) rank */
2945: const PetscInt neighbor = 10;
2946: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2947: PetscInt i,dLocal;
2948: i = stag->n[0];
2949: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2950: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2951: }
2952: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2953: idxGlobalAll[count] = -1;
2954: }
2955: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
2956: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2957: }
2958: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
2959: idxGlobalAll[count] = -1;
2960: }
2961: ++i;
2962: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2963: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2964: idxGlobalAll[count] = -1;
2965: }
2966: }
2967: } else {
2968: PetscInt i;
2969: /* Right Down Front dummies */
2970: for(i=0; i<ghostOffsetEnd[0]; ++i) {
2971: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2972: idxGlobalAll[count] = -1;
2973: }
2974: }
2975: }
2976: }
2977: } else {
2978: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2979: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2980: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2981: idxGlobalAll[count] = -1;
2982: }
2983: }
2984: }
2985: }
2987: /* Front partial ghost layer: rows 2/3: (Middle) Front, on (Middle) (Middle) */
2988: {
2989: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2990: const PetscInt j = jghost - ghostOffsetStart[1];
2992: /* (Middle) Front partial dummy row: columns 1/3: Left (Middle) Front, on Left (Middle) (Middle) */
2993: if (!dummyStart[0]) {
2994: const PetscInt neighbor = 12;
2995: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
2996: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2997: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2998: PetscInt dLocal;
2999: 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 */
3000: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3001: }
3002: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3003: idxGlobalAll[count] = -1;
3004: }
3005: }
3006: } else {
3007: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3008: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3009: idxGlobalAll[count] = -1;
3010: }
3011: }
3012: }
3014: /* (Middle) Front partial dummy row: columns 2/3: (Middle) (Middle) Front, on (Middle) (Middle) (Middle) */
3015: {
3016: const PetscInt neighbor = 13;
3017: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3018: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3019: const PetscInt i = ighost - ghostOffsetStart[0];
3020: PetscInt dLocal;
3021: 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 */
3022: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3023: }
3024: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3025: idxGlobalAll[count] = -1;
3026: }
3027: }
3028: }
3030: if (!dummyEnd[0]) {
3031: /* (Middle) Front partial dummy row: columns 3/3: Right (Middle) Front, on Right (Middle) (Middle) */
3032: const PetscInt neighbor = 14;
3033: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3034: PetscInt i;
3035: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3036: PetscInt dLocal;
3037: for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3038: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3039: }
3040: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3041: idxGlobalAll[count] = -1;
3042: }
3043: }
3044: } else {
3045: /* Right (Middle) Front partial dummy element, on (Middle) (Middle) (Middle) rank */
3046: const PetscInt neighbor = 13;
3047: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3048: PetscInt i,dLocal;
3049: i = stag->n[0];
3050: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3051: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3052: }
3053: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3054: idxGlobalAll[count] = -1;
3055: }
3056: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3057: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3058: }
3059: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3060: idxGlobalAll[count] = -1;
3061: }
3062: ++i;
3063: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3064: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3065: idxGlobalAll[count] = -1;
3066: }
3067: }
3068: }
3069: }
3070: }
3072: /* Front partial ghost layer: rows 3/3: Up Front, on Up (Middle) */
3073: if (!dummyEnd[1]) {
3074: PetscInt j;
3075: for (j=0; j<ghostOffsetEnd[1]; ++j) {
3077: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left Up (Middle) */
3078: if (!star && !dummyStart[0]) {
3079: const PetscInt neighbor = 15;
3080: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3081: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3082: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3083: PetscInt dLocal;
3084: 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 */
3085: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3086: }
3087: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3088: idxGlobalAll[count] = -1;
3089: }
3090: }
3091: } else {
3092: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3093: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3094: idxGlobalAll[count] = -1;
3095: }
3096: }
3097: }
3099: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) Up (Middle) */
3100: {
3101: const PetscInt neighbor = 16;
3102: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3103: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3104: const PetscInt i = ighost - ghostOffsetStart[0];
3105: PetscInt dLocal;
3106: 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 */
3107: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3108: }
3109: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3110: idxGlobalAll[count] = -1;
3111: }
3112: }
3113: }
3115: if (!star && !dummyEnd[0]) {
3116: /* Up Front partial dummy row: columsn 3/3: Right Up Front, on Right Up (Middle) */
3117: const PetscInt neighbor = 17;
3118: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3119: PetscInt i;
3120: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3121: PetscInt dLocal;
3122: for (d=0,dLocal=0; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Vertex, back down edge, back face, back left edge */
3123: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3124: }
3125: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3126: idxGlobalAll[count] = -1;
3127: }
3128: }
3129: } else if (dummyEnd[0]) {
3130: /* Right Up Front partial dummy element, on (Middle) Up (Middle) */
3131: const PetscInt neighbor = 16;
3132: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3133: PetscInt i,dLocal;
3134: i = stag->n[0];
3135: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3136: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3137: }
3138: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3139: idxGlobalAll[count] = -1;
3140: }
3141: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3142: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3143: }
3144: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3145: idxGlobalAll[count] = -1;
3146: }
3147: ++i;
3148: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3149: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3150: idxGlobalAll[count] = -1;
3151: }
3152: }
3153: } else {
3154: /* Right Up Front dummies */
3155: PetscInt i;
3156: /* Right Down Front dummies */
3157: for(i=0; i<ghostOffsetEnd[0]; ++i) {
3158: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3159: idxGlobalAll[count] = -1;
3160: }
3161: }
3162: }
3163: }
3164: } else {
3165: /* Up Front partial dummy row */
3166: PetscInt j = stag->n[1];
3168: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) (Middle) */
3169: if (!dummyStart[0]) {
3170: const PetscInt neighbor = 12;
3171: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3172: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3173: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3174: PetscInt dLocal;
3175: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3176: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3177: }
3178: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3179: idxGlobalAll[count] = -1;
3180: }
3181: }
3182: } else {
3183: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3184: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3185: idxGlobalAll[count] = -1;
3186: }
3187: }
3188: }
3190: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) (Middle) */
3191: {
3192: const PetscInt neighbor = 13;
3193: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3194: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3195: const PetscInt i = ighost - ghostOffsetStart[0];
3196: PetscInt dLocal;
3197: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3198: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3199: }
3200: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3201: idxGlobalAll[count] = -1;
3202: }
3203: }
3204: }
3206: if (!dummyEnd[0]) {
3207: /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) (Middle) */
3208: const PetscInt neighbor = 14;
3209: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3210: PetscInt i;
3211: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3212: PetscInt dLocal;
3213: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3214: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3215: }
3216: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3217: idxGlobalAll[count] = -1;
3218: }
3219: }
3220: } else {
3221: /* Right Up Front partial dummy element, on this rank */
3222: const PetscInt neighbor = 13;
3223: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3224: PetscInt i,dLocal;
3225: i = stag->n[0];
3226: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3227: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3228: }
3229: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3230: idxGlobalAll[count] = -1;
3231: }
3232: ++i;
3233: for(; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3234: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3235: idxGlobalAll[count] = -1;
3236: }
3237: }
3238: }
3239: ++j;
3240: /* Up Front additional dummy rows */
3241: for(; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
3242: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3243: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3244: idxGlobalAll[count] = -1;
3245: }
3246: }
3247: }
3248: }
3249: /* Front additional dummy layers */
3250: ++k;
3251: for (; k<stag->n[2] + ghostOffsetEnd[2]; ++k) {
3252: for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
3253: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3254: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3255: idxGlobalAll[count] = -1;
3256: }
3257: }
3258: }
3259: }
3260: }
3262: /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
3263: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
3264: PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
3265: return(0);
3266: }
3268: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM dm)
3269: {
3270: PetscErrorCode ierr;
3271: DM_Stag * const stag = (DM_Stag*)dm->data;
3272: const PetscInt epe = stag->entriesPerElement;
3273: const PetscInt epr = stag->nGhost[0]*epe;
3274: const PetscInt epl = stag->nGhost[1]*epr;
3277: PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
3278: stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] = 0;
3279: stag->locationOffsets[DMSTAG_BACK_DOWN] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + stag->dof[0];
3280: stag->locationOffsets[DMSTAG_BACK_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epe;
3281: stag->locationOffsets[DMSTAG_BACK_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN] + stag->dof[1];
3282: stag->locationOffsets[DMSTAG_BACK] = stag->locationOffsets[DMSTAG_BACK_LEFT] + stag->dof[1];
3283: stag->locationOffsets[DMSTAG_BACK_RIGHT] = stag->locationOffsets[DMSTAG_BACK_LEFT] + epe;
3284: stag->locationOffsets[DMSTAG_BACK_UP_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epr;
3285: stag->locationOffsets[DMSTAG_BACK_UP] = stag->locationOffsets[DMSTAG_BACK_DOWN] + epr;
3286: stag->locationOffsets[DMSTAG_BACK_UP_RIGHT] = stag->locationOffsets[DMSTAG_BACK_UP_LEFT] + epe;
3287: stag->locationOffsets[DMSTAG_DOWN_LEFT] = stag->locationOffsets[DMSTAG_BACK] + stag->dof[2];
3288: stag->locationOffsets[DMSTAG_DOWN] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + stag->dof[1];
3289: stag->locationOffsets[DMSTAG_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epe;
3290: stag->locationOffsets[DMSTAG_LEFT] = stag->locationOffsets[DMSTAG_DOWN] + stag->dof[2];
3291: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[2];
3292: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
3293: stag->locationOffsets[DMSTAG_UP_LEFT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epr;
3294: stag->locationOffsets[DMSTAG_UP] = stag->locationOffsets[DMSTAG_DOWN] + epr;
3295: stag->locationOffsets[DMSTAG_UP_RIGHT] = stag->locationOffsets[DMSTAG_UP_LEFT] + epe;
3296: stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epl;
3297: stag->locationOffsets[DMSTAG_FRONT_DOWN] = stag->locationOffsets[DMSTAG_BACK_DOWN] + epl;
3298: stag->locationOffsets[DMSTAG_FRONT_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epe;
3299: stag->locationOffsets[DMSTAG_FRONT_LEFT] = stag->locationOffsets[DMSTAG_BACK_LEFT] + epl;
3300: stag->locationOffsets[DMSTAG_FRONT] = stag->locationOffsets[DMSTAG_BACK] + epl;
3301: stag->locationOffsets[DMSTAG_FRONT_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_LEFT] + epe;
3302: stag->locationOffsets[DMSTAG_FRONT_UP_LEFT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epr;
3303: stag->locationOffsets[DMSTAG_FRONT_UP] = stag->locationOffsets[DMSTAG_FRONT_DOWN] + epr;
3304: stag->locationOffsets[DMSTAG_FRONT_UP_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_UP_LEFT] + epe;
3305: return(0);
3306: }
3308: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_3d(DM dm)
3309: {
3310: PetscErrorCode ierr;
3311: DM_Stag * const stag = (DM_Stag*)dm->data;
3312: PetscInt *idxLocal,*idxGlobal,*globalOffsetsRecomputed;
3313: const PetscInt *globalOffsets;
3314: PetscInt count,d,entriesPerEdge,entriesPerFace,eprGhost,eplGhost,ghostOffsetStart[3],ghostOffsetEnd[3];
3315: IS isLocal,isGlobal;
3316: PetscBool dummyEnd[3];
3319: DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsetsRecomputed); /* note that we don't actually use all of these. An available optimization is to pass them, when available */
3320: globalOffsets = globalOffsetsRecomputed;
3321: PetscMalloc1(stag->entries,&idxLocal);
3322: PetscMalloc1(stag->entries,&idxGlobal);
3323: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
3324: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
3325: entriesPerEdge = stag->dof[0] + stag->dof[1];
3326: eprGhost = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row */
3327: eplGhost = stag->nGhost[1]*eprGhost; /* epl = entries per (element) layer */
3328: count = 0;
3329: for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
3330: for (d=0; d<3; ++d) ghostOffsetEnd[d] = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
3331: {
3332: const PetscInt neighbor = 13;
3333: const PetscInt epr = stag->entriesPerElement * stag->n[0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
3334: const PetscInt epl = epr * stag->n[1] + (dummyEnd[1] ? stag->n[0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
3335: const PetscInt epFaceRow = entriesPerFace * stag->n[0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3336: const PetscInt start0 = 0;
3337: const PetscInt start1 = 0;
3338: const PetscInt start2 = 0;
3339: const PetscInt startGhost0 = ghostOffsetStart[0];
3340: const PetscInt startGhost1 = ghostOffsetStart[1];
3341: const PetscInt startGhost2 = ghostOffsetStart[2];
3342: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
3343: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
3344: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
3345: const PetscBool extra0 = dummyEnd[0];
3346: const PetscBool extra1 = dummyEnd[1];
3347: const PetscBool extra2 = dummyEnd[2];
3348: 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);
3349: }
3350: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
3351: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
3352: {
3353: Vec local,global;
3354: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
3355: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
3356: VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
3357: VecDestroy(&global);
3358: VecDestroy(&local);
3359: }
3360: ISDestroy(&isLocal);
3361: ISDestroy(&isGlobal);
3362: if (globalOffsetsRecomputed) {
3363: PetscFree(globalOffsetsRecomputed);
3364: }
3365: return(0);
3366: }