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 count = x*y*z; 461 ierr = PetscMalloc(x*y*z*sizeof(PetscInt),&idx);CHKERRQ(ierr); 462 left = xs - Xs; right = left + x; 463 bottom = ys - Ys; top = bottom + y; 464 down = zs - Zs; up = down + z; 465 count = 0; 466 for (i=down; i<up; i++) { 467 for (j=bottom; j<top; j++) { 468 for (k=left; k<right; k++) { 469 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 470 } 471 } 472 } 473 474 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 475 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 476 ierr = PetscLogObjectParent(da,ltog);CHKERRQ(ierr); 477 ierr = ISDestroy(&from);CHKERRQ(ierr); 478 ierr = ISDestroy(&to);CHKERRQ(ierr); 479 480 /* global to local must include ghost points within the domain, 481 but not ghost points outside the domain that aren't periodic */ 482 if (stencil_type == DMDA_STENCIL_BOX || o > 0) { 483 count = (IXe-IXs)*(IYe-IYs)*(IZe-IZs); 484 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 485 486 left = IXs - Xs; right = left + (IXe-IXs); 487 bottom = IYs - Ys; top = bottom + (IYe-IYs); 488 down = IZs - Zs; up = down + (IZe-IZs); 489 count = 0; 490 for (i=down; i<up; i++) { 491 for (j=bottom; j<top; j++) { 492 for (k=left; k<right; k++) { 493 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 494 } 495 } 496 } 497 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 498 499 } else { 500 /* This is way ugly! We need to list the funny cross type region */ 501 count = ((ys-IYs) + (IYe-ye))*x*z + ((xs-IXs) + (IXe-xe))*y*z + ((zs-IZs) + (IZe-ze))*x*y + x*y*z; 502 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 503 504 left = xs - Xs; right = left + x; 505 bottom = ys - Ys; top = bottom + y; 506 down = zs - Zs; up = down + z; 507 count = 0; 508 /* the bottom chunck */ 509 for (i=(IZs-Zs); i<down; i++) { 510 for (j=bottom; j<top; j++) { 511 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 512 } 513 } 514 /* the middle piece */ 515 for (i=down; i<up; i++) { 516 /* front */ 517 for (j=(IYs-Ys); j<bottom; j++) { 518 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 519 } 520 /* middle */ 521 for (j=bottom; j<top; j++) { 522 for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 523 } 524 /* back */ 525 for (j=top; j<top+IYe-ye; j++) { 526 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 527 } 528 } 529 /* the top piece */ 530 for (i=up; i<up+IZe-ze; i++) { 531 for (j=bottom; j<top; j++) { 532 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 533 } 534 } 535 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 536 } 537 538 /* determine who lies on each side of use stored in n24 n25 n26 539 n21 n22 n23 540 n18 n19 n20 541 542 n15 n16 n17 543 n12 n14 544 n9 n10 n11 545 546 n6 n7 n8 547 n3 n4 n5 548 n0 n1 n2 549 */ 550 551 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ 552 /* Assume Nodes are Internal to the Cube */ 553 n0 = rank - m*n - m - 1; 554 n1 = rank - m*n - m; 555 n2 = rank - m*n - m + 1; 556 n3 = rank - m*n -1; 557 n4 = rank - m*n; 558 n5 = rank - m*n + 1; 559 n6 = rank - m*n + m - 1; 560 n7 = rank - m*n + m; 561 n8 = rank - m*n + m + 1; 562 563 n9 = rank - m - 1; 564 n10 = rank - m; 565 n11 = rank - m + 1; 566 n12 = rank - 1; 567 n14 = rank + 1; 568 n15 = rank + m - 1; 569 n16 = rank + m; 570 n17 = rank + m + 1; 571 572 n18 = rank + m*n - m - 1; 573 n19 = rank + m*n - m; 574 n20 = rank + m*n - m + 1; 575 n21 = rank + m*n - 1; 576 n22 = rank + m*n; 577 n23 = rank + m*n + 1; 578 n24 = rank + m*n + m - 1; 579 n25 = rank + m*n + m; 580 n26 = rank + m*n + m + 1; 581 582 /* Assume Pieces are on Faces of Cube */ 583 584 if (xs == 0) { /* First assume not corner or edge */ 585 n0 = rank -1 - (m*n); 586 n3 = rank + m -1 - (m*n); 587 n6 = rank + 2*m -1 - (m*n); 588 n9 = rank -1; 589 n12 = rank + m -1; 590 n15 = rank + 2*m -1; 591 n18 = rank -1 + (m*n); 592 n21 = rank + m -1 + (m*n); 593 n24 = rank + 2*m -1 + (m*n); 594 } 595 596 if (xe == M) { /* First assume not corner or edge */ 597 n2 = rank -2*m +1 - (m*n); 598 n5 = rank - m +1 - (m*n); 599 n8 = rank +1 - (m*n); 600 n11 = rank -2*m +1; 601 n14 = rank - m +1; 602 n17 = rank +1; 603 n20 = rank -2*m +1 + (m*n); 604 n23 = rank - m +1 + (m*n); 605 n26 = rank +1 + (m*n); 606 } 607 608 if (ys==0) { /* First assume not corner or edge */ 609 n0 = rank + m * (n-1) -1 - (m*n); 610 n1 = rank + m * (n-1) - (m*n); 611 n2 = rank + m * (n-1) +1 - (m*n); 612 n9 = rank + m * (n-1) -1; 613 n10 = rank + m * (n-1); 614 n11 = rank + m * (n-1) +1; 615 n18 = rank + m * (n-1) -1 + (m*n); 616 n19 = rank + m * (n-1) + (m*n); 617 n20 = rank + m * (n-1) +1 + (m*n); 618 } 619 620 if (ye == N) { /* First assume not corner or edge */ 621 n6 = rank - m * (n-1) -1 - (m*n); 622 n7 = rank - m * (n-1) - (m*n); 623 n8 = rank - m * (n-1) +1 - (m*n); 624 n15 = rank - m * (n-1) -1; 625 n16 = rank - m * (n-1); 626 n17 = rank - m * (n-1) +1; 627 n24 = rank - m * (n-1) -1 + (m*n); 628 n25 = rank - m * (n-1) + (m*n); 629 n26 = rank - m * (n-1) +1 + (m*n); 630 } 631 632 if (zs == 0) { /* First assume not corner or edge */ 633 n0 = size - (m*n) + rank - m - 1; 634 n1 = size - (m*n) + rank - m; 635 n2 = size - (m*n) + rank - m + 1; 636 n3 = size - (m*n) + rank - 1; 637 n4 = size - (m*n) + rank; 638 n5 = size - (m*n) + rank + 1; 639 n6 = size - (m*n) + rank + m - 1; 640 n7 = size - (m*n) + rank + m ; 641 n8 = size - (m*n) + rank + m + 1; 642 } 643 644 if (ze == P) { /* First assume not corner or edge */ 645 n18 = (m*n) - (size-rank) - m - 1; 646 n19 = (m*n) - (size-rank) - m; 647 n20 = (m*n) - (size-rank) - m + 1; 648 n21 = (m*n) - (size-rank) - 1; 649 n22 = (m*n) - (size-rank); 650 n23 = (m*n) - (size-rank) + 1; 651 n24 = (m*n) - (size-rank) + m - 1; 652 n25 = (m*n) - (size-rank) + m; 653 n26 = (m*n) - (size-rank) + m + 1; 654 } 655 656 if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */ 657 n0 = size - m*n + rank + m-1 - m; 658 n3 = size - m*n + rank + m-1; 659 n6 = size - m*n + rank + m-1 + m; 660 } 661 662 if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */ 663 n18 = m*n - (size - rank) + m-1 - m; 664 n21 = m*n - (size - rank) + m-1; 665 n24 = m*n - (size - rank) + m-1 + m; 666 } 667 668 if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */ 669 n0 = rank + m*n -1 - m*n; 670 n9 = rank + m*n -1; 671 n18 = rank + m*n -1 + m*n; 672 } 673 674 if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */ 675 n6 = rank - m*(n-1) + m-1 - m*n; 676 n15 = rank - m*(n-1) + m-1; 677 n24 = rank - m*(n-1) + m-1 + m*n; 678 } 679 680 if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */ 681 n2 = size - (m*n-rank) - (m-1) - m; 682 n5 = size - (m*n-rank) - (m-1); 683 n8 = size - (m*n-rank) - (m-1) + m; 684 } 685 686 if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */ 687 n20 = m*n - (size - rank) - (m-1) - m; 688 n23 = m*n - (size - rank) - (m-1); 689 n26 = m*n - (size - rank) - (m-1) + m; 690 } 691 692 if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */ 693 n2 = rank + m*(n-1) - (m-1) - m*n; 694 n11 = rank + m*(n-1) - (m-1); 695 n20 = rank + m*(n-1) - (m-1) + m*n; 696 } 697 698 if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */ 699 n8 = rank - m*n +1 - m*n; 700 n17 = rank - m*n +1; 701 n26 = rank - m*n +1 + m*n; 702 } 703 704 if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */ 705 n0 = size - m + rank -1; 706 n1 = size - m + rank; 707 n2 = size - m + rank +1; 708 } 709 710 if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */ 711 n18 = m*n - (size - rank) + m*(n-1) -1; 712 n19 = m*n - (size - rank) + m*(n-1); 713 n20 = m*n - (size - rank) + m*(n-1) +1; 714 } 715 716 if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */ 717 n6 = size - (m*n-rank) - m * (n-1) -1; 718 n7 = size - (m*n-rank) - m * (n-1); 719 n8 = size - (m*n-rank) - m * (n-1) +1; 720 } 721 722 if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */ 723 n24 = rank - (size-m) -1; 724 n25 = rank - (size-m); 725 n26 = rank - (size-m) +1; 726 } 727 728 /* Check for Corners */ 729 if ((xs==0) && (ys==0) && (zs==0)) { n0 = size -1;} 730 if ((xs==0) && (ys==0) && (ze==P)) { n18 = m*n-1;} 731 if ((xs==0) && (ye==N) && (zs==0)) { n6 = (size-1)-m*(n-1);} 732 if ((xs==0) && (ye==N) && (ze==P)) { n24 = m-1;} 733 if ((xe==M) && (ys==0) && (zs==0)) { n2 = size-m;} 734 if ((xe==M) && (ys==0) && (ze==P)) { n20 = m*n-m;} 735 if ((xe==M) && (ye==N) && (zs==0)) { n8 = size-m*n;} 736 if ((xe==M) && (ye==N) && (ze==P)) { n26 = 0;} 737 738 /* Check for when not X,Y, and Z Periodic */ 739 740 /* If not X periodic */ 741 if (bx != DMDA_BOUNDARY_PERIODIC) { 742 if (xs==0) {n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;} 743 if (xe==M) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;} 744 } 745 746 /* If not Y periodic */ 747 if (by != DMDA_BOUNDARY_PERIODIC) { 748 if (ys==0) {n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;} 749 if (ye==N) {n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;} 750 } 751 752 /* If not Z periodic */ 753 if (bz != DMDA_BOUNDARY_PERIODIC) { 754 if (zs==0) {n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;} 755 if (ze==P) {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;} 756 } 757 758 ierr = PetscMalloc(27*sizeof(PetscInt),&dd->neighbors);CHKERRQ(ierr); 759 dd->neighbors[0] = n0; 760 dd->neighbors[1] = n1; 761 dd->neighbors[2] = n2; 762 dd->neighbors[3] = n3; 763 dd->neighbors[4] = n4; 764 dd->neighbors[5] = n5; 765 dd->neighbors[6] = n6; 766 dd->neighbors[7] = n7; 767 dd->neighbors[8] = n8; 768 dd->neighbors[9] = n9; 769 dd->neighbors[10] = n10; 770 dd->neighbors[11] = n11; 771 dd->neighbors[12] = n12; 772 dd->neighbors[13] = rank; 773 dd->neighbors[14] = n14; 774 dd->neighbors[15] = n15; 775 dd->neighbors[16] = n16; 776 dd->neighbors[17] = n17; 777 dd->neighbors[18] = n18; 778 dd->neighbors[19] = n19; 779 dd->neighbors[20] = n20; 780 dd->neighbors[21] = n21; 781 dd->neighbors[22] = n22; 782 dd->neighbors[23] = n23; 783 dd->neighbors[24] = n24; 784 dd->neighbors[25] = n25; 785 dd->neighbors[26] = n26; 786 787 /* If star stencil then delete the corner neighbors */ 788 if (stencil_type == DMDA_STENCIL_STAR && o == 0) { 789 /* save information about corner neighbors */ 790 sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7; 791 sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18; 792 sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25; 793 sn26 = n26; 794 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = 795 n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; 796 } 797 798 799 ierr = PetscMalloc((Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt),&idx);CHKERRQ(ierr); 800 ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt));CHKERRQ(ierr); 801 802 nn = 0; 803 /* Bottom Level */ 804 for (k=0; k<s_z; k++) { 805 for (i=1; i<=s_y; i++) { 806 if (n0 >= 0) { /* left below */ 807 x_t = lx[n0 % m]; 808 y_t = ly[(n0 % (m*n))/m]; 809 z_t = lz[n0 / (m*n)]; 810 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; 811 if (twod && (s_t < 0)) {s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x;} /* 2D case */ 812 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 813 } 814 if (n1 >= 0) { /* directly below */ 815 x_t = x; 816 y_t = ly[(n1 % (m*n))/m]; 817 z_t = lz[n1 / (m*n)]; 818 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 819 if (twod && (s_t < 0)) {s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t;} /* 2D case */ 820 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 821 } 822 if (n2 >= 0) { /* right below */ 823 x_t = lx[n2 % m]; 824 y_t = ly[(n2 % (m*n))/m]; 825 z_t = lz[n2 / (m*n)]; 826 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 827 if (twod && (s_t < 0)) {s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t;} /* 2D case */ 828 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 829 } 830 } 831 832 for (i=0; i<y; i++) { 833 if (n3 >= 0) { /* directly left */ 834 x_t = lx[n3 % m]; 835 y_t = y; 836 z_t = lz[n3 / (m*n)]; 837 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 838 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 */ 839 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 840 } 841 842 if (n4 >= 0) { /* middle */ 843 x_t = x; 844 y_t = y; 845 z_t = lz[n4 / (m*n)]; 846 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 847 if (twod && (s_t < 0)) {s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 848 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 849 } else if (bz == DMDA_BOUNDARY_MIRROR) { 850 for (j=0; j<x; j++) { idx[nn++] = 0;} 851 } 852 853 if (n5 >= 0) { /* directly right */ 854 x_t = lx[n5 % m]; 855 y_t = y; 856 z_t = lz[n5 / (m*n)]; 857 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 858 if (twod && (s_t < 0)) {s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 859 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 860 } 861 } 862 863 for (i=1; i<=s_y; i++) { 864 if (n6 >= 0) { /* left above */ 865 x_t = lx[n6 % m]; 866 y_t = ly[(n6 % (m*n))/m]; 867 z_t = lz[n6 / (m*n)]; 868 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 869 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 */ 870 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 871 } 872 if (n7 >= 0) { /* directly above */ 873 x_t = x; 874 y_t = ly[(n7 % (m*n))/m]; 875 z_t = lz[n7 / (m*n)]; 876 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 877 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 */ 878 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 879 } 880 if (n8 >= 0) { /* right above */ 881 x_t = lx[n8 % m]; 882 y_t = ly[(n8 % (m*n))/m]; 883 z_t = lz[n8 / (m*n)]; 884 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 885 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 */ 886 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 887 } 888 } 889 } 890 891 /* Middle Level */ 892 for (k=0; k<z; k++) { 893 for (i=1; i<=s_y; i++) { 894 if (n9 >= 0) { /* left below */ 895 x_t = lx[n9 % m]; 896 y_t = ly[(n9 % (m*n))/m]; 897 /* z_t = z; */ 898 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 899 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 900 } 901 if (n10 >= 0) { /* directly below */ 902 x_t = x; 903 y_t = ly[(n10 % (m*n))/m]; 904 /* z_t = z; */ 905 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 906 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 907 } else if (by == DMDA_BOUNDARY_MIRROR) { 908 for (j=0; j<x; j++) { idx[nn++] = 0;} 909 } 910 if (n11 >= 0) { /* right below */ 911 x_t = lx[n11 % m]; 912 y_t = ly[(n11 % (m*n))/m]; 913 /* z_t = z; */ 914 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 915 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 916 } 917 } 918 919 for (i=0; i<y; i++) { 920 if (n12 >= 0) { /* directly left */ 921 x_t = lx[n12 % m]; 922 y_t = y; 923 /* z_t = z; */ 924 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 925 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 926 } else if (bx == DMDA_BOUNDARY_MIRROR) { 927 for (j=0; j<s_x; j++) { idx[nn++] = 0;} 928 } 929 930 /* Interior */ 931 s_t = bases[rank] + i*x + k*x*y; 932 for (j=0; j<x; j++) { idx[nn++] = s_t++;} 933 934 if (n14 >= 0) { /* directly right */ 935 x_t = lx[n14 % m]; 936 y_t = y; 937 /* z_t = z; */ 938 s_t = bases[n14] + i*x_t + k*x_t*y_t; 939 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 940 } else if (bx == DMDA_BOUNDARY_MIRROR) { 941 for (j=0; j<s_x; j++) { idx[nn++] = 0;} 942 } 943 } 944 945 for (i=1; i<=s_y; i++) { 946 if (n15 >= 0) { /* left above */ 947 x_t = lx[n15 % m]; 948 y_t = ly[(n15 % (m*n))/m]; 949 /* z_t = z; */ 950 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 951 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 952 } 953 if (n16 >= 0) { /* directly above */ 954 x_t = x; 955 y_t = ly[(n16 % (m*n))/m]; 956 /* z_t = z; */ 957 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 958 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 959 } else if (by == DMDA_BOUNDARY_MIRROR) { 960 for (j=0; j<x; j++) { idx[nn++] = 0;} 961 } 962 if (n17 >= 0) { /* right above */ 963 x_t = lx[n17 % m]; 964 y_t = ly[(n17 % (m*n))/m]; 965 /* z_t = z; */ 966 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 967 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 968 } 969 } 970 } 971 972 /* Upper Level */ 973 for (k=0; k<s_z; k++) { 974 for (i=1; i<=s_y; i++) { 975 if (n18 >= 0) { /* left below */ 976 x_t = lx[n18 % m]; 977 y_t = ly[(n18 % (m*n))/m]; 978 /* z_t = lz[n18 / (m*n)]; */ 979 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 980 if (twod && (s_t >= M*N*P)) {s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t;} /* 2d case */ 981 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 982 } 983 if (n19 >= 0) { /* directly below */ 984 x_t = x; 985 y_t = ly[(n19 % (m*n))/m]; 986 /* z_t = lz[n19 / (m*n)]; */ 987 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 988 if (twod && (s_t >= M*N*P)) {s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t;} /* 2d case */ 989 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 990 } 991 if (n20 >= 0) { /* right below */ 992 x_t = lx[n20 % m]; 993 y_t = ly[(n20 % (m*n))/m]; 994 /* z_t = lz[n20 / (m*n)]; */ 995 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 996 if (twod && (s_t >= M*N*P)) {s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t;} /* 2d case */ 997 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 998 } 999 } 1000 1001 for (i=0; i<y; i++) { 1002 if (n21 >= 0) { /* directly left */ 1003 x_t = lx[n21 % m]; 1004 y_t = y; 1005 /* z_t = lz[n21 / (m*n)]; */ 1006 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1007 if (twod && (s_t >= M*N*P)) {s_t = bases[n21] + (i+1)*x_t - s_x;} /* 2d case */ 1008 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1009 } 1010 1011 if (n22 >= 0) { /* middle */ 1012 x_t = x; 1013 y_t = y; 1014 /* z_t = lz[n22 / (m*n)]; */ 1015 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1016 if (twod && (s_t >= M*N*P)) {s_t = bases[n22] + i*x_t;} /* 2d case */ 1017 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1018 } else if (bz == DMDA_BOUNDARY_MIRROR) { 1019 for (j=0; j<x; j++) { idx[nn++] = 0;} 1020 } 1021 1022 if (n23 >= 0) { /* directly right */ 1023 x_t = lx[n23 % m]; 1024 y_t = y; 1025 /* z_t = lz[n23 / (m*n)]; */ 1026 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1027 if (twod && (s_t >= M*N*P)) {s_t = bases[n23] + i*x_t;} /* 2d case */ 1028 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1029 } 1030 } 1031 1032 for (i=1; i<=s_y; i++) { 1033 if (n24 >= 0) { /* left above */ 1034 x_t = lx[n24 % m]; 1035 y_t = ly[(n24 % (m*n))/m]; 1036 /* z_t = lz[n24 / (m*n)]; */ 1037 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1038 if (twod && (s_t >= M*N*P)) {s_t = bases[n24] + i*x_t - s_x;} /* 2d case */ 1039 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1040 } 1041 if (n25 >= 0) { /* directly above */ 1042 x_t = x; 1043 y_t = ly[(n25 % (m*n))/m]; 1044 /* z_t = lz[n25 / (m*n)]; */ 1045 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1046 if (twod && (s_t >= M*N*P)) {s_t = bases[n25] + (i-1)*x_t;} /* 2d case */ 1047 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1048 } 1049 if (n26 >= 0) { /* right above */ 1050 x_t = lx[n26 % m]; 1051 y_t = ly[(n26 % (m*n))/m]; 1052 /* z_t = lz[n26 / (m*n)]; */ 1053 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1054 if (twod && (s_t >= M*N*P)) {s_t = bases[n26] + (i-1)*x_t;} /* 2d case */ 1055 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1056 } 1057 } 1058 } 1059 1060 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 1061 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1062 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 1063 ierr = ISDestroy(&to);CHKERRQ(ierr); 1064 ierr = ISDestroy(&from);CHKERRQ(ierr); 1065 1066 if (stencil_type == DMDA_STENCIL_STAR && o == 0) { 1067 n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; 1068 n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; 1069 n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; 1070 n26 = sn26; 1071 } 1072 1073 if (((stencil_type == DMDA_STENCIL_STAR) || 1074 (bx != DMDA_BOUNDARY_PERIODIC && bx) || 1075 (by != DMDA_BOUNDARY_PERIODIC && by) || 1076 (bz != DMDA_BOUNDARY_PERIODIC && bz)) && o == 0) { 1077 /* 1078 Recompute the local to global mappings, this time keeping the 1079 information about the cross corner processor numbers. 1080 */ 1081 nn = 0; 1082 /* Bottom Level */ 1083 for (k=0; k<s_z; k++) { 1084 for (i=1; i<=s_y; i++) { 1085 if (n0 >= 0) { /* left below */ 1086 x_t = lx[n0 % m]; 1087 y_t = ly[(n0 % (m*n))/m]; 1088 z_t = lz[n0 / (m*n)]; 1089 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; 1090 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1091 } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) { 1092 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1093 } 1094 if (n1 >= 0) { /* directly below */ 1095 x_t = x; 1096 y_t = ly[(n1 % (m*n))/m]; 1097 z_t = lz[n1 / (m*n)]; 1098 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1099 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1100 } else if (Ys-ys < 0 && Zs-zs < 0) { 1101 for (j=0; j<x; j++) { idx[nn++] = -1;} 1102 } 1103 if (n2 >= 0) { /* right below */ 1104 x_t = lx[n2 % m]; 1105 y_t = ly[(n2 % (m*n))/m]; 1106 z_t = lz[n2 / (m*n)]; 1107 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1108 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1109 } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) { 1110 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1111 } 1112 } 1113 1114 for (i=0; i<y; i++) { 1115 if (n3 >= 0) { /* directly left */ 1116 x_t = lx[n3 % m]; 1117 y_t = y; 1118 z_t = lz[n3 / (m*n)]; 1119 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1120 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1121 } else if (Xs-xs < 0 && Zs-zs < 0) { 1122 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1123 } 1124 1125 if (n4 >= 0) { /* middle */ 1126 x_t = x; 1127 y_t = y; 1128 z_t = lz[n4 / (m*n)]; 1129 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1130 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1131 } else if (Zs-zs < 0) { 1132 if (bz == DMDA_BOUNDARY_MIRROR) { 1133 for (j=0; j<x; j++) { idx[nn++] = 0;} 1134 } else { 1135 for (j=0; j<x; j++) { idx[nn++] = -1;} 1136 } 1137 } 1138 1139 if (n5 >= 0) { /* directly right */ 1140 x_t = lx[n5 % m]; 1141 y_t = y; 1142 z_t = lz[n5 / (m*n)]; 1143 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1144 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1145 } else if (xe-Xe < 0 && Zs-zs < 0) { 1146 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1147 } 1148 } 1149 1150 for (i=1; i<=s_y; i++) { 1151 if (n6 >= 0) { /* left above */ 1152 x_t = lx[n6 % m]; 1153 y_t = ly[(n6 % (m*n))/m]; 1154 z_t = lz[n6 / (m*n)]; 1155 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1156 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1157 } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) { 1158 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1159 } 1160 if (n7 >= 0) { /* directly above */ 1161 x_t = x; 1162 y_t = ly[(n7 % (m*n))/m]; 1163 z_t = lz[n7 / (m*n)]; 1164 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1165 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1166 } else if (ye-Ye < 0 && Zs-zs < 0) { 1167 for (j=0; j<x; j++) { idx[nn++] = -1;} 1168 } 1169 if (n8 >= 0) { /* right above */ 1170 x_t = lx[n8 % m]; 1171 y_t = ly[(n8 % (m*n))/m]; 1172 z_t = lz[n8 / (m*n)]; 1173 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1174 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1175 } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) { 1176 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1177 } 1178 } 1179 } 1180 1181 /* Middle Level */ 1182 for (k=0; k<z; k++) { 1183 for (i=1; i<=s_y; i++) { 1184 if (n9 >= 0) { /* left below */ 1185 x_t = lx[n9 % m]; 1186 y_t = ly[(n9 % (m*n))/m]; 1187 /* z_t = z; */ 1188 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1189 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1190 } else if (Xs-xs < 0 && Ys-ys < 0) { 1191 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1192 } 1193 if (n10 >= 0) { /* directly below */ 1194 x_t = x; 1195 y_t = ly[(n10 % (m*n))/m]; 1196 /* z_t = z; */ 1197 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1198 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1199 } else if (Ys-ys < 0) { 1200 if (by == DMDA_BOUNDARY_MIRROR) { 1201 for (j=0; j<x; j++) { idx[nn++] = -1;} 1202 } else { 1203 for (j=0; j<x; j++) { idx[nn++] = -1;} 1204 } 1205 } 1206 if (n11 >= 0) { /* right below */ 1207 x_t = lx[n11 % m]; 1208 y_t = ly[(n11 % (m*n))/m]; 1209 /* z_t = z; */ 1210 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1211 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1212 } else if (xe-Xe < 0 && Ys-ys < 0) { 1213 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1214 } 1215 } 1216 1217 for (i=0; i<y; i++) { 1218 if (n12 >= 0) { /* directly left */ 1219 x_t = lx[n12 % m]; 1220 y_t = y; 1221 /* z_t = z; */ 1222 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 1223 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1224 } else if (Xs-xs < 0) { 1225 if (bx == DMDA_BOUNDARY_MIRROR) { 1226 for (j=0; j<s_x; j++) { idx[nn++] = 0;} 1227 } else { 1228 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1229 } 1230 } 1231 1232 /* Interior */ 1233 s_t = bases[rank] + i*x + k*x*y; 1234 for (j=0; j<x; j++) { idx[nn++] = s_t++;} 1235 1236 if (n14 >= 0) { /* directly right */ 1237 x_t = lx[n14 % m]; 1238 y_t = y; 1239 /* z_t = z; */ 1240 s_t = bases[n14] + i*x_t + k*x_t*y_t; 1241 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1242 } else if (xe-Xe < 0) { 1243 if (bx == DMDA_BOUNDARY_MIRROR) { 1244 for (j=0; j<s_x; j++) { idx[nn++] = 0;} 1245 } else { 1246 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1247 } 1248 } 1249 } 1250 1251 for (i=1; i<=s_y; i++) { 1252 if (n15 >= 0) { /* left above */ 1253 x_t = lx[n15 % m]; 1254 y_t = ly[(n15 % (m*n))/m]; 1255 /* z_t = z; */ 1256 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 1257 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1258 } else if (Xs-xs < 0 && ye-Ye < 0) { 1259 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1260 } 1261 if (n16 >= 0) { /* directly above */ 1262 x_t = x; 1263 y_t = ly[(n16 % (m*n))/m]; 1264 /* z_t = z; */ 1265 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 1266 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1267 } else if (ye-Ye < 0) { 1268 if (by == DMDA_BOUNDARY_MIRROR) { 1269 for (j=0; j<x; j++) { idx[nn++] = 0;} 1270 } else { 1271 for (j=0; j<x; j++) { idx[nn++] = -1;} 1272 } 1273 } 1274 if (n17 >= 0) { /* right above */ 1275 x_t = lx[n17 % m]; 1276 y_t = ly[(n17 % (m*n))/m]; 1277 /* z_t = z; */ 1278 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 1279 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1280 } else if (xe-Xe < 0 && ye-Ye < 0) { 1281 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1282 } 1283 } 1284 } 1285 1286 /* Upper Level */ 1287 for (k=0; k<s_z; k++) { 1288 for (i=1; i<=s_y; i++) { 1289 if (n18 >= 0) { /* left below */ 1290 x_t = lx[n18 % m]; 1291 y_t = ly[(n18 % (m*n))/m]; 1292 /* z_t = lz[n18 / (m*n)]; */ 1293 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1294 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1295 } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) { 1296 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1297 } 1298 if (n19 >= 0) { /* directly below */ 1299 x_t = x; 1300 y_t = ly[(n19 % (m*n))/m]; 1301 /* z_t = lz[n19 / (m*n)]; */ 1302 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1303 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1304 } else if (Ys-ys < 0 && ze-Ze < 0) { 1305 for (j=0; j<x; j++) { idx[nn++] = -1;} 1306 } 1307 if (n20 >= 0) { /* right below */ 1308 x_t = lx[n20 % m]; 1309 y_t = ly[(n20 % (m*n))/m]; 1310 /* z_t = lz[n20 / (m*n)]; */ 1311 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1312 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1313 } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) { 1314 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1315 } 1316 } 1317 1318 for (i=0; i<y; i++) { 1319 if (n21 >= 0) { /* directly left */ 1320 x_t = lx[n21 % m]; 1321 y_t = y; 1322 /* z_t = lz[n21 / (m*n)]; */ 1323 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1324 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1325 } else if (Xs-xs < 0 && ze-Ze < 0) { 1326 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1327 } 1328 1329 if (n22 >= 0) { /* middle */ 1330 x_t = x; 1331 y_t = y; 1332 /* z_t = lz[n22 / (m*n)]; */ 1333 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1334 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1335 } else if (ze-Ze < 0) { 1336 if (bz == DMDA_BOUNDARY_MIRROR) { 1337 for (j=0; j<x; j++) { idx[nn++] = 0;} 1338 } else { 1339 for (j=0; j<x; j++) { idx[nn++] = -1;} 1340 } 1341 } 1342 1343 if (n23 >= 0) { /* directly right */ 1344 x_t = lx[n23 % m]; 1345 y_t = y; 1346 /* z_t = lz[n23 / (m*n)]; */ 1347 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1348 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1349 } else if (xe-Xe < 0 && ze-Ze < 0) { 1350 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1351 } 1352 } 1353 1354 for (i=1; i<=s_y; i++) { 1355 if (n24 >= 0) { /* left above */ 1356 x_t = lx[n24 % m]; 1357 y_t = ly[(n24 % (m*n))/m]; 1358 /* z_t = lz[n24 / (m*n)]; */ 1359 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1360 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1361 } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) { 1362 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1363 } 1364 if (n25 >= 0) { /* directly above */ 1365 x_t = x; 1366 y_t = ly[(n25 % (m*n))/m]; 1367 /* z_t = lz[n25 / (m*n)]; */ 1368 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1369 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1370 } else if (ye-Ye < 0 && ze-Ze < 0) { 1371 for (j=0; j<x; j++) { idx[nn++] = -1;} 1372 } 1373 if (n26 >= 0) { /* right above */ 1374 x_t = lx[n26 % m]; 1375 y_t = ly[(n26 % (m*n))/m]; 1376 /* z_t = lz[n26 / (m*n)]; */ 1377 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1378 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1379 } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) { 1380 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1381 } 1382 } 1383 } 1384 } 1385 /* 1386 Set the local to global ordering in the global vector, this allows use 1387 of VecSetValuesLocal(). 1388 */ 1389 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr); 1390 ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr); 1391 ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1392 ierr = ISGetIndices(ltogis, &idx_full); 1393 ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1394 ierr = ISRestoreIndices(ltogis, &idx_full); 1395 ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr); 1396 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1397 ierr = ISDestroy(<ogis);CHKERRQ(ierr); 1398 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 1399 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1400 1401 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1402 dd->m = m; dd->n = n; dd->p = p; 1403 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1404 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; 1405 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; 1406 1407 ierr = VecDestroy(&local);CHKERRQ(ierr); 1408 ierr = VecDestroy(&global);CHKERRQ(ierr); 1409 1410 dd->gtol = gtol; 1411 dd->ltog = ltog; 1412 dd->idx = idx_cpy; 1413 dd->Nl = nn*dof; 1414 dd->base = base; 1415 da->ops->view = DMView_DA_3d; 1416 dd->ltol = PETSC_NULL; 1417 dd->ao = PETSC_NULL; 1418 1419 PetscFunctionReturn(0); 1420 } 1421 1422 1423 #undef __FUNCT__ 1424 #define __FUNCT__ "DMDACreate3d" 1425 /*@C 1426 DMDACreate3d - Creates an object that will manage the communication of three-dimensional 1427 regular array data that is distributed across some processors. 1428 1429 Collective on MPI_Comm 1430 1431 Input Parameters: 1432 + comm - MPI communicator 1433 . bx,by,bz - type of ghost nodes the array have. 1434 Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC. 1435 . stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX) 1436 . 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 1437 from the command line with -da_grid_x <M> -da_grid_y <N> -da_grid_z <P>) 1438 . m,n,p - corresponding number of processors in each dimension 1439 (or PETSC_DECIDE to have calculated) 1440 . dof - number of degrees of freedom per node 1441 . s - stencil width 1442 - lx, ly, lz - arrays containing the number of nodes in each cell along 1443 the x, y, and z coordinates, or PETSC_NULL. If non-null, these 1444 must be of length as m,n,p and the corresponding 1445 m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of 1446 the ly[] must N, sum of the lz[] must be P 1447 1448 Output Parameter: 1449 . da - the resulting distributed array object 1450 1451 Options Database Key: 1452 + -dm_view - Calls DMView() at the conclusion of DMDACreate3d() 1453 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 1454 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 1455 . -da_grid_z <nz> - number of grid points in z direction, if P < 0 1456 . -da_processors_x <MX> - number of processors in x direction 1457 . -da_processors_y <MY> - number of processors in y direction 1458 . -da_processors_z <MZ> - number of processors in z direction 1459 . -da_refine_x <rx> - refinement ratio in x direction 1460 . -da_refine_y <ry> - refinement ratio in y direction 1461 . -da_refine_z <rz>- refinement ratio in z directio 1462 - -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0 1463 1464 Level: beginner 1465 1466 Notes: 1467 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 1468 standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 1469 the standard 27-pt stencil. 1470 1471 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 1472 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 1473 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 1474 1475 .keywords: distributed array, create, three-dimensional 1476 1477 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 1478 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(), 1479 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 1480 1481 @*/ 1482 PetscErrorCode DMDACreate3d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz,DMDAStencilType stencil_type,PetscInt M, 1483 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) 1484 { 1485 PetscErrorCode ierr; 1486 1487 PetscFunctionBegin; 1488 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 1489 ierr = DMDASetDim(*da, 3);CHKERRQ(ierr); 1490 ierr = DMDASetSizes(*da, M, N, P);CHKERRQ(ierr); 1491 ierr = DMDASetNumProcs(*da, m, n, p);CHKERRQ(ierr); 1492 ierr = DMDASetBoundaryType(*da, bx, by, bz);CHKERRQ(ierr); 1493 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 1494 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 1495 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 1496 ierr = DMDASetOwnershipRanges(*da, lx, ly, lz);CHKERRQ(ierr); 1497 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 1498 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 1499 ierr = DMSetUp(*da);CHKERRQ(ierr); 1500 ierr = DMViewFromOptions(*da,"-dm_view");CHKERRQ(ierr); 1501 PetscFunctionReturn(0); 1502 } 1503