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