1 /*
2 Code for manipulating distributed regular 3d arrays in parallel.
3 File created by Peter Mell 7/14/95
4 */
5
6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
7
8 #include <petscdraw.h>
DMView_DA_3d(DM da,PetscViewer viewer)9 static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer)
10 {
11 PetscMPIInt rank;
12 PetscBool isascii, isdraw, isglvis, isbinary;
13 DM_DA *dd = (DM_DA *)da->data;
14 Vec coordinates;
15 #if defined(PETSC_HAVE_MATLAB)
16 PetscBool ismatlab;
17 #endif
18
19 PetscFunctionBegin;
20 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));
21
22 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
23 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
24 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
25 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
26 #if defined(PETSC_HAVE_MATLAB)
27 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
28 #endif
29 if (isascii) {
30 PetscViewerFormat format;
31
32 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
33 PetscCall(PetscViewerGetFormat(viewer, &format));
34 if (format == PETSC_VIEWER_LOAD_BALANCE) {
35 PetscInt i, nmax = 0, nmin = PETSC_INT_MAX, navg = 0, *nz, nzlocal;
36 DMDALocalInfo info;
37 PetscMPIInt size;
38 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
39 PetscCall(DMDAGetLocalInfo(da, &info));
40 nzlocal = info.xm * info.ym * info.zm;
41 PetscCall(PetscMalloc1(size, &nz));
42 PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da)));
43 for (i = 0; i < size; i++) {
44 nmax = PetscMax(nmax, nz[i]);
45 nmin = PetscMin(nmin, nz[i]);
46 navg += nz[i];
47 }
48 PetscCall(PetscFree(nz));
49 navg = navg / size;
50 PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax));
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53 if (format != PETSC_VIEWER_ASCII_GLVIS) {
54 DMDALocalInfo info;
55 PetscCall(DMDAGetLocalInfo(da, &info));
56 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " P %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " p %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->P, dd->m, dd->n, dd->p, dd->w, dd->s));
57 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Z range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs,
58 info.xs + info.xm, info.ys, info.ys + info.ym, info.zs, info.zs + info.zm));
59 PetscCall(DMGetCoordinates(da, &coordinates));
60 #if !defined(PETSC_USE_COMPLEX)
61 if (coordinates) {
62 PetscInt last;
63 const PetscReal *coors;
64 PetscCall(VecGetArrayRead(coordinates, &coors));
65 PetscCall(VecGetLocalSize(coordinates, &last));
66 last = last - 3;
67 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Lower left corner %g %g %g : Upper right %g %g %g\n", (double)coors[0], (double)coors[1], (double)coors[2], (double)coors[last], (double)coors[last + 1], (double)coors[last + 2]));
68 PetscCall(VecRestoreArrayRead(coordinates, &coors));
69 }
70 #endif
71 PetscCall(PetscViewerFlush(viewer));
72 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
73 } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer));
74 } else if (isdraw) {
75 PetscDraw draw;
76 PetscReal ymin = -1.0, ymax = (PetscReal)dd->N;
77 PetscReal xmin = -1.0, xmax = (PetscReal)((dd->M + 2) * dd->P), x, y, ycoord, xcoord;
78 PetscInt k, plane, base;
79 const PetscInt *idx;
80 char node[10];
81 PetscBool isnull;
82
83 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
84 PetscCall(PetscDrawIsNull(draw, &isnull));
85 if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
86
87 PetscCall(PetscDrawCheckResizedWindow(draw));
88 PetscCall(PetscDrawClear(draw));
89 PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));
90
91 PetscDrawCollectiveBegin(draw);
92 /* first processor draw all node lines */
93 if (rank == 0) {
94 for (k = 0; k < dd->P; k++) {
95 ymin = 0.0;
96 ymax = (PetscReal)(dd->N - 1);
97 for (xmin = (PetscReal)(k * (dd->M + 1)); xmin < (PetscReal)(dd->M + (k * (dd->M + 1))); xmin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK));
98 xmin = (PetscReal)(k * (dd->M + 1));
99 xmax = xmin + (PetscReal)(dd->M - 1);
100 for (ymin = 0; ymin < (PetscReal)dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
101 }
102 }
103 PetscDrawCollectiveEnd(draw);
104 PetscCall(PetscDrawFlush(draw));
105 PetscCall(PetscDrawPause(draw));
106
107 PetscDrawCollectiveBegin(draw);
108 /*Go through and draw for each plane*/
109 for (k = 0; k < dd->P; k++) {
110 if ((k >= dd->zs) && (k < dd->ze)) {
111 /* draw my box */
112 ymin = dd->ys;
113 ymax = dd->ye - 1;
114 xmin = dd->xs / dd->w + (dd->M + 1) * k;
115 xmax = (dd->xe - 1) / dd->w + (dd->M + 1) * k;
116
117 PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
118 PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
119 PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
120 PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED));
121
122 xmin = dd->xs / dd->w;
123 xmax = (dd->xe - 1) / dd->w;
124
125 /* identify which processor owns the box */
126 PetscCall(PetscSNPrintf(node, sizeof(node), "%d", rank));
127 PetscCall(PetscDrawString(draw, xmin + (dd->M + 1) * k + .2, ymin + .3, PETSC_DRAW_RED, node));
128 /* put in numbers*/
129 base = (dd->base + (dd->xe - dd->xs) * (dd->ye - dd->ys) * (k - dd->zs)) / dd->w;
130 for (y = ymin; y <= ymax; y++) {
131 for (x = xmin + (dd->M + 1) * k; x <= xmax + (dd->M + 1) * k; x++) {
132 PetscCall(PetscSNPrintf(node, sizeof(node), "%" PetscInt_FMT, base++));
133 PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
134 }
135 }
136 }
137 }
138 PetscDrawCollectiveEnd(draw);
139 PetscCall(PetscDrawFlush(draw));
140 PetscCall(PetscDrawPause(draw));
141
142 PetscDrawCollectiveBegin(draw);
143 for (k = 0 - dd->s; k < dd->P + dd->s; k++) {
144 /* Go through and draw for each plane */
145 if ((k >= dd->Zs) && (k < dd->Ze)) {
146 /* overlay ghost numbers, useful for error checking */
147 base = (dd->Xe - dd->Xs) * (dd->Ye - dd->Ys) * (k - dd->Zs) / dd->w;
148 PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx));
149 plane = k;
150 /* Keep z wrap around points on the drawing */
151 if (k < 0) plane = dd->P + k;
152 if (k >= dd->P) plane = k - dd->P;
153 ymin = dd->Ys;
154 ymax = dd->Ye;
155 xmin = (dd->M + 1) * plane * dd->w;
156 xmax = (dd->M + 1) * plane * dd->w + dd->M * dd->w;
157 for (y = ymin; y < ymax; y++) {
158 for (x = xmin + dd->Xs; x < xmin + dd->Xe; x += dd->w) {
159 PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, idx[base]));
160 ycoord = y;
161 /*Keep y wrap around points on drawing */
162 if (y < 0) ycoord = dd->N + y;
163 if (y >= dd->N) ycoord = y - dd->N;
164 xcoord = x; /* Keep x wrap points on drawing */
165 if (x < xmin) xcoord = xmax - (xmin - x);
166 if (x >= xmax) xcoord = xmin + (x - xmax);
167 PetscCall(PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node));
168 base++;
169 }
170 }
171 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
172 }
173 }
174 PetscDrawCollectiveEnd(draw);
175 PetscCall(PetscDrawFlush(draw));
176 PetscCall(PetscDrawPause(draw));
177 PetscCall(PetscDrawSave(draw));
178 } else if (isglvis) {
179 PetscCall(DMView_DA_GLVis(da, viewer));
180 } else if (isbinary) {
181 PetscCall(DMView_DA_Binary(da, viewer));
182 #if defined(PETSC_HAVE_MATLAB)
183 } else if (ismatlab) {
184 PetscCall(DMView_DA_Matlab(da, viewer));
185 #endif
186 }
187 PetscFunctionReturn(PETSC_SUCCESS);
188 }
189
DMSetUp_DA_3D(DM da)190 PetscErrorCode DMSetUp_DA_3D(DM da)
191 {
192 DM_DA *dd = (DM_DA *)da->data;
193 const PetscInt M = dd->M;
194 const PetscInt N = dd->N;
195 const PetscInt P = dd->P;
196 PetscMPIInt m, n, p;
197 const PetscInt dof = dd->w;
198 const PetscInt s = dd->s;
199 DMBoundaryType bx = dd->bx;
200 DMBoundaryType by = dd->by;
201 DMBoundaryType bz = dd->bz;
202 DMDAStencilType stencil_type = dd->stencil_type;
203 PetscInt *lx = dd->lx;
204 PetscInt *ly = dd->ly;
205 PetscInt *lz = dd->lz;
206 MPI_Comm comm;
207 PetscMPIInt rank, size;
208 PetscInt xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0;
209 PetscInt Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm;
210 PetscInt left, right, up, down, bottom, top, i, j, k, *idx, nn;
211 PetscMPIInt n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14;
212 PetscMPIInt n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26;
213 PetscInt *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z;
214 PetscMPIInt sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0;
215 PetscMPIInt sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0;
216 PetscMPIInt sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0;
217 Vec local, global;
218 VecScatter gtol;
219 IS to, from;
220 PetscBool twod;
221
222 PetscFunctionBegin;
223 PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && by != DM_BOUNDARY_MIRROR && bz != DM_BOUNDARY_MIRROR), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box stencil");
224 PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
225 #if !defined(PETSC_USE_64BIT_INDICES)
226 PetscCheck(((PetscInt64)M) * ((PetscInt64)N) * ((PetscInt64)P) * ((PetscInt64)dof) <= (PetscInt64)PETSC_MPI_INT_MAX, comm, PETSC_ERR_INT_OVERFLOW, "Mesh of %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " (dof) is too large for 32-bit indices", M, N, P, dof);
227 #endif
228 PetscCall(PetscMPIIntCast(dd->m, &m));
229 PetscCall(PetscMPIIntCast(dd->n, &n));
230 PetscCall(PetscMPIIntCast(dd->p, &p));
231
232 PetscCallMPI(MPI_Comm_size(comm, &size));
233 PetscCallMPI(MPI_Comm_rank(comm, &rank));
234
235 if (m != PETSC_DECIDE) {
236 PetscCheck(m >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %d", m);
237 PetscCheck(m <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %d %d", m, size);
238 }
239 if (n != PETSC_DECIDE) {
240 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %d", n);
241 PetscCheck(n <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %d %d", n, size);
242 }
243 if (p != PETSC_DECIDE) {
244 PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Z direction: %d", p);
245 PetscCheck(p <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Z direction: %d %d", p, size);
246 }
247 PetscCheck(m <= 0 || n <= 0 || p <= 0 || m * n * p == size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m %d * n %d * p %d != size %d", m, n, p, size);
248
249 /* Partition the array among the processors */
250 if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
251 m = size / (n * p);
252 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
253 n = size / (m * p);
254 } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
255 p = size / (m * n);
256 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
257 /* try for squarish distribution */
258 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N * p)));
259 if (!m) m = 1;
260 while (m > 0) {
261 n = size / (m * p);
262 if (m * n * p == size) break;
263 m--;
264 }
265 PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad p value: p = %d", p);
266 if (M > N && m < n) {
267 PetscMPIInt _m = m;
268 m = n;
269 n = _m;
270 }
271 } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
272 /* try for squarish distribution */
273 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
274 if (!m) m = 1;
275 while (m > 0) {
276 p = size / (m * n);
277 if (m * n * p == size) break;
278 m--;
279 }
280 PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad n value: n = %d", n);
281 if (M > P && m < p) {
282 PetscMPIInt _m = m;
283 m = p;
284 p = _m;
285 }
286 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
287 /* try for squarish distribution */
288 n = (int)(0.5 + PetscSqrtReal(((PetscReal)N) * ((PetscReal)size) / ((PetscReal)P * m)));
289 if (!n) n = 1;
290 while (n > 0) {
291 p = size / (m * n);
292 if (m * n * p == size) break;
293 n--;
294 }
295 PetscCheck(n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad m value: m = %d", n);
296 if (N > P && n < p) {
297 PetscMPIInt _n = n;
298 n = p;
299 p = _n;
300 }
301 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
302 /* try for squarish distribution */
303 n = (PetscMPIInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.)));
304 if (!n) n = 1;
305 while (n > 0) {
306 pm = size / n;
307 if (n * pm == size) break;
308 n--;
309 }
310 if (!n) n = 1;
311 m = (PetscMPIInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
312 if (!m) m = 1;
313 while (m > 0) {
314 p = size / (m * n);
315 if (m * n * p == size) break;
316 m--;
317 }
318 if (M > P && m < p) {
319 PetscMPIInt _m = m;
320 m = p;
321 p = _m;
322 }
323 } else PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition");
324
325 PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_PLIB, "Could not find good partition");
326 PetscCheck(M >= m, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %d", M, m);
327 PetscCheck(N >= n, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %d", N, n);
328 PetscCheck(P >= p, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in z direction is too fine! %" PetscInt_FMT " %d", P, p);
329
330 /*
331 Determine locally owned region
332 [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
333 */
334
335 if (!lx) {
336 PetscCall(PetscMalloc1(m, &dd->lx));
337 lx = dd->lx;
338 for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m));
339 }
340 x = lx[rank % m];
341 xs = 0;
342 for (i = 0; i < (rank % m); i++) xs += lx[i];
343 PetscCheck(x >= s || (m <= 1 && bx != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, x, s);
344
345 if (!ly) {
346 PetscCall(PetscMalloc1(n, &dd->ly));
347 ly = dd->ly;
348 for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n));
349 }
350 y = ly[(rank % (m * n)) / m];
351 PetscCheck(y >= s || (n <= 1 && by != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, y, s);
352
353 ys = 0;
354 for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i];
355
356 if (!lz) {
357 PetscCall(PetscMalloc1(p, &dd->lz));
358 lz = dd->lz;
359 for (i = 0; i < p; i++) lz[i] = P / p + ((P % p) > (i % p));
360 }
361 z = lz[rank / (m * n)];
362
363 PetscCheck((x > s) || (bx != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", x, s);
364 PetscCheck((y > s) || (by != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", y, s);
365 PetscCheck((z > s) || (bz != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local z-width of domain z %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", z, s);
366
367 /* note this is different than x- and y-, as we will handle as an important special
368 case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems
369 in a 3D code. Additional code for this case is noted with "2d case" comments */
370 twod = PETSC_FALSE;
371 if (P == 1) twod = PETSC_TRUE;
372 else PetscCheck(z >= s || (p <= 1 && bz != DM_BOUNDARY_PERIODIC), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local z-width of domain z %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, z, s);
373 zs = 0;
374 for (i = 0; i < (rank / (m * n)); i++) zs += lz[i];
375 ye = ys + y;
376 xe = xs + x;
377 ze = zs + z;
378
379 /* determine ghost region (Xs) and region scattered into (IXs) */
380 if (xs - s > 0) {
381 Xs = xs - s;
382 IXs = xs - s;
383 } else {
384 if (bx) Xs = xs - s;
385 else Xs = 0;
386 IXs = 0;
387 }
388 if (xe + s <= M) {
389 Xe = xe + s;
390 IXe = xe + s;
391 } else {
392 if (bx) {
393 Xs = xs - s;
394 Xe = xe + s;
395 } else Xe = M;
396 IXe = M;
397 }
398
399 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
400 IXs = xs - s;
401 IXe = xe + s;
402 Xs = xs - s;
403 Xe = xe + s;
404 }
405
406 if (ys - s > 0) {
407 Ys = ys - s;
408 IYs = ys - s;
409 } else {
410 if (by) Ys = ys - s;
411 else Ys = 0;
412 IYs = 0;
413 }
414 if (ye + s <= N) {
415 Ye = ye + s;
416 IYe = ye + s;
417 } else {
418 if (by) Ye = ye + s;
419 else Ye = N;
420 IYe = N;
421 }
422
423 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
424 IYs = ys - s;
425 IYe = ye + s;
426 Ys = ys - s;
427 Ye = ye + s;
428 }
429
430 if (zs - s > 0) {
431 Zs = zs - s;
432 IZs = zs - s;
433 } else {
434 if (bz) Zs = zs - s;
435 else Zs = 0;
436 IZs = 0;
437 }
438 if (ze + s <= P) {
439 Ze = ze + s;
440 IZe = ze + s;
441 } else {
442 if (bz) Ze = ze + s;
443 else Ze = P;
444 IZe = P;
445 }
446
447 if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
448 IZs = zs - s;
449 IZe = ze + s;
450 Zs = zs - s;
451 Ze = ze + s;
452 }
453
454 /* Resize all X parameters to reflect w */
455 s_x = s;
456 s_y = s;
457 s_z = s;
458
459 /* determine starting point of each processor */
460 nn = x * y * z;
461 PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims));
462 PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm));
463 bases[0] = 0;
464 for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
465 for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
466 base = bases[rank] * dof;
467
468 /* allocate the base parallel and sequential vectors */
469 dd->Nlocal = x * y * z * dof;
470 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global));
471 dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof;
472 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local));
473
474 /* generate global to local vector scatter and local to global mapping*/
475
476 /* global to local must include ghost points within the domain,
477 but not ghost points outside the domain that aren't periodic */
478 PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx));
479 if (stencil_type == DMDA_STENCIL_BOX) {
480 left = IXs - Xs;
481 right = left + (IXe - IXs);
482 bottom = IYs - Ys;
483 top = bottom + (IYe - IYs);
484 down = IZs - Zs;
485 up = down + (IZe - IZs);
486 count = 0;
487 for (i = down; i < up; i++) {
488 for (j = bottom; j < top; j++) {
489 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
490 }
491 }
492 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
493 } else {
494 /* This is way ugly! We need to list the funny cross type region */
495 left = xs - Xs;
496 right = left + x;
497 bottom = ys - Ys;
498 top = bottom + y;
499 down = zs - Zs;
500 up = down + z;
501 count = 0;
502 /* the bottom chunk */
503 for (i = (IZs - Zs); i < down; i++) {
504 for (j = bottom; j < top; j++) {
505 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
506 }
507 }
508 /* the middle piece */
509 for (i = down; i < up; i++) {
510 /* front */
511 for (j = (IYs - Ys); j < bottom; j++) {
512 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
513 }
514 /* middle */
515 for (j = bottom; j < top; j++) {
516 for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
517 }
518 /* back */
519 for (j = top; j < top + IYe - ye; j++) {
520 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
521 }
522 }
523 /* the top piece */
524 for (i = up; i < up + IZe - ze; i++) {
525 for (j = bottom; j < top; j++) {
526 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
527 }
528 }
529 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
530 }
531
532 /* determine who lies on each side of use stored in n24 n25 n26
533 n21 n22 n23
534 n18 n19 n20
535
536 n15 n16 n17
537 n12 n14
538 n9 n10 n11
539
540 n6 n7 n8
541 n3 n4 n5
542 n0 n1 n2
543 */
544
545 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
546 /* Assume Nodes are Internal to the Cube */
547 n0 = rank - m * n - m - 1;
548 n1 = rank - m * n - m;
549 n2 = rank - m * n - m + 1;
550 n3 = rank - m * n - 1;
551 n4 = rank - m * n;
552 n5 = rank - m * n + 1;
553 n6 = rank - m * n + m - 1;
554 n7 = rank - m * n + m;
555 n8 = rank - m * n + m + 1;
556
557 n9 = rank - m - 1;
558 n10 = rank - m;
559 n11 = rank - m + 1;
560 n12 = rank - 1;
561 n14 = rank + 1;
562 n15 = rank + m - 1;
563 n16 = rank + m;
564 n17 = rank + m + 1;
565
566 n18 = rank + m * n - m - 1;
567 n19 = rank + m * n - m;
568 n20 = rank + m * n - m + 1;
569 n21 = rank + m * n - 1;
570 n22 = rank + m * n;
571 n23 = rank + m * n + 1;
572 n24 = rank + m * n + m - 1;
573 n25 = rank + m * n + m;
574 n26 = rank + m * n + m + 1;
575
576 /* Assume Pieces are on Faces of Cube */
577
578 if (xs == 0) { /* First assume not corner or edge */
579 n0 = rank - 1 - (m * n);
580 n3 = rank + m - 1 - (m * n);
581 n6 = rank + 2 * m - 1 - (m * n);
582 n9 = rank - 1;
583 n12 = rank + m - 1;
584 n15 = rank + 2 * m - 1;
585 n18 = rank - 1 + (m * n);
586 n21 = rank + m - 1 + (m * n);
587 n24 = rank + 2 * m - 1 + (m * n);
588 }
589
590 if (xe == M) { /* First assume not corner or edge */
591 n2 = rank - 2 * m + 1 - (m * n);
592 n5 = rank - m + 1 - (m * n);
593 n8 = rank + 1 - (m * n);
594 n11 = rank - 2 * m + 1;
595 n14 = rank - m + 1;
596 n17 = rank + 1;
597 n20 = rank - 2 * m + 1 + (m * n);
598 n23 = rank - m + 1 + (m * n);
599 n26 = rank + 1 + (m * n);
600 }
601
602 if (ys == 0) { /* First assume not corner or edge */
603 n0 = rank + m * (n - 1) - 1 - (m * n);
604 n1 = rank + m * (n - 1) - (m * n);
605 n2 = rank + m * (n - 1) + 1 - (m * n);
606 n9 = rank + m * (n - 1) - 1;
607 n10 = rank + m * (n - 1);
608 n11 = rank + m * (n - 1) + 1;
609 n18 = rank + m * (n - 1) - 1 + (m * n);
610 n19 = rank + m * (n - 1) + (m * n);
611 n20 = rank + m * (n - 1) + 1 + (m * n);
612 }
613
614 if (ye == N) { /* First assume not corner or edge */
615 n6 = rank - m * (n - 1) - 1 - (m * n);
616 n7 = rank - m * (n - 1) - (m * n);
617 n8 = rank - m * (n - 1) + 1 - (m * n);
618 n15 = rank - m * (n - 1) - 1;
619 n16 = rank - m * (n - 1);
620 n17 = rank - m * (n - 1) + 1;
621 n24 = rank - m * (n - 1) - 1 + (m * n);
622 n25 = rank - m * (n - 1) + (m * n);
623 n26 = rank - m * (n - 1) + 1 + (m * n);
624 }
625
626 if (zs == 0) { /* First assume not corner or edge */
627 n0 = size - (m * n) + rank - m - 1;
628 n1 = size - (m * n) + rank - m;
629 n2 = size - (m * n) + rank - m + 1;
630 n3 = size - (m * n) + rank - 1;
631 n4 = size - (m * n) + rank;
632 n5 = size - (m * n) + rank + 1;
633 n6 = size - (m * n) + rank + m - 1;
634 n7 = size - (m * n) + rank + m;
635 n8 = size - (m * n) + rank + m + 1;
636 }
637
638 if (ze == P) { /* First assume not corner or edge */
639 n18 = (m * n) - (size - rank) - m - 1;
640 n19 = (m * n) - (size - rank) - m;
641 n20 = (m * n) - (size - rank) - m + 1;
642 n21 = (m * n) - (size - rank) - 1;
643 n22 = (m * n) - (size - rank);
644 n23 = (m * n) - (size - rank) + 1;
645 n24 = (m * n) - (size - rank) + m - 1;
646 n25 = (m * n) - (size - rank) + m;
647 n26 = (m * n) - (size - rank) + m + 1;
648 }
649
650 if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */
651 n0 = size - m * n + rank + m - 1 - m;
652 n3 = size - m * n + rank + m - 1;
653 n6 = size - m * n + rank + m - 1 + m;
654 }
655
656 if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */
657 n18 = m * n - (size - rank) + m - 1 - m;
658 n21 = m * n - (size - rank) + m - 1;
659 n24 = m * n - (size - rank) + m - 1 + m;
660 }
661
662 if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */
663 n0 = rank + m * n - 1 - m * n;
664 n9 = rank + m * n - 1;
665 n18 = rank + m * n - 1 + m * n;
666 }
667
668 if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */
669 n6 = rank - m * (n - 1) + m - 1 - m * n;
670 n15 = rank - m * (n - 1) + m - 1;
671 n24 = rank - m * (n - 1) + m - 1 + m * n;
672 }
673
674 if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */
675 n2 = size - (m * n - rank) - (m - 1) - m;
676 n5 = size - (m * n - rank) - (m - 1);
677 n8 = size - (m * n - rank) - (m - 1) + m;
678 }
679
680 if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */
681 n20 = m * n - (size - rank) - (m - 1) - m;
682 n23 = m * n - (size - rank) - (m - 1);
683 n26 = m * n - (size - rank) - (m - 1) + m;
684 }
685
686 if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */
687 n2 = rank + m * (n - 1) - (m - 1) - m * n;
688 n11 = rank + m * (n - 1) - (m - 1);
689 n20 = rank + m * (n - 1) - (m - 1) + m * n;
690 }
691
692 if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */
693 n8 = rank - m * n + 1 - m * n;
694 n17 = rank - m * n + 1;
695 n26 = rank - m * n + 1 + m * n;
696 }
697
698 if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */
699 n0 = size - m + rank - 1;
700 n1 = size - m + rank;
701 n2 = size - m + rank + 1;
702 }
703
704 if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */
705 n18 = m * n - (size - rank) + m * (n - 1) - 1;
706 n19 = m * n - (size - rank) + m * (n - 1);
707 n20 = m * n - (size - rank) + m * (n - 1) + 1;
708 }
709
710 if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */
711 n6 = size - (m * n - rank) - m * (n - 1) - 1;
712 n7 = size - (m * n - rank) - m * (n - 1);
713 n8 = size - (m * n - rank) - m * (n - 1) + 1;
714 }
715
716 if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */
717 n24 = rank - (size - m) - 1;
718 n25 = rank - (size - m);
719 n26 = rank - (size - m) + 1;
720 }
721
722 /* Check for Corners */
723 if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1;
724 if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1;
725 if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1);
726 if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1;
727 if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m;
728 if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m;
729 if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n;
730 if ((xe == M) && (ye == N) && (ze == P)) n26 = 0;
731
732 /* Check for when not X,Y, and Z Periodic */
733
734 /* If not X periodic */
735 if (bx != DM_BOUNDARY_PERIODIC) {
736 if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;
737 if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;
738 }
739
740 /* If not Y periodic */
741 if (by != DM_BOUNDARY_PERIODIC) {
742 if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;
743 if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;
744 }
745
746 /* If not Z periodic */
747 if (bz != DM_BOUNDARY_PERIODIC) {
748 if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;
749 if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;
750 }
751
752 PetscCall(PetscMalloc1(27, &dd->neighbors));
753
754 dd->neighbors[0] = n0;
755 dd->neighbors[1] = n1;
756 dd->neighbors[2] = n2;
757 dd->neighbors[3] = n3;
758 dd->neighbors[4] = n4;
759 dd->neighbors[5] = n5;
760 dd->neighbors[6] = n6;
761 dd->neighbors[7] = n7;
762 dd->neighbors[8] = n8;
763 dd->neighbors[9] = n9;
764 dd->neighbors[10] = n10;
765 dd->neighbors[11] = n11;
766 dd->neighbors[12] = n12;
767 dd->neighbors[13] = rank;
768 dd->neighbors[14] = n14;
769 dd->neighbors[15] = n15;
770 dd->neighbors[16] = n16;
771 dd->neighbors[17] = n17;
772 dd->neighbors[18] = n18;
773 dd->neighbors[19] = n19;
774 dd->neighbors[20] = n20;
775 dd->neighbors[21] = n21;
776 dd->neighbors[22] = n22;
777 dd->neighbors[23] = n23;
778 dd->neighbors[24] = n24;
779 dd->neighbors[25] = n25;
780 dd->neighbors[26] = n26;
781
782 /* If star stencil then delete the corner neighbors */
783 if (stencil_type == DMDA_STENCIL_STAR) {
784 /* save information about corner neighbors */
785 sn0 = n0;
786 sn1 = n1;
787 sn2 = n2;
788 sn3 = n3;
789 sn5 = n5;
790 sn6 = n6;
791 sn7 = n7;
792 sn8 = n8;
793 sn9 = n9;
794 sn11 = n11;
795 sn15 = n15;
796 sn17 = n17;
797 sn18 = n18;
798 sn19 = n19;
799 sn20 = n20;
800 sn21 = n21;
801 sn23 = n23;
802 sn24 = n24;
803 sn25 = n25;
804 sn26 = n26;
805 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
806 }
807
808 PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx));
809
810 nn = 0;
811 /* Bottom Level */
812 for (k = 0; k < s_z; k++) {
813 for (i = 1; i <= s_y; i++) {
814 if (n0 >= 0) { /* left below */
815 x_t = lx[n0 % m];
816 y_t = ly[(n0 % (m * n)) / m];
817 z_t = lz[n0 / (m * n)];
818 s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x - (s_z - k - 1) * x_t * y_t;
819 if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */
820 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
821 }
822 if (n1 >= 0) { /* directly below */
823 x_t = x;
824 y_t = ly[(n1 % (m * n)) / m];
825 z_t = lz[n1 / (m * n)];
826 s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
827 if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
828 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
829 }
830 if (n2 >= 0) { /* right below */
831 x_t = lx[n2 % m];
832 y_t = ly[(n2 % (m * n)) / m];
833 z_t = lz[n2 / (m * n)];
834 s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
835 if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
836 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
837 }
838 }
839
840 for (i = 0; i < y; i++) {
841 if (n3 >= 0) { /* directly left */
842 x_t = lx[n3 % m];
843 y_t = y;
844 z_t = lz[n3 / (m * n)];
845 s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
846 if (twod && (s_t < 0)) s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - x_t * y_t; /* 2D case */
847 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
848 }
849
850 if (n4 >= 0) { /* middle */
851 x_t = x;
852 y_t = y;
853 z_t = lz[n4 / (m * n)];
854 s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
855 if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
856 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
857 } else if (bz == DM_BOUNDARY_MIRROR) {
858 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
859 }
860
861 if (n5 >= 0) { /* directly right */
862 x_t = lx[n5 % m];
863 y_t = y;
864 z_t = lz[n5 / (m * n)];
865 s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
866 if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
867 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
868 }
869 }
870
871 for (i = 1; i <= s_y; i++) {
872 if (n6 >= 0) { /* left above */
873 x_t = lx[n6 % m];
874 y_t = ly[(n6 % (m * n)) / m];
875 z_t = lz[n6 / (m * n)];
876 s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
877 if (twod && (s_t < 0)) s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - x_t * y_t; /* 2D case */
878 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
879 }
880 if (n7 >= 0) { /* directly above */
881 x_t = x;
882 y_t = ly[(n7 % (m * n)) / m];
883 z_t = lz[n7 / (m * n)];
884 s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
885 if (twod && (s_t < 0)) s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
886 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
887 }
888 if (n8 >= 0) { /* right above */
889 x_t = lx[n8 % m];
890 y_t = ly[(n8 % (m * n)) / m];
891 z_t = lz[n8 / (m * n)];
892 s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
893 if (twod && (s_t < 0)) s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
894 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
895 }
896 }
897 }
898
899 /* Middle Level */
900 for (k = 0; k < z; k++) {
901 for (i = 1; i <= s_y; i++) {
902 if (n9 >= 0) { /* left below */
903 x_t = lx[n9 % m];
904 y_t = ly[(n9 % (m * n)) / m];
905 /* z_t = z; */
906 s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
907 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
908 }
909 if (n10 >= 0) { /* directly below */
910 x_t = x;
911 y_t = ly[(n10 % (m * n)) / m];
912 /* z_t = z; */
913 s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
914 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
915 } else if (by == DM_BOUNDARY_MIRROR) {
916 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
917 }
918 if (n11 >= 0) { /* right below */
919 x_t = lx[n11 % m];
920 y_t = ly[(n11 % (m * n)) / m];
921 /* z_t = z; */
922 s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
923 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
924 }
925 }
926
927 for (i = 0; i < y; i++) {
928 if (n12 >= 0) { /* directly left */
929 x_t = lx[n12 % m];
930 y_t = y;
931 /* z_t = z; */
932 s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
933 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
934 } else if (bx == DM_BOUNDARY_MIRROR) {
935 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
936 }
937
938 /* Interior */
939 s_t = bases[rank] + i * x + k * x * y;
940 for (j = 0; j < x; j++) idx[nn++] = s_t++;
941
942 if (n14 >= 0) { /* directly right */
943 x_t = lx[n14 % m];
944 y_t = y;
945 /* z_t = z; */
946 s_t = bases[n14] + i * x_t + k * x_t * y_t;
947 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
948 } else if (bx == DM_BOUNDARY_MIRROR) {
949 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
950 }
951 }
952
953 for (i = 1; i <= s_y; i++) {
954 if (n15 >= 0) { /* left above */
955 x_t = lx[n15 % m];
956 y_t = ly[(n15 % (m * n)) / m];
957 /* z_t = z; */
958 s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
959 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
960 }
961 if (n16 >= 0) { /* directly above */
962 x_t = x;
963 y_t = ly[(n16 % (m * n)) / m];
964 /* z_t = z; */
965 s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
966 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
967 } else if (by == DM_BOUNDARY_MIRROR) {
968 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
969 }
970 if (n17 >= 0) { /* right above */
971 x_t = lx[n17 % m];
972 y_t = ly[(n17 % (m * n)) / m];
973 /* z_t = z; */
974 s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
975 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
976 }
977 }
978 }
979
980 /* Upper Level */
981 for (k = 0; k < s_z; k++) {
982 for (i = 1; i <= s_y; i++) {
983 if (n18 >= 0) { /* left below */
984 x_t = lx[n18 % m];
985 y_t = ly[(n18 % (m * n)) / m];
986 /* z_t = lz[n18 / (m*n)]; */
987 s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
988 if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */
989 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
990 }
991 if (n19 >= 0) { /* directly below */
992 x_t = x;
993 y_t = ly[(n19 % (m * n)) / m];
994 /* z_t = lz[n19 / (m*n)]; */
995 s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
996 if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
997 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
998 }
999 if (n20 >= 0) { /* right below */
1000 x_t = lx[n20 % m];
1001 y_t = ly[(n20 % (m * n)) / m];
1002 /* z_t = lz[n20 / (m*n)]; */
1003 s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1004 if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
1005 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1006 }
1007 }
1008
1009 for (i = 0; i < y; i++) {
1010 if (n21 >= 0) { /* directly left */
1011 x_t = lx[n21 % m];
1012 y_t = y;
1013 /* z_t = lz[n21 / (m*n)]; */
1014 s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1015 if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */
1016 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1017 }
1018
1019 if (n22 >= 0) { /* middle */
1020 x_t = x;
1021 y_t = y;
1022 /* z_t = lz[n22 / (m*n)]; */
1023 s_t = bases[n22] + i * x_t + k * x_t * y_t;
1024 if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */
1025 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1026 } else if (bz == DM_BOUNDARY_MIRROR) {
1027 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1028 }
1029
1030 if (n23 >= 0) { /* directly right */
1031 x_t = lx[n23 % m];
1032 y_t = y;
1033 /* z_t = lz[n23 / (m*n)]; */
1034 s_t = bases[n23] + i * x_t + k * x_t * y_t;
1035 if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */
1036 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1037 }
1038 }
1039
1040 for (i = 1; i <= s_y; i++) {
1041 if (n24 >= 0) { /* left above */
1042 x_t = lx[n24 % m];
1043 y_t = ly[(n24 % (m * n)) / m];
1044 /* z_t = lz[n24 / (m*n)]; */
1045 s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1046 if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */
1047 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1048 }
1049 if (n25 >= 0) { /* directly above */
1050 x_t = x;
1051 y_t = ly[(n25 % (m * n)) / m];
1052 /* z_t = lz[n25 / (m*n)]; */
1053 s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1054 if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */
1055 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1056 }
1057 if (n26 >= 0) { /* right above */
1058 x_t = lx[n26 % m];
1059 y_t = ly[(n26 % (m * n)) / m];
1060 /* z_t = lz[n26 / (m*n)]; */
1061 s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1062 if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */
1063 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1064 }
1065 }
1066 }
1067
1068 PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
1069 PetscCall(VecScatterCreate(global, from, local, to, >ol));
1070 PetscCall(ISDestroy(&to));
1071 PetscCall(ISDestroy(&from));
1072
1073 if (stencil_type == DMDA_STENCIL_STAR) {
1074 n0 = sn0;
1075 n1 = sn1;
1076 n2 = sn2;
1077 n3 = sn3;
1078 n5 = sn5;
1079 n6 = sn6;
1080 n7 = sn7;
1081 n8 = sn8;
1082 n9 = sn9;
1083 n11 = sn11;
1084 n15 = sn15;
1085 n17 = sn17;
1086 n18 = sn18;
1087 n19 = sn19;
1088 n20 = sn20;
1089 n21 = sn21;
1090 n23 = sn23;
1091 n24 = sn24;
1092 n25 = sn25;
1093 n26 = sn26;
1094 }
1095
1096 if ((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz)) {
1097 /*
1098 Recompute the local to global mappings, this time keeping the
1099 information about the cross corner processor numbers.
1100 */
1101 nn = 0;
1102 /* Bottom Level */
1103 for (k = 0; k < s_z; k++) {
1104 for (i = 1; i <= s_y; i++) {
1105 if (n0 >= 0) { /* left below */
1106 x_t = lx[n0 % m];
1107 y_t = ly[(n0 % (m * n)) / m];
1108 z_t = lz[n0 / (m * n)];
1109 s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x - (s_z - k - 1) * x_t * y_t;
1110 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1111 } else if (Xs - xs < 0 && Ys - ys < 0 && Zs - zs < 0) {
1112 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1113 }
1114 if (n1 >= 0) { /* directly below */
1115 x_t = x;
1116 y_t = ly[(n1 % (m * n)) / m];
1117 z_t = lz[n1 / (m * n)];
1118 s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1119 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1120 } else if (Ys - ys < 0 && Zs - zs < 0) {
1121 for (j = 0; j < x; j++) idx[nn++] = -1;
1122 }
1123 if (n2 >= 0) { /* right below */
1124 x_t = lx[n2 % m];
1125 y_t = ly[(n2 % (m * n)) / m];
1126 z_t = lz[n2 / (m * n)];
1127 s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1128 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1129 } else if (xe - Xe < 0 && Ys - ys < 0 && Zs - zs < 0) {
1130 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1131 }
1132 }
1133
1134 for (i = 0; i < y; i++) {
1135 if (n3 >= 0) { /* directly left */
1136 x_t = lx[n3 % m];
1137 y_t = y;
1138 z_t = lz[n3 / (m * n)];
1139 s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1140 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1141 } else if (Xs - xs < 0 && Zs - zs < 0) {
1142 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1143 }
1144
1145 if (n4 >= 0) { /* middle */
1146 x_t = x;
1147 y_t = y;
1148 z_t = lz[n4 / (m * n)];
1149 s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1150 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1151 } else if (Zs - zs < 0) {
1152 if (bz == DM_BOUNDARY_MIRROR) {
1153 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k) * x * y;
1154 } else {
1155 for (j = 0; j < x; j++) idx[nn++] = -1;
1156 }
1157 }
1158
1159 if (n5 >= 0) { /* directly right */
1160 x_t = lx[n5 % m];
1161 y_t = y;
1162 z_t = lz[n5 / (m * n)];
1163 s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1164 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1165 } else if (xe - Xe < 0 && Zs - zs < 0) {
1166 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1167 }
1168 }
1169
1170 for (i = 1; i <= s_y; i++) {
1171 if (n6 >= 0) { /* left above */
1172 x_t = lx[n6 % m];
1173 y_t = ly[(n6 % (m * n)) / m];
1174 z_t = lz[n6 / (m * n)];
1175 s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1176 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1177 } else if (Xs - xs < 0 && ye - Ye < 0 && Zs - zs < 0) {
1178 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1179 }
1180 if (n7 >= 0) { /* directly above */
1181 x_t = x;
1182 y_t = ly[(n7 % (m * n)) / m];
1183 z_t = lz[n7 / (m * n)];
1184 s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1185 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1186 } else if (ye - Ye < 0 && Zs - zs < 0) {
1187 for (j = 0; j < x; j++) idx[nn++] = -1;
1188 }
1189 if (n8 >= 0) { /* right above */
1190 x_t = lx[n8 % m];
1191 y_t = ly[(n8 % (m * n)) / m];
1192 z_t = lz[n8 / (m * n)];
1193 s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1194 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1195 } else if (xe - Xe < 0 && ye - Ye < 0 && Zs - zs < 0) {
1196 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1197 }
1198 }
1199 }
1200
1201 /* Middle Level */
1202 for (k = 0; k < z; k++) {
1203 for (i = 1; i <= s_y; i++) {
1204 if (n9 >= 0) { /* left below */
1205 x_t = lx[n9 % m];
1206 y_t = ly[(n9 % (m * n)) / m];
1207 /* z_t = z; */
1208 s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1209 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1210 } else if (Xs - xs < 0 && Ys - ys < 0) {
1211 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1212 }
1213 if (n10 >= 0) { /* directly below */
1214 x_t = x;
1215 y_t = ly[(n10 % (m * n)) / m];
1216 /* z_t = z; */
1217 s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1218 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1219 } else if (Ys - ys < 0) {
1220 if (by == DM_BOUNDARY_MIRROR) {
1221 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i + 1) * x;
1222 } else {
1223 for (j = 0; j < x; j++) idx[nn++] = -1;
1224 }
1225 }
1226 if (n11 >= 0) { /* right below */
1227 x_t = lx[n11 % m];
1228 y_t = ly[(n11 % (m * n)) / m];
1229 /* z_t = z; */
1230 s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1231 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1232 } else if (xe - Xe < 0 && Ys - ys < 0) {
1233 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1234 }
1235 }
1236
1237 for (i = 0; i < y; i++) {
1238 if (n12 >= 0) { /* directly left */
1239 x_t = lx[n12 % m];
1240 y_t = y;
1241 /* z_t = z; */
1242 s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
1243 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1244 } else if (Xs - xs < 0) {
1245 if (bx == DM_BOUNDARY_MIRROR) {
1246 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x + 1;
1247 } else {
1248 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1249 }
1250 }
1251
1252 /* Interior */
1253 s_t = bases[rank] + i * x + k * x * y;
1254 for (j = 0; j < x; j++) idx[nn++] = s_t++;
1255
1256 if (n14 >= 0) { /* directly right */
1257 x_t = lx[n14 % m];
1258 y_t = y;
1259 /* z_t = z; */
1260 s_t = bases[n14] + i * x_t + k * x_t * y_t;
1261 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1262 } else if (xe - Xe < 0) {
1263 if (bx == DM_BOUNDARY_MIRROR) {
1264 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x - 1;
1265 } else {
1266 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1267 }
1268 }
1269 }
1270
1271 for (i = 1; i <= s_y; i++) {
1272 if (n15 >= 0) { /* left above */
1273 x_t = lx[n15 % m];
1274 y_t = ly[(n15 % (m * n)) / m];
1275 /* z_t = z; */
1276 s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
1277 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1278 } else if (Xs - xs < 0 && ye - Ye < 0) {
1279 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1280 }
1281 if (n16 >= 0) { /* directly above */
1282 x_t = x;
1283 y_t = ly[(n16 % (m * n)) / m];
1284 /* z_t = z; */
1285 s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
1286 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1287 } else if (ye - Ye < 0) {
1288 if (by == DM_BOUNDARY_MIRROR) {
1289 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i - 1) * x;
1290 } else {
1291 for (j = 0; j < x; j++) idx[nn++] = -1;
1292 }
1293 }
1294 if (n17 >= 0) { /* right above */
1295 x_t = lx[n17 % m];
1296 y_t = ly[(n17 % (m * n)) / m];
1297 /* z_t = z; */
1298 s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
1299 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1300 } else if (xe - Xe < 0 && ye - Ye < 0) {
1301 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1302 }
1303 }
1304 }
1305
1306 /* Upper Level */
1307 for (k = 0; k < s_z; k++) {
1308 for (i = 1; i <= s_y; i++) {
1309 if (n18 >= 0) { /* left below */
1310 x_t = lx[n18 % m];
1311 y_t = ly[(n18 % (m * n)) / m];
1312 /* z_t = lz[n18 / (m*n)]; */
1313 s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1314 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1315 } else if (Xs - xs < 0 && Ys - ys < 0 && ze - Ze < 0) {
1316 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1317 }
1318 if (n19 >= 0) { /* directly below */
1319 x_t = x;
1320 y_t = ly[(n19 % (m * n)) / m];
1321 /* z_t = lz[n19 / (m*n)]; */
1322 s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1323 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1324 } else if (Ys - ys < 0 && ze - Ze < 0) {
1325 for (j = 0; j < x; j++) idx[nn++] = -1;
1326 }
1327 if (n20 >= 0) { /* right below */
1328 x_t = lx[n20 % m];
1329 y_t = ly[(n20 % (m * n)) / m];
1330 /* z_t = lz[n20 / (m*n)]; */
1331 s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1332 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1333 } else if (xe - Xe < 0 && Ys - ys < 0 && ze - Ze < 0) {
1334 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1335 }
1336 }
1337
1338 for (i = 0; i < y; i++) {
1339 if (n21 >= 0) { /* directly left */
1340 x_t = lx[n21 % m];
1341 y_t = y;
1342 /* z_t = lz[n21 / (m*n)]; */
1343 s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1344 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1345 } else if (Xs - xs < 0 && ze - Ze < 0) {
1346 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1347 }
1348
1349 if (n22 >= 0) { /* middle */
1350 x_t = x;
1351 y_t = y;
1352 /* z_t = lz[n22 / (m*n)]; */
1353 s_t = bases[n22] + i * x_t + k * x_t * y_t;
1354 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1355 } else if (ze - Ze < 0) {
1356 if (bz == DM_BOUNDARY_MIRROR) {
1357 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 2) * x * y + i * x;
1358 } else {
1359 for (j = 0; j < x; j++) idx[nn++] = -1;
1360 }
1361 }
1362
1363 if (n23 >= 0) { /* directly right */
1364 x_t = lx[n23 % m];
1365 y_t = y;
1366 /* z_t = lz[n23 / (m*n)]; */
1367 s_t = bases[n23] + i * x_t + k * x_t * y_t;
1368 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1369 } else if (xe - Xe < 0 && ze - Ze < 0) {
1370 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1371 }
1372 }
1373
1374 for (i = 1; i <= s_y; i++) {
1375 if (n24 >= 0) { /* left above */
1376 x_t = lx[n24 % m];
1377 y_t = ly[(n24 % (m * n)) / m];
1378 /* z_t = lz[n24 / (m*n)]; */
1379 s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1380 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1381 } else if (Xs - xs < 0 && ye - Ye < 0 && ze - Ze < 0) {
1382 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1383 }
1384 if (n25 >= 0) { /* directly above */
1385 x_t = x;
1386 y_t = ly[(n25 % (m * n)) / m];
1387 /* z_t = lz[n25 / (m*n)]; */
1388 s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1389 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1390 } else if (ye - Ye < 0 && ze - Ze < 0) {
1391 for (j = 0; j < x; j++) idx[nn++] = -1;
1392 }
1393 if (n26 >= 0) { /* right above */
1394 x_t = lx[n26 % m];
1395 y_t = ly[(n26 % (m * n)) / m];
1396 /* z_t = lz[n26 / (m*n)]; */
1397 s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1398 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1399 } else if (xe - Xe < 0 && ye - Ye < 0 && ze - Ze < 0) {
1400 for (j = 0; j < s_x; j++) idx[nn++] = -1;
1401 }
1402 }
1403 }
1404 }
1405 /*
1406 Set the local to global ordering in the global vector, this allows use
1407 of VecSetValuesLocal().
1408 */
1409 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap));
1410
1411 PetscCall(PetscFree2(bases, ldims));
1412 dd->m = m;
1413 dd->n = n;
1414 dd->p = p;
1415 /* note PETSc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1416 dd->xs = xs * dof;
1417 dd->xe = xe * dof;
1418 dd->ys = ys;
1419 dd->ye = ye;
1420 dd->zs = zs;
1421 dd->ze = ze;
1422 dd->Xs = Xs * dof;
1423 dd->Xe = Xe * dof;
1424 dd->Ys = Ys;
1425 dd->Ye = Ye;
1426 dd->Zs = Zs;
1427 dd->Ze = Ze;
1428
1429 PetscCall(VecDestroy(&local));
1430 PetscCall(VecDestroy(&global));
1431
1432 dd->gtol = gtol;
1433 dd->base = base;
1434 da->ops->view = DMView_DA_3d;
1435 dd->ltol = NULL;
1436 dd->ao = NULL;
1437 PetscFunctionReturn(PETSC_SUCCESS);
1438 }
1439
1440 /*@
1441 DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1442 regular array data that is distributed across one or more MPI processes.
1443
1444 Collective
1445
1446 Input Parameters:
1447 + comm - MPI communicator
1448 . bx - type of x ghost nodes the array have.
1449 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1450 . by - type of y ghost nodes the array have.
1451 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1452 . bz - type of z ghost nodes the array have.
1453 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1454 . stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`)
1455 . M - global dimension in x direction of the array
1456 . N - global dimension in y direction of the array
1457 . P - global dimension in z direction of the array
1458 . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated)
1459 . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated)
1460 . p - corresponding number of processors in z dimension (or `PETSC_DECIDE` to have calculated)
1461 . dof - number of degrees of freedom per node
1462 . s - stencil width
1463 . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`.
1464 . ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`.
1465 - lz - arrays containing the number of nodes in each cell along the z coordinates, or `NULL`.
1466
1467 Output Parameter:
1468 . da - the resulting distributed array object
1469
1470 Options Database Keys:
1471 + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()`
1472 . -da_grid_x <nx> - number of grid points in x direction
1473 . -da_grid_y <ny> - number of grid points in y direction
1474 . -da_grid_z <nz> - number of grid points in z direction
1475 . -da_processors_x <MX> - number of processors in x direction
1476 . -da_processors_y <MY> - number of processors in y direction
1477 . -da_processors_z <MZ> - number of processors in z direction
1478 . -da_bd_x <bx> - boundary type in x direction
1479 . -da_bd_y <by> - boundary type in y direction
1480 . -da_bd_z <bz> - boundary type in x direction
1481 . -da_bd_all <bt> - boundary type in all directions
1482 . -da_refine_x <rx> - refinement ratio in x direction
1483 . -da_refine_y <ry> - refinement ratio in y direction
1484 . -da_refine_z <rz> - refinement ratio in z directio
1485 - -da_refine <n> - refine the `DMDA` n times before creating it
1486
1487 Level: beginner
1488
1489 Notes:
1490 If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the
1491 corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`,
1492 sum of the `ly` must `N`, sum of the `lz` must be `P`.
1493
1494 The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1495 standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1496 the standard 27-pt stencil.
1497
1498 The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
1499 The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
1500 and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed.
1501
1502 You must call `DMSetUp()` after this call before using this `DM`.
1503
1504 To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1505 but before `DMSetUp()`.
1506
1507 .seealso: [](sec_struct), `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
1508 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
1509 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
1510 `DMStagCreate3d()`, `DMBoundaryType`
1511 @*/
DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M,PetscInt N,PetscInt P,PetscInt m,PetscInt n,PetscInt p,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[],DM * da)1512 PetscErrorCode DMDACreate3d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt P, PetscInt m, PetscInt n, PetscInt p, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[], DM *da)
1513 {
1514 PetscFunctionBegin;
1515 PetscCall(DMDACreate(comm, da));
1516 PetscCall(DMSetDimension(*da, 3));
1517 PetscCall(DMDASetSizes(*da, M, N, P));
1518 PetscCall(DMDASetNumProcs(*da, m, n, p));
1519 PetscCall(DMDASetBoundaryType(*da, bx, by, bz));
1520 PetscCall(DMDASetDof(*da, dof));
1521 PetscCall(DMDASetStencilType(*da, stencil_type));
1522 PetscCall(DMDASetStencilWidth(*da, s));
1523 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz));
1524 PetscFunctionReturn(PETSC_SUCCESS);
1525 }
1526