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