Actual source code: da2.c
2: #include <petsc/private/dmdaimpl.h>
3: #include <petscdraw.h>
5: static PetscErrorCode DMView_DA_2d(DM da, PetscViewer viewer)
6: {
7: PetscMPIInt rank;
8: PetscBool iascii, isdraw, isglvis, isbinary;
9: DM_DA *dd = (DM_DA *)da->data;
10: #if defined(PETSC_HAVE_MATLAB)
11: PetscBool ismatlab;
12: #endif
14: MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank);
16: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
17: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
18: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis);
19: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
20: #if defined(PETSC_HAVE_MATLAB)
21: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab);
22: #endif
23: if (iascii) {
24: PetscViewerFormat format;
26: PetscViewerGetFormat(viewer, &format);
27: if (format == PETSC_VIEWER_LOAD_BALANCE) {
28: PetscInt i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal;
29: DMDALocalInfo info;
30: PetscMPIInt size;
31: MPI_Comm_size(PetscObjectComm((PetscObject)da), &size);
32: DMDAGetLocalInfo(da, &info);
33: nzlocal = info.xm * info.ym;
34: PetscMalloc1(size, &nz);
35: MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da));
36: for (i = 0; i < (PetscInt)size; i++) {
37: nmax = PetscMax(nmax, nz[i]);
38: nmin = PetscMin(nmin, nz[i]);
39: navg += nz[i];
40: }
41: PetscFree(nz);
42: navg = navg / size;
43: PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax);
44: return 0;
45: }
46: if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) {
47: DMDALocalInfo info;
48: DMDAGetLocalInfo(da, &info);
49: PetscViewerASCIIPushSynchronized(viewer);
50: PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->m, dd->n, dd->w, dd->s);
51: PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs, info.xs + info.xm, info.ys, info.ys + info.ym);
52: PetscViewerFlush(viewer);
53: PetscViewerASCIIPopSynchronized(viewer);
54: } else if (format == PETSC_VIEWER_ASCII_GLVIS) DMView_DA_GLVis(da, viewer);
55: else DMView_DA_VTK(da, viewer);
56: } else if (isdraw) {
57: PetscDraw draw;
58: double ymin = -1 * dd->s - 1, ymax = dd->N + dd->s;
59: double xmin = -1 * dd->s - 1, xmax = dd->M + dd->s;
60: double x, y;
61: PetscInt base;
62: const PetscInt *idx;
63: char node[10];
64: PetscBool isnull;
66: PetscViewerDrawGetDraw(viewer, 0, &draw);
67: PetscDrawIsNull(draw, &isnull);
68: if (isnull) return 0;
70: PetscDrawCheckResizedWindow(draw);
71: PetscDrawClear(draw);
72: PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax);
74: PetscDrawCollectiveBegin(draw);
75: /* first processor draw all node lines */
76: if (rank == 0) {
77: ymin = 0.0;
78: ymax = dd->N - 1;
79: for (xmin = 0; xmin < dd->M; xmin++) PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK);
80: xmin = 0.0;
81: xmax = dd->M - 1;
82: for (ymin = 0; ymin < dd->N; ymin++) PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK);
83: }
84: PetscDrawCollectiveEnd(draw);
85: PetscDrawFlush(draw);
86: PetscDrawPause(draw);
88: PetscDrawCollectiveBegin(draw);
89: /* draw my box */
90: xmin = dd->xs / dd->w;
91: xmax = (dd->xe - 1) / dd->w;
92: ymin = dd->ys;
93: ymax = dd->ye - 1;
94: PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED);
95: PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED);
96: PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED);
97: PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED);
98: /* put in numbers */
99: base = (dd->base) / dd->w;
100: for (y = ymin; y <= ymax; y++) {
101: for (x = xmin; x <= xmax; x++) {
102: PetscSNPrintf(node, sizeof(node), "%d", (int)base++);
103: PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node);
104: }
105: }
106: PetscDrawCollectiveEnd(draw);
107: PetscDrawFlush(draw);
108: PetscDrawPause(draw);
110: PetscDrawCollectiveBegin(draw);
111: /* overlay ghost numbers, useful for error checking */
112: ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx);
113: base = 0;
114: xmin = dd->Xs;
115: xmax = dd->Xe;
116: ymin = dd->Ys;
117: ymax = dd->Ye;
118: for (y = ymin; y < ymax; y++) {
119: for (x = xmin; x < xmax; x++) {
120: if ((base % dd->w) == 0) {
121: PetscSNPrintf(node, sizeof(node), "%d", (int)(idx[base / dd->w]));
122: PetscDrawString(draw, x / dd->w, y, PETSC_DRAW_BLUE, node);
123: }
124: base++;
125: }
126: }
127: ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx);
128: PetscDrawCollectiveEnd(draw);
129: PetscDrawFlush(draw);
130: PetscDrawPause(draw);
131: PetscDrawSave(draw);
132: } else if (isglvis) {
133: DMView_DA_GLVis(da, viewer);
134: } else if (isbinary) {
135: DMView_DA_Binary(da, viewer);
136: #if defined(PETSC_HAVE_MATLAB)
137: } else if (ismatlab) {
138: DMView_DA_Matlab(da, viewer);
139: #endif
140: }
141: return 0;
142: }
144: #if defined(new)
145: /*
146: DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
147: function lives on a DMDA
149: y ~= (F(u + ha) - F(u))/h,
150: where F = nonlinear function, as set by SNESSetFunction()
151: u = current iterate
152: h = difference interval
153: */
154: PetscErrorCode DMDAGetDiagonal_MFFD(DM da, Vec U, Vec a)
155: {
156: PetscScalar h, *aa, *ww, v;
157: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
158: PetscInt gI, nI;
159: MatStencil stencil;
160: DMDALocalInfo info;
162: (*ctx->func)(0, U, a, ctx->funcctx);
163: (*ctx->funcisetbase)(U, ctx->funcctx);
165: VecGetArray(U, &ww);
166: VecGetArray(a, &aa);
168: nI = 0;
169: h = ww[gI];
170: if (h == 0.0) h = 1.0;
171: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
172: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
173: h *= epsilon;
175: ww[gI] += h;
176: (*ctx->funci)(i, w, &v, ctx->funcctx);
177: aa[nI] = (v - aa[nI]) / h;
178: ww[gI] -= h;
179: nI++;
181: VecRestoreArray(U, &ww);
182: VecRestoreArray(a, &aa);
183: return 0;
184: }
185: #endif
187: PetscErrorCode DMSetUp_DA_2D(DM da)
188: {
189: DM_DA *dd = (DM_DA *)da->data;
190: const PetscInt M = dd->M;
191: const PetscInt N = dd->N;
192: PetscInt m = dd->m;
193: PetscInt n = dd->n;
194: const PetscInt dof = dd->w;
195: const PetscInt s = dd->s;
196: DMBoundaryType bx = dd->bx;
197: DMBoundaryType by = dd->by;
198: DMDAStencilType stencil_type = dd->stencil_type;
199: PetscInt *lx = dd->lx;
200: PetscInt *ly = dd->ly;
201: MPI_Comm comm;
202: PetscMPIInt rank, size;
203: PetscInt xs, xe, ys, ye, x, y, Xs, Xe, Ys, Ye, IXs, IXe, IYs, IYe;
204: PetscInt up, down, left, right, i, n0, n1, n2, n3, n5, n6, n7, n8, *idx, nn;
205: PetscInt xbase, *bases, *ldims, j, x_t, y_t, s_t, base, count;
206: PetscInt s_x, s_y; /* s proportionalized to w */
207: PetscInt sn0 = 0, sn2 = 0, sn6 = 0, sn8 = 0;
208: Vec local, global;
209: VecScatter gtol;
210: IS to, from;
213: PetscObjectGetComm((PetscObject)da, &comm);
214: #if !defined(PETSC_USE_64BIT_INDICES)
216: #endif
218: MPI_Comm_size(comm, &size);
219: MPI_Comm_rank(comm, &rank);
221: dd->p = 1;
222: if (m != PETSC_DECIDE) {
225: }
226: if (n != PETSC_DECIDE) {
229: }
231: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
232: if (n != PETSC_DECIDE) {
233: m = size / n;
234: } else if (m != PETSC_DECIDE) {
235: n = size / m;
236: } else {
237: /* try for squarish distribution */
238: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
239: if (!m) m = 1;
240: while (m > 0) {
241: n = size / m;
242: if (m * n == size) break;
243: m--;
244: }
245: if (M > N && m < n) {
246: PetscInt _m = m;
247: m = n;
248: n = _m;
249: }
250: }
257: /*
258: Determine locally owned region
259: xs is the first local node number, x is the number of local nodes
260: */
261: if (!lx) {
262: PetscMalloc1(m, &dd->lx);
263: lx = dd->lx;
264: for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i);
265: }
266: x = lx[rank % m];
267: xs = 0;
268: for (i = 0; i < (rank % m); i++) xs += lx[i];
269: if (PetscDefined(USE_DEBUG)) {
270: left = xs;
271: for (i = (rank % m); i < m; i++) left += lx[i];
273: }
275: /*
276: Determine locally owned region
277: ys is the first local node number, y is the number of local nodes
278: */
279: if (!ly) {
280: PetscMalloc1(n, &dd->ly);
281: ly = dd->ly;
282: for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i);
283: }
284: y = ly[rank / m];
285: ys = 0;
286: for (i = 0; i < (rank / m); i++) ys += ly[i];
287: if (PetscDefined(USE_DEBUG)) {
288: left = ys;
289: for (i = (rank / m); i < n; i++) left += ly[i];
291: }
293: /*
294: check if the scatter requires more than one process neighbor or wraps around
295: the domain more than once
296: */
299: xe = xs + x;
300: ye = ys + y;
302: /* determine ghost region (Xs) and region scattered into (IXs) */
303: if (xs - s > 0) {
304: Xs = xs - s;
305: IXs = xs - s;
306: } else {
307: if (bx) {
308: Xs = xs - s;
309: } else {
310: Xs = 0;
311: }
312: IXs = 0;
313: }
314: if (xe + s <= M) {
315: Xe = xe + s;
316: IXe = xe + s;
317: } else {
318: if (bx) {
319: Xs = xs - s;
320: Xe = xe + s;
321: } else {
322: Xe = M;
323: }
324: IXe = M;
325: }
327: if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
328: IXs = xs - s;
329: IXe = xe + s;
330: Xs = xs - s;
331: Xe = xe + s;
332: }
334: if (ys - s > 0) {
335: Ys = ys - s;
336: IYs = ys - s;
337: } else {
338: if (by) {
339: Ys = ys - s;
340: } else {
341: Ys = 0;
342: }
343: IYs = 0;
344: }
345: if (ye + s <= N) {
346: Ye = ye + s;
347: IYe = ye + s;
348: } else {
349: if (by) {
350: Ye = ye + s;
351: } else {
352: Ye = N;
353: }
354: IYe = N;
355: }
357: if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
358: IYs = ys - s;
359: IYe = ye + s;
360: Ys = ys - s;
361: Ye = ye + s;
362: }
364: /* stencil length in each direction */
365: s_x = s;
366: s_y = s;
368: /* determine starting point of each processor */
369: nn = x * y;
370: PetscMalloc2(size + 1, &bases, size, &ldims);
371: MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm);
372: bases[0] = 0;
373: for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
374: for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
375: base = bases[rank] * dof;
377: /* allocate the base parallel and sequential vectors */
378: dd->Nlocal = x * y * dof;
379: VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global);
380: dd->nlocal = (Xe - Xs) * (Ye - Ys) * dof;
381: VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local);
383: /* generate global to local vector scatter and local to global mapping*/
385: /* global to local must include ghost points within the domain,
386: but not ghost points outside the domain that aren't periodic */
387: PetscMalloc1((IXe - IXs) * (IYe - IYs), &idx);
388: if (stencil_type == DMDA_STENCIL_BOX) {
389: left = IXs - Xs;
390: right = left + (IXe - IXs);
391: down = IYs - Ys;
392: up = down + (IYe - IYs);
393: count = 0;
394: for (i = down; i < up; i++) {
395: for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
396: }
397: ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to);
399: } else {
400: /* must drop into cross shape region */
401: /* ---------|
402: | top |
403: |--- ---| up
404: | middle |
405: | |
406: ---- ---- down
407: | bottom |
408: -----------
409: Xs xs xe Xe */
410: left = xs - Xs;
411: right = left + x;
412: down = ys - Ys;
413: up = down + y;
414: count = 0;
415: /* bottom */
416: for (i = (IYs - Ys); i < down; i++) {
417: for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
418: }
419: /* middle */
420: for (i = down; i < up; i++) {
421: for (j = (IXs - Xs); j < (IXe - Xs); j++) idx[count++] = j + i * (Xe - Xs);
422: }
423: /* top */
424: for (i = up; i < up + IYe - ye; i++) {
425: for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
426: }
427: ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to);
428: }
430: /* determine who lies on each side of us stored in n6 n7 n8
431: n3 n5
432: n0 n1 n2
433: */
435: /* Assume the Non-Periodic Case */
436: n1 = rank - m;
437: if (rank % m) {
438: n0 = n1 - 1;
439: } else {
440: n0 = -1;
441: }
442: if ((rank + 1) % m) {
443: n2 = n1 + 1;
444: n5 = rank + 1;
445: n8 = rank + m + 1;
446: if (n8 >= m * n) n8 = -1;
447: } else {
448: n2 = -1;
449: n5 = -1;
450: n8 = -1;
451: }
452: if (rank % m) {
453: n3 = rank - 1;
454: n6 = n3 + m;
455: if (n6 >= m * n) n6 = -1;
456: } else {
457: n3 = -1;
458: n6 = -1;
459: }
460: n7 = rank + m;
461: if (n7 >= m * n) n7 = -1;
463: if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) {
464: /* Modify for Periodic Cases */
465: /* Handle all four corners */
466: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m - 1;
467: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
468: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size - m;
469: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size - 1;
471: /* Handle Top and Bottom Sides */
472: if (n1 < 0) n1 = rank + m * (n - 1);
473: if (n7 < 0) n7 = rank - m * (n - 1);
474: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
475: if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
476: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
477: if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;
479: /* Handle Left and Right Sides */
480: if (n3 < 0) n3 = rank + (m - 1);
481: if (n5 < 0) n5 = rank - (m - 1);
482: if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
483: if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
484: if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
485: if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
486: } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
487: if (n1 < 0) n1 = rank + m * (n - 1);
488: if (n7 < 0) n7 = rank - m * (n - 1);
489: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
490: if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
491: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
492: if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;
493: } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
494: if (n3 < 0) n3 = rank + (m - 1);
495: if (n5 < 0) n5 = rank - (m - 1);
496: if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
497: if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
498: if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
499: if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
500: }
502: PetscMalloc1(9, &dd->neighbors);
504: dd->neighbors[0] = n0;
505: dd->neighbors[1] = n1;
506: dd->neighbors[2] = n2;
507: dd->neighbors[3] = n3;
508: dd->neighbors[4] = rank;
509: dd->neighbors[5] = n5;
510: dd->neighbors[6] = n6;
511: dd->neighbors[7] = n7;
512: dd->neighbors[8] = n8;
514: if (stencil_type == DMDA_STENCIL_STAR) {
515: /* save corner processor numbers */
516: sn0 = n0;
517: sn2 = n2;
518: sn6 = n6;
519: sn8 = n8;
520: n0 = n2 = n6 = n8 = -1;
521: }
523: PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx);
525: nn = 0;
526: xbase = bases[rank];
527: for (i = 1; i <= s_y; i++) {
528: if (n0 >= 0) { /* left below */
529: x_t = lx[n0 % m];
530: y_t = ly[(n0 / m)];
531: s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
532: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
533: }
535: if (n1 >= 0) { /* directly below */
536: x_t = x;
537: y_t = ly[(n1 / m)];
538: s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
539: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
540: } else if (by == DM_BOUNDARY_MIRROR) {
541: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
542: }
544: if (n2 >= 0) { /* right below */
545: x_t = lx[n2 % m];
546: y_t = ly[(n2 / m)];
547: s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
548: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
549: }
550: }
552: for (i = 0; i < y; i++) {
553: if (n3 >= 0) { /* directly left */
554: x_t = lx[n3 % m];
555: /* y_t = y; */
556: s_t = bases[n3] + (i + 1) * x_t - s_x;
557: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
558: } else if (bx == DM_BOUNDARY_MIRROR) {
559: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
560: }
562: for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */
564: if (n5 >= 0) { /* directly right */
565: x_t = lx[n5 % m];
566: /* y_t = y; */
567: s_t = bases[n5] + (i)*x_t;
568: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
569: } else if (bx == DM_BOUNDARY_MIRROR) {
570: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
571: }
572: }
574: for (i = 1; i <= s_y; i++) {
575: if (n6 >= 0) { /* left above */
576: x_t = lx[n6 % m];
577: /* y_t = ly[(n6/m)]; */
578: s_t = bases[n6] + (i)*x_t - s_x;
579: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
580: }
582: if (n7 >= 0) { /* directly above */
583: x_t = x;
584: /* y_t = ly[(n7/m)]; */
585: s_t = bases[n7] + (i - 1) * x_t;
586: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
587: } else if (by == DM_BOUNDARY_MIRROR) {
588: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
589: }
591: if (n8 >= 0) { /* right above */
592: x_t = lx[n8 % m];
593: /* y_t = ly[(n8/m)]; */
594: s_t = bases[n8] + (i - 1) * x_t;
595: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
596: }
597: }
599: ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from);
600: VecScatterCreate(global, from, local, to, >ol);
601: ISDestroy(&to);
602: ISDestroy(&from);
604: if (stencil_type == DMDA_STENCIL_STAR) {
605: n0 = sn0;
606: n2 = sn2;
607: n6 = sn6;
608: n8 = sn8;
609: }
611: if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) {
612: /*
613: Recompute the local to global mappings, this time keeping the
614: information about the cross corner processor numbers and any ghosted
615: but not periodic indices.
616: */
617: nn = 0;
618: xbase = bases[rank];
619: for (i = 1; i <= s_y; i++) {
620: if (n0 >= 0) { /* left below */
621: x_t = lx[n0 % m];
622: y_t = ly[(n0 / m)];
623: s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
624: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
625: } else if (xs - Xs > 0 && ys - Ys > 0) {
626: for (j = 0; j < s_x; j++) idx[nn++] = -1;
627: }
628: if (n1 >= 0) { /* directly below */
629: x_t = x;
630: y_t = ly[(n1 / m)];
631: s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
632: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
633: } else if (ys - Ys > 0) {
634: if (by == DM_BOUNDARY_MIRROR) {
635: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
636: } else {
637: for (j = 0; j < x; j++) idx[nn++] = -1;
638: }
639: }
640: if (n2 >= 0) { /* right below */
641: x_t = lx[n2 % m];
642: y_t = ly[(n2 / m)];
643: s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
644: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
645: } else if (Xe - xe > 0 && ys - Ys > 0) {
646: for (j = 0; j < s_x; j++) idx[nn++] = -1;
647: }
648: }
650: for (i = 0; i < y; i++) {
651: if (n3 >= 0) { /* directly left */
652: x_t = lx[n3 % m];
653: /* y_t = y; */
654: s_t = bases[n3] + (i + 1) * x_t - s_x;
655: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
656: } else if (xs - Xs > 0) {
657: if (bx == DM_BOUNDARY_MIRROR) {
658: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
659: } else {
660: for (j = 0; j < s_x; j++) idx[nn++] = -1;
661: }
662: }
664: for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */
666: if (n5 >= 0) { /* directly right */
667: x_t = lx[n5 % m];
668: /* y_t = y; */
669: s_t = bases[n5] + (i)*x_t;
670: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
671: } else if (Xe - xe > 0) {
672: if (bx == DM_BOUNDARY_MIRROR) {
673: for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
674: } else {
675: for (j = 0; j < s_x; j++) idx[nn++] = -1;
676: }
677: }
678: }
680: for (i = 1; i <= s_y; i++) {
681: if (n6 >= 0) { /* left above */
682: x_t = lx[n6 % m];
683: /* y_t = ly[(n6/m)]; */
684: s_t = bases[n6] + (i)*x_t - s_x;
685: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
686: } else if (xs - Xs > 0 && Ye - ye > 0) {
687: for (j = 0; j < s_x; j++) idx[nn++] = -1;
688: }
689: if (n7 >= 0) { /* directly above */
690: x_t = x;
691: /* y_t = ly[(n7/m)]; */
692: s_t = bases[n7] + (i - 1) * x_t;
693: for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
694: } else if (Ye - ye > 0) {
695: if (by == DM_BOUNDARY_MIRROR) {
696: for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
697: } else {
698: for (j = 0; j < x; j++) idx[nn++] = -1;
699: }
700: }
701: if (n8 >= 0) { /* right above */
702: x_t = lx[n8 % m];
703: /* y_t = ly[(n8/m)]; */
704: s_t = bases[n8] + (i - 1) * x_t;
705: for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
706: } else if (Xe - xe > 0 && Ye - ye > 0) {
707: for (j = 0; j < s_x; j++) idx[nn++] = -1;
708: }
709: }
710: }
711: /*
712: Set the local to global ordering in the global vector, this allows use
713: of VecSetValuesLocal().
714: */
715: ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap);
717: PetscFree2(bases, ldims);
718: dd->m = m;
719: dd->n = n;
720: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
721: dd->xs = xs * dof;
722: dd->xe = xe * dof;
723: dd->ys = ys;
724: dd->ye = ye;
725: dd->zs = 0;
726: dd->ze = 1;
727: dd->Xs = Xs * dof;
728: dd->Xe = Xe * dof;
729: dd->Ys = Ys;
730: dd->Ye = Ye;
731: dd->Zs = 0;
732: dd->Ze = 1;
734: VecDestroy(&local);
735: VecDestroy(&global);
737: dd->gtol = gtol;
738: dd->base = base;
739: da->ops->view = DMView_DA_2d;
740: dd->ltol = NULL;
741: dd->ao = NULL;
742: return 0;
743: }
745: /*@C
746: DMDACreate2d - Creates an object that will manage the communication of two-dimensional
747: regular array data that is distributed across some processors.
749: Collective
751: Input Parameters:
752: + comm - MPI communicator
753: . bx,by - type of ghost nodes the array have.
754: Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
755: . stencil_type - stencil type. Use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
756: . M,N - global dimension in each direction of the array
757: . m,n - corresponding number of processors in each dimension
758: (or `PETSC_DECIDE` to have calculated)
759: . dof - number of degrees of freedom per node
760: . s - stencil width
761: - lx, ly - arrays containing the number of nodes in each cell along
762: the x and y coordinates, or NULL. If non-null, these
763: must be of length as m and n, and the corresponding
764: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
765: must be M, and the sum of the ly[] entries must be N.
767: Output Parameter:
768: . da - the resulting distributed array object
770: Options Database Keys:
771: + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate2d()`
772: . -da_grid_x <nx> - number of grid points in x direction
773: . -da_grid_y <ny> - number of grid points in y direction
774: . -da_processors_x <nx> - number of processors in x direction
775: . -da_processors_y <ny> - number of processors in y direction
776: . -da_refine_x <rx> - refinement ratio in x direction
777: . -da_refine_y <ry> - refinement ratio in y direction
778: - -da_refine <n> - refine the DMDA n times before creating
780: Level: beginner
782: Notes:
783: The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
784: standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
785: the standard 9-pt stencil.
787: The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
788: The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
789: and DMCreateLocalVector() and calls to `VecDuplicate()` if more are needed.
791: You must call `DMSetUp()` after this call before using this `DM`.
793: If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
794: but before `DMSetUp()`.
796: .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
797: `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
798: `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
799: `DMStagCreate2d()`
800: @*/
802: PetscErrorCode DMDACreate2d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt m, PetscInt n, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], DM *da)
803: {
804: DMDACreate(comm, da);
805: DMSetDimension(*da, 2);
806: DMDASetSizes(*da, M, N, 1);
807: DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
808: DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);
809: DMDASetDof(*da, dof);
810: DMDASetStencilType(*da, stencil_type);
811: DMDASetStencilWidth(*da, s);
812: DMDASetOwnershipRanges(*da, lx, ly, NULL);
813: return 0;
814: }