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