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