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