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