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