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