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