Actual source code: stag3d.c
petsc-3.14.6 2021-03-30
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: }
93: DMStagGetCorners(dmCoord,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
94: min[0] = xmin; min[1]= ymin; min[2] = zmin;
95: h[0] = (xmax-xmin)/stagCoord->N[0];
96: h[1] = (ymax-ymin)/stagCoord->N[1];
97: h[2] = (zmax-zmin)/stagCoord->N[2];
99: for (ind[2]=start[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 (PetscDefined(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: }
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 (!PetscDefined(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: }
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])) * stag->entriesPerElement
925: + (dummyEnd[0] ? (nNonDummyGhost[1] * stag->n[2] + stag->n[1] * nNonDummyGhost[2] - stag->n[1] * stag->n[2]) * entriesPerFace : 0)
926: + (dummyEnd[1] ? (nNonDummyGhost[0] * stag->n[2] + stag->n[0] * nNonDummyGhost[2] - stag->n[0] * stag->n[2]) * entriesPerFace : 0)
927: + (dummyEnd[2] ? (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * entriesPerFace : 0)
928: + (dummyEnd[0] && dummyEnd[1] ? nNonDummyGhost[2] * entriesPerEdge : 0)
929: + (dummyEnd[2] && dummyEnd[0] ? nNonDummyGhost[1] * entriesPerEdge : 0)
930: + (dummyEnd[1] && dummyEnd[2] ? nNonDummyGhost[0] * entriesPerEdge : 0)
931: + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ? entriesPerCorner : 0);
932: } else {
933: entriesToTransferTotal = nNonDummyGhost[0] * nNonDummyGhost[1] * nNonDummyGhost[2] * stag->entriesPerElement
934: + (dummyEnd[0] ? nNonDummyGhost[1] * nNonDummyGhost[2] * entriesPerFace : 0)
935: + (dummyEnd[1] ? nNonDummyGhost[2] * nNonDummyGhost[0] * entriesPerFace : 0)
936: + (dummyEnd[2] ? nNonDummyGhost[0] * nNonDummyGhost[1] * entriesPerFace : 0)
937: + (dummyEnd[0] && dummyEnd[1] ? nNonDummyGhost[2] * entriesPerEdge : 0)
938: + (dummyEnd[2] && dummyEnd[0] ? nNonDummyGhost[1] * entriesPerEdge : 0)
939: + (dummyEnd[1] && dummyEnd[2] ? nNonDummyGhost[0] * entriesPerEdge : 0)
940: + (dummyEnd[0] && dummyEnd[1] && dummyEnd[2] ? entriesPerCorner : 0);
941: }
942: }
944: /* Allocate arrays to populate */
945: PetscMalloc1(entriesToTransferTotal,&idxLocal);
946: PetscMalloc1(entriesToTransferTotal,&idxGlobal);
948: /* Counts into idxLocal/idxGlobal */
949: count = 0;
951: /* Loop over each of the 27 neighbor, populating idxLocal and idxGlobal. These
952: cases are principally distinguished by
954: 1. The loop bounds (i/ighost,j/jghost,k/kghost)
955: 2. the strides used in the loop body, which depend on whether the current
956: rank and/or its neighbor are a non-periodic right/top/front boundary, which has additional
957: points in the global representation.
958: - If the neighboring rank is a right/top/front boundary,
959: then eprNeighbor (entries per element row on the neighbor) and/or eplNeighbor (entries per element layer on the neighbor)
960: are different.
961: - If this rank is a non-periodic right/top/front boundary (checking entries of stag->lastRank),
962: there is an extra loop over 1-3 boundary faces)
964: Here, we do not include "dummy" dof (points in the local representation which
965: do not correspond to any global dof). This, and the fact that we iterate over points in terms of
966: increasing global ordering, are the main two differences from the construction of
967: the Local-to-global map, which iterates over points in local order, and does include dummy points. */
969: /* LEFT DOWN BACK */
970: if (!star && !dummyStart[0] && !dummyStart[1] && !dummyStart[2]) {
971: const PetscInt neighbor = 0;
972: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
973: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
974: const PetscInt epFaceRow = -1;
975: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
976: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
977: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
978: const PetscInt startGhost0 = 0;
979: const PetscInt startGhost1 = 0;
980: const PetscInt startGhost2 = 0;
981: const PetscInt endGhost0 = ghostOffsetStart[0];
982: const PetscInt endGhost1 = ghostOffsetStart[1];
983: const PetscInt endGhost2 = ghostOffsetStart[2];
984: const PetscBool extra0 = PETSC_FALSE;
985: const PetscBool extra1 = PETSC_FALSE;
986: const PetscBool extra2 = PETSC_FALSE;
987: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
988: }
990: /* DOWN BACK */
991: if (!star && !dummyStart[1] && !dummyStart[2]) {
992: const PetscInt neighbor = 1;
993: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
994: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
995: const PetscInt epFaceRow = -1;
996: const PetscInt start0 = 0;
997: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
998: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
999: const PetscInt startGhost0 = ghostOffsetStart[0];
1000: const PetscInt startGhost1 = 0;
1001: const PetscInt startGhost2 = 0;
1002: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1003: const PetscInt endGhost1 = ghostOffsetStart[1];
1004: const PetscInt endGhost2 = ghostOffsetStart[2];
1005: const PetscBool extra0 = dummyEnd[0];
1006: const PetscBool extra1 = PETSC_FALSE;
1007: const PetscBool extra2 = PETSC_FALSE;
1008: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1009: }
1011: /* RIGHT DOWN BACK */
1012: if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyStart[2]) {
1013: const PetscInt neighbor = 2;
1014: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1015: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1016: const PetscInt epFaceRow = -1;
1017: const PetscInt start0 = 0;
1018: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1019: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1020: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1021: const PetscInt startGhost1 = 0;
1022: const PetscInt startGhost2 = 0;
1023: const PetscInt endGhost0 = stag->nGhost[0];
1024: const PetscInt endGhost1 = ghostOffsetStart[1];
1025: const PetscInt endGhost2 = ghostOffsetStart[2];
1026: const PetscBool extra0 = PETSC_FALSE;
1027: const PetscBool extra1 = PETSC_FALSE;
1028: const PetscBool extra2 = PETSC_FALSE;
1029: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1030: }
1032: /* LEFT BACK */
1033: if (!star && !dummyStart[0] && !dummyStart[2]) {
1034: const PetscInt neighbor = 3;
1035: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1036: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* May be a top boundary */
1037: const PetscInt epFaceRow = -1;
1038: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1039: const PetscInt start1 = 0;
1040: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1041: const PetscInt startGhost0 = 0;
1042: const PetscInt startGhost1 = ghostOffsetStart[1];
1043: const PetscInt startGhost2 = 0;
1044: const PetscInt endGhost0 = ghostOffsetStart[0];
1045: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1046: const PetscInt endGhost2 = ghostOffsetStart[2];
1047: const PetscBool extra0 = PETSC_FALSE;
1048: const PetscBool extra1 = dummyEnd[1];
1049: const PetscBool extra2 = PETSC_FALSE;
1050: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1051: }
1053: /* BACK */
1054: if (!dummyStart[2]) {
1055: const PetscInt neighbor = 4;
1056: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1057: 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 */
1058: const PetscInt epFaceRow = -1;
1059: const PetscInt start0 = 0;
1060: const PetscInt start1 = 0;
1061: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1062: const PetscInt startGhost0 = ghostOffsetStart[0];
1063: const PetscInt startGhost1 = ghostOffsetStart[1];
1064: const PetscInt startGhost2 = 0;
1065: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1066: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1067: const PetscInt endGhost2 = ghostOffsetStart[2];
1068: const PetscBool extra0 = dummyEnd[0];
1069: const PetscBool extra1 = dummyEnd[1];
1070: const PetscBool extra2 = PETSC_FALSE;
1071: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1072: }
1074: /* RIGHT BACK */
1075: if (!star && !dummyEnd[0] && !dummyStart[2]) {
1076: const PetscInt neighbor = 5;
1077: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1078: 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 */
1079: const PetscInt epFaceRow = -1;
1080: const PetscInt start0 = 0;
1081: const PetscInt start1 = 0;
1082: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1083: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1084: const PetscInt startGhost1 = ghostOffsetStart[1];
1085: const PetscInt startGhost2 = 0;
1086: const PetscInt endGhost0 = stag->nGhost[0];
1087: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1088: const PetscInt endGhost2 = ghostOffsetStart[2];
1089: const PetscBool extra0 = PETSC_FALSE;
1090: const PetscBool extra1 = dummyEnd[1];
1091: const PetscBool extra2 = PETSC_FALSE;
1092: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1093: }
1095: /* LEFT UP BACK */
1096: if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyStart[2]) {
1097: const PetscInt neighbor = 6;
1098: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1099: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1100: const PetscInt epFaceRow = -1;
1101: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1102: const PetscInt start1 = 0;
1103: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1104: const PetscInt startGhost0 = 0;
1105: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1106: const PetscInt startGhost2 = 0;
1107: const PetscInt endGhost0 = ghostOffsetStart[0];
1108: const PetscInt endGhost1 = stag->nGhost[1];
1109: const PetscInt endGhost2 = ghostOffsetStart[2];
1110: const PetscBool extra0 = PETSC_FALSE;
1111: const PetscBool extra1 = PETSC_FALSE;
1112: const PetscBool extra2 = PETSC_FALSE;
1113: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1114: }
1116: /* UP BACK */
1117: if (!star && !dummyEnd[1] && !dummyStart[2]) {
1118: const PetscInt neighbor = 7;
1119: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1120: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1121: const PetscInt epFaceRow = -1;
1122: const PetscInt start0 = 0;
1123: const PetscInt start1 = 0;
1124: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1125: const PetscInt startGhost0 = ghostOffsetStart[0];
1126: const PetscInt startGhost1 = stag->nGhost[1]-ghostOffsetEnd[1];
1127: const PetscInt startGhost2 = 0;
1128: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1129: const PetscInt endGhost1 = stag->nGhost[1];
1130: const PetscInt endGhost2 = ghostOffsetStart[2];
1131: const PetscBool extra0 = dummyEnd[0];
1132: const PetscBool extra1 = PETSC_FALSE;
1133: const PetscBool extra2 = PETSC_FALSE;
1134: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1135: }
1137: /* RIGHT UP BACK */
1138: if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyStart[2]) {
1139: const PetscInt neighbor = 8;
1140: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1141: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1142: const PetscInt epFaceRow = -1;
1143: const PetscInt start0 = 0;
1144: const PetscInt start1 = 0;
1145: const PetscInt start2 = nNeighbors[neighbor][2] - ghostOffsetStart[2];
1146: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1147: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1148: const PetscInt startGhost2 = 0;
1149: const PetscInt endGhost0 = stag->nGhost[0];
1150: const PetscInt endGhost1 = stag->nGhost[1];
1151: const PetscInt endGhost2 = ghostOffsetStart[2];
1152: const PetscBool extra0 = PETSC_FALSE;
1153: const PetscBool extra1 = PETSC_FALSE;
1154: const PetscBool extra2 = PETSC_FALSE;
1155: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1156: }
1158: /* LEFT DOWN */
1159: if (!star && !dummyStart[0] && !dummyStart[1]) {
1160: const PetscInt neighbor = 9;
1161: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1162: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1163: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
1164: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1165: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1166: const PetscInt start2 = 0;
1167: const PetscInt startGhost0 = 0;
1168: const PetscInt startGhost1 = 0;
1169: const PetscInt startGhost2 = ghostOffsetStart[2];
1170: const PetscInt endGhost0 = ghostOffsetStart[0];
1171: const PetscInt endGhost1 = ghostOffsetStart[1];
1172: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1173: const PetscBool extra0 = PETSC_FALSE;
1174: const PetscBool extra1 = PETSC_FALSE;
1175: const PetscBool extra2 = dummyEnd[2];
1176: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1177: }
1179: /* DOWN */
1180: if (!dummyStart[1]) {
1181: const PetscInt neighbor = 10;
1182: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1183: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1184: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1185: const PetscInt start0 = 0;
1186: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1187: const PetscInt start2 = 0;
1188: const PetscInt startGhost0 = ghostOffsetStart[0];
1189: const PetscInt startGhost1 = 0;
1190: const PetscInt startGhost2 = ghostOffsetStart[2];
1191: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1192: const PetscInt endGhost1 = ghostOffsetStart[1];
1193: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1194: const PetscBool extra0 = dummyEnd[0];
1195: const PetscBool extra1 = PETSC_FALSE;
1196: const PetscBool extra2 = dummyEnd[2];
1197: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1198: }
1200: /* RIGHT DOWN */
1201: if (!star && !dummyEnd[0] && !dummyStart[1]) {
1202: const PetscInt neighbor = 11;
1203: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1204: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1205: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1206: const PetscInt start0 = 0;
1207: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1208: const PetscInt start2 = 0;
1209: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1210: const PetscInt startGhost1 = 0;
1211: const PetscInt startGhost2 = ghostOffsetStart[2];
1212: const PetscInt endGhost0 = stag->nGhost[0];
1213: const PetscInt endGhost1 = ghostOffsetStart[1];
1214: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1215: const PetscBool extra0 = PETSC_FALSE;
1216: const PetscBool extra1 = PETSC_FALSE;
1217: const PetscBool extra2 = dummyEnd[2];
1218: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1219: }
1221: /* LEFT */
1222: if (!dummyStart[0]) {
1223: const PetscInt neighbor = 12;
1224: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1225: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1226: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
1227: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1228: const PetscInt start1 = 0;
1229: const PetscInt start2 = 0;
1230: const PetscInt startGhost0 = 0;
1231: const PetscInt startGhost1 = ghostOffsetStart[1];
1232: const PetscInt startGhost2 = ghostOffsetStart[2];
1233: const PetscInt endGhost0 = ghostOffsetStart[0];
1234: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1235: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1236: const PetscBool extra0 = PETSC_FALSE;
1237: const PetscBool extra1 = dummyEnd[1];
1238: const PetscBool extra2 = dummyEnd[2];
1239: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1240: }
1242: /* (HERE) */
1243: {
1244: const PetscInt neighbor = 13;
1245: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1246: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1247: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
1248: const PetscInt start0 = 0;
1249: const PetscInt start1 = 0;
1250: const PetscInt start2 = 0;
1251: const PetscInt startGhost0 = ghostOffsetStart[0];
1252: const PetscInt startGhost1 = ghostOffsetStart[1];
1253: const PetscInt startGhost2 = ghostOffsetStart[2];
1254: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1255: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1256: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1257: const PetscBool extra0 = dummyEnd[0];
1258: const PetscBool extra1 = dummyEnd[1];
1259: const PetscBool extra2 = dummyEnd[2];
1260: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1261: }
1263: /* RIGHT */
1264: if (!dummyEnd[0]) {
1265: const PetscInt neighbor = 14;
1266: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1267: 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 */
1268: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1269: const PetscInt start0 = 0;
1270: const PetscInt start1 = 0;
1271: const PetscInt start2 = 0;
1272: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1273: const PetscInt startGhost1 = ghostOffsetStart[1];
1274: const PetscInt startGhost2 = ghostOffsetStart[2];
1275: const PetscInt endGhost0 = stag->nGhost[0];
1276: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1277: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1278: const PetscBool extra0 = PETSC_FALSE;
1279: const PetscBool extra1 = dummyEnd[1];
1280: const PetscBool extra2 = dummyEnd[2];
1281: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1282: }
1284: /* LEFT UP */
1285: if (!star && !dummyStart[0] && !dummyEnd[1]) {
1286: const PetscInt neighbor = 15;
1287: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1288: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1289: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
1290: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1291: const PetscInt start1 = 0;
1292: const PetscInt start2 = 0;
1293: const PetscInt startGhost0 = 0;
1294: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1295: const PetscInt startGhost2 = ghostOffsetStart[2];
1296: const PetscInt endGhost0 = ghostOffsetStart[0];
1297: const PetscInt endGhost1 = stag->nGhost[1];
1298: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1299: const PetscBool extra0 = PETSC_FALSE;
1300: const PetscBool extra1 = PETSC_FALSE;
1301: const PetscBool extra2 = dummyEnd[2];
1302: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1303: }
1305: /* UP */
1306: if (!dummyEnd[1]) {
1307: const PetscInt neighbor = 16;
1308: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1309: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1310: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
1311: const PetscInt start0 = 0;
1312: const PetscInt start1 = 0;
1313: const PetscInt start2 = 0;
1314: const PetscInt startGhost0 = ghostOffsetStart[0];
1315: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1316: const PetscInt startGhost2 = ghostOffsetStart[2];
1317: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1318: const PetscInt endGhost1 = stag->nGhost[1];
1319: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1320: const PetscBool extra0 = dummyEnd[0];
1321: const PetscBool extra1 = PETSC_FALSE;
1322: const PetscBool extra2 = dummyEnd[2];
1323: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1324: }
1326: /* RIGHT UP */
1327: if (!star && !dummyEnd[0] && !dummyEnd[1]) {
1328: const PetscInt neighbor = 17;
1329: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1330: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1331: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
1332: const PetscInt start0 = 0;
1333: const PetscInt start1 = 0;
1334: const PetscInt start2 = 0;
1335: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1336: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1337: const PetscInt startGhost2 = ghostOffsetStart[2];
1338: const PetscInt endGhost0 = stag->nGhost[0];
1339: const PetscInt endGhost1 = stag->nGhost[1];
1340: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1341: const PetscBool extra0 = PETSC_FALSE;
1342: const PetscBool extra1 = PETSC_FALSE;
1343: const PetscBool extra2 = dummyEnd[2];
1344: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1345: }
1347: /* LEFT DOWN FRONT */
1348: if (!star && !dummyStart[0] && !dummyStart[1] && !dummyEnd[2]) {
1349: const PetscInt neighbor = 18;
1350: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1351: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1352: const PetscInt epFaceRow = -1;
1353: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1354: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1355: const PetscInt start2 = 0;
1356: const PetscInt startGhost0 = 0;
1357: const PetscInt startGhost1 = 0;
1358: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1359: const PetscInt endGhost0 = ghostOffsetStart[0];
1360: const PetscInt endGhost1 = ghostOffsetStart[1];
1361: const PetscInt endGhost2 = stag->nGhost[2];
1362: const PetscBool extra0 = PETSC_FALSE;
1363: const PetscBool extra1 = PETSC_FALSE;
1364: const PetscBool extra2 = PETSC_FALSE;
1365: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1366: }
1368: /* DOWN FRONT */
1369: if (!star && !dummyStart[1] && !dummyEnd[2]) {
1370: const PetscInt neighbor = 19;
1371: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1372: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1373: const PetscInt epFaceRow = -1;
1374: const PetscInt start0 = 0;
1375: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1376: const PetscInt start2 = 0;
1377: const PetscInt startGhost0 = ghostOffsetStart[0];
1378: const PetscInt startGhost1 = 0;
1379: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1380: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1381: const PetscInt endGhost1 = ghostOffsetStart[1];
1382: const PetscInt endGhost2 = stag->nGhost[2];
1383: const PetscBool extra0 = dummyEnd[0];
1384: const PetscBool extra1 = PETSC_FALSE;
1385: const PetscBool extra2 = PETSC_FALSE;
1386: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1387: }
1389: /* RIGHT DOWN FRONT */
1390: if (!star && !dummyEnd[0] && !dummyStart[1] && !dummyEnd[2]) {
1391: const PetscInt neighbor = 20;
1392: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1393: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1];
1394: const PetscInt epFaceRow = -1;
1395: const PetscInt start0 = 0;
1396: const PetscInt start1 = nNeighbors[neighbor][1] - ghostOffsetStart[1];
1397: const PetscInt start2 = 0;
1398: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1399: const PetscInt startGhost1 = 0;
1400: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1401: const PetscInt endGhost0 = stag->nGhost[0];
1402: const PetscInt endGhost1 = ghostOffsetStart[1];
1403: const PetscInt endGhost2 = stag->nGhost[2];
1404: const PetscBool extra0 = PETSC_FALSE;
1405: const PetscBool extra1 = PETSC_FALSE;
1406: const PetscBool extra2 = PETSC_FALSE;
1407: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1408: }
1410: /* LEFT FRONT */
1411: if (!star && !dummyStart[0] && !dummyEnd[2]) {
1412: const PetscInt neighbor = 21;
1413: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1414: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1415: const PetscInt epFaceRow = -1;
1416: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1417: const PetscInt start1 = 0;
1418: const PetscInt start2 = 0;
1419: const PetscInt startGhost0 = 0;
1420: const PetscInt startGhost1 = ghostOffsetStart[1];
1421: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1422: const PetscInt endGhost0 = ghostOffsetStart[0];
1423: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1424: const PetscInt endGhost2 = stag->nGhost[2];
1425: const PetscBool extra0 = PETSC_FALSE;
1426: const PetscBool extra1 = dummyEnd[1];
1427: const PetscBool extra2 = PETSC_FALSE;
1428: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1429: }
1431: /* FRONT */
1432: if (!dummyEnd[2]) {
1433: const PetscInt neighbor = 22;
1434: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1435: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (dummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1436: const PetscInt epFaceRow = -1;
1437: const PetscInt start0 = 0;
1438: const PetscInt start1 = 0;
1439: const PetscInt start2 = 0;
1440: const PetscInt startGhost0 = ghostOffsetStart[0];
1441: const PetscInt startGhost1 = ghostOffsetStart[1];
1442: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1443: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1444: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1445: const PetscInt endGhost2 = stag->nGhost[2];
1446: const PetscBool extra0 = dummyEnd[0];
1447: const PetscBool extra1 = dummyEnd[1];
1448: const PetscBool extra2 = PETSC_FALSE;
1449: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1450: }
1452: /* RIGHT FRONT */
1453: if (!star && !dummyEnd[0] && !dummyEnd[2]) {
1454: const PetscInt neighbor = 23;
1455: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1456: 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 */
1457: const PetscInt epFaceRow = -1;
1458: const PetscInt start0 = 0;
1459: const PetscInt start1 = 0;
1460: const PetscInt start2 = 0;
1461: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1462: const PetscInt startGhost1 = ghostOffsetStart[1];
1463: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1464: const PetscInt endGhost0 = stag->nGhost[0];
1465: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1466: const PetscInt endGhost2 = stag->nGhost[2];
1467: const PetscBool extra0 = PETSC_FALSE;
1468: const PetscBool extra1 = dummyEnd[1];
1469: const PetscBool extra2 = PETSC_FALSE;
1470: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1471: }
1473: /* LEFT UP FRONT */
1474: if (!star && !dummyStart[0] && !dummyEnd[1] && !dummyEnd[2]) {
1475: const PetscInt neighbor = 24;
1476: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0];
1477: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1478: const PetscInt epFaceRow = -1;
1479: const PetscInt start0 = nNeighbors[neighbor][0] - ghostOffsetStart[0];
1480: const PetscInt start1 = 0;
1481: const PetscInt start2 = 0;
1482: const PetscInt startGhost0 = 0;
1483: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1484: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1485: const PetscInt endGhost0 = ghostOffsetStart[0];
1486: const PetscInt endGhost1 = stag->nGhost[1];
1487: const PetscInt endGhost2 = stag->nGhost[2];
1488: const PetscBool extra0 = PETSC_FALSE;
1489: const PetscBool extra1 = PETSC_FALSE;
1490: const PetscBool extra2 = PETSC_FALSE;
1491: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1492: }
1494: /* UP FRONT */
1495: if (!star && !dummyEnd[1] && !dummyEnd[2]) {
1496: const PetscInt neighbor = 25;
1497: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1498: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1499: const PetscInt epFaceRow = -1;
1500: const PetscInt start0 = 0;
1501: const PetscInt start1 = 0;
1502: const PetscInt start2 = 0;
1503: const PetscInt startGhost0 = ghostOffsetStart[0];
1504: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1505: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1506: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1507: const PetscInt endGhost1 = stag->nGhost[1];
1508: const PetscInt endGhost2 = stag->nGhost[2];
1509: const PetscBool extra0 = dummyEnd[0];
1510: const PetscBool extra1 = PETSC_FALSE;
1511: const PetscBool extra2 = PETSC_FALSE;
1512: DMStagSetUpBuildScatterPopulateIdx_3d(stag,&count,idxLocal,idxGlobal,entriesPerEdge,entriesPerFace,eprNeighbor,eplNeighbor,eprGhost,eplGhost,epFaceRow,globalOffsets[stag->neighbors[neighbor]],start0,start1,start2,startGhost0,startGhost1,startGhost2,endGhost0,endGhost1,endGhost2,extra0,extra1,extra2);
1513: }
1515: /* RIGHT UP FRONT */
1516: if (!star && !dummyEnd[0] && !dummyEnd[1] && !dummyEnd[2]) {
1517: const PetscInt neighbor = 26;
1518: const PetscInt eprNeighbor = stag->entriesPerElement * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1519: const PetscInt eplNeighbor = eprNeighbor * nNeighbors[neighbor][1] + (nextToDummyEnd[1] ? nNeighbors[neighbor][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1520: const PetscInt epFaceRow = -1;
1521: const PetscInt start0 = 0;
1522: const PetscInt start1 = 0;
1523: const PetscInt start2 = 0;
1524: const PetscInt startGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
1525: const PetscInt startGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
1526: const PetscInt startGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
1527: const PetscInt endGhost0 = stag->nGhost[0];
1528: const PetscInt endGhost1 = stag->nGhost[1];
1529: const PetscInt endGhost2 = stag->nGhost[2];
1530: const PetscBool extra0 = PETSC_FALSE;
1531: const PetscBool extra1 = PETSC_FALSE;
1532: const PetscBool extra2 = PETSC_FALSE;
1533: 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);
1534: }
1536: if (count != entriesToTransferTotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries computed in gtol (%d) is not as expected (%d)",count,entriesToTransferTotal);
1538: /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
1539: {
1540: Vec local,global;
1541: IS isLocal,isGlobal;
1542: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
1543: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
1544: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
1545: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
1546: VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
1547: VecDestroy(&global);
1548: VecDestroy(&local);
1549: ISDestroy(&isLocal); /* frees idxLocal */
1550: ISDestroy(&isGlobal); /* free idxGlobal */
1551: }
1553: return(0);
1554: }
1556: /* Note: this function assumes that DMBoundary types of none, ghosted, and periodic are the only ones of interest.
1557: Adding support for others should be done very carefully. */
1558: static PetscErrorCode DMStagSetUpBuildL2G_3d(DM dm,const PetscInt *globalOffsets)
1559: {
1560: PetscErrorCode ierr;
1561: const DM_Stag * const stag = (DM_Stag*)dm->data;
1562: PetscInt *idxGlobalAll;
1563: PetscInt d,count,ighost,jghost,kghost,ghostOffsetStart[3],ghostOffsetEnd[3],entriesPerFace,entriesPerEdge;
1564: PetscInt nNeighbors[27][3],eprNeighbor[27],eplNeighbor[27],globalOffsetNeighbor[27];
1565: PetscBool nextToDummyEnd[3],dummyStart[3],dummyEnd[3],star;
1569: /* Check stencil type */
1570: 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]);
1571: 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");
1572: star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR);
1574: /* Convenience variables */
1575: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
1576: entriesPerEdge = stag->dof[0] + stag->dof[1];
1577: for (d=0;d<3;++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d]-2);
1579: /* Ghost offsets (may not be the ghost width, since we always have at least one ghost element on the right/top/front) */
1580: for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
1581: for (d=0; d<3; ++d) ghostOffsetEnd[d] = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
1583: /* Whether the ghost region includes dummies. Currently equivalent to being a non-periodic boundary. */
1584: for (d=0; d<3; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1585: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1587: /* Compute numbers of elements on each neighbor and offset*/
1588: {
1589: PetscInt i;
1590: for (i=0; i<27; ++i) {
1591: const PetscInt neighborRank = stag->neighbors[i];
1592: if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 13) */
1593: nNeighbors[i][0] = stag->l[0][neighborRank % stag->nRanks[0]];
1594: nNeighbors[i][1] = stag->l[1][neighborRank % (stag->nRanks[0] * stag->nRanks[1]) / stag->nRanks[0]];
1595: nNeighbors[i][2] = stag->l[2][neighborRank / (stag->nRanks[0] * stag->nRanks[1])];
1596: globalOffsetNeighbor[i] = globalOffsets[stag->neighbors[i]];
1597: } /* else leave uninitialized - error if accessed */
1598: }
1599: }
1601: /* Precompute elements per row and layer on neighbor (zero unused) */
1602: PetscMemzero(eprNeighbor,sizeof(eprNeighbor));
1603: PetscMemzero(eplNeighbor,sizeof(eplNeighbor));
1604: if (stag->neighbors[0] >= 0) {
1605: eprNeighbor[0] = stag->entriesPerElement * nNeighbors[0][0];
1606: eplNeighbor[0] = eprNeighbor[0] * nNeighbors[0][1];
1607: }
1608: if (stag->neighbors[1] >= 0) {
1609: eprNeighbor[1] = stag->entriesPerElement * nNeighbors[1][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1610: eplNeighbor[1] = eprNeighbor[1] * nNeighbors[1][1];
1611: }
1612: if (stag->neighbors[2] >= 0) {
1613: eprNeighbor[2] = stag->entriesPerElement * nNeighbors[2][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1614: eplNeighbor[2] = eprNeighbor[2] * nNeighbors[2][1];
1615: }
1616: if (stag->neighbors[3] >= 0) {
1617: eprNeighbor[3] = stag->entriesPerElement * nNeighbors[3][0];
1618: eplNeighbor[3] = eprNeighbor[3] * nNeighbors[3][1] + (dummyEnd[1] ? nNeighbors[3][0] * entriesPerFace : 0); /* May be a top boundary */
1619: }
1620: if (stag->neighbors[4] >= 0) {
1621: eprNeighbor[4] = stag->entriesPerElement * nNeighbors[4][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1622: 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 */
1623: }
1624: if (stag->neighbors[5] >= 0) {
1625: eprNeighbor[5] = stag->entriesPerElement * nNeighbors[5][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1626: 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 */
1627: }
1628: if (stag->neighbors[6] >= 0) {
1629: eprNeighbor[6] = stag->entriesPerElement * nNeighbors[6][0];
1630: eplNeighbor[6] = eprNeighbor[6] * nNeighbors[6][1] + (nextToDummyEnd[1] ? nNeighbors[6][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1631: }
1632: if (stag->neighbors[7] >= 0) {
1633: eprNeighbor[7] = stag->entriesPerElement * nNeighbors[7][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1634: eplNeighbor[7] = eprNeighbor[7] * nNeighbors[7][1] + (nextToDummyEnd[1] ? nNeighbors[7][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1635: }
1636: if (stag->neighbors[8] >= 0) {
1637: eprNeighbor[8] = stag->entriesPerElement * nNeighbors[8][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1638: eplNeighbor[8] = eprNeighbor[8] * nNeighbors[8][1] + (nextToDummyEnd[1] ? nNeighbors[8][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1639: }
1640: if (stag->neighbors[9] >= 0) {
1641: eprNeighbor[9] = stag->entriesPerElement * nNeighbors[9][0];
1642: eplNeighbor[9] = eprNeighbor[9] * nNeighbors[9][1];
1643: }
1644: if (stag->neighbors[10] >= 0) {
1645: eprNeighbor[10] = stag->entriesPerElement * nNeighbors[10][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1646: eplNeighbor[10] = eprNeighbor[10] * nNeighbors[10][1];
1647: }
1648: if (stag->neighbors[11] >= 0) {
1649: eprNeighbor[11] = stag->entriesPerElement * nNeighbors[11][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1650: eplNeighbor[11] = eprNeighbor[11] * nNeighbors[11][1];
1651: }
1652: if (stag->neighbors[12] >= 0) {
1653: eprNeighbor[12] = stag->entriesPerElement * nNeighbors[12][0];
1654: eplNeighbor[12] = eprNeighbor[12] * nNeighbors[12][1] + (dummyEnd[1] ? nNeighbors[12][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1655: }
1656: if (stag->neighbors[13] >= 0) {
1657: eprNeighbor[13] = stag->entriesPerElement * nNeighbors[13][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
1658: eplNeighbor[13] = eprNeighbor[13] * nNeighbors[13][1] + (dummyEnd[1] ? nNeighbors[13][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
1659: }
1660: if (stag->neighbors[14] >= 0) {
1661: eprNeighbor[14] = stag->entriesPerElement * nNeighbors[14][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1662: 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 */
1663: }
1664: if (stag->neighbors[15] >= 0) {
1665: eprNeighbor[15] = stag->entriesPerElement * nNeighbors[15][0];
1666: eplNeighbor[15] = eprNeighbor[15] * nNeighbors[15][1] + (nextToDummyEnd[1] ? nNeighbors[15][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1667: }
1668: if (stag->neighbors[16] >= 0) {
1669: eprNeighbor[16] = stag->entriesPerElement * nNeighbors[16][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1670: eplNeighbor[16] = eprNeighbor[16] * nNeighbors[16][1] + (nextToDummyEnd[1] ? nNeighbors[16][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1671: }
1672: if (stag->neighbors[17] >= 0) {
1673: eprNeighbor[17] = stag->entriesPerElement * nNeighbors[17][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1674: eplNeighbor[17] = eprNeighbor[17] * nNeighbors[17][1] + (nextToDummyEnd[1] ? nNeighbors[17][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1675: }
1676: if (stag->neighbors[18] >= 0) {
1677: eprNeighbor[18] = stag->entriesPerElement * nNeighbors[18][0];
1678: eplNeighbor[18] = eprNeighbor[18] * nNeighbors[18][1];
1679: }
1680: if (stag->neighbors[19] >= 0) {
1681: eprNeighbor[19] = stag->entriesPerElement * nNeighbors[19][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1682: eplNeighbor[19] = eprNeighbor[19] * nNeighbors[19][1];
1683: }
1684: if (stag->neighbors[20] >= 0) {
1685: eprNeighbor[20] = stag->entriesPerElement * nNeighbors[20][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1686: eplNeighbor[20] = eprNeighbor[20] * nNeighbors[20][1];
1687: }
1688: if (stag->neighbors[21] >= 0) {
1689: eprNeighbor[21] = stag->entriesPerElement * nNeighbors[21][0];
1690: eplNeighbor[21] = eprNeighbor[21] * nNeighbors[21][1] + (dummyEnd[1] ? nNeighbors[21][0] * entriesPerFace : 0); /* We+neighbor may be a top boundary */
1691: }
1692: if (stag->neighbors[22] >= 0) {
1693: eprNeighbor[22] = stag->entriesPerElement * nNeighbors[22][0] + (dummyEnd[0] ? entriesPerFace : 0); /* neighbor is a right boundary if we are*/
1694: eplNeighbor[22] = eprNeighbor[22] * nNeighbors[22][1] + (dummyEnd[1] ? nNeighbors[22][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* May be a top boundary */
1695: }
1696: if (stag->neighbors[23] >= 0) {
1697: eprNeighbor[23] = stag->entriesPerElement * nNeighbors[23][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1698: 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 */
1699: }
1700: if (stag->neighbors[24] >= 0) {
1701: eprNeighbor[24] = stag->entriesPerElement * nNeighbors[24][0];
1702: eplNeighbor[24] = eprNeighbor[24] * nNeighbors[24][1] + (nextToDummyEnd[1] ? nNeighbors[24][0] * entriesPerFace : 0); /* Neighbor may be a top boundary */
1703: }
1704: if (stag->neighbors[25] >= 0) {
1705: eprNeighbor[25] = stag->entriesPerElement * nNeighbors[25][0] + (dummyEnd[0] ? entriesPerFace : 0); /* We+neighbor may be a right boundary */
1706: eplNeighbor[25] = eprNeighbor[25] * nNeighbors[25][1] + (nextToDummyEnd[1] ? nNeighbors[25][0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1707: }
1708: if (stag->neighbors[26] >= 0) {
1709: eprNeighbor[26] = stag->entriesPerElement * nNeighbors[26][0] + (nextToDummyEnd[0] ? entriesPerFace : 0); /* Neighbor may be a right boundary */
1710: eplNeighbor[26] = eprNeighbor[26] * nNeighbors[26][1] + (nextToDummyEnd[1] ? nNeighbors[26][0] * entriesPerFace + (nextToDummyEnd[0] ? entriesPerEdge : 0) : 0); /* Neighbor may be a top boundary */
1711: }
1713: /* Populate idxGlobalAll */
1714: PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
1715: count = 0;
1717: /* Loop over layers 1/3 : Back */
1718: if (!dummyStart[2]) {
1719: for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
1720: 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) */
1722: /* Loop over rows 1/3: Down Back*/
1723: if (!star && !dummyStart[1]) {
1724: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1725: 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)*/
1727: /* Loop over columns 1/3: Left Back Down*/
1728: if (!dummyStart[0]) {
1729: const PetscInt neighbor = 0;
1730: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1731: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1732: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1733: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1734: }
1735: }
1736: } else {
1737: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1738: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1739: idxGlobalAll[count] = -1;
1740: }
1741: }
1742: }
1744: /* Loop over columns 2/3: (Middle) Down Back */
1745: {
1746: const PetscInt neighbor = 1;
1747: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1748: const PetscInt i = ighost - ghostOffsetStart[0];
1749: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1750: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1751: }
1752: }
1753: }
1755: /* Loop over columns 3/3: Right Down Back */
1756: if (!dummyEnd[0]) {
1757: const PetscInt neighbor = 2;
1758: PetscInt i;
1759: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1760: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1761: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1762: }
1763: }
1764: } else {
1765: /* Partial dummy entries on (Middle) Down Back neighbor */
1766: const PetscInt neighbor = 1;
1767: PetscInt i,dLocal;
1768: i = stag->n[0];
1769: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1770: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1771: }
1772: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1773: idxGlobalAll[count] = -1;
1774: }
1775: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1776: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1777: }
1778: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1779: idxGlobalAll[count] = -1;
1780: }
1781: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1782: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1783: }
1784: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1785: idxGlobalAll[count] = -1;
1786: }
1787: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1788: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1789: }
1790: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1791: idxGlobalAll[count] = -1;
1792: }
1793: ++i;
1794: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1795: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1796: idxGlobalAll[count] = -1;
1797: }
1798: }
1799: }
1800: }
1801: } else {
1802: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
1803: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
1804: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1805: idxGlobalAll[count] = -1;
1806: }
1807: }
1808: }
1809: }
1811: /* Loop over rows 2/3: (Middle) Back */
1812: {
1813: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
1814: const PetscInt j = jghost - ghostOffsetStart[1];
1816: /* Loop over columns 1/3: Left (Middle) Back */
1817: if (!star && !dummyStart[0]) {
1818: const PetscInt neighbor = 3;
1819: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1820: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1821: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1822: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1823: }
1824: }
1825: } else {
1826: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1827: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1828: idxGlobalAll[count] = -1;
1829: }
1830: }
1831: }
1833: /* Loop over columns 2/3: (Middle) (Middle) Back */
1834: {
1835: const PetscInt neighbor = 4;
1836: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1837: const PetscInt i = ighost - ghostOffsetStart[0];
1838: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1839: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1840: }
1841: }
1842: }
1844: /* Loop over columns 3/3: Right (Middle) Back */
1845: if (!star && !dummyEnd[0]) {
1846: const PetscInt neighbor = 5;
1847: PetscInt i;
1848: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1849: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1850: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1851: }
1852: }
1853: } else if (dummyEnd[0]) {
1854: /* Partial dummy entries on (Middle) (Middle) Back rank */
1855: const PetscInt neighbor = 4;
1856: PetscInt i,dLocal;
1857: i = stag->n[0];
1858: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1859: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1860: }
1861: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1862: idxGlobalAll[count] = -1;
1863: }
1864: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1865: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1866: }
1867: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1868: idxGlobalAll[count] = -1;
1869: }
1870: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1871: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1872: }
1873: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1874: idxGlobalAll[count] = -1;
1875: }
1876: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1877: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1878: }
1879: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1880: idxGlobalAll[count] = -1;
1881: }
1882: ++i;
1883: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1884: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1885: idxGlobalAll[count] = -1;
1886: }
1887: }
1888: } else {
1889: /* Right (Middle) Back dummies */
1890: PetscInt i;
1891: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1892: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1893: idxGlobalAll[count] = -1;
1894: }
1895: }
1896: }
1897: }
1898: }
1900: /* Loop over rows 3/3: Up Back */
1901: if (!star && !dummyEnd[1]) {
1902: PetscInt j;
1903: for (j=0; j<ghostOffsetEnd[1]; ++j) {
1904: /* Loop over columns 1/3: Left Up Back*/
1905: if (!dummyStart[0]) {
1906: const PetscInt neighbor = 6;
1907: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1908: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1909: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1910: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1911: }
1912: }
1913: } else {
1914: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1915: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1916: idxGlobalAll[count] = -1;
1917: }
1918: }
1919: }
1921: /* Loop over columns 2/3: (Middle) Up Back*/
1922: {
1923: const PetscInt neighbor = 7;
1924: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
1925: const PetscInt i = ighost - ghostOffsetStart[0];
1926: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1927: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1928: }
1929: }
1930: }
1932: /* Loop over columns 3/3: Right Up Back*/
1933: if (!dummyEnd[0]) {
1934: const PetscInt neighbor = 8;
1935: PetscInt i;
1936: for (i=0; i<ghostOffsetEnd[0]; ++i) {
1937: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1938: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1939: }
1940: }
1941: } else {
1942: /* Partial dummies on (Middle) Up Back */
1943: const PetscInt neighbor = 7;
1944: PetscInt i,dLocal;
1945: i = stag->n[0];
1946: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
1947: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1948: }
1949: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
1950: idxGlobalAll[count] = -1;
1951: }
1952: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
1953: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1954: }
1955: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
1956: idxGlobalAll[count] = -1;
1957: }
1958: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
1959: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1960: }
1961: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
1962: idxGlobalAll[count] = -1;
1963: }
1964: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
1965: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
1966: }
1967: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
1968: idxGlobalAll[count] = -1;
1969: }
1970: ++i;
1971: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
1972: for (d=0; d<stag->entriesPerElement; ++d,++count) {
1973: idxGlobalAll[count] = -1;
1974: }
1975: }
1976: }
1977: }
1978: } else if (dummyEnd[1]) {
1979: /* Up Back partial dummy row */
1980: PetscInt j = stag->n[1];
1982: if (!star && !dummyStart[0]) {
1983: /* Up Back partial dummy row: Loop over columns 1/3: Left Up Back, on Left (Middle) Back rank */
1984: const PetscInt neighbor = 3;
1985: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
1986: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
1987: PetscInt dLocal;
1988: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
1989: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1990: }
1992: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
1993: idxGlobalAll[count] = -1;
1994: }
1995: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
1996: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
1997: }
1998: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
1999: idxGlobalAll[count] = -1;
2000: }
2001: }
2002: } else {
2003: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2004: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2005: idxGlobalAll[count] = -1;
2006: }
2007: }
2008: }
2010: /* Up Back partial dummy row: Loop over columns 2/3: (Middle) Up Back, on (Middle) (Middle) Back rank */
2011: {
2012: const PetscInt neighbor = 4;
2013: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2014: const PetscInt i = ighost - ghostOffsetStart[0];
2015: PetscInt dLocal;
2016: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2017: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2018: }
2019: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2020: idxGlobalAll[count] = -1;
2021: }
2022: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2023: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2024: }
2025: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2026: idxGlobalAll[count] = -1;
2027: }
2028: }
2029: }
2031: /* Up Back partial dummy row: Loop over columns 3/3: Right Up Back, on Right (Middle) Back rank */
2032: if (!star && !dummyEnd[0]) {
2033: const PetscInt neighbor = 5;
2034: PetscInt i;
2035: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2036: PetscInt dLocal;
2037: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2038: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2039: }
2040: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2041: idxGlobalAll[count] = -1;
2042: }
2043: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2044: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2045: }
2046: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2047: idxGlobalAll[count] = -1;
2048: }
2049: }
2050: } else if (dummyEnd[0]) {
2051: /* Up Right Back partial dummy element, on (Middle) (Middle) Back rank */
2052: const PetscInt neighbor = 4;
2053: PetscInt i,dLocal;
2054: i = stag->n[0];
2055: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2056: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2057: }
2058: 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 */
2059: idxGlobalAll[count] = -1;
2060: }
2061: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2062: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2063: }
2064: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2065: idxGlobalAll[count] = -1;
2066: }
2067: ++i;
2068: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2069: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2070: idxGlobalAll[count] = -1;
2071: }
2072: }
2073: } else {
2074: /* Up Right Back dummies */
2075: PetscInt i;
2076: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2077: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2078: idxGlobalAll[count] = -1;
2079: }
2080: }
2081: }
2082: ++j;
2083: /* Up Back additional dummy rows */
2084: for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2085: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2086: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2087: idxGlobalAll[count] = -1;
2088: }
2089: }
2090: }
2091: } else {
2092: /* Up Back dummy rows */
2093: PetscInt j;
2094: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2095: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2096: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2097: idxGlobalAll[count] = -1;
2098: }
2099: }
2100: }
2101: }
2102: }
2103: } else {
2104: for (kghost = 0; kghost<ghostOffsetStart[2]; ++kghost) {
2105: for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
2106: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2107: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2108: idxGlobalAll[count] = -1;
2109: }
2110: }
2111: }
2112: }
2113: }
2115: /* Loop over layers 2/3 : (Middle) */
2116: {
2117: for (kghost = ghostOffsetStart[2]; kghost<stag->nGhost[2]-ghostOffsetEnd[2]; ++kghost) {
2118: const PetscInt k = kghost - ghostOffsetStart[2];
2120: /* Loop over rows 1/3: Down (Middle) */
2121: if (!dummyStart[1]) {
2122: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2123: const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* down neighbor (10) always exists */
2125: /* Loop over columns 1/3: Left Down (Middle) */
2126: if (!star && !dummyStart[0]) {
2127: const PetscInt neighbor = 9;
2128: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2129: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2130: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2131: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2132: }
2133: }
2134: } else {
2135: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2136: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2137: idxGlobalAll[count] = -1;
2138: }
2139: }
2140: }
2142: /* Loop over columns 2/3: (Middle) Down (Middle) */
2143: {
2144: const PetscInt neighbor = 10;
2145: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2146: const PetscInt i = ighost - ghostOffsetStart[0];
2147: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2148: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2149: }
2150: }
2151: }
2153: /* Loop over columns 3/3: Right Down (Middle) */
2154: if (!star && !dummyEnd[0]) {
2155: const PetscInt neighbor = 11;
2156: PetscInt i;
2157: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2158: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2159: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2160: }
2161: }
2162: } else if (dummyEnd[0]) {
2163: /* Partial dummies on (Middle) Down (Middle) neighbor */
2164: const PetscInt neighbor = 10;
2165: PetscInt i,dLocal;
2166: i = stag->n[0];
2167: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2168: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2169: }
2170: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2171: idxGlobalAll[count] = -1;
2172: }
2173: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2174: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2175: }
2176: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2177: idxGlobalAll[count] = -1;
2178: }
2179: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2180: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2181: }
2182: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2183: idxGlobalAll[count] = -1;
2184: }
2185: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2186: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2187: }
2188: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2189: idxGlobalAll[count] = -1;
2190: }
2191: ++i;
2192: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2193: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2194: idxGlobalAll[count] = -1;
2195: }
2196: }
2197: } else {
2198: /* Right Down (Middle) dummies */
2199: PetscInt i;
2200: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2201: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2202: idxGlobalAll[count] = -1;
2203: }
2204: }
2205: }
2206: }
2207: } else {
2208: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2209: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2210: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2211: idxGlobalAll[count] = -1;
2212: }
2213: }
2214: }
2215: }
2217: /* Loop over rows 2/3: (Middle) (Middle) */
2218: {
2219: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2220: const PetscInt j = jghost - ghostOffsetStart[1];
2222: /* Loop over columns 1/3: Left (Middle) (Middle) */
2223: if (!dummyStart[0]) {
2224: const PetscInt neighbor = 12;
2225: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2226: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2227: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2228: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2229: }
2230: }
2231: } else {
2232: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2233: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2234: idxGlobalAll[count] = -1;
2235: }
2236: }
2237: }
2239: /* Loop over columns 2/3: (Middle) (Middle) (Middle) aka (Here) */
2240: {
2241: const PetscInt neighbor = 13;
2242: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2243: const PetscInt i = ighost - ghostOffsetStart[0];
2244: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2245: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2246: }
2247: }
2248: }
2250: /* Loop over columns 3/3: Right (Middle) (Middle) */
2251: if (!dummyEnd[0]) {
2252: const PetscInt neighbor = 14;
2253: PetscInt i;
2254: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2255: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2256: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2257: }
2258: }
2259: } else {
2260: /* Partial dummies on this rank */
2261: const PetscInt neighbor = 13;
2262: PetscInt i,dLocal;
2263: i = stag->n[0];
2264: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2265: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2266: }
2267: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2268: idxGlobalAll[count] = -1;
2269: }
2270: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2271: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor]+ j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2272: }
2273: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2274: idxGlobalAll[count] = -1;
2275: }
2276: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2277: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2278: }
2279: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2280: idxGlobalAll[count] = -1;
2281: }
2282: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2283: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2284: }
2285: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2286: idxGlobalAll[count] = -1;
2287: }
2288: ++i;
2289: for (;i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2290: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2291: idxGlobalAll[count] = -1;
2292: }
2293: }
2294: }
2295: }
2296: }
2298: /* Loop over rows 3/3: Up (Middle) */
2299: if (!dummyEnd[1]) {
2300: PetscInt j;
2301: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2303: /* Loop over columns 1/3: Left Up (Middle) */
2304: if (!star && !dummyStart[0]) {
2305: const PetscInt neighbor = 15;
2306: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2307: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2308: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2309: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2310: }
2311: }
2312: } else {
2313: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2314: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2315: idxGlobalAll[count] = -1;
2316: }
2317: }
2318: }
2320: /* Loop over columns 2/3: (Middle) Up (Middle) */
2321: {
2322: const PetscInt neighbor = 16;
2323: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2324: const PetscInt i = ighost - ghostOffsetStart[0];
2325: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2326: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2327: }
2328: }
2329: }
2331: /* Loop over columns 3/3: Right Up (Middle) */
2332: if (!star && !dummyEnd[0]) {
2333: const PetscInt neighbor = 17;
2334: PetscInt i;
2335: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2336: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2337: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2338: }
2339: }
2340: } else if (dummyEnd[0]) {
2341: /* Partial dummies on (Middle) Up (Middle) rank */
2342: const PetscInt neighbor = 16;
2343: PetscInt i,dLocal;
2344: i = stag->n[0];
2345: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2346: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2347: }
2348: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2349: idxGlobalAll[count] = -1;
2350: }
2351: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2352: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2353: }
2354: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2355: idxGlobalAll[count] = -1;
2356: }
2357: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2358: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2359: }
2360: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2361: idxGlobalAll[count] = -1;
2362: }
2363: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2364: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2365: }
2366: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2367: idxGlobalAll[count] = -1;
2368: }
2369: ++i;
2370: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2371: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2372: idxGlobalAll[count] = -1;
2373: }
2374: }
2375: } else {
2376: /* Right Up (Middle) dumies */
2377: PetscInt i;
2378: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2379: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2380: idxGlobalAll[count] = -1;
2381: }
2382: }
2383: }
2384: }
2385: } else {
2386: /* Up (Middle) partial dummy row */
2387: PetscInt j = stag->n[1];
2389: /* Up (Middle) partial dummy row: colums 1/3: Left Up (Middle), on Left (Middle) (Middle) rank */
2390: if (!dummyStart[0]) {
2391: const PetscInt neighbor = 12;
2392: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2393: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2394: PetscInt dLocal;
2395: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2396: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2397: }
2398: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2399: idxGlobalAll[count] = -1;
2400: }
2401: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2402: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2403: }
2404: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2405: idxGlobalAll[count] = -1;
2406: }
2407: }
2408: } else {
2409: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2410: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2411: idxGlobalAll[count] = -1;
2412: }
2413: }
2414: }
2416: /* Up (Middle) partial dummy row: columns 2/3: (Middle) Up (Middle), on (Middle) (Middle) (Middle) = here rank */
2417: {
2418: const PetscInt neighbor = 13;
2419: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2420: const PetscInt i = ighost - ghostOffsetStart[0];
2421: PetscInt dLocal;
2422: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2423: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2424: }
2425: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2426: idxGlobalAll[count] = -1;
2427: }
2428: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2429: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2430: }
2431: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2432: idxGlobalAll[count] = -1;
2433: }
2434: }
2435: }
2437: if (!dummyEnd[0]) {
2438: /* Up (Middle) partial dummy row: columns 3/3: Right Up (Middle), on Right (Middle) (Middle) rank */
2439: const PetscInt neighbor = 14;
2440: PetscInt i;
2441: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2442: PetscInt dLocal;
2443: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2444: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2445: }
2446: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2447: idxGlobalAll[count] = -1;
2448: }
2449: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2450: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2451: }
2452: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2453: idxGlobalAll[count] = -1;
2454: }
2455: }
2456: } else {
2457: /* Up (Middle) Right partial dummy element, on (Middle) (Middle) (Middle) = here rank */
2458: const PetscInt neighbor = 13;
2459: PetscInt i,dLocal;
2460: i = stag->n[0];
2461: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2462: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2463: }
2464: 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 */
2465: idxGlobalAll[count] = -1;
2466: }
2467: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2468: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2469: }
2470: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2471: idxGlobalAll[count] = -1;
2472: }
2473: ++i;
2474: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2475: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2476: idxGlobalAll[count] = -1;
2477: }
2478: }
2479: }
2480: ++j;
2481: /* Up (Middle) additional dummy rows */
2482: for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2483: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2484: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2485: idxGlobalAll[count] = -1;
2486: }
2487: }
2488: }
2489: }
2490: }
2491: }
2493: /* Loop over layers 3/3 : Front */
2494: if (!dummyEnd[2]) {
2495: PetscInt k;
2496: for (k=0; k<ghostOffsetEnd[2]; ++k) {
2498: /* Loop over rows 1/3: Down Front */
2499: if (!star && !dummyStart[1]) {
2500: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2501: const PetscInt j = nNeighbors[19][1] - ghostOffsetStart[1] + jghost; /* Constant for whole row (use down front neighbor) */
2503: /* Loop over columns 1/3: Left Down Front */
2504: if (!dummyStart[0]) {
2505: const PetscInt neighbor = 18;
2506: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2507: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2508: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2509: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2510: }
2511: }
2512: } else {
2513: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2514: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2515: idxGlobalAll[count] = -1;
2516: }
2517: }
2518: }
2520: /* Loop over columns 2/3: (Middle) Down Front */
2521: {
2522: const PetscInt neighbor = 19;
2523: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2524: const PetscInt i = ighost - ghostOffsetStart[0];
2525: for (d=0; d<stag->entriesPerElement; ++d,++count) { /* vertex, 2 edges, and face associated with back face */
2526: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2527: }
2528: }
2529: }
2531: /* Loop over columns 3/3: Right Down Front */
2532: if (!dummyEnd[0]) {
2533: const PetscInt neighbor = 20;
2534: PetscInt i;
2535: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2536: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2537: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2538: }
2539: }
2540: } else {
2541: /* Partial dummy element on (Middle) Down Front rank */
2542: const PetscInt neighbor = 19;
2543: PetscInt i,dLocal;
2544: i = stag->n[0];
2545: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2546: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2547: }
2548: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2549: idxGlobalAll[count] = -1;
2550: }
2551: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2552: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2553: }
2554: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2555: idxGlobalAll[count] = -1;
2556: }
2557: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2558: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2559: }
2560: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2561: idxGlobalAll[count] = -1;
2562: }
2563: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2564: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2565: }
2566: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2567: idxGlobalAll[count] = -1;
2568: }
2569: ++i;
2570: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2571: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2572: idxGlobalAll[count] = -1;
2573: }
2574: }
2575: }
2576: }
2577: } else {
2578: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2579: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2580: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2581: idxGlobalAll[count] = -1;
2582: }
2583: }
2584: }
2585: }
2587: /* Loop over rows 2/3: (Middle) Front */
2588: {
2589: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2590: const PetscInt j = jghost - ghostOffsetStart[1];
2592: /* Loop over columns 1/3: Left (Middle) Front*/
2593: if (!star && !dummyStart[0]) {
2594: const PetscInt neighbor = 21;
2595: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2596: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2597: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2598: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2599: }
2600: }
2601: } else {
2602: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2603: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2604: idxGlobalAll[count] = -1;
2605: }
2606: }
2607: }
2609: /* Loop over columns 2/3: (Middle) (Middle) Front*/
2610: {
2611: const PetscInt neighbor = 22;
2612: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2613: const PetscInt i = ighost - ghostOffsetStart[0];
2614: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2615: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2616: }
2617: }
2618: }
2620: /* Loop over columns 3/3: Right (Middle) Front*/
2621: if (!star && !dummyEnd[0]) {
2622: const PetscInt neighbor = 23;
2623: PetscInt i;
2624: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2625: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2626: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2627: }
2628: }
2629: } else if (dummyEnd[0]) {
2630: /* Partial dummy element on (Middle) (Middle) Front element */
2631: const PetscInt neighbor = 22;
2632: PetscInt i,dLocal;
2633: i = stag->n[0];
2634: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2635: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2636: }
2637: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2638: idxGlobalAll[count] = -1;
2639: }
2640: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2641: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2642: }
2643: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2644: idxGlobalAll[count] = -1;
2645: }
2646: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2647: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2648: }
2649: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2650: idxGlobalAll[count] = -1;
2651: }
2652: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2653: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2654: }
2655: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2656: idxGlobalAll[count] = -1;
2657: }
2658: ++i;
2659: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2660: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2661: idxGlobalAll[count] = -1;
2662: }
2663: }
2664: } else {
2665: /* Right (Middle) Front dummies */
2666: PetscInt i;
2667: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2668: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2669: idxGlobalAll[count] = -1;
2670: }
2671: }
2672: }
2673: }
2674: }
2676: /* Loop over rows 3/3: Up Front */
2677: if (!star && !dummyEnd[1]) {
2678: PetscInt j;
2679: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2681: /* Loop over columns 1/3: Left Up Front */
2682: if (!dummyStart[0]) {
2683: const PetscInt neighbor = 24;
2684: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2685: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2686: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2687: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2688: }
2689: }
2690: } else {
2691: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2692: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2693: idxGlobalAll[count] = -1;
2694: }
2695: }
2696: }
2698: /* Loop over columns 2/3: (Middle) Up Front */
2699: {
2700: const PetscInt neighbor = 25;
2701: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2702: const PetscInt i = ighost - ghostOffsetStart[0];
2703: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2704: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2705: }
2706: }
2707: }
2709: /* Loop over columns 3/3: Right Up Front */
2710: if (!dummyEnd[0]) {
2711: const PetscInt neighbor = 26;
2712: PetscInt i;
2713: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2714: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2715: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2716: }
2717: }
2718: } else {
2719: /* Partial dummy element on (Middle) Up Front rank */
2720: const PetscInt neighbor = 25;
2721: PetscInt i,dLocal;
2722: i = stag->n[0];
2723: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2724: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2725: }
2726: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2727: idxGlobalAll[count] = -1;
2728: }
2729: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* Back left edge */
2730: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2731: }
2732: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face */
2733: idxGlobalAll[count] = -1;
2734: }
2735: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2736: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2737: }
2738: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++dLocal, ++count) { /* Dummies on down face */
2739: idxGlobalAll[count] = -1;
2740: }
2741: for (; dLocal<stag->dof[0]+3*stag->dof[1]+3*stag->dof[2]; ++d, ++dLocal, ++count) { /* Left face */
2742: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*stag->entriesPerElement + d;
2743: }
2744: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on element */
2745: idxGlobalAll[count] = -1;
2746: }
2747: ++i;
2748: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2749: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2750: idxGlobalAll[count] = -1;
2751: }
2752: }
2753: }
2754: }
2755: } else if (dummyEnd[1]) {
2756: /* Up Front partial dummy row */
2757: PetscInt j = stag->n[1];
2759: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) Front rank */
2760: if (!star && !dummyStart[0]) {
2761: const PetscInt neighbor = 21;
2762: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2763: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2764: PetscInt dLocal;
2765: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2766: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2767: }
2768: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2769: idxGlobalAll[count] = -1;
2770: }
2771: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2772: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2773: }
2774: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2775: idxGlobalAll[count] = -1;
2776: }
2777: }
2778: } else {
2779: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2780: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2781: idxGlobalAll[count] = -1;
2782: }
2783: }
2784: }
2786: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) Front rank */
2787: {
2788: const PetscInt neighbor = 22;
2789: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2790: const PetscInt i = ighost - ghostOffsetStart[0];
2791: PetscInt dLocal;
2792: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2793: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2794: }
2795: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2796: idxGlobalAll[count] = -1;
2797: }
2798: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2799: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2800: }
2801: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2802: idxGlobalAll[count] = -1;
2803: }
2804: }
2805: }
2807: if (!star && !dummyEnd[0]) {
2808: /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) Front rank */
2809: const PetscInt neighbor = 23;
2810: PetscInt i;
2811: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2812: PetscInt dLocal;
2813: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and back down edge */
2814: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2815: }
2816: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back face and back left edge */
2817: idxGlobalAll[count] = -1;
2818: }
2819: for (; dLocal<stag->dof[0]+3*stag->dof[1]+2*stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge and down face */
2820: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2821: }
2822: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies left face and element */
2823: idxGlobalAll[count] = -1;
2824: }
2825: }
2827: } else if (dummyEnd[0]) {
2828: /* Partial Right Up Front dummy element, on (Middle) (Middle) Front rank */
2829: const PetscInt neighbor = 22;
2830: PetscInt i,dLocal;
2831: i = stag->n[0];
2832: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2833: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2834: }
2835: for (; dLocal<stag->dof[0]+2*stag->dof[1]+stag->dof[2]; ++dLocal,++count) { /* Dummies on back down edge, back face and back left edge */
2836: idxGlobalAll[count] = -1;
2837: }
2838: for (; dLocal<stag->dof[0]+3*stag->dof[1]+stag->dof[2]; ++d, ++dLocal, ++count) { /* Down left edge */
2839: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*eprNeighbor[neighbor] + i*entriesPerFace + d; /* Note i increment by face */
2840: }
2841: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies down face, left face and element */
2842: idxGlobalAll[count] = -1;
2843: }
2844: ++i;
2845: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2846: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2847: idxGlobalAll[count] = -1;
2848: }
2849: }
2850: } else {
2851: /* Right Up Front dummies */
2852: PetscInt i;
2853: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2854: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2855: idxGlobalAll[count] = -1;
2856: }
2857: }
2858: }
2859: ++j;
2860: /* Up Front additional dummy rows */
2861: for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
2862: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2863: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2864: idxGlobalAll[count] = -1;
2865: }
2866: }
2867: }
2868: } else {
2869: /* Up Front dummies */
2870: PetscInt j;
2871: for (j=0; j<ghostOffsetEnd[1]; ++j) {
2872: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2873: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2874: idxGlobalAll[count] = -1;
2875: }
2876: }
2877: }
2878: }
2879: }
2880: } else {
2881: PetscInt k = stag->n[2];
2883: /* Front partial ghost layer: rows 1/3: Down Front, on Down (Middle) */
2884: if (!dummyStart[1]) {
2885: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2886: const PetscInt j = nNeighbors[10][1] - ghostOffsetStart[1] + jghost; /* wrt down neighbor (10) */
2888: /* Down Front partial ghost row: colums 1/3: Left Down Front, on Left Down (Middle) */
2889: if (!star && !dummyStart[0]) {
2890: const PetscInt neighbor = 9;
2891: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0]; /* Note that we can't be a right boundary */
2892: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2893: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2894: PetscInt dLocal;
2895: 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 */
2896: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2897: }
2898: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2899: idxGlobalAll[count] = -1;
2900: }
2901: }
2902: } else {
2903: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2904: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2905: idxGlobalAll[count] = -1;
2906: }
2907: }
2908: }
2910: /* Down Front partial ghost row: columns 2/3: (Middle) Down Front, on (Middle) Down (Middle) */
2911: {
2912: const PetscInt neighbor = 10;
2913: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2914: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
2915: const PetscInt i = ighost - ghostOffsetStart[0];
2916: PetscInt dLocal;
2917: 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 */
2918: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2919: }
2920: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2921: idxGlobalAll[count] = -1;
2922: }
2923: }
2924: }
2926: if (!star && !dummyEnd[0]) {
2927: /* Down Front partial dummy row: columns 3/3: Right Down Front, on Right Down (Middle) */
2928: const PetscInt neighbor = 11;
2929: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
2930: PetscInt i;
2931: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2932: PetscInt dLocal;
2933: 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 */
2934: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2935: }
2936: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
2937: idxGlobalAll[count] = -1;
2938: }
2939: }
2940: } else if (dummyEnd[0]) {
2941: /* Right Down Front partial dummy element, living on (Middle) Down (Middle) rank */
2942: const PetscInt neighbor = 10;
2943: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
2944: PetscInt i,dLocal;
2945: i = stag->n[0];
2946: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
2947: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2948: }
2949: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
2950: idxGlobalAll[count] = -1;
2951: }
2952: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
2953: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
2954: }
2955: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
2956: idxGlobalAll[count] = -1;
2957: }
2958: ++i;
2959: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
2960: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2961: idxGlobalAll[count] = -1;
2962: }
2963: }
2964: } else {
2965: PetscInt i;
2966: /* Right Down Front dummies */
2967: for (i=0; i<ghostOffsetEnd[0]; ++i) {
2968: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2969: idxGlobalAll[count] = -1;
2970: }
2971: }
2972: }
2973: }
2974: } else {
2975: for (jghost = 0; jghost<ghostOffsetStart[1]; ++jghost) {
2976: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
2977: for (d=0; d<stag->entriesPerElement; ++d,++count) {
2978: idxGlobalAll[count] = -1;
2979: }
2980: }
2981: }
2982: }
2984: /* Front partial ghost layer: rows 2/3: (Middle) Front, on (Middle) (Middle) */
2985: {
2986: for (jghost = ghostOffsetStart[1]; jghost<stag->nGhost[1]-ghostOffsetEnd[1]; ++jghost) {
2987: const PetscInt j = jghost - ghostOffsetStart[1];
2989: /* (Middle) Front partial dummy row: columns 1/3: Left (Middle) Front, on Left (Middle) (Middle) */
2990: if (!dummyStart[0]) {
2991: const PetscInt neighbor = 12;
2992: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
2993: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
2994: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
2995: PetscInt dLocal;
2996: 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 */
2997: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
2998: }
2999: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3000: idxGlobalAll[count] = -1;
3001: }
3002: }
3003: } else {
3004: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3005: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3006: idxGlobalAll[count] = -1;
3007: }
3008: }
3009: }
3011: /* (Middle) Front partial dummy row: columns 2/3: (Middle) (Middle) Front, on (Middle) (Middle) (Middle) */
3012: {
3013: const PetscInt neighbor = 13;
3014: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3015: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3016: const PetscInt i = ighost - ghostOffsetStart[0];
3017: PetscInt dLocal;
3018: 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 */
3019: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3020: }
3021: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3022: idxGlobalAll[count] = -1;
3023: }
3024: }
3025: }
3027: if (!dummyEnd[0]) {
3028: /* (Middle) Front partial dummy row: columns 3/3: Right (Middle) Front, on Right (Middle) (Middle) */
3029: const PetscInt neighbor = 14;
3030: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3031: PetscInt i;
3032: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3033: PetscInt dLocal;
3034: 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 */
3035: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3036: }
3037: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3038: idxGlobalAll[count] = -1;
3039: }
3040: }
3041: } else {
3042: /* Right (Middle) Front partial dummy element, on (Middle) (Middle) (Middle) rank */
3043: const PetscInt neighbor = 13;
3044: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3045: PetscInt i,dLocal;
3046: i = stag->n[0];
3047: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3048: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3049: }
3050: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3051: idxGlobalAll[count] = -1;
3052: }
3053: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3054: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3055: }
3056: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3057: idxGlobalAll[count] = -1;
3058: }
3059: ++i;
3060: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3061: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3062: idxGlobalAll[count] = -1;
3063: }
3064: }
3065: }
3066: }
3067: }
3069: /* Front partial ghost layer: rows 3/3: Up Front, on Up (Middle) */
3070: if (!dummyEnd[1]) {
3071: PetscInt j;
3072: for (j=0; j<ghostOffsetEnd[1]; ++j) {
3074: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left Up (Middle) */
3075: if (!star && !dummyStart[0]) {
3076: const PetscInt neighbor = 15;
3077: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3078: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3079: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3080: PetscInt dLocal;
3081: 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 */
3082: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3083: }
3084: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3085: idxGlobalAll[count] = -1;
3086: }
3087: }
3088: } else {
3089: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3090: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3091: idxGlobalAll[count] = -1;
3092: }
3093: }
3094: }
3096: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) Up (Middle) */
3097: {
3098: const PetscInt neighbor = 16;
3099: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3100: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3101: const PetscInt i = ighost - ghostOffsetStart[0];
3102: PetscInt dLocal;
3103: 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 */
3104: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3105: }
3106: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3107: idxGlobalAll[count] = -1;
3108: }
3109: }
3110: }
3112: if (!star && !dummyEnd[0]) {
3113: /* Up Front partial dummy row: columsn 3/3: Right Up Front, on Right Up (Middle) */
3114: const PetscInt neighbor = 17;
3115: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3116: PetscInt i;
3117: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3118: PetscInt dLocal;
3119: 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 */
3120: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note j increment by face */
3121: }
3122: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining points */
3123: idxGlobalAll[count] = -1;
3124: }
3125: }
3126: } else if (dummyEnd[0]) {
3127: /* Right Up Front partial dummy element, on (Middle) Up (Middle) */
3128: const PetscInt neighbor = 16;
3129: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We+neighbor may be a right boundary */
3130: PetscInt i,dLocal;
3131: i = stag->n[0];
3132: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3133: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3134: }
3135: for (; dLocal<stag->dof[0]+stag->dof[1]; ++dLocal,++count) { /* Dummies on back down edge */
3136: idxGlobalAll[count] = -1;
3137: }
3138: for (; dLocal<stag->dof[0]+2*stag->dof[1]; ++d, ++dLocal, ++count) { /* left back edge */
3139: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerFace + d; /* Note i increment by face */
3140: }
3141: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3142: idxGlobalAll[count] = -1;
3143: }
3144: ++i;
3145: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3146: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3147: idxGlobalAll[count] = -1;
3148: }
3149: }
3150: } else {
3151: /* Right Up Front dummies */
3152: PetscInt i;
3153: /* Right Down Front dummies */
3154: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3155: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3156: idxGlobalAll[count] = -1;
3157: }
3158: }
3159: }
3160: }
3161: } else {
3162: /* Up Front partial dummy row */
3163: PetscInt j = stag->n[1];
3165: /* Up Front partial dummy row: columns 1/3: Left Up Front, on Left (Middle) (Middle) */
3166: if (!dummyStart[0]) {
3167: const PetscInt neighbor = 12;
3168: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0];
3169: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3170: const PetscInt i = nNeighbors[neighbor][0] - ghostOffsetStart[0] + ighost;
3171: PetscInt dLocal;
3172: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3173: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3174: }
3175: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3176: idxGlobalAll[count] = -1;
3177: }
3178: }
3179: } else {
3180: for (ighost = 0; ighost<ghostOffsetStart[0]; ++ighost) {
3181: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3182: idxGlobalAll[count] = -1;
3183: }
3184: }
3185: }
3187: /* Up Front partial dummy row: columns 2/3: (Middle) Up Front, on (Middle) (Middle) (Middle) */
3188: {
3189: const PetscInt neighbor = 13;
3190: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3191: for (ighost = ghostOffsetStart[0]; ighost<stag->nGhost[0]-ghostOffsetEnd[0]; ++ighost) {
3192: const PetscInt i = ighost - ghostOffsetStart[0];
3193: PetscInt dLocal;
3194: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3195: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3196: }
3197: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3198: idxGlobalAll[count] = -1;
3199: }
3200: }
3201: }
3203: if (!dummyEnd[0]) {
3204: /* Up Front partial dummy row: columns 3/3: Right Up Front, on Right (Middle) (Middle) */
3205: const PetscInt neighbor = 14;
3206: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (nextToDummyEnd[0] ? entriesPerEdge : 0); /* Neighbor may be a right boundary */
3207: PetscInt i;
3208: for (i=0; i<ghostOffsetEnd[0]; ++i) {
3209: PetscInt dLocal;
3210: for (d=0,dLocal=0; dLocal<stag->dof[0]+stag->dof[1]; ++d, ++dLocal, ++count) { /* Vertex and down back edge */
3211: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3212: }
3213: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3214: idxGlobalAll[count] = -1;
3215: }
3216: }
3217: } else {
3218: /* Right Up Front partial dummy element, on this rank */
3219: const PetscInt neighbor = 13;
3220: const PetscInt epFaceRow = entriesPerFace * nNeighbors[neighbor][0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3221: PetscInt i,dLocal;
3222: i = stag->n[0];
3223: for (d=0,dLocal=0; dLocal<stag->dof[0]; ++d, ++dLocal, ++count) { /* Vertex */
3224: idxGlobalAll[count] = globalOffsetNeighbor[neighbor] + k*eplNeighbor[neighbor] + j*epFaceRow + i*entriesPerEdge + d; /* Note increments */
3225: }
3226: for (; dLocal<stag->entriesPerElement; ++dLocal, ++count) { /* Dummies on remaining */
3227: idxGlobalAll[count] = -1;
3228: }
3229: ++i;
3230: for (; i<stag->n[0] + ghostOffsetEnd[0]; ++i) {
3231: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3232: idxGlobalAll[count] = -1;
3233: }
3234: }
3235: }
3236: ++j;
3237: /* Up Front additional dummy rows */
3238: for (; j<stag->n[1] + ghostOffsetEnd[1]; ++j) {
3239: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3240: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3241: idxGlobalAll[count] = -1;
3242: }
3243: }
3244: }
3245: }
3246: /* Front additional dummy layers */
3247: ++k;
3248: for (; k<stag->n[2] + ghostOffsetEnd[2]; ++k) {
3249: for (jghost = 0; jghost<stag->nGhost[1]; ++jghost) {
3250: for (ighost = 0; ighost<stag->nGhost[0]; ++ighost) {
3251: for (d=0; d<stag->entriesPerElement; ++d,++count) {
3252: idxGlobalAll[count] = -1;
3253: }
3254: }
3255: }
3256: }
3257: }
3259: /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
3260: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
3261: PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
3262: return(0);
3263: }
3265: static PetscErrorCode DMStagComputeLocationOffsets_3d(DM dm)
3266: {
3267: PetscErrorCode ierr;
3268: DM_Stag * const stag = (DM_Stag*)dm->data;
3269: const PetscInt epe = stag->entriesPerElement;
3270: const PetscInt epr = stag->nGhost[0]*epe;
3271: const PetscInt epl = stag->nGhost[1]*epr;
3274: PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
3275: stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] = 0;
3276: stag->locationOffsets[DMSTAG_BACK_DOWN] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + stag->dof[0];
3277: stag->locationOffsets[DMSTAG_BACK_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epe;
3278: stag->locationOffsets[DMSTAG_BACK_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN] + stag->dof[1];
3279: stag->locationOffsets[DMSTAG_BACK] = stag->locationOffsets[DMSTAG_BACK_LEFT] + stag->dof[1];
3280: stag->locationOffsets[DMSTAG_BACK_RIGHT] = stag->locationOffsets[DMSTAG_BACK_LEFT] + epe;
3281: stag->locationOffsets[DMSTAG_BACK_UP_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epr;
3282: stag->locationOffsets[DMSTAG_BACK_UP] = stag->locationOffsets[DMSTAG_BACK_DOWN] + epr;
3283: stag->locationOffsets[DMSTAG_BACK_UP_RIGHT] = stag->locationOffsets[DMSTAG_BACK_UP_LEFT] + epe;
3284: stag->locationOffsets[DMSTAG_DOWN_LEFT] = stag->locationOffsets[DMSTAG_BACK] + stag->dof[2];
3285: stag->locationOffsets[DMSTAG_DOWN] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + stag->dof[1];
3286: stag->locationOffsets[DMSTAG_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epe;
3287: stag->locationOffsets[DMSTAG_LEFT] = stag->locationOffsets[DMSTAG_DOWN] + stag->dof[2];
3288: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[2];
3289: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
3290: stag->locationOffsets[DMSTAG_UP_LEFT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epr;
3291: stag->locationOffsets[DMSTAG_UP] = stag->locationOffsets[DMSTAG_DOWN] + epr;
3292: stag->locationOffsets[DMSTAG_UP_RIGHT] = stag->locationOffsets[DMSTAG_UP_LEFT] + epe;
3293: stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] = stag->locationOffsets[DMSTAG_BACK_DOWN_LEFT] + epl;
3294: stag->locationOffsets[DMSTAG_FRONT_DOWN] = stag->locationOffsets[DMSTAG_BACK_DOWN] + epl;
3295: stag->locationOffsets[DMSTAG_FRONT_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epe;
3296: stag->locationOffsets[DMSTAG_FRONT_LEFT] = stag->locationOffsets[DMSTAG_BACK_LEFT] + epl;
3297: stag->locationOffsets[DMSTAG_FRONT] = stag->locationOffsets[DMSTAG_BACK] + epl;
3298: stag->locationOffsets[DMSTAG_FRONT_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_LEFT] + epe;
3299: stag->locationOffsets[DMSTAG_FRONT_UP_LEFT] = stag->locationOffsets[DMSTAG_FRONT_DOWN_LEFT] + epr;
3300: stag->locationOffsets[DMSTAG_FRONT_UP] = stag->locationOffsets[DMSTAG_FRONT_DOWN] + epr;
3301: stag->locationOffsets[DMSTAG_FRONT_UP_RIGHT] = stag->locationOffsets[DMSTAG_FRONT_UP_LEFT] + epe;
3302: return(0);
3303: }
3305: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_3d(DM dm)
3306: {
3307: PetscErrorCode ierr;
3308: DM_Stag * const stag = (DM_Stag*)dm->data;
3309: PetscInt *idxLocal,*idxGlobal,*globalOffsetsRecomputed;
3310: const PetscInt *globalOffsets;
3311: PetscInt count,d,entriesPerEdge,entriesPerFace,eprGhost,eplGhost,ghostOffsetStart[3],ghostOffsetEnd[3];
3312: IS isLocal,isGlobal;
3313: PetscBool dummyEnd[3];
3316: DMStagSetUpBuildGlobalOffsets_3d(dm,&globalOffsetsRecomputed); /* note that we don't actually use all of these. An available optimization is to pass them, when available */
3317: globalOffsets = globalOffsetsRecomputed;
3318: PetscMalloc1(stag->entries,&idxLocal);
3319: PetscMalloc1(stag->entries,&idxGlobal);
3320: for (d=0; d<3; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
3321: entriesPerFace = stag->dof[0] + 2*stag->dof[1] + stag->dof[2];
3322: entriesPerEdge = stag->dof[0] + stag->dof[1];
3323: eprGhost = stag->nGhost[0]*stag->entriesPerElement; /* epr = entries per (element) row */
3324: eplGhost = stag->nGhost[1]*eprGhost; /* epl = entries per (element) layer */
3325: count = 0;
3326: for (d=0; d<3; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
3327: for (d=0; d<3; ++d) ghostOffsetEnd[d] = stag->startGhost[d]+stag->nGhost[d] - (stag->start[d]+stag->n[d]);
3328: {
3329: const PetscInt neighbor = 13;
3330: const PetscInt epr = stag->entriesPerElement * stag->n[0] + (dummyEnd[0] ? entriesPerFace : 0); /* We may be a right boundary */
3331: const PetscInt epl = epr * stag->n[1] + (dummyEnd[1] ? stag->n[0] * entriesPerFace + (dummyEnd[0] ? entriesPerEdge : 0) : 0); /* We may be a top boundary */
3332: const PetscInt epFaceRow = entriesPerFace * stag->n[0] + (dummyEnd[0] ? entriesPerEdge : 0); /* We may be a right boundary */
3333: const PetscInt start0 = 0;
3334: const PetscInt start1 = 0;
3335: const PetscInt start2 = 0;
3336: const PetscInt startGhost0 = ghostOffsetStart[0];
3337: const PetscInt startGhost1 = ghostOffsetStart[1];
3338: const PetscInt startGhost2 = ghostOffsetStart[2];
3339: const PetscInt endGhost0 = stag->nGhost[0] - ghostOffsetEnd[0];
3340: const PetscInt endGhost1 = stag->nGhost[1] - ghostOffsetEnd[1];
3341: const PetscInt endGhost2 = stag->nGhost[2] - ghostOffsetEnd[2];
3342: const PetscBool extra0 = dummyEnd[0];
3343: const PetscBool extra1 = dummyEnd[1];
3344: const PetscBool extra2 = dummyEnd[2];
3345: 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);
3346: }
3347: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
3348: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
3349: {
3350: Vec local,global;
3351: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
3352: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
3353: VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
3354: VecDestroy(&global);
3355: VecDestroy(&local);
3356: }
3357: ISDestroy(&isLocal);
3358: ISDestroy(&isGlobal);
3359: if (globalOffsetsRecomputed) {
3360: PetscFree(globalOffsetsRecomputed);
3361: }
3362: return(0);
3363: }