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