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