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