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