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