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