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