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