Actual source code: da3.c
1: /*$Id: da3.c,v 1.136 2001/09/07 20:12:17 bsmith Exp $*/
3: /*
4: Code for manipulating distributed regular 3d arrays in parallel.
5: File created by Peter Mell 7/14/95
6: */
8: #include src/dm/da/daimpl.h
10: #if defined (PETSC_HAVE_AMS)
11: EXTERN_C_BEGIN
12: EXTERN int AMSSetFieldBlock_DA(AMS_Memory,char *,Vec);
13: EXTERN_C_END
14: #endif
16: #undef __FUNCT__
18: int DAView_3d(DA da,PetscViewer viewer)
19: {
20: int rank,ierr;
21: PetscTruth isascii,isdraw,isbinary;
24: MPI_Comm_rank(da->comm,&rank);
26: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
27: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
28: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
29: if (isascii) {
30: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %d N %d P %d m %d n %d p %d w %d s %dn",
31: rank,da->M,da->N,da->P,da->m,da->n,da->p,da->w,da->s);
32: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %d %d, Y range of indices: %d %d, Z range of indices: %d %dn",
33: da->xs,da->xe,da->ys,da->ye,da->zs,da->ze);
34: #if !defined(PETSC_USE_COMPLEX)
35: if (da->coordinates) {
36: int last;
37: PetscReal *coors;
38: VecGetArray(da->coordinates,&coors);
39: VecGetLocalSize(da->coordinates,&last);
40: last = last - 3;
41: PetscViewerASCIISynchronizedPrintf(viewer,"Lower left corner %g %g %g : Upper right %g %g %gn",
42: coors[0],coors[1],coors[2],coors[last],coors[last+1],coors[last+2]);
43: VecRestoreArray(da->coordinates,&coors);
44: }
45: #endif
46: PetscViewerFlush(viewer);
47: } else if (isdraw) {
48: PetscDraw draw;
49: PetscReal ymin = -1.0,ymax = (PetscReal)da->N;
50: PetscReal xmin = -1.0,xmax = (PetscReal)((da->M+2)*da->P),x,y,ycoord,xcoord;
51: int k,plane,base,*idx;
52: char node[10];
53: PetscTruth isnull;
55: PetscViewerDrawGetDraw(viewer,0,&draw);
56: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
57: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
58: PetscDrawSynchronizedClear(draw);
60: /* first processor draw all node lines */
61: if (!rank) {
62: for (k=0; k<da->P; k++) {
63: ymin = 0.0; ymax = (PetscReal)(da->N - 1);
64: for (xmin=(PetscReal)(k*(da->M+1)); xmin<(PetscReal)(da->M+(k*(da->M+1))); xmin++) {
65: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
66: }
67:
68: xmin = (PetscReal)(k*(da->M+1)); xmax = xmin + (PetscReal)(da->M - 1);
69: for (ymin=0; ymin<(PetscReal)da->N; ymin++) {
70: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
71: }
72: }
73: }
74: PetscDrawSynchronizedFlush(draw);
75: PetscDrawPause(draw);
77: for (k=0; k<da->P; k++) { /*Go through and draw for each plane*/
78: if ((k >= da->zs) && (k < da->ze)) {
79: /* draw my box */
80: ymin = da->ys;
81: ymax = da->ye - 1;
82: xmin = da->xs/da->w + (da->M+1)*k;
83: xmax =(da->xe-1)/da->w + (da->M+1)*k;
85: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
86: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
87: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
88: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
90: xmin = da->xs/da->w;
91: xmax =(da->xe-1)/da->w;
93: /* put in numbers*/
94: base = (da->base+(da->xe-da->xs)*(da->ye-da->ys)*(k-da->zs))/da->w;
96: /* Identify which processor owns the box */
97: sprintf(node,"%d",rank);
98: PetscDrawString(draw,xmin+(da->M+1)*k+.2,ymin+.3,PETSC_DRAW_RED,node);
100: for (y=ymin; y<=ymax; y++) {
101: for (x=xmin+(da->M+1)*k; x<=xmax+(da->M+1)*k; x++) {
102: sprintf(node,"%d",base++);
103: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
104: }
105: }
106:
107: }
108: }
109: PetscDrawSynchronizedFlush(draw);
110: PetscDrawPause(draw);
112: for (k=0-da->s; k<da->P+da->s; k++) {
113: /* Go through and draw for each plane */
114: if ((k >= da->Zs) && (k < da->Ze)) {
115:
116: /* overlay ghost numbers, useful for error checking */
117: base = (da->Xe-da->Xs)*(da->Ye-da->Ys)*(k-da->Zs); idx = da->idx;
118: plane=k;
119: /* Keep z wrap around points on the dradrawg */
120: if (k<0) { plane=da->P+k; }
121: if (k>=da->P) { plane=k-da->P; }
122: ymin = da->Ys; ymax = da->Ye;
123: xmin = (da->M+1)*plane*da->w;
124: xmax = (da->M+1)*plane*da->w+da->M*da->w;
125: for (y=ymin; y<ymax; y++) {
126: for (x=xmin+da->Xs; x<xmin+da->Xe; x+=da->w) {
127: sprintf(node,"%d",idx[base]/da->w);
128: ycoord = y;
129: /*Keep y wrap around points on drawing */
130: if (y<0) { ycoord = da->N+y; }
132: if (y>=da->N) { ycoord = y-da->N; }
133: xcoord = x; /* Keep x wrap points on drawing */
135: if (x<xmin) { xcoord = xmax - (xmin-x); }
136: if (x>=xmax) { xcoord = xmin + (x-xmax); }
137: PetscDrawString(draw,xcoord/da->w,ycoord,PETSC_DRAW_BLUE,node);
138: base+=da->w;
139: }
140: }
141: }
142: }
143: PetscDrawSynchronizedFlush(draw);
144: PetscDrawPause(draw);
145: } else if (isbinary) {
146: DAView_Binary(da,viewer);
147: } else {
148: SETERRQ1(1,"Viewer type %s not supported for DA 3d",((PetscObject)viewer)->type_name);
149: }
150: return(0);
151: }
153: EXTERN int DAPublish_Petsc(PetscObject);
155: #undef __FUNCT__
157: /*@C
158: DACreate3d - Creates an object that will manage the communication of three-dimensional
159: regular array data that is distributed across some processors.
161: Collective on MPI_Comm
163: Input Parameters:
164: + comm - MPI communicator
165: . wrap - type of periodicity the array should have, if any. Use one
166: of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, DA_XYPERIODIC, DA_XYZPERIODIC, DA_XZPERIODIC, or DA_YZPERIODIC.
167: . stencil_type - Type of stencil (DA_STENCIL_STAR or DA_STENCIL_BOX)
168: . M,N,P - global dimension in each direction of the array
169: . m,n,p - corresponding number of processors in each dimension
170: (or PETSC_DECIDE to have calculated)
171: . dof - number of degrees of freedom per node
172: . lx, ly, lz - arrays containing the number of nodes in each cell along
173: the x, y, and z coordinates, or PETSC_NULL. If non-null, these
174: must be of length as m,n,p and the corresponding
175: m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of
176: the ly[] must n, sum of the lz[] must be P
177: - s - stencil width
179: Output Parameter:
180: . inra - the resulting distributed array object
182: Options Database Key:
183: + -da_view - Calls DAView() at the conclusion of DACreate3d()
184: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
185: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
186: - -da_grid_z <nz> - number of grid points in z direction, if P < 0
188: Level: beginner
190: Notes:
191: The stencil type DA_STENCIL_STAR with width 1 corresponds to the
192: standard 7-pt stencil, while DA_STENCIL_BOX with width 1 denotes
193: the standard 27-pt stencil.
195: The array data itself is NOT stored in the DA, it is stored in Vec objects;
196: The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
197: and DACreateLocalVector() and calls to VecDuplicate() if more are needed.
199: .keywords: distributed array, create, three-dimensional
201: .seealso: DADestroy(), DAView(), DACreate1d(), DACreate2d(), DAGlobalToLocalBegin(),
202: DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(),
203: DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()
205: @*/
206: int DACreate3d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type,int M,
207: int N,int P,int m,int n,int p,int dof,int s,int *lx,int *ly,int *lz,DA *inra)
208: {
209: int rank,size,ierr,start,end,pm;
210: int xs,xe,ys,ye,zs,ze,x,y,z,Xs,Xe,Ys,Ye,Zs,Ze;
211: int left,up,down,bottom,top,i,j,k,*idx,nn,*flx = 0,*fly = 0,*flz = 0;
212: int n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n14;
213: int n15,n16,n17,n18,n19,n20,n21,n22,n23,n24,n25,n26;
214: int *bases,*ldims,x_t,y_t,z_t,s_t,base,count,s_x,s_y,s_z;
215: int tM = M,tN = N,tP = P;
216: int sn0 = 0,sn1 = 0,sn2 = 0,sn3 = 0,sn5 = 0,sn6 = 0,sn7 = 0;
217: int sn8 = 0,sn9 = 0,sn11 = 0,sn15 = 0,sn24 = 0,sn25 = 0,sn26 = 0;
218: int sn17 = 0,sn18 = 0,sn19 = 0,sn20 = 0,sn21 = 0,sn23 = 0,refine_x = 2, refine_y = 2, refine_z = 2;
219: PetscTruth flg1,flg2;
220: DA da;
221: Vec local,global;
222: VecScatter ltog,gtol;
223: IS to,from;
227: *inra = 0;
228: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
229: DMInitializePackage(PETSC_NULL);
230: #endif
232: if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %d",dof);
233: if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %d",s);
235: PetscOptionsBegin(comm,PETSC_NULL,"3d DA Options","DA");
236: if (M < 0){
237: tM = -M;
238: PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate3d",tM,&tM,PETSC_NULL);
239: }
240: if (N < 0){
241: tN = -N;
242: PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DACreate3d",tN,&tN,PETSC_NULL);
243: }
244: if (P < 0){
245: tP = -P;
246: PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DACreate3d",tP,&tP,PETSC_NULL);
247: }
248: PetscOptionsInt("-da_processors_x","Number of processors in x direction","DACreate3d",m,&m,PETSC_NULL);
249: PetscOptionsInt("-da_processors_y","Number of processors in y direction","DACreate3d",n,&n,PETSC_NULL);
250: PetscOptionsInt("-da_processors_z","Number of processors in z direction","DACreate3d",p,&p,PETSC_NULL);
251: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DACreate3d",refine_x,&refine_x,PETSC_NULL);
252: PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DACreate3d",refine_y,&refine_y,PETSC_NULL);
253: PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DACreate3d",refine_z,&refine_z,PETSC_NULL);
254: PetscOptionsEnd();
255: M = tM; N = tN; P = tP;
257: PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
258: da->bops->publish = DAPublish_Petsc;
259: da->ops->createglobalvector = DACreateGlobalVector;
260: da->ops->getinterpolation = DAGetInterpolation;
261: da->ops->getcoloring = DAGetColoring;
262: da->ops->getmatrix = DAGetMatrix;
263: da->ops->refine = DARefine;
265: PetscLogObjectCreate(da);
266: PetscLogObjectMemory(da,sizeof(struct _p_DA));
267: da->dim = 3;
268: da->interptype = DA_Q1;
269: da->refine_x = refine_x;
270: da->refine_y = refine_y;
271: da->refine_z = refine_z;
272: PetscMalloc(dof*sizeof(char*),&da->fieldname);
273: PetscMemzero(da->fieldname,dof*sizeof(char*));
275: MPI_Comm_size(comm,&size);
276: MPI_Comm_rank(comm,&rank);
278: if (m != PETSC_DECIDE) {
279: if (m < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %d",m);}
280: else if (m > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %d %d",m,size);}
281: }
282: if (n != PETSC_DECIDE) {
283: if (n < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %d",n);}
284: else if (n > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %d %d",n,size);}
285: }
286: if (p != PETSC_DECIDE) {
287: if (p < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %d",p);}
288: else if (p > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %d %d",p,size);}
289: }
291: /* Partition the array among the processors */
292: if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) {
293: m = size/(n*p);
294: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
295: n = size/(m*p);
296: } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
297: p = size/(m*n);
298: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) {
299: /* try for squarish distribution */
300: m = (int)(0.5 + sqrt(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p)));
301: if (!m) m = 1;
302: while (m > 0) {
303: n = size/(m*p);
304: if (m*n*p == size) break;
305: m--;
306: }
307: if (!m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %d",p);
308: if (M > N && m < n) {int _m = m; m = n; n = _m;}
309: } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) {
310: /* try for squarish distribution */
311: m = (int)(0.5 + sqrt(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
312: if (!m) m = 1;
313: while (m > 0) {
314: p = size/(m*n);
315: if (m*n*p == size) break;
316: m--;
317: }
318: if (!m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %d",n);
319: if (M > P && m < p) {int _m = m; m = p; p = _m;}
320: } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
321: /* try for squarish distribution */
322: n = (int)(0.5 + sqrt(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m)));
323: if (!n) n = 1;
324: while (n > 0) {
325: p = size/(m*n);
326: if (m*n*p == size) break;
327: n--;
328: }
329: if (!n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %d",n);
330: if (N > P && n < p) {int _n = n; n = p; p = _n;}
331: } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) {
332: /* try for squarish distribution */
333: n = (int)(0.5 + pow(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),1./3.));
334: if (!n) n = 1;
335: while (n > 0) {
336: pm = size/n;
337: if (n*pm == size) break;
338: n--;
339: }
340: if (!n) n = 1;
341: m = (int)(0.5 + sqrt(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
342: if (!m) m = 1;
343: while (m > 0) {
344: p = size/(m*n);
345: if (m*n*p == size) break;
346: m--;
347: }
348: if (M > P && m < p) {int _m = m; m = p; p = _m;}
349: } else if (m*n*p != size) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
351: if (m*n*p != size) SETERRQ(PETSC_ERR_PLIB,"Could not find good partition");
352: if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %d %d",M,m);
353: if (N < n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %d %d",N,n);
354: if (P < p) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in z direction is too fine! %d %d",P,p);
356: PetscOptionsHasName(PETSC_NULL,"-da_partition_nodes_at_end",&flg2);
357: /*
358: Determine locally owned region
359: [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes
360: */
361: if (lx) { /* user decided distribution */
362: x = lx[rank % m];
363: xs = 0;
364: for (i=0; i<(rank%m); i++) { xs += lx[i];}
365: if (x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column width is too thin for stencil! %d %d",x,s);
366: } else if (flg2) {
367: SETERRQ(PETSC_ERR_SUP,"-da_partition_nodes_at_end not supported");
368: } else { /* Normal PETSc distribution */
369: x = M/m + ((M % m) > (rank % m));
370: if (x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column width is too thin for stencil! %d %d",x,s);
371: if ((M % m) > (rank % m)) { xs = (rank % m)*x; }
372: else { xs = (M % m)*(x+1) + ((rank % m)-(M % m))*x; }
373: PetscMalloc(m*sizeof(int),&lx);
374: flx = lx;
375: for (i=0; i<m; i++) {
376: lx[i] = M/m + ((M % m) > (i % m));
377: }
378: }
379: if (ly) { /* user decided distribution */
380: y = ly[(rank % (m*n))/m];
381: if (y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row width is too thin for stencil! %d %d",y,s);
382: ys = 0;
383: for (i=0; i<(rank % (m*n))/m; i++) { ys += ly[i];}
384: } else if (flg2) {
385: SETERRQ(PETSC_ERR_SUP,"-da_partition_nodes_at_end not supported");
386: } else { /* Normal PETSc distribution */
387: y = N/n + ((N % n) > ((rank % (m*n)) /m));
388: if (y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row width is too thin for stencil! %d %d",y,s);
389: if ((N % n) > ((rank % (m*n)) /m)) {ys = ((rank % (m*n))/m)*y;}
390: else {ys = (N % n)*(y+1) + (((rank % (m*n))/m)-(N % n))*y;}
391: PetscMalloc(n*sizeof(int),&ly);
392: fly = ly;
393: for (i=0; i<n; i++) {
394: ly[i] = N/n + ((N % n) > (i % n));
395: }
396: }
397: if (lz) { /* user decided distribution */
398: z = lz[rank/(m*n)];
399: if (z < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Plane width is too thin for stencil! %d %d",z,s);
400: zs = 0;
401: for (i=0; i<(rank/(m*n)); i++) { zs += lz[i];}
402: } else if (flg2) {
403: SETERRQ(PETSC_ERR_SUP,"-da_partition_nodes_at_end not supported");
404: } else { /* Normal PETSc distribution */
405: z = P/p + ((P % p) > (rank / (m*n)));
406: if (z < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Plane width is too thin for stencil! %d %d",z,s);
407: if ((P % p) > (rank / (m*n))) {zs = (rank/(m*n))*z;}
408: else {zs = (P % p)*(z+1) + ((rank/(m*n))-(P % p))*z;}
409: PetscMalloc(p*sizeof(int),&lz);
410: flz = lz;
411: for (i=0; i<p; i++) {
412: lz[i] = P/p + ((P % p) > (i % p));
413: }
414: }
415: ye = ys + y;
416: xe = xs + x;
417: ze = zs + z;
419: /* determine ghost region */
420: /* Assume No Periodicity */
421: if (xs-s > 0) Xs = xs - s; else Xs = 0;
422: if (ys-s > 0) Ys = ys - s; else Ys = 0;
423: if (zs-s > 0) Zs = zs - s; else Zs = 0;
424: if (xe+s <= M) Xe = xe + s; else Xe = M;
425: if (ye+s <= N) Ye = ye + s; else Ye = N;
426: if (ze+s <= P) Ze = ze + s; else Ze = P;
428: /* X Periodic */
429: if (DAXPeriodic(wrap)){
430: Xs = xs - s;
431: Xe = xe + s;
432: }
434: /* Y Periodic */
435: if (DAYPeriodic(wrap)){
436: Ys = ys - s;
437: Ye = ye + s;
438: }
440: /* Z Periodic */
441: if (DAZPeriodic(wrap)){
442: Zs = zs - s;
443: Ze = ze + s;
444: }
446: /* Resize all X parameters to reflect w */
447: x *= dof;
448: xs *= dof;
449: xe *= dof;
450: Xs *= dof;
451: Xe *= dof;
452: s_x = s*dof;
453: s_y = s;
454: s_z = s;
456: /* determine starting point of each processor */
457: nn = x*y*z;
458: ierr = PetscMalloc((2*size+1)*sizeof(int),&bases);
459: ldims = (int*)(bases+size+1);
460: ierr = MPI_Allgather(&nn,1,MPI_INT,ldims,1,MPI_INT,comm);
461: bases[0] = 0;
462: for (i=1; i<=size; i++) {
463: bases[i] = ldims[i-1];
464: }
465: for (i=1; i<=size; i++) {
466: bases[i] += bases[i-1];
467: }
469: /* allocate the base parallel and sequential vectors */
470: da->Nlocal = x*y*z;
471: VecCreateMPIWithArray(comm,da->Nlocal,PETSC_DECIDE,0,&global);
472: VecSetBlockSize(global,dof);
473: da->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs);
474: VecCreateSeqWithArray(MPI_COMM_SELF,da->nlocal,0,&local);
475: VecSetBlockSize(local,dof);
477: /* generate appropriate vector scatters */
478: /* local to global inserts non-ghost point region into global */
479: VecGetOwnershipRange(global,&start,&end);
480: ISCreateStride(comm,x*y*z,start,1,&to);
482: left = xs - Xs;
483: bottom = ys - Ys; top = bottom + y;
484: down = zs - Zs; up = down + z;
485: count = x*(top-bottom)*(up-down);
486: PetscMalloc(count*sizeof(int),&idx);
487: count = 0;
488: for (i=down; i<up; i++) {
489: for (j=bottom; j<top; j++) {
490: for (k=0; k<x; k++) {
491: idx[count++] = (left+j*(Xe-Xs))+i*(Xe-Xs)*(Ye-Ys) + k;
492: }
493: }
494: }
495: ISCreateGeneral(comm,count,idx,&from);
496: PetscFree(idx);
498: VecScatterCreate(local,from,global,to,<og);
499: PetscLogObjectParent(da,to);
500: PetscLogObjectParent(da,from);
501: PetscLogObjectParent(da,ltog);
502: ISDestroy(from);
503: ISDestroy(to);
505: /* global to local must include ghost points */
506: if (stencil_type == DA_STENCIL_BOX) {
507: ISCreateStride(comm,(Xe-Xs)*(Ye-Ys)*(Ze-Zs),0,1,&to);
508: } else {
509: /* This is way ugly! We need to list the funny cross type region */
510: /* the bottom chunck */
511: left = xs - Xs;
512: bottom = ys - Ys; top = bottom + y;
513: down = zs - Zs; up = down + z;
514: count = down*(top-bottom)*x +
515: (up-down)*(bottom*x + (top-bottom)*(Xe-Xs) + (Ye-Ys-top)*x) +
516: (Ze-Zs-up)*(top-bottom)*x;
517: PetscMalloc(count*sizeof(int),&idx);
518: count = 0;
519: for (i=0; i<down; i++) {
520: for (j=bottom; j<top; j++) {
521: for (k=0; k<x; k++) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
522: }
523: }
524: /* the middle piece */
525: for (i=down; i<up; i++) {
526: /* front */
527: for (j=0; j<bottom; j++) {
528: for (k=0; k<x; k++) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
529: }
530: /* middle */
531: for (j=bottom; j<top; j++) {
532: for (k=0; k<Xe-Xs; k++) idx[count++] = j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
533: }
534: /* back */
535: for (j=top; j<Ye-Ys; j++) {
536: for (k=0; k<x; k++) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
537: }
538: }
539: /* the top piece */
540: for (i=up; i<Ze-Zs; i++) {
541: for (j=bottom; j<top; j++) {
542: for (k=0; k<x; k++) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
543: }
544: }
545: ISCreateGeneral(comm,count,idx,&to);
546: PetscFree(idx);
547: }
549: /* determine who lies on each side of use stored in n24 n25 n26
550: n21 n22 n23
551: n18 n19 n20
553: n15 n16 n17
554: n12 n14
555: n9 n10 n11
557: n6 n7 n8
558: n3 n4 n5
559: n0 n1 n2
560: */
561:
562: /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
563:
564: /* Assume Nodes are Internal to the Cube */
565:
566: n0 = rank - m*n - m - 1;
567: n1 = rank - m*n - m;
568: n2 = rank - m*n - m + 1;
569: n3 = rank - m*n -1;
570: n4 = rank - m*n;
571: n5 = rank - m*n + 1;
572: n6 = rank - m*n + m - 1;
573: n7 = rank - m*n + m;
574: n8 = rank - m*n + m + 1;
576: n9 = rank - m - 1;
577: n10 = rank - m;
578: n11 = rank - m + 1;
579: n12 = rank - 1;
580: n14 = rank + 1;
581: n15 = rank + m - 1;
582: n16 = rank + m;
583: n17 = rank + m + 1;
585: n18 = rank + m*n - m - 1;
586: n19 = rank + m*n - m;
587: n20 = rank + m*n - m + 1;
588: n21 = rank + m*n - 1;
589: n22 = rank + m*n;
590: n23 = rank + m*n + 1;
591: n24 = rank + m*n + m - 1;
592: n25 = rank + m*n + m;
593: n26 = rank + m*n + m + 1;
595: /* Assume Pieces are on Faces of Cube */
597: if (xs == 0) { /* First assume not corner or edge */
598: n0 = rank -1 - (m*n);
599: n3 = rank + m -1 - (m*n);
600: n6 = rank + 2*m -1 - (m*n);
601: n9 = rank -1;
602: n12 = rank + m -1;
603: n15 = rank + 2*m -1;
604: n18 = rank -1 + (m*n);
605: n21 = rank + m -1 + (m*n);
606: n24 = rank + 2*m -1 + (m*n);
607: }
609: if (xe == M*dof) { /* First assume not corner or edge */
610: n2 = rank -2*m +1 - (m*n);
611: n5 = rank - m +1 - (m*n);
612: n8 = rank +1 - (m*n);
613: n11 = rank -2*m +1;
614: n14 = rank - m +1;
615: n17 = rank +1;
616: n20 = rank -2*m +1 + (m*n);
617: n23 = rank - m +1 + (m*n);
618: n26 = rank +1 + (m*n);
619: }
621: if (ys==0) { /* First assume not corner or edge */
622: n0 = rank + m * (n-1) -1 - (m*n);
623: n1 = rank + m * (n-1) - (m*n);
624: n2 = rank + m * (n-1) +1 - (m*n);
625: n9 = rank + m * (n-1) -1;
626: n10 = rank + m * (n-1);
627: n11 = rank + m * (n-1) +1;
628: n18 = rank + m * (n-1) -1 + (m*n);
629: n19 = rank + m * (n-1) + (m*n);
630: n20 = rank + m * (n-1) +1 + (m*n);
631: }
633: if (ye == N) { /* First assume not corner or edge */
634: n6 = rank - m * (n-1) -1 - (m*n);
635: n7 = rank - m * (n-1) - (m*n);
636: n8 = rank - m * (n-1) +1 - (m*n);
637: n15 = rank - m * (n-1) -1;
638: n16 = rank - m * (n-1);
639: n17 = rank - m * (n-1) +1;
640: n24 = rank - m * (n-1) -1 + (m*n);
641: n25 = rank - m * (n-1) + (m*n);
642: n26 = rank - m * (n-1) +1 + (m*n);
643: }
644:
645: if (zs == 0) { /* First assume not corner or edge */
646: n0 = size - (m*n) + rank - m - 1;
647: n1 = size - (m*n) + rank - m;
648: n2 = size - (m*n) + rank - m + 1;
649: n3 = size - (m*n) + rank - 1;
650: n4 = size - (m*n) + rank;
651: n5 = size - (m*n) + rank + 1;
652: n6 = size - (m*n) + rank + m - 1;
653: n7 = size - (m*n) + rank + m ;
654: n8 = size - (m*n) + rank + m + 1;
655: }
657: if (ze == P) { /* First assume not corner or edge */
658: n18 = (m*n) - (size-rank) - m - 1;
659: n19 = (m*n) - (size-rank) - m;
660: n20 = (m*n) - (size-rank) - m + 1;
661: n21 = (m*n) - (size-rank) - 1;
662: n22 = (m*n) - (size-rank);
663: n23 = (m*n) - (size-rank) + 1;
664: n24 = (m*n) - (size-rank) + m - 1;
665: n25 = (m*n) - (size-rank) + m;
666: n26 = (m*n) - (size-rank) + m + 1;
667: }
669: if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
670: n0 = size - m*n + rank + m-1 - m;
671: n3 = size - m*n + rank + m-1;
672: n6 = size - m*n + rank + m-1 + m;
673: }
674:
675: if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
676: n18 = m*n - (size - rank) + m-1 - m;
677: n21 = m*n - (size - rank) + m-1;
678: n24 = m*n - (size - rank) + m-1 + m;
679: }
681: if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
682: n0 = rank + m*n -1 - m*n;
683: n9 = rank + m*n -1;
684: n18 = rank + m*n -1 + m*n;
685: }
687: if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
688: n6 = rank - m*(n-1) + m-1 - m*n;
689: n15 = rank - m*(n-1) + m-1;
690: n24 = rank - m*(n-1) + m-1 + m*n;
691: }
693: if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
694: n2 = size - (m*n-rank) - (m-1) - m;
695: n5 = size - (m*n-rank) - (m-1);
696: n8 = size - (m*n-rank) - (m-1) + m;
697: }
699: if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
700: n20 = m*n - (size - rank) - (m-1) - m;
701: n23 = m*n - (size - rank) - (m-1);
702: n26 = m*n - (size - rank) - (m-1) + m;
703: }
705: if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
706: n2 = rank + m*(n-1) - (m-1) - m*n;
707: n11 = rank + m*(n-1) - (m-1);
708: n20 = rank + m*(n-1) - (m-1) + m*n;
709: }
711: if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
712: n8 = rank - m*n +1 - m*n;
713: n17 = rank - m*n +1;
714: n26 = rank - m*n +1 + m*n;
715: }
717: if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
718: n0 = size - m + rank -1;
719: n1 = size - m + rank;
720: n2 = size - m + rank +1;
721: }
723: if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
724: n18 = m*n - (size - rank) + m*(n-1) -1;
725: n19 = m*n - (size - rank) + m*(n-1);
726: n20 = m*n - (size - rank) + m*(n-1) +1;
727: }
729: if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
730: n6 = size - (m*n-rank) - m * (n-1) -1;
731: n7 = size - (m*n-rank) - m * (n-1);
732: n8 = size - (m*n-rank) - m * (n-1) +1;
733: }
735: if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
736: n24 = rank - (size-m) -1;
737: n25 = rank - (size-m);
738: n26 = rank - (size-m) +1;
739: }
741: /* Check for Corners */
742: if ((xs==0) && (ys==0) && (zs==0)) { n0 = size -1;}
743: if ((xs==0) && (ys==0) && (ze==P)) { n18 = m*n-1;}
744: if ((xs==0) && (ye==N) && (zs==0)) { n6 = (size-1)-m*(n-1);}
745: if ((xs==0) && (ye==N) && (ze==P)) { n24 = m-1;}
746: if ((xe==M*dof) && (ys==0) && (zs==0)) { n2 = size-m;}
747: if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
748: if ((xe==M*dof) && (ye==N) && (zs==0)) { n8 = size-m*n;}
749: if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}
751: /* Check for when not X,Y, and Z Periodic */
753: /* If not X periodic */
754: if ((wrap != DA_XPERIODIC) && (wrap != DA_XYPERIODIC) &&
755: (wrap != DA_XZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
756: if (xs==0) {n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;}
757: if (xe==M*dof) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
758: }
760: /* If not Y periodic */
761: if ((wrap != DA_YPERIODIC) && (wrap != DA_XYPERIODIC) &&
762: (wrap != DA_YZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
763: if (ys==0) {n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;}
764: if (ye==N) {n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
765: }
767: /* If not Z periodic */
768: if ((wrap != DA_ZPERIODIC) && (wrap != DA_XZPERIODIC) &&
769: (wrap != DA_YZPERIODIC) && (wrap != DA_XYZPERIODIC)) {
770: if (zs==0) {n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;}
771: if (ze==P) {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
772: }
774: /* If star stencil then delete the corner neighbors */
775: if (stencil_type == DA_STENCIL_STAR) {
776: /* save information about corner neighbors */
777: sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7;
778: sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18;
779: sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25;
780: sn26 = n26;
781: n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 =
782: n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1;
783: }
786: PetscMalloc((Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(int),&idx);
787: PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(int));
789: nn = 0;
791: /* Bottom Level */
792: for (k=0; k<s_z; k++) {
793: for (i=1; i<=s_y; i++) {
794: if (n0 >= 0) { /* left below */
795: x_t = lx[n0 % m]*dof;
796: y_t = ly[(n0 % (m*n))/m];
797: z_t = lz[n0 / (m*n)];
798: 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;
799: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
800: }
801: if (n1 >= 0) { /* directly below */
802: x_t = x;
803: y_t = ly[(n1 % (m*n))/m];
804: z_t = lz[n1 / (m*n)];
805: s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
806: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
807: }
808: if (n2 >= 0) { /* right below */
809: x_t = lx[n2 % m]*dof;
810: y_t = ly[(n2 % (m*n))/m];
811: z_t = lz[n2 / (m*n)];
812: s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
813: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
814: }
815: }
817: for (i=0; i<y; i++) {
818: if (n3 >= 0) { /* directly left */
819: x_t = lx[n3 % m]*dof;
820: y_t = y;
821: z_t = lz[n3 / (m*n)];
822: s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
823: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
824: }
826: if (n4 >= 0) { /* middle */
827: x_t = x;
828: y_t = y;
829: z_t = lz[n4 / (m*n)];
830: s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
831: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
832: }
834: if (n5 >= 0) { /* directly right */
835: x_t = lx[n5 % m]*dof;
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: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
840: }
841: }
843: for (i=1; i<=s_y; i++) {
844: if (n6 >= 0) { /* left above */
845: x_t = lx[n6 % m]*dof;
846: y_t = ly[(n6 % (m*n))/m];
847: z_t = lz[n6 / (m*n)];
848: s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
849: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
850: }
851: if (n7 >= 0) { /* directly above */
852: x_t = x;
853: y_t = ly[(n7 % (m*n))/m];
854: z_t = lz[n7 / (m*n)];
855: s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
856: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
857: }
858: if (n8 >= 0) { /* right above */
859: x_t = lx[n8 % m]*dof;
860: y_t = ly[(n8 % (m*n))/m];
861: z_t = lz[n8 / (m*n)];
862: s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
863: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
864: }
865: }
866: }
868: /* Middle Level */
869: for (k=0; k<z; k++) {
870: for (i=1; i<=s_y; i++) {
871: if (n9 >= 0) { /* left below */
872: x_t = lx[n9 % m]*dof;
873: y_t = ly[(n9 % (m*n))/m];
874: /* z_t = z; */
875: s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
876: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
877: }
878: if (n10 >= 0) { /* directly below */
879: x_t = x;
880: y_t = ly[(n10 % (m*n))/m];
881: /* z_t = z; */
882: s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
883: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
884: }
885: if (n11 >= 0) { /* right below */
886: x_t = lx[n11 % m]*dof;
887: y_t = ly[(n11 % (m*n))/m];
888: /* z_t = z; */
889: s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
890: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
891: }
892: }
894: for (i=0; i<y; i++) {
895: if (n12 >= 0) { /* directly left */
896: x_t = lx[n12 % m]*dof;
897: y_t = y;
898: /* z_t = z; */
899: s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
900: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
901: }
903: /* Interior */
904: s_t = bases[rank] + i*x + k*x*y;
905: for (j=0; j<x; j++) { idx[nn++] = s_t++;}
907: if (n14 >= 0) { /* directly right */
908: x_t = lx[n14 % m]*dof;
909: y_t = y;
910: /* z_t = z; */
911: s_t = bases[n14] + i*x_t + k*x_t*y_t;
912: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
913: }
914: }
916: for (i=1; i<=s_y; i++) {
917: if (n15 >= 0) { /* left above */
918: x_t = lx[n15 % m]*dof;
919: y_t = ly[(n15 % (m*n))/m];
920: /* z_t = z; */
921: s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
922: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
923: }
924: if (n16 >= 0) { /* directly above */
925: x_t = x;
926: y_t = ly[(n16 % (m*n))/m];
927: /* z_t = z; */
928: s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
929: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
930: }
931: if (n17 >= 0) { /* right above */
932: x_t = lx[n17 % m]*dof;
933: y_t = ly[(n17 % (m*n))/m];
934: /* z_t = z; */
935: s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
936: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
937: }
938: }
939: }
940:
941: /* Upper Level */
942: for (k=0; k<s_z; k++) {
943: for (i=1; i<=s_y; i++) {
944: if (n18 >= 0) { /* left below */
945: x_t = lx[n18 % m]*dof;
946: y_t = ly[(n18 % (m*n))/m];
947: /* z_t = lz[n18 / (m*n)]; */
948: s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
949: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
950: }
951: if (n19 >= 0) { /* directly below */
952: x_t = x;
953: y_t = ly[(n19 % (m*n))/m];
954: /* z_t = lz[n19 / (m*n)]; */
955: s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
956: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
957: }
958: if (n20 >= 0) { /* right below */
959: x_t = lx[n20 % m]*dof;
960: y_t = ly[(n20 % (m*n))/m];
961: /* z_t = lz[n20 / (m*n)]; */
962: s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
963: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
964: }
965: }
967: for (i=0; i<y; i++) {
968: if (n21 >= 0) { /* directly left */
969: x_t = lx[n21 % m]*dof;
970: y_t = y;
971: /* z_t = lz[n21 / (m*n)]; */
972: s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
973: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
974: }
976: if (n22 >= 0) { /* middle */
977: x_t = x;
978: y_t = y;
979: /* z_t = lz[n22 / (m*n)]; */
980: s_t = bases[n22] + i*x_t + k*x_t*y_t;
981: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
982: }
984: if (n23 >= 0) { /* directly right */
985: x_t = lx[n23 % m]*dof;
986: y_t = y;
987: /* z_t = lz[n23 / (m*n)]; */
988: s_t = bases[n23] + i*x_t + k*x_t*y_t;
989: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
990: }
991: }
993: for (i=1; i<=s_y; i++) {
994: if (n24 >= 0) { /* left above */
995: x_t = lx[n24 % m]*dof;
996: y_t = ly[(n24 % (m*n))/m];
997: /* z_t = lz[n24 / (m*n)]; */
998: s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
999: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1000: }
1001: if (n25 >= 0) { /* directly above */
1002: x_t = x;
1003: y_t = ly[(n25 % (m*n))/m];
1004: /* z_t = lz[n25 / (m*n)]; */
1005: s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1006: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1007: }
1008: if (n26 >= 0) { /* right above */
1009: x_t = lx[n26 % m]*dof;
1010: y_t = ly[(n26 % (m*n))/m];
1011: /* z_t = lz[n26 / (m*n)]; */
1012: s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1013: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1014: }
1015: }
1016: }
1017: base = bases[rank];
1018: ISCreateGeneral(comm,nn,idx,&from);
1019: VecScatterCreate(global,from,local,to,>ol);
1020: PetscLogObjectParent(da,gtol);
1021: PetscLogObjectParent(da,to);
1022: PetscLogObjectParent(da,from);
1023: ISDestroy(to);
1024: ISDestroy(from);
1025: da->stencil_type = stencil_type;
1026: da->M = M; da->N = N; da->P = P;
1027: da->m = m; da->n = n; da->p = p;
1028: da->w = dof; da->s = s;
1029: da->xs = xs; da->xe = xe; da->ys = ys; da->ye = ye; da->zs = zs; da->ze = ze;
1030: da->Xs = Xs; da->Xe = Xe; da->Ys = Ys; da->Ye = Ye; da->Zs = Zs; da->Ze = Ze;
1032: VecDestroy(local);
1033: VecDestroy(global);
1035: if (stencil_type == DA_STENCIL_STAR) {
1036: /*
1037: Recompute the local to global mappings, this time keeping the
1038: information about the cross corner processor numbers.
1039: */
1040: n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7;
1041: n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18;
1042: n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25;
1043: n26 = sn26;
1045: nn = 0;
1047: /* Bottom Level */
1048: for (k=0; k<s_z; k++) {
1049: for (i=1; i<=s_y; i++) {
1050: if (n0 >= 0) { /* left below */
1051: x_t = lx[n0 % m]*dof;
1052: y_t = ly[(n0 % (m*n))/m];
1053: z_t = lz[n0 / (m*n)];
1054: 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;
1055: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1056: }
1057: if (n1 >= 0) { /* directly below */
1058: x_t = x;
1059: y_t = ly[(n1 % (m*n))/m];
1060: z_t = lz[n1 / (m*n)];
1061: s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1062: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1063: }
1064: if (n2 >= 0) { /* right below */
1065: x_t = lx[n2 % m]*dof;
1066: y_t = ly[(n2 % (m*n))/m];
1067: z_t = lz[n2 / (m*n)];
1068: s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1069: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1070: }
1071: }
1073: for (i=0; i<y; i++) {
1074: if (n3 >= 0) { /* directly left */
1075: x_t = lx[n3 % m]*dof;
1076: y_t = y;
1077: z_t = lz[n3 / (m*n)];
1078: s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1079: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1080: }
1082: if (n4 >= 0) { /* middle */
1083: x_t = x;
1084: y_t = y;
1085: z_t = lz[n4 / (m*n)];
1086: s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1087: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1088: }
1090: if (n5 >= 0) { /* directly right */
1091: x_t = lx[n5 % m]*dof;
1092: y_t = y;
1093: z_t = lz[n5 / (m*n)];
1094: s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1095: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1096: }
1097: }
1099: for (i=1; i<=s_y; i++) {
1100: if (n6 >= 0) { /* left above */
1101: x_t = lx[n6 % m]*dof;
1102: y_t = ly[(n6 % (m*n))/m];
1103: z_t = lz[n6 / (m*n)];
1104: s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1105: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1106: }
1107: if (n7 >= 0) { /* directly above */
1108: x_t = x;
1109: y_t = ly[(n7 % (m*n))/m];
1110: z_t = lz[n7 / (m*n)];
1111: s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1112: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1113: }
1114: if (n8 >= 0) { /* right above */
1115: x_t = lx[n8 % m]*dof;
1116: y_t = ly[(n8 % (m*n))/m];
1117: z_t = lz[n8 / (m*n)];
1118: s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1119: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1120: }
1121: }
1122: }
1124: /* Middle Level */
1125: for (k=0; k<z; k++) {
1126: for (i=1; i<=s_y; i++) {
1127: if (n9 >= 0) { /* left below */
1128: x_t = lx[n9 % m]*dof;
1129: y_t = ly[(n9 % (m*n))/m];
1130: /* z_t = z; */
1131: s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1132: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1133: }
1134: if (n10 >= 0) { /* directly below */
1135: x_t = x;
1136: y_t = ly[(n10 % (m*n))/m];
1137: /* z_t = z; */
1138: s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1139: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1140: }
1141: if (n11 >= 0) { /* right below */
1142: x_t = lx[n11 % m]*dof;
1143: y_t = ly[(n11 % (m*n))/m];
1144: /* z_t = z; */
1145: s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1146: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1147: }
1148: }
1150: for (i=0; i<y; i++) {
1151: if (n12 >= 0) { /* directly left */
1152: x_t = lx[n12 % m]*dof;
1153: y_t = y;
1154: /* z_t = z; */
1155: s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
1156: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1157: }
1159: /* Interior */
1160: s_t = bases[rank] + i*x + k*x*y;
1161: for (j=0; j<x; j++) { idx[nn++] = s_t++;}
1163: if (n14 >= 0) { /* directly right */
1164: x_t = lx[n14 % m]*dof;
1165: y_t = y;
1166: /* z_t = z; */
1167: s_t = bases[n14] + i*x_t + k*x_t*y_t;
1168: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1169: }
1170: }
1172: for (i=1; i<=s_y; i++) {
1173: if (n15 >= 0) { /* left above */
1174: x_t = lx[n15 % m]*dof;
1175: y_t = ly[(n15 % (m*n))/m];
1176: /* z_t = z; */
1177: s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
1178: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1179: }
1180: if (n16 >= 0) { /* directly above */
1181: x_t = x;
1182: y_t = ly[(n16 % (m*n))/m];
1183: /* z_t = z; */
1184: s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
1185: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1186: }
1187: if (n17 >= 0) { /* right above */
1188: x_t = lx[n17 % m]*dof;
1189: y_t = ly[(n17 % (m*n))/m];
1190: /* z_t = z; */
1191: s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
1192: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1193: }
1194: }
1195: }
1196:
1197: /* Upper Level */
1198: for (k=0; k<s_z; k++) {
1199: for (i=1; i<=s_y; i++) {
1200: if (n18 >= 0) { /* left below */
1201: x_t = lx[n18 % m]*dof;
1202: y_t = ly[(n18 % (m*n))/m];
1203: /* z_t = lz[n18 / (m*n)]; */
1204: s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1205: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1206: }
1207: if (n19 >= 0) { /* directly below */
1208: x_t = x;
1209: y_t = ly[(n19 % (m*n))/m];
1210: /* z_t = lz[n19 / (m*n)]; */
1211: s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1212: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1213: }
1214: if (n20 >= 0) { /* right below */
1215: x_t = lx[n20 % m]*dof;
1216: y_t = ly[(n20 % (m*n))/m];
1217: /* z_t = lz[n20 / (m*n)]; */
1218: s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1219: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1220: }
1221: }
1223: for (i=0; i<y; i++) {
1224: if (n21 >= 0) { /* directly left */
1225: x_t = lx[n21 % m]*dof;
1226: y_t = y;
1227: /* z_t = lz[n21 / (m*n)]; */
1228: s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
1229: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1230: }
1232: if (n22 >= 0) { /* middle */
1233: x_t = x;
1234: y_t = y;
1235: /* z_t = lz[n22 / (m*n)]; */
1236: s_t = bases[n22] + i*x_t + k*x_t*y_t;
1237: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1238: }
1240: if (n23 >= 0) { /* directly right */
1241: x_t = lx[n23 % m]*dof;
1242: y_t = y;
1243: /* z_t = lz[n23 / (m*n)]; */
1244: s_t = bases[n23] + i*x_t + k*x_t*y_t;
1245: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1246: }
1247: }
1249: for (i=1; i<=s_y; i++) {
1250: if (n24 >= 0) { /* left above */
1251: x_t = lx[n24 % m]*dof;
1252: y_t = ly[(n24 % (m*n))/m];
1253: /* z_t = lz[n24 / (m*n)]; */
1254: s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
1255: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1256: }
1257: if (n25 >= 0) { /* directly above */
1258: x_t = x;
1259: y_t = ly[(n25 % (m*n))/m];
1260: /* z_t = lz[n25 / (m*n)]; */
1261: s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1262: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1263: }
1264: if (n26 >= 0) { /* right above */
1265: x_t = lx[n26 % m]*dof;
1266: y_t = ly[(n26 % (m*n))/m];
1267: /* z_t = lz[n26 / (m*n)]; */
1268: s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1269: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1270: }
1271: }
1272: }
1273: }
1274: /* redo idx to include "missing" ghost points */
1275: /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
1276:
1277: /* Assume Nodes are Internal to the Cube */
1278:
1279: n0 = rank - m*n - m - 1;
1280: n1 = rank - m*n - m;
1281: n2 = rank - m*n - m + 1;
1282: n3 = rank - m*n -1;
1283: n4 = rank - m*n;
1284: n5 = rank - m*n + 1;
1285: n6 = rank - m*n + m - 1;
1286: n7 = rank - m*n + m;
1287: n8 = rank - m*n + m + 1;
1289: n9 = rank - m - 1;
1290: n10 = rank - m;
1291: n11 = rank - m + 1;
1292: n12 = rank - 1;
1293: n14 = rank + 1;
1294: n15 = rank + m - 1;
1295: n16 = rank + m;
1296: n17 = rank + m + 1;
1298: n18 = rank + m*n - m - 1;
1299: n19 = rank + m*n - m;
1300: n20 = rank + m*n - m + 1;
1301: n21 = rank + m*n - 1;
1302: n22 = rank + m*n;
1303: n23 = rank + m*n + 1;
1304: n24 = rank + m*n + m - 1;
1305: n25 = rank + m*n + m;
1306: n26 = rank + m*n + m + 1;
1308: /* Assume Pieces are on Faces of Cube */
1310: if (xs == 0) { /* First assume not corner or edge */
1311: n0 = rank -1 - (m*n);
1312: n3 = rank + m -1 - (m*n);
1313: n6 = rank + 2*m -1 - (m*n);
1314: n9 = rank -1;
1315: n12 = rank + m -1;
1316: n15 = rank + 2*m -1;
1317: n18 = rank -1 + (m*n);
1318: n21 = rank + m -1 + (m*n);
1319: n24 = rank + 2*m -1 + (m*n);
1320: }
1322: if (xe == M*dof) { /* First assume not corner or edge */
1323: n2 = rank -2*m +1 - (m*n);
1324: n5 = rank - m +1 - (m*n);
1325: n8 = rank +1 - (m*n);
1326: n11 = rank -2*m +1;
1327: n14 = rank - m +1;
1328: n17 = rank +1;
1329: n20 = rank -2*m +1 + (m*n);
1330: n23 = rank - m +1 + (m*n);
1331: n26 = rank +1 + (m*n);
1332: }
1334: if (ys==0) { /* First assume not corner or edge */
1335: n0 = rank + m * (n-1) -1 - (m*n);
1336: n1 = rank + m * (n-1) - (m*n);
1337: n2 = rank + m * (n-1) +1 - (m*n);
1338: n9 = rank + m * (n-1) -1;
1339: n10 = rank + m * (n-1);
1340: n11 = rank + m * (n-1) +1;
1341: n18 = rank + m * (n-1) -1 + (m*n);
1342: n19 = rank + m * (n-1) + (m*n);
1343: n20 = rank + m * (n-1) +1 + (m*n);
1344: }
1346: if (ye == N) { /* First assume not corner or edge */
1347: n6 = rank - m * (n-1) -1 - (m*n);
1348: n7 = rank - m * (n-1) - (m*n);
1349: n8 = rank - m * (n-1) +1 - (m*n);
1350: n15 = rank - m * (n-1) -1;
1351: n16 = rank - m * (n-1);
1352: n17 = rank - m * (n-1) +1;
1353: n24 = rank - m * (n-1) -1 + (m*n);
1354: n25 = rank - m * (n-1) + (m*n);
1355: n26 = rank - m * (n-1) +1 + (m*n);
1356: }
1357:
1358: if (zs == 0) { /* First assume not corner or edge */
1359: n0 = size - (m*n) + rank - m - 1;
1360: n1 = size - (m*n) + rank - m;
1361: n2 = size - (m*n) + rank - m + 1;
1362: n3 = size - (m*n) + rank - 1;
1363: n4 = size - (m*n) + rank;
1364: n5 = size - (m*n) + rank + 1;
1365: n6 = size - (m*n) + rank + m - 1;
1366: n7 = size - (m*n) + rank + m ;
1367: n8 = size - (m*n) + rank + m + 1;
1368: }
1370: if (ze == P) { /* First assume not corner or edge */
1371: n18 = (m*n) - (size-rank) - m - 1;
1372: n19 = (m*n) - (size-rank) - m;
1373: n20 = (m*n) - (size-rank) - m + 1;
1374: n21 = (m*n) - (size-rank) - 1;
1375: n22 = (m*n) - (size-rank);
1376: n23 = (m*n) - (size-rank) + 1;
1377: n24 = (m*n) - (size-rank) + m - 1;
1378: n25 = (m*n) - (size-rank) + m;
1379: n26 = (m*n) - (size-rank) + m + 1;
1380: }
1382: if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
1383: n0 = size - m*n + rank + m-1 - m;
1384: n3 = size - m*n + rank + m-1;
1385: n6 = size - m*n + rank + m-1 + m;
1386: }
1387:
1388: if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
1389: n18 = m*n - (size - rank) + m-1 - m;
1390: n21 = m*n - (size - rank) + m-1;
1391: n24 = m*n - (size - rank) + m-1 + m;
1392: }
1394: if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
1395: n0 = rank + m*n -1 - m*n;
1396: n9 = rank + m*n -1;
1397: n18 = rank + m*n -1 + m*n;
1398: }
1400: if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
1401: n6 = rank - m*(n-1) + m-1 - m*n;
1402: n15 = rank - m*(n-1) + m-1;
1403: n24 = rank - m*(n-1) + m-1 + m*n;
1404: }
1406: if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
1407: n2 = size - (m*n-rank) - (m-1) - m;
1408: n5 = size - (m*n-rank) - (m-1);
1409: n8 = size - (m*n-rank) - (m-1) + m;
1410: }
1412: if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
1413: n20 = m*n - (size - rank) - (m-1) - m;
1414: n23 = m*n - (size - rank) - (m-1);
1415: n26 = m*n - (size - rank) - (m-1) + m;
1416: }
1418: if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
1419: n2 = rank + m*(n-1) - (m-1) - m*n;
1420: n11 = rank + m*(n-1) - (m-1);
1421: n20 = rank + m*(n-1) - (m-1) + m*n;
1422: }
1424: if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
1425: n8 = rank - m*n +1 - m*n;
1426: n17 = rank - m*n +1;
1427: n26 = rank - m*n +1 + m*n;
1428: }
1430: if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
1431: n0 = size - m + rank -1;
1432: n1 = size - m + rank;
1433: n2 = size - m + rank +1;
1434: }
1436: if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
1437: n18 = m*n - (size - rank) + m*(n-1) -1;
1438: n19 = m*n - (size - rank) + m*(n-1);
1439: n20 = m*n - (size - rank) + m*(n-1) +1;
1440: }
1442: if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
1443: n6 = size - (m*n-rank) - m * (n-1) -1;
1444: n7 = size - (m*n-rank) - m * (n-1);
1445: n8 = size - (m*n-rank) - m * (n-1) +1;
1446: }
1448: if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
1449: n24 = rank - (size-m) -1;
1450: n25 = rank - (size-m);
1451: n26 = rank - (size-m) +1;
1452: }
1454: /* Check for Corners */
1455: if ((xs==0) && (ys==0) && (zs==0)) { n0 = size -1;}
1456: if ((xs==0) && (ys==0) && (ze==P)) { n18 = m*n-1;}
1457: if ((xs==0) && (ye==N) && (zs==0)) { n6 = (size-1)-m*(n-1);}
1458: if ((xs==0) && (ye==N) && (ze==P)) { n24 = m-1;}
1459: if ((xe==M*dof) && (ys==0) && (zs==0)) { n2 = size-m;}
1460: if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
1461: if ((xe==M*dof) && (ye==N) && (zs==0)) { n8 = size-m*n;}
1462: if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}
1464: /* Check for when not X,Y, and Z Periodic */
1466: /* If not X periodic */
1467: if (!DAXPeriodic(wrap)){
1468: if (xs==0) {n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;}
1469: if (xe==M*dof) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
1470: }
1472: /* If not Y periodic */
1473: if (!DAYPeriodic(wrap)){
1474: if (ys==0) {n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;}
1475: if (ye==N) {n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
1476: }
1478: /* If not Z periodic */
1479: if (!DAZPeriodic(wrap)){
1480: if (zs==0) {n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;}
1481: if (ze==P) {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
1482: }
1484: nn = 0;
1486: /* Bottom Level */
1487: for (k=0; k<s_z; k++) {
1488: for (i=1; i<=s_y; i++) {
1489: if (n0 >= 0) { /* left below */
1490: x_t = lx[n0 % m]*dof;
1491: y_t = ly[(n0 % (m*n))/m];
1492: z_t = lz[n0 / (m*n)];
1493: 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;
1494: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1495: }
1496: if (n1 >= 0) { /* directly below */
1497: x_t = x;
1498: y_t = ly[(n1 % (m*n))/m];
1499: z_t = lz[n1 / (m*n)];
1500: s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1501: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1502: }
1503: if (n2 >= 0) { /* right below */
1504: x_t = lx[n2 % m]*dof;
1505: y_t = ly[(n2 % (m*n))/m];
1506: z_t = lz[n2 / (m*n)];
1507: s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
1508: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1509: }
1510: }
1512: for (i=0; i<y; i++) {
1513: if (n3 >= 0) { /* directly left */
1514: x_t = lx[n3 % m]*dof;
1515: y_t = y;
1516: z_t = lz[n3 / (m*n)];
1517: s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1518: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1519: }
1521: if (n4 >= 0) { /* middle */
1522: x_t = x;
1523: y_t = y;
1524: z_t = lz[n4 / (m*n)];
1525: s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1526: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1527: }
1529: if (n5 >= 0) { /* directly right */
1530: x_t = lx[n5 % m]*dof;
1531: y_t = y;
1532: z_t = lz[n5 / (m*n)];
1533: s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1534: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1535: }
1536: }
1538: for (i=1; i<=s_y; i++) {
1539: if (n6 >= 0) { /* left above */
1540: x_t = lx[n6 % m]*dof;
1541: y_t = ly[(n6 % (m*n))/m];
1542: z_t = lz[n6 / (m*n)];
1543: s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1544: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1545: }
1546: if (n7 >= 0) { /* directly above */
1547: x_t = x;
1548: y_t = ly[(n7 % (m*n))/m];
1549: z_t = lz[n7 / (m*n)];
1550: s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1551: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1552: }
1553: if (n8 >= 0) { /* right above */
1554: x_t = lx[n8 % m]*dof;
1555: y_t = ly[(n8 % (m*n))/m];
1556: z_t = lz[n8 / (m*n)];
1557: s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
1558: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1559: }
1560: }
1561: }
1563: /* Middle Level */
1564: for (k=0; k<z; k++) {
1565: for (i=1; i<=s_y; i++) {
1566: if (n9 >= 0) { /* left below */
1567: x_t = lx[n9 % m]*dof;
1568: y_t = ly[(n9 % (m*n))/m];
1569: /* z_t = z; */
1570: s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1571: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1572: }
1573: if (n10 >= 0) { /* directly below */
1574: x_t = x;
1575: y_t = ly[(n10 % (m*n))/m];
1576: /* z_t = z; */
1577: s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1578: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1579: }
1580: if (n11 >= 0) { /* right below */
1581: x_t = lx[n11 % m]*dof;
1582: y_t = ly[(n11 % (m*n))/m];
1583: /* z_t = z; */
1584: s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1585: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1586: }
1587: }
1589: for (i=0; i<y; i++) {
1590: if (n12 >= 0) { /* directly left */
1591: x_t = lx[n12 % m]*dof;
1592: y_t = y;
1593: /* z_t = z; */
1594: s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
1595: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1596: }
1598: /* Interior */
1599: s_t = bases[rank] + i*x + k*x*y;
1600: for (j=0; j<x; j++) { idx[nn++] = s_t++;}
1602: if (n14 >= 0) { /* directly right */
1603: x_t = lx[n14 % m]*dof;
1604: y_t = y;
1605: /* z_t = z; */
1606: s_t = bases[n14] + i*x_t + k*x_t*y_t;
1607: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1608: }
1609: }
1611: for (i=1; i<=s_y; i++) {
1612: if (n15 >= 0) { /* left above */
1613: x_t = lx[n15 % m]*dof;
1614: y_t = ly[(n15 % (m*n))/m];
1615: /* z_t = z; */
1616: s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
1617: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1618: }
1619: if (n16 >= 0) { /* directly above */
1620: x_t = x;
1621: y_t = ly[(n16 % (m*n))/m];
1622: /* z_t = z; */
1623: s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
1624: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1625: }
1626: if (n17 >= 0) { /* right above */
1627: x_t = lx[n17 % m]*dof;
1628: y_t = ly[(n17 % (m*n))/m];
1629: /* z_t = z; */
1630: s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
1631: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1632: }
1633: }
1634: }
1635:
1636: /* Upper Level */
1637: for (k=0; k<s_z; k++) {
1638: for (i=1; i<=s_y; i++) {
1639: if (n18 >= 0) { /* left below */
1640: x_t = lx[n18 % m]*dof;
1641: y_t = ly[(n18 % (m*n))/m];
1642: /* z_t = lz[n18 / (m*n)]; */
1643: s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
1644: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1645: }
1646: if (n19 >= 0) { /* directly below */
1647: x_t = x;
1648: y_t = ly[(n19 % (m*n))/m];
1649: /* z_t = lz[n19 / (m*n)]; */
1650: s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1651: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1652: }
1653: if (n20 >= 0) { /* right belodof */
1654: x_t = lx[n20 % m]*dof;
1655: y_t = ly[(n20 % (m*n))/m];
1656: /* z_t = lz[n20 / (m*n)]; */
1657: s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
1658: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1659: }
1660: }
1662: for (i=0; i<y; i++) {
1663: if (n21 >= 0) { /* directly left */
1664: x_t = lx[n21 % m]*dof;
1665: y_t = y;
1666: /* z_t = lz[n21 / (m*n)]; */
1667: s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
1668: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1669: }
1671: if (n22 >= 0) { /* middle */
1672: x_t = x;
1673: y_t = y;
1674: /* z_t = lz[n22 / (m*n)]; */
1675: s_t = bases[n22] + i*x_t + k*x_t*y_t;
1676: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1677: }
1679: if (n23 >= 0) { /* directly right */
1680: x_t = lx[n23 % m]*dof;
1681: y_t = y;
1682: /* z_t = lz[n23 / (m*n)]; */
1683: s_t = bases[n23] + i*x_t + k*x_t*y_t;
1684: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1685: }
1686: }
1688: for (i=1; i<=s_y; i++) {
1689: if (n24 >= 0) { /* left above */
1690: x_t = lx[n24 % m]*dof;
1691: y_t = ly[(n24 % (m*n))/m];
1692: /* z_t = lz[n24 / (m*n)]; */
1693: s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
1694: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1695: }
1696: if (n25 >= 0) { /* directly above */
1697: x_t = x;
1698: y_t = ly[(n25 % (m*n))/m];
1699: /* z_t = lz[n25 / (m*n)]; */
1700: s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
1701: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1702: }
1703: if (n26 >= 0) { /* right above */
1704: x_t = lx[n26 % m]*dof;
1705: y_t = ly[(n26 % (m*n))/m];
1706: /* z_t = lz[n26 / (m*n)]; */
1707: s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
1708: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1709: }
1710: }
1711: }
1712: PetscFree(bases);
1713: da->gtol = gtol;
1714: da->ltog = ltog;
1715: da->idx = idx;
1716: da->Nl = nn;
1717: da->base = base;
1718: da->ops->view = DAView_3d;
1719: da->wrap = wrap;
1720: *inra = da;
1722: /*
1723: Set the local to global ordering in the global vector, this allows use
1724: of VecSetValuesLocal().
1725: */
1726: ierr = ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
1727: ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
1728: PetscLogObjectParent(da,da->ltogmap);
1730: da->ltol = PETSC_NULL;
1731: da->ao = PETSC_NULL;
1733: if (!flx) {
1734: PetscMalloc(m*sizeof(int),&flx);
1735: PetscMemcpy(flx,lx,m*sizeof(int));
1736: }
1737: if (!fly) {
1738: PetscMalloc(n*sizeof(int),&fly);
1739: PetscMemcpy(fly,ly,n*sizeof(int));
1740: }
1741: if (!flz) {
1742: PetscMalloc(p*sizeof(int),&flz);
1743: PetscMemcpy(flz,lz,p*sizeof(int));
1744: }
1745: da->lx = flx;
1746: da->ly = fly;
1747: da->lz = flz;
1749: PetscOptionsHasName(PETSC_NULL,"-da_view",&flg1);
1750: if (flg1) {DAView(da,PETSC_VIEWER_STDOUT_(da->comm));}
1751: PetscOptionsHasName(PETSC_NULL,"-da_view_draw",&flg1);
1752: if (flg1) {DAView(da,PETSC_VIEWER_DRAW_(da->comm));}
1753: PetscOptionsHasName(PETSC_NULL,"-help",&flg1);
1754: if (flg1) {DAPrintHelp(da);}
1755: PetscPublishAll(da);
1757: return(0);
1758: }