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