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