Actual source code: da3.c
1: /*
2: Code for manipulating distributed regular 3d arrays in parallel.
3: File created by Peter Mell 7/14/95
4: */
6: #include <petsc/private/dmdaimpl.h>
8: #include <petscdraw.h>
9: static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer)
10: {
11: PetscMPIInt rank;
12: PetscBool iascii, 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
19: PetscFunctionBegin;
20: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));
22: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
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 (iascii) {
30: PetscViewerFormat format;
32: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
33: PetscCall(PetscViewerGetFormat(viewer, &format));
34: if (format == PETSC_VIEWER_LOAD_BALANCE) {
35: PetscInt i, nmax = 0, nmin = PETSC_MAX_INT, 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 < (PetscInt)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_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && 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 PetscCall(DMView_DA_VTK(da, viewer));
75: } else if (isdraw) {
76: PetscDraw draw;
77: PetscReal ymin = -1.0, ymax = (PetscReal)dd->N;
78: PetscReal xmin = -1.0, xmax = (PetscReal)((dd->M + 2) * dd->P), x, y, ycoord, xcoord;
79: PetscInt k, plane, base;
80: const PetscInt *idx;
81: char node[10];
82: PetscBool isnull;
84: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
85: PetscCall(PetscDrawIsNull(draw, &isnull));
86: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
88: PetscCall(PetscDrawCheckResizedWindow(draw));
89: PetscCall(PetscDrawClear(draw));
90: PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));
92: PetscDrawCollectiveBegin(draw);
93: /* first processor draw all node lines */
94: if (rank == 0) {
95: for (k = 0; k < dd->P; k++) {
96: ymin = 0.0;
97: ymax = (PetscReal)(dd->N - 1);
98: 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));
99: xmin = (PetscReal)(k * (dd->M + 1));
100: xmax = xmin + (PetscReal)(dd->M - 1);
101: for (ymin = 0; ymin < (PetscReal)dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
102: }
103: }
104: PetscDrawCollectiveEnd(draw);
105: PetscCall(PetscDrawFlush(draw));
106: PetscCall(PetscDrawPause(draw));
108: PetscDrawCollectiveBegin(draw);
109: /*Go through and draw for each plane*/
110: for (k = 0; k < dd->P; k++) {
111: if ((k >= dd->zs) && (k < dd->ze)) {
112: /* draw my box */
113: ymin = dd->ys;
114: ymax = dd->ye - 1;
115: xmin = dd->xs / dd->w + (dd->M + 1) * k;
116: xmax = (dd->xe - 1) / dd->w + (dd->M + 1) * k;
118: PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
119: PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
120: PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
121: PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED));
123: xmin = dd->xs / dd->w;
124: xmax = (dd->xe - 1) / dd->w;
126: /* identify which processor owns the box */
127: PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)rank));
128: PetscCall(PetscDrawString(draw, xmin + (dd->M + 1) * k + .2, ymin + .3, PETSC_DRAW_RED, node));
129: /* put in numbers*/
130: base = (dd->base + (dd->xe - dd->xs) * (dd->ye - dd->ys) * (k - dd->zs)) / dd->w;
131: for (y = ymin; y <= ymax; y++) {
132: for (x = xmin + (dd->M + 1) * k; x <= xmax + (dd->M + 1) * k; x++) {
133: PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++));
134: PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
135: }
136: }
137: }
138: }
139: PetscDrawCollectiveEnd(draw);
140: PetscCall(PetscDrawFlush(draw));
141: PetscCall(PetscDrawPause(draw));
143: PetscDrawCollectiveBegin(draw);
144: for (k = 0 - dd->s; k < dd->P + dd->s; k++) {
145: /* Go through and draw for each plane */
146: if ((k >= dd->Zs) && (k < dd->Ze)) {
147: /* overlay ghost numbers, useful for error checking */
148: base = (dd->Xe - dd->Xs) * (dd->Ye - dd->Ys) * (k - dd->Zs) / dd->w;
149: PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx));
150: plane = k;
151: /* Keep z wrap around points on the drawing */
152: if (k < 0) plane = dd->P + k;
153: if (k >= dd->P) plane = k - dd->P;
154: ymin = dd->Ys;
155: ymax = dd->Ye;
156: xmin = (dd->M + 1) * plane * dd->w;
157: xmax = (dd->M + 1) * plane * dd->w + dd->M * dd->w;
158: for (y = ymin; y < ymax; y++) {
159: for (x = xmin + dd->Xs; x < xmin + dd->Xe; x += dd->w) {
160: PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, idx[base]));
161: ycoord = y;
162: /*Keep y wrap around points on drawing */
163: if (y < 0) ycoord = dd->N + y;
164: if (y >= dd->N) ycoord = y - dd->N;
165: xcoord = x; /* Keep x wrap points on drawing */
166: if (x < xmin) xcoord = xmax - (xmin - x);
167: if (x >= xmax) xcoord = xmin + (x - xmax);
168: PetscCall(PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node));
169: base++;
170: }
171: }
172: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
173: }
174: }
175: PetscDrawCollectiveEnd(draw);
176: PetscCall(PetscDrawFlush(draw));
177: PetscCall(PetscDrawPause(draw));
178: PetscCall(PetscDrawSave(draw));
179: } else if (isglvis) {
180: PetscCall(DMView_DA_GLVis(da, viewer));
181: } else if (isbinary) {
182: PetscCall(DMView_DA_Binary(da, viewer));
183: #if defined(PETSC_HAVE_MATLAB)
184: } else if (ismatlab) {
185: PetscCall(DMView_DA_Matlab(da, viewer));
186: #endif
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode DMSetUp_DA_3D(DM da)
192: {
193: DM_DA *dd = (DM_DA *)da->data;
194: const PetscInt M = dd->M;
195: const PetscInt N = dd->N;
196: const PetscInt P = dd->P;
197: PetscInt m = dd->m;
198: PetscInt n = dd->n;
199: PetscInt p = dd->p;
200: const PetscInt dof = dd->w;
201: const PetscInt s = dd->s;
202: DMBoundaryType bx = dd->bx;
203: DMBoundaryType by = dd->by;
204: DMBoundaryType bz = dd->bz;
205: DMDAStencilType stencil_type = dd->stencil_type;
206: PetscInt *lx = dd->lx;
207: PetscInt *ly = dd->ly;
208: PetscInt *lz = dd->lz;
209: MPI_Comm comm;
210: PetscMPIInt rank, size;
211: PetscInt xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0;
212: PetscInt Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm;
213: PetscInt left, right, up, down, bottom, top, i, j, k, *idx, nn;
214: PetscInt n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14;
215: PetscInt n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26;
216: PetscInt *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z;
217: PetscInt sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0;
218: PetscInt sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0;
219: PetscInt sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0;
220: Vec local, global;
221: VecScatter gtol;
222: IS to, from;
223: PetscBool twod;
225: PetscFunctionBegin;
226: 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");
227: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
228: #if !defined(PETSC_USE_64BIT_INDICES)
229: 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);
230: #endif
232: PetscCallMPI(MPI_Comm_size(comm, &size));
233: PetscCallMPI(MPI_Comm_rank(comm, &rank));
235: if (m != PETSC_DECIDE) {
236: PetscCheck(m >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m);
237: PetscCheck(m <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %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: %" PetscInt_FMT, n);
241: PetscCheck(n <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %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: %" PetscInt_FMT, p);
245: PetscCheck(p <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Z direction: %" PetscInt_FMT " %d", p, size);
246: }
247: 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);
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 = %" PetscInt_FMT, p);
266: if (M > N && m < n) {
267: PetscInt _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 = %" PetscInt_FMT, n);
281: if (M > P && m < p) {
282: PetscInt _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 = %" PetscInt_FMT, n);
296: if (N > P && n < p) {
297: PetscInt _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 = (PetscInt)(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 = (PetscInt)(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: PetscInt _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");
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 " %" PetscInt_FMT, M, m);
327: PetscCheck(N >= n, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n);
328: PetscCheck(P >= p, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in z direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, P, p);
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: */
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);
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);
353: ys = 0;
354: for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i];
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)];
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);
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;
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: }
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: }
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: }
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: }
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: }
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: }
454: /* Resize all X parameters to reflect w */
455: s_x = s;
456: s_y = s;
457: s_z = s;
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;
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));
474: /* generate global to local vector scatter and local to global mapping*/
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: }
532: /* determine who lies on each side of use stored in n24 n25 n26
533: n21 n22 n23
534: n18 n19 n20
536: n15 n16 n17
537: n12 n14
538: n9 n10 n11
540: n6 n7 n8
541: n3 n4 n5
542: n0 n1 n2
543: */
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;
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;
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;
576: /* Assume Pieces are on Faces of Cube */
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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;
732: /* Check for when not X,Y, and Z Periodic */
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: }
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: }
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: }
752: PetscCall(PetscMalloc1(27, &dd->neighbors));
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;
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: }
808: PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx));
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: }
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: }
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: }
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: }
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: }
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: }
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: }
938: /* Interior */
939: s_t = bases[rank] + i * x + k * x * y;
940: for (j = 0; j < x; j++) idx[nn++] = s_t++;
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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));
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }
1252: /* Interior */
1253: s_t = bases[rank] + i * x + k * x * y;
1254: for (j = 0; j < x; j++) idx[nn++] = s_t++;
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: }
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: }
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: }
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: }
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: }
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: }
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));
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;
1429: PetscCall(VecDestroy(&local));
1430: PetscCall(VecDestroy(&global));
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: }
1440: /*@C
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.
1444: Collective
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`.
1467: Output Parameter:
1468: . da - the resulting distributed array object
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_refine_x <rx> - refinement ratio in x direction
1479: . -da_refine_y <ry> - refinement ratio in y direction
1480: . -da_refine_z <rz> - refinement ratio in z directio
1481: - -da_refine <n> - refine the `DMDA` n times before creating it
1483: Level: beginner
1485: Notes:
1486: If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the
1487: corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`,
1488: sum of the `ly` must `N`, sum of the `lz` must be `P`.
1490: The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1491: standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1492: the standard 27-pt stencil.
1494: The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
1495: The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
1496: and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed.
1498: You must call `DMSetUp()` after this call before using this `DM`.
1500: To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1501: but before `DMSetUp()`.
1503: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
1504: `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
1505: `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
1506: `DMStagCreate3d()`, `DMBoundaryType`
1507: @*/
1508: 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)
1509: {
1510: PetscFunctionBegin;
1511: PetscCall(DMDACreate(comm, da));
1512: PetscCall(DMSetDimension(*da, 3));
1513: PetscCall(DMDASetSizes(*da, M, N, P));
1514: PetscCall(DMDASetNumProcs(*da, m, n, p));
1515: PetscCall(DMDASetBoundaryType(*da, bx, by, bz));
1516: PetscCall(DMDASetDof(*da, dof));
1517: PetscCall(DMDASetStencilType(*da, stencil_type));
1518: PetscCall(DMDASetStencilWidth(*da, s));
1519: PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz));
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }