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: /* note this is different than x- and y-, as we will handle as an important special
364: case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems
365: in a 3D code. Additional code for this case is noted with "2d case" comments */
366: twod = PETSC_FALSE;
367: if (P == 1) twod = PETSC_TRUE;
368: 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);
369: zs = 0;
370: for (i = 0; i < (rank / (m * n)); i++) zs += lz[i];
371: ye = ys + y;
372: xe = xs + x;
373: ze = zs + z;
375: /* determine ghost region (Xs) and region scattered into (IXs) */
376: if (xs - s > 0) {
377: Xs = xs - s;
378: IXs = xs - s;
379: } else {
380: if (bx) Xs = xs - s;
381: else Xs = 0;
382: IXs = 0;
383: }
384: if (xe + s <= M) {
385: Xe = xe + s;
386: IXe = xe + s;
387: } else {
388: if (bx) {
389: Xs = xs - s;
390: Xe = xe + s;
391: } else Xe = M;
392: IXe = M;
393: }
395: if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
396: IXs = xs - s;
397: IXe = xe + s;
398: Xs = xs - s;
399: Xe = xe + s;
400: }
402: if (ys - s > 0) {
403: Ys = ys - s;
404: IYs = ys - s;
405: } else {
406: if (by) Ys = ys - s;
407: else Ys = 0;
408: IYs = 0;
409: }
410: if (ye + s <= N) {
411: Ye = ye + s;
412: IYe = ye + s;
413: } else {
414: if (by) Ye = ye + s;
415: else Ye = N;
416: IYe = N;
417: }
419: if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
420: IYs = ys - s;
421: IYe = ye + s;
422: Ys = ys - s;
423: Ye = ye + s;
424: }
426: if (zs - s > 0) {
427: Zs = zs - s;
428: IZs = zs - s;
429: } else {
430: if (bz) Zs = zs - s;
431: else Zs = 0;
432: IZs = 0;
433: }
434: if (ze + s <= P) {
435: Ze = ze + s;
436: IZe = ze + s;
437: } else {
438: if (bz) Ze = ze + s;
439: else Ze = P;
440: IZe = P;
441: }
443: if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) {
444: IZs = zs - s;
445: IZe = ze + s;
446: Zs = zs - s;
447: Ze = ze + s;
448: }
450: /* Resize all X parameters to reflect w */
451: s_x = s;
452: s_y = s;
453: s_z = s;
455: /* determine starting point of each processor */
456: nn = x * y * z;
457: PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims));
458: PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm));
459: bases[0] = 0;
460: for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
461: for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
462: base = bases[rank] * dof;
464: /* allocate the base parallel and sequential vectors */
465: dd->Nlocal = x * y * z * dof;
466: PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global));
467: dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof;
468: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local));
470: /* generate global to local vector scatter and local to global mapping*/
472: /* global to local must include ghost points within the domain,
473: but not ghost points outside the domain that aren't periodic */
474: PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx));
475: if (stencil_type == DMDA_STENCIL_BOX) {
476: left = IXs - Xs;
477: right = left + (IXe - IXs);
478: bottom = IYs - Ys;
479: top = bottom + (IYe - IYs);
480: down = IZs - Zs;
481: up = down + (IZe - IZs);
482: count = 0;
483: for (i = down; i < up; i++) {
484: for (j = bottom; j < top; j++) {
485: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
486: }
487: }
488: PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
489: } else {
490: /* This is way ugly! We need to list the funny cross type region */
491: left = xs - Xs;
492: right = left + x;
493: bottom = ys - Ys;
494: top = bottom + y;
495: down = zs - Zs;
496: up = down + z;
497: count = 0;
498: /* the bottom chunk */
499: for (i = (IZs - Zs); i < down; i++) {
500: for (j = bottom; j < top; j++) {
501: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
502: }
503: }
504: /* the middle piece */
505: for (i = down; i < up; i++) {
506: /* front */
507: for (j = (IYs - Ys); j < bottom; j++) {
508: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
509: }
510: /* middle */
511: for (j = bottom; j < top; j++) {
512: for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
513: }
514: /* back */
515: for (j = top; j < top + IYe - ye; j++) {
516: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
517: }
518: }
519: /* the top piece */
520: for (i = up; i < up + IZe - ze; i++) {
521: for (j = bottom; j < top; j++) {
522: for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k;
523: }
524: }
525: PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
526: }
528: /* determine who lies on each side of use stored in n24 n25 n26
529: n21 n22 n23
530: n18 n19 n20
532: n15 n16 n17
533: n12 n14
534: n9 n10 n11
536: n6 n7 n8
537: n3 n4 n5
538: n0 n1 n2
539: */
541: /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
542: /* Assume Nodes are Internal to the Cube */
543: n0 = rank - m * n - m - 1;
544: n1 = rank - m * n - m;
545: n2 = rank - m * n - m + 1;
546: n3 = rank - m * n - 1;
547: n4 = rank - m * n;
548: n5 = rank - m * n + 1;
549: n6 = rank - m * n + m - 1;
550: n7 = rank - m * n + m;
551: n8 = rank - m * n + m + 1;
553: n9 = rank - m - 1;
554: n10 = rank - m;
555: n11 = rank - m + 1;
556: n12 = rank - 1;
557: n14 = rank + 1;
558: n15 = rank + m - 1;
559: n16 = rank + m;
560: n17 = rank + m + 1;
562: n18 = rank + m * n - m - 1;
563: n19 = rank + m * n - m;
564: n20 = rank + m * n - m + 1;
565: n21 = rank + m * n - 1;
566: n22 = rank + m * n;
567: n23 = rank + m * n + 1;
568: n24 = rank + m * n + m - 1;
569: n25 = rank + m * n + m;
570: n26 = rank + m * n + m + 1;
572: /* Assume Pieces are on Faces of Cube */
574: if (xs == 0) { /* First assume not corner or edge */
575: n0 = rank - 1 - (m * n);
576: n3 = rank + m - 1 - (m * n);
577: n6 = rank + 2 * m - 1 - (m * n);
578: n9 = rank - 1;
579: n12 = rank + m - 1;
580: n15 = rank + 2 * m - 1;
581: n18 = rank - 1 + (m * n);
582: n21 = rank + m - 1 + (m * n);
583: n24 = rank + 2 * m - 1 + (m * n);
584: }
586: if (xe == M) { /* First assume not corner or edge */
587: n2 = rank - 2 * m + 1 - (m * n);
588: n5 = rank - m + 1 - (m * n);
589: n8 = rank + 1 - (m * n);
590: n11 = rank - 2 * m + 1;
591: n14 = rank - m + 1;
592: n17 = rank + 1;
593: n20 = rank - 2 * m + 1 + (m * n);
594: n23 = rank - m + 1 + (m * n);
595: n26 = rank + 1 + (m * n);
596: }
598: if (ys == 0) { /* First assume not corner or edge */
599: n0 = rank + m * (n - 1) - 1 - (m * n);
600: n1 = rank + m * (n - 1) - (m * n);
601: n2 = rank + m * (n - 1) + 1 - (m * n);
602: n9 = rank + m * (n - 1) - 1;
603: n10 = rank + m * (n - 1);
604: n11 = rank + m * (n - 1) + 1;
605: n18 = rank + m * (n - 1) - 1 + (m * n);
606: n19 = rank + m * (n - 1) + (m * n);
607: n20 = rank + m * (n - 1) + 1 + (m * n);
608: }
610: if (ye == N) { /* First assume not corner or edge */
611: n6 = rank - m * (n - 1) - 1 - (m * n);
612: n7 = rank - m * (n - 1) - (m * n);
613: n8 = rank - m * (n - 1) + 1 - (m * n);
614: n15 = rank - m * (n - 1) - 1;
615: n16 = rank - m * (n - 1);
616: n17 = rank - m * (n - 1) + 1;
617: n24 = rank - m * (n - 1) - 1 + (m * n);
618: n25 = rank - m * (n - 1) + (m * n);
619: n26 = rank - m * (n - 1) + 1 + (m * n);
620: }
622: if (zs == 0) { /* First assume not corner or edge */
623: n0 = size - (m * n) + rank - m - 1;
624: n1 = size - (m * n) + rank - m;
625: n2 = size - (m * n) + rank - m + 1;
626: n3 = size - (m * n) + rank - 1;
627: n4 = size - (m * n) + rank;
628: n5 = size - (m * n) + rank + 1;
629: n6 = size - (m * n) + rank + m - 1;
630: n7 = size - (m * n) + rank + m;
631: n8 = size - (m * n) + rank + m + 1;
632: }
634: if (ze == P) { /* First assume not corner or edge */
635: n18 = (m * n) - (size - rank) - m - 1;
636: n19 = (m * n) - (size - rank) - m;
637: n20 = (m * n) - (size - rank) - m + 1;
638: n21 = (m * n) - (size - rank) - 1;
639: n22 = (m * n) - (size - rank);
640: n23 = (m * n) - (size - rank) + 1;
641: n24 = (m * n) - (size - rank) + m - 1;
642: n25 = (m * n) - (size - rank) + m;
643: n26 = (m * n) - (size - rank) + m + 1;
644: }
646: if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */
647: n0 = size - m * n + rank + m - 1 - m;
648: n3 = size - m * n + rank + m - 1;
649: n6 = size - m * n + rank + m - 1 + m;
650: }
652: if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */
653: n18 = m * n - (size - rank) + m - 1 - m;
654: n21 = m * n - (size - rank) + m - 1;
655: n24 = m * n - (size - rank) + m - 1 + m;
656: }
658: if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */
659: n0 = rank + m * n - 1 - m * n;
660: n9 = rank + m * n - 1;
661: n18 = rank + m * n - 1 + m * n;
662: }
664: if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */
665: n6 = rank - m * (n - 1) + m - 1 - m * n;
666: n15 = rank - m * (n - 1) + m - 1;
667: n24 = rank - m * (n - 1) + m - 1 + m * n;
668: }
670: if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */
671: n2 = size - (m * n - rank) - (m - 1) - m;
672: n5 = size - (m * n - rank) - (m - 1);
673: n8 = size - (m * n - rank) - (m - 1) + m;
674: }
676: if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */
677: n20 = m * n - (size - rank) - (m - 1) - m;
678: n23 = m * n - (size - rank) - (m - 1);
679: n26 = m * n - (size - rank) - (m - 1) + m;
680: }
682: if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */
683: n2 = rank + m * (n - 1) - (m - 1) - m * n;
684: n11 = rank + m * (n - 1) - (m - 1);
685: n20 = rank + m * (n - 1) - (m - 1) + m * n;
686: }
688: if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */
689: n8 = rank - m * n + 1 - m * n;
690: n17 = rank - m * n + 1;
691: n26 = rank - m * n + 1 + m * n;
692: }
694: if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */
695: n0 = size - m + rank - 1;
696: n1 = size - m + rank;
697: n2 = size - m + rank + 1;
698: }
700: if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */
701: n18 = m * n - (size - rank) + m * (n - 1) - 1;
702: n19 = m * n - (size - rank) + m * (n - 1);
703: n20 = m * n - (size - rank) + m * (n - 1) + 1;
704: }
706: if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */
707: n6 = size - (m * n - rank) - m * (n - 1) - 1;
708: n7 = size - (m * n - rank) - m * (n - 1);
709: n8 = size - (m * n - rank) - m * (n - 1) + 1;
710: }
712: if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */
713: n24 = rank - (size - m) - 1;
714: n25 = rank - (size - m);
715: n26 = rank - (size - m) + 1;
716: }
718: /* Check for Corners */
719: if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1;
720: if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1;
721: if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1);
722: if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1;
723: if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m;
724: if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m;
725: if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n;
726: if ((xe == M) && (ye == N) && (ze == P)) n26 = 0;
728: /* Check for when not X,Y, and Z Periodic */
730: /* If not X periodic */
731: if (bx != DM_BOUNDARY_PERIODIC) {
732: if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;
733: if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;
734: }
736: /* If not Y periodic */
737: if (by != DM_BOUNDARY_PERIODIC) {
738: if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;
739: if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;
740: }
742: /* If not Z periodic */
743: if (bz != DM_BOUNDARY_PERIODIC) {
744: if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;
745: if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;
746: }
748: PetscCall(PetscMalloc1(27, &dd->neighbors));
750: dd->neighbors[0] = n0;
751: dd->neighbors[1] = n1;
752: dd->neighbors[2] = n2;
753: dd->neighbors[3] = n3;
754: dd->neighbors[4] = n4;
755: dd->neighbors[5] = n5;
756: dd->neighbors[6] = n6;
757: dd->neighbors[7] = n7;
758: dd->neighbors[8] = n8;
759: dd->neighbors[9] = n9;
760: dd->neighbors[10] = n10;
761: dd->neighbors[11] = n11;
762: dd->neighbors[12] = n12;
763: dd->neighbors[13] = rank;
764: dd->neighbors[14] = n14;
765: dd->neighbors[15] = n15;
766: dd->neighbors[16] = n16;
767: dd->neighbors[17] = n17;
768: dd->neighbors[18] = n18;
769: dd->neighbors[19] = n19;
770: dd->neighbors[20] = n20;
771: dd->neighbors[21] = n21;
772: dd->neighbors[22] = n22;
773: dd->neighbors[23] = n23;
774: dd->neighbors[24] = n24;
775: dd->neighbors[25] = n25;
776: dd->neighbors[26] = n26;
778: /* If star stencil then delete the corner neighbors */
779: if (stencil_type == DMDA_STENCIL_STAR) {
780: /* save information about corner neighbors */
781: sn0 = n0;
782: sn1 = n1;
783: sn2 = n2;
784: sn3 = n3;
785: sn5 = n5;
786: sn6 = n6;
787: sn7 = n7;
788: sn8 = n8;
789: sn9 = n9;
790: sn11 = n11;
791: sn15 = n15;
792: sn17 = n17;
793: sn18 = n18;
794: sn19 = n19;
795: sn20 = n20;
796: sn21 = n21;
797: sn23 = n23;
798: sn24 = n24;
799: sn25 = n25;
800: sn26 = n26;
801: n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
802: }
804: PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx));
806: nn = 0;
807: /* Bottom Level */
808: for (k = 0; k < s_z; k++) {
809: for (i = 1; i <= s_y; i++) {
810: if (n0 >= 0) { /* left below */
811: x_t = lx[n0 % m];
812: y_t = ly[(n0 % (m * n)) / m];
813: z_t = lz[n0 / (m * n)];
814: 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;
815: if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */
816: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
817: }
818: if (n1 >= 0) { /* directly below */
819: x_t = x;
820: y_t = ly[(n1 % (m * n)) / m];
821: z_t = lz[n1 / (m * n)];
822: s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
823: if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
824: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
825: }
826: if (n2 >= 0) { /* right below */
827: x_t = lx[n2 % m];
828: y_t = ly[(n2 % (m * n)) / m];
829: z_t = lz[n2 / (m * n)];
830: s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
831: if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */
832: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
833: }
834: }
836: for (i = 0; i < y; i++) {
837: if (n3 >= 0) { /* directly left */
838: x_t = lx[n3 % m];
839: y_t = y;
840: z_t = lz[n3 / (m * n)];
841: s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
842: 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 */
843: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
844: }
846: if (n4 >= 0) { /* middle */
847: x_t = x;
848: y_t = y;
849: z_t = lz[n4 / (m * n)];
850: s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
851: if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
852: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
853: } else if (bz == DM_BOUNDARY_MIRROR) {
854: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
855: }
857: if (n5 >= 0) { /* directly right */
858: x_t = lx[n5 % m];
859: y_t = y;
860: z_t = lz[n5 / (m * n)];
861: s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
862: if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */
863: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
864: }
865: }
867: for (i = 1; i <= s_y; i++) {
868: if (n6 >= 0) { /* left above */
869: x_t = lx[n6 % m];
870: y_t = ly[(n6 % (m * n)) / m];
871: z_t = lz[n6 / (m * n)];
872: s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
873: 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 */
874: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
875: }
876: if (n7 >= 0) { /* directly above */
877: x_t = x;
878: y_t = ly[(n7 % (m * n)) / m];
879: z_t = lz[n7 / (m * n)];
880: s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
881: 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 */
882: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
883: }
884: if (n8 >= 0) { /* right above */
885: x_t = lx[n8 % m];
886: y_t = ly[(n8 % (m * n)) / m];
887: z_t = lz[n8 / (m * n)];
888: s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
889: 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 */
890: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
891: }
892: }
893: }
895: /* Middle Level */
896: for (k = 0; k < z; k++) {
897: for (i = 1; i <= s_y; i++) {
898: if (n9 >= 0) { /* left below */
899: x_t = lx[n9 % m];
900: y_t = ly[(n9 % (m * n)) / m];
901: /* z_t = z; */
902: s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
903: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
904: }
905: if (n10 >= 0) { /* directly below */
906: x_t = x;
907: y_t = ly[(n10 % (m * n)) / m];
908: /* z_t = z; */
909: s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
910: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
911: } else if (by == DM_BOUNDARY_MIRROR) {
912: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
913: }
914: if (n11 >= 0) { /* right below */
915: x_t = lx[n11 % m];
916: y_t = ly[(n11 % (m * n)) / m];
917: /* z_t = z; */
918: s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
919: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
920: }
921: }
923: for (i = 0; i < y; i++) {
924: if (n12 >= 0) { /* directly left */
925: x_t = lx[n12 % m];
926: y_t = y;
927: /* z_t = z; */
928: s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
929: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
930: } else if (bx == DM_BOUNDARY_MIRROR) {
931: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
932: }
934: /* Interior */
935: s_t = bases[rank] + i * x + k * x * y;
936: for (j = 0; j < x; j++) idx[nn++] = s_t++;
938: if (n14 >= 0) { /* directly right */
939: x_t = lx[n14 % m];
940: y_t = y;
941: /* z_t = z; */
942: s_t = bases[n14] + i * x_t + k * x_t * y_t;
943: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
944: } else if (bx == DM_BOUNDARY_MIRROR) {
945: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
946: }
947: }
949: for (i = 1; i <= s_y; i++) {
950: if (n15 >= 0) { /* left above */
951: x_t = lx[n15 % m];
952: y_t = ly[(n15 % (m * n)) / m];
953: /* z_t = z; */
954: s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
955: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
956: }
957: if (n16 >= 0) { /* directly above */
958: x_t = x;
959: y_t = ly[(n16 % (m * n)) / m];
960: /* z_t = z; */
961: s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
962: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
963: } else if (by == DM_BOUNDARY_MIRROR) {
964: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
965: }
966: if (n17 >= 0) { /* right above */
967: x_t = lx[n17 % m];
968: y_t = ly[(n17 % (m * n)) / m];
969: /* z_t = z; */
970: s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
971: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
972: }
973: }
974: }
976: /* Upper Level */
977: for (k = 0; k < s_z; k++) {
978: for (i = 1; i <= s_y; i++) {
979: if (n18 >= 0) { /* left below */
980: x_t = lx[n18 % m];
981: y_t = ly[(n18 % (m * n)) / m];
982: /* z_t = lz[n18 / (m*n)]; */
983: s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
984: if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */
985: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
986: }
987: if (n19 >= 0) { /* directly below */
988: x_t = x;
989: y_t = ly[(n19 % (m * n)) / m];
990: /* z_t = lz[n19 / (m*n)]; */
991: s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
992: if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
993: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
994: }
995: if (n20 >= 0) { /* right below */
996: x_t = lx[n20 % m];
997: y_t = ly[(n20 % (m * n)) / m];
998: /* z_t = lz[n20 / (m*n)]; */
999: s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1000: if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */
1001: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1002: }
1003: }
1005: for (i = 0; i < y; i++) {
1006: if (n21 >= 0) { /* directly left */
1007: x_t = lx[n21 % m];
1008: y_t = y;
1009: /* z_t = lz[n21 / (m*n)]; */
1010: s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1011: if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */
1012: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1013: }
1015: if (n22 >= 0) { /* middle */
1016: x_t = x;
1017: y_t = y;
1018: /* z_t = lz[n22 / (m*n)]; */
1019: s_t = bases[n22] + i * x_t + k * x_t * y_t;
1020: if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */
1021: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1022: } else if (bz == DM_BOUNDARY_MIRROR) {
1023: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1024: }
1026: if (n23 >= 0) { /* directly right */
1027: x_t = lx[n23 % m];
1028: y_t = y;
1029: /* z_t = lz[n23 / (m*n)]; */
1030: s_t = bases[n23] + i * x_t + k * x_t * y_t;
1031: if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */
1032: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1033: }
1034: }
1036: for (i = 1; i <= s_y; i++) {
1037: if (n24 >= 0) { /* left above */
1038: x_t = lx[n24 % m];
1039: y_t = ly[(n24 % (m * n)) / m];
1040: /* z_t = lz[n24 / (m*n)]; */
1041: s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1042: if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */
1043: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1044: }
1045: if (n25 >= 0) { /* directly above */
1046: x_t = x;
1047: y_t = ly[(n25 % (m * n)) / m];
1048: /* z_t = lz[n25 / (m*n)]; */
1049: s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1050: if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */
1051: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1052: }
1053: if (n26 >= 0) { /* right above */
1054: x_t = lx[n26 % m];
1055: y_t = ly[(n26 % (m * n)) / m];
1056: /* z_t = lz[n26 / (m*n)]; */
1057: s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1058: if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */
1059: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1060: }
1061: }
1062: }
1064: PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
1065: PetscCall(VecScatterCreate(global, from, local, to, >ol));
1066: PetscCall(ISDestroy(&to));
1067: PetscCall(ISDestroy(&from));
1069: if (stencil_type == DMDA_STENCIL_STAR) {
1070: n0 = sn0;
1071: n1 = sn1;
1072: n2 = sn2;
1073: n3 = sn3;
1074: n5 = sn5;
1075: n6 = sn6;
1076: n7 = sn7;
1077: n8 = sn8;
1078: n9 = sn9;
1079: n11 = sn11;
1080: n15 = sn15;
1081: n17 = sn17;
1082: n18 = sn18;
1083: n19 = sn19;
1084: n20 = sn20;
1085: n21 = sn21;
1086: n23 = sn23;
1087: n24 = sn24;
1088: n25 = sn25;
1089: n26 = sn26;
1090: }
1092: if (((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz))) {
1093: /*
1094: Recompute the local to global mappings, this time keeping the
1095: information about the cross corner processor numbers.
1096: */
1097: nn = 0;
1098: /* Bottom Level */
1099: for (k = 0; k < s_z; k++) {
1100: for (i = 1; i <= s_y; i++) {
1101: if (n0 >= 0) { /* left below */
1102: x_t = lx[n0 % m];
1103: y_t = ly[(n0 % (m * n)) / m];
1104: z_t = lz[n0 / (m * n)];
1105: 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;
1106: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1107: } else if (Xs - xs < 0 && Ys - ys < 0 && Zs - zs < 0) {
1108: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1109: }
1110: if (n1 >= 0) { /* directly below */
1111: x_t = x;
1112: y_t = ly[(n1 % (m * n)) / m];
1113: z_t = lz[n1 / (m * n)];
1114: s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1115: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1116: } else if (Ys - ys < 0 && Zs - zs < 0) {
1117: for (j = 0; j < x; j++) idx[nn++] = -1;
1118: }
1119: if (n2 >= 0) { /* right below */
1120: x_t = lx[n2 % m];
1121: y_t = ly[(n2 % (m * n)) / m];
1122: z_t = lz[n2 / (m * n)];
1123: s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t;
1124: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1125: } else if (xe - Xe < 0 && Ys - ys < 0 && Zs - zs < 0) {
1126: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1127: }
1128: }
1130: for (i = 0; i < y; i++) {
1131: if (n3 >= 0) { /* directly left */
1132: x_t = lx[n3 % m];
1133: y_t = y;
1134: z_t = lz[n3 / (m * n)];
1135: s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1136: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1137: } else if (Xs - xs < 0 && Zs - zs < 0) {
1138: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1139: }
1141: if (n4 >= 0) { /* middle */
1142: x_t = x;
1143: y_t = y;
1144: z_t = lz[n4 / (m * n)];
1145: s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1146: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1147: } else if (Zs - zs < 0) {
1148: if (bz == DM_BOUNDARY_MIRROR) {
1149: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y;
1150: } else {
1151: for (j = 0; j < x; j++) idx[nn++] = -1;
1152: }
1153: }
1155: if (n5 >= 0) { /* directly right */
1156: x_t = lx[n5 % m];
1157: y_t = y;
1158: z_t = lz[n5 / (m * n)];
1159: s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1160: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1161: } else if (xe - Xe < 0 && Zs - zs < 0) {
1162: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1163: }
1164: }
1166: for (i = 1; i <= s_y; i++) {
1167: if (n6 >= 0) { /* left above */
1168: x_t = lx[n6 % m];
1169: y_t = ly[(n6 % (m * n)) / m];
1170: z_t = lz[n6 / (m * n)];
1171: s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1172: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1173: } else if (Xs - xs < 0 && ye - Ye < 0 && Zs - zs < 0) {
1174: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1175: }
1176: if (n7 >= 0) { /* directly above */
1177: x_t = x;
1178: y_t = ly[(n7 % (m * n)) / m];
1179: z_t = lz[n7 / (m * n)];
1180: s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1181: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1182: } else if (ye - Ye < 0 && Zs - zs < 0) {
1183: for (j = 0; j < x; j++) idx[nn++] = -1;
1184: }
1185: if (n8 >= 0) { /* right above */
1186: x_t = lx[n8 % m];
1187: y_t = ly[(n8 % (m * n)) / m];
1188: z_t = lz[n8 / (m * n)];
1189: s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t;
1190: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1191: } else if (xe - Xe < 0 && ye - Ye < 0 && Zs - zs < 0) {
1192: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1193: }
1194: }
1195: }
1197: /* Middle Level */
1198: for (k = 0; k < z; k++) {
1199: for (i = 1; i <= s_y; i++) {
1200: if (n9 >= 0) { /* left below */
1201: x_t = lx[n9 % m];
1202: y_t = ly[(n9 % (m * n)) / m];
1203: /* z_t = z; */
1204: s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1205: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1206: } else if (Xs - xs < 0 && Ys - ys < 0) {
1207: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1208: }
1209: if (n10 >= 0) { /* directly below */
1210: x_t = x;
1211: y_t = ly[(n10 % (m * n)) / m];
1212: /* z_t = z; */
1213: s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1214: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1215: } else if (Ys - ys < 0) {
1216: if (by == DM_BOUNDARY_MIRROR) {
1217: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x;
1218: } else {
1219: for (j = 0; j < x; j++) idx[nn++] = -1;
1220: }
1221: }
1222: if (n11 >= 0) { /* right below */
1223: x_t = lx[n11 % m];
1224: y_t = ly[(n11 % (m * n)) / m];
1225: /* z_t = z; */
1226: s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1227: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1228: } else if (xe - Xe < 0 && Ys - ys < 0) {
1229: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1230: }
1231: }
1233: for (i = 0; i < y; i++) {
1234: if (n12 >= 0) { /* directly left */
1235: x_t = lx[n12 % m];
1236: y_t = y;
1237: /* z_t = z; */
1238: s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t;
1239: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1240: } else if (Xs - xs < 0) {
1241: if (bx == DM_BOUNDARY_MIRROR) {
1242: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x;
1243: } else {
1244: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1245: }
1246: }
1248: /* Interior */
1249: s_t = bases[rank] + i * x + k * x * y;
1250: for (j = 0; j < x; j++) idx[nn++] = s_t++;
1252: if (n14 >= 0) { /* directly right */
1253: x_t = lx[n14 % m];
1254: y_t = y;
1255: /* z_t = z; */
1256: s_t = bases[n14] + i * x_t + k * x_t * y_t;
1257: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1258: } else if (xe - Xe < 0) {
1259: if (bx == DM_BOUNDARY_MIRROR) {
1260: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x;
1261: } else {
1262: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1263: }
1264: }
1265: }
1267: for (i = 1; i <= s_y; i++) {
1268: if (n15 >= 0) { /* left above */
1269: x_t = lx[n15 % m];
1270: y_t = ly[(n15 % (m * n)) / m];
1271: /* z_t = z; */
1272: s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t;
1273: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1274: } else if (Xs - xs < 0 && ye - Ye < 0) {
1275: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1276: }
1277: if (n16 >= 0) { /* directly above */
1278: x_t = x;
1279: y_t = ly[(n16 % (m * n)) / m];
1280: /* z_t = z; */
1281: s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t;
1282: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1283: } else if (ye - Ye < 0) {
1284: if (by == DM_BOUNDARY_MIRROR) {
1285: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x;
1286: } else {
1287: for (j = 0; j < x; j++) idx[nn++] = -1;
1288: }
1289: }
1290: if (n17 >= 0) { /* right above */
1291: x_t = lx[n17 % m];
1292: y_t = ly[(n17 % (m * n)) / m];
1293: /* z_t = z; */
1294: s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t;
1295: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1296: } else if (xe - Xe < 0 && ye - Ye < 0) {
1297: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1298: }
1299: }
1300: }
1302: /* Upper Level */
1303: for (k = 0; k < s_z; k++) {
1304: for (i = 1; i <= s_y; i++) {
1305: if (n18 >= 0) { /* left below */
1306: x_t = lx[n18 % m];
1307: y_t = ly[(n18 % (m * n)) / m];
1308: /* z_t = lz[n18 / (m*n)]; */
1309: s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t;
1310: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1311: } else if (Xs - xs < 0 && Ys - ys < 0 && ze - Ze < 0) {
1312: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1313: }
1314: if (n19 >= 0) { /* directly below */
1315: x_t = x;
1316: y_t = ly[(n19 % (m * n)) / m];
1317: /* z_t = lz[n19 / (m*n)]; */
1318: s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1319: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1320: } else if (Ys - ys < 0 && ze - Ze < 0) {
1321: for (j = 0; j < x; j++) idx[nn++] = -1;
1322: }
1323: if (n20 >= 0) { /* right below */
1324: x_t = lx[n20 % m];
1325: y_t = ly[(n20 % (m * n)) / m];
1326: /* z_t = lz[n20 / (m*n)]; */
1327: s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t;
1328: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1329: } else if (xe - Xe < 0 && Ys - ys < 0 && ze - Ze < 0) {
1330: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1331: }
1332: }
1334: for (i = 0; i < y; i++) {
1335: if (n21 >= 0) { /* directly left */
1336: x_t = lx[n21 % m];
1337: y_t = y;
1338: /* z_t = lz[n21 / (m*n)]; */
1339: s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t;
1340: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1341: } else if (Xs - xs < 0 && ze - Ze < 0) {
1342: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1343: }
1345: if (n22 >= 0) { /* middle */
1346: x_t = x;
1347: y_t = y;
1348: /* z_t = lz[n22 / (m*n)]; */
1349: s_t = bases[n22] + i * x_t + k * x_t * y_t;
1350: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1351: } else if (ze - Ze < 0) {
1352: if (bz == DM_BOUNDARY_MIRROR) {
1353: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x;
1354: } else {
1355: for (j = 0; j < x; j++) idx[nn++] = -1;
1356: }
1357: }
1359: if (n23 >= 0) { /* directly right */
1360: x_t = lx[n23 % m];
1361: y_t = y;
1362: /* z_t = lz[n23 / (m*n)]; */
1363: s_t = bases[n23] + i * x_t + k * x_t * y_t;
1364: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1365: } else if (xe - Xe < 0 && ze - Ze < 0) {
1366: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1367: }
1368: }
1370: for (i = 1; i <= s_y; i++) {
1371: if (n24 >= 0) { /* left above */
1372: x_t = lx[n24 % m];
1373: y_t = ly[(n24 % (m * n)) / m];
1374: /* z_t = lz[n24 / (m*n)]; */
1375: s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t;
1376: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1377: } else if (Xs - xs < 0 && ye - Ye < 0 && ze - Ze < 0) {
1378: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1379: }
1380: if (n25 >= 0) { /* directly above */
1381: x_t = x;
1382: y_t = ly[(n25 % (m * n)) / m];
1383: /* z_t = lz[n25 / (m*n)]; */
1384: s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t;
1385: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
1386: } else if (ye - Ye < 0 && ze - Ze < 0) {
1387: for (j = 0; j < x; j++) idx[nn++] = -1;
1388: }
1389: if (n26 >= 0) { /* right above */
1390: x_t = lx[n26 % m];
1391: y_t = ly[(n26 % (m * n)) / m];
1392: /* z_t = lz[n26 / (m*n)]; */
1393: s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t;
1394: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
1395: } else if (xe - Xe < 0 && ye - Ye < 0 && ze - Ze < 0) {
1396: for (j = 0; j < s_x; j++) idx[nn++] = -1;
1397: }
1398: }
1399: }
1400: }
1401: /*
1402: Set the local to global ordering in the global vector, this allows use
1403: of VecSetValuesLocal().
1404: */
1405: PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap));
1407: PetscCall(PetscFree2(bases, ldims));
1408: dd->m = m;
1409: dd->n = n;
1410: dd->p = p;
1411: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1412: dd->xs = xs * dof;
1413: dd->xe = xe * dof;
1414: dd->ys = ys;
1415: dd->ye = ye;
1416: dd->zs = zs;
1417: dd->ze = ze;
1418: dd->Xs = Xs * dof;
1419: dd->Xe = Xe * dof;
1420: dd->Ys = Ys;
1421: dd->Ye = Ye;
1422: dd->Zs = Zs;
1423: dd->Ze = Ze;
1425: PetscCall(VecDestroy(&local));
1426: PetscCall(VecDestroy(&global));
1428: dd->gtol = gtol;
1429: dd->base = base;
1430: da->ops->view = DMView_DA_3d;
1431: dd->ltol = NULL;
1432: dd->ao = NULL;
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@C
1437: DMDACreate3d - Creates an object that will manage the communication of three-dimensional
1438: regular array data that is distributed across some processors.
1440: Collective
1442: Input Parameters:
1443: + comm - MPI communicator
1444: . bx - type of x ghost nodes the array have.
1445: Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1446: . by - type of y ghost nodes the array have.
1447: Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1448: . bz - type of z ghost nodes the array have.
1449: Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
1450: . stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`)
1451: . M - global dimension in x direction of the array
1452: . N - global dimension in y direction of the array
1453: . P - global dimension in z direction of the array
1454: . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated)
1455: . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated)
1456: . p - corresponding number of processors in z dimension (or `PETSC_DECIDE` to have calculated)
1457: . dof - number of degrees of freedom per node
1458: . s - stencil width
1459: . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`.
1460: . ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`.
1461: - lz - arrays containing the number of nodes in each cell along the z coordinates, or `NULL`.
1463: Output Parameter:
1464: . da - the resulting distributed array object
1466: Options Database Keys:
1467: + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()`
1468: . -da_grid_x <nx> - number of grid points in x direction
1469: . -da_grid_y <ny> - number of grid points in y direction
1470: . -da_grid_z <nz> - number of grid points in z direction
1471: . -da_processors_x <MX> - number of processors in x direction
1472: . -da_processors_y <MY> - number of processors in y direction
1473: . -da_processors_z <MZ> - number of processors in z direction
1474: . -da_refine_x <rx> - refinement ratio in x direction
1475: . -da_refine_y <ry> - refinement ratio in y direction
1476: . -da_refine_z <rz> - refinement ratio in z directio
1477: - -da_refine <n> - refine the DMDA n times before creating it
1479: Level: beginner
1481: Notes:
1482: If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the
1483: corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`,
1484: sum of the `ly` must `N`, sum of the `lz` must be `P`.
1486: The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
1487: standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
1488: the standard 27-pt stencil.
1490: The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
1491: The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
1492: and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed.
1494: You must call `DMSetUp()` after this call before using this `DM`.
1496: If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
1497: but before `DMSetUp()`.
1499: .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
1500: `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
1501: `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
1502: `DMStagCreate3d()`
1503: @*/
1504: 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)
1505: {
1506: PetscFunctionBegin;
1507: PetscCall(DMDACreate(comm, da));
1508: PetscCall(DMSetDimension(*da, 3));
1509: PetscCall(DMDASetSizes(*da, M, N, P));
1510: PetscCall(DMDASetNumProcs(*da, m, n, p));
1511: PetscCall(DMDASetBoundaryType(*da, bx, by, bz));
1512: PetscCall(DMDASetDof(*da, dof));
1513: PetscCall(DMDASetStencilType(*da, stencil_type));
1514: PetscCall(DMDASetStencilWidth(*da, s));
1515: PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz));
1516: PetscFunctionReturn(PETSC_SUCCESS);
1517: }