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