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