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