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/dmdaimpl.h> /*I "petscdmda.h" I*/ 8 9 #include <petscdraw.h> 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMView_DA_3d" 12 PetscErrorCode DMView_DA_3d(DM da,PetscViewer viewer) 13 { 14 PetscErrorCode ierr; 15 PetscMPIInt rank; 16 PetscBool iascii,isdraw,isbinary; 17 DM_DA *dd = (DM_DA*)da->data; 18 #if defined(PETSC_HAVE_MATLAB_ENGINE) 19 PetscBool ismatlab; 20 #endif 21 22 PetscFunctionBegin; 23 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 24 25 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 26 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 27 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 28 #if defined(PETSC_HAVE_MATLAB_ENGINE) 29 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);CHKERRQ(ierr); 30 #endif 31 if (iascii) { 32 PetscViewerFormat format; 33 34 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 35 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 36 if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) { 37 DMDALocalInfo info; 38 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 39 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); 40 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D, Z range of indices: %D %D\n", 41 info.xs,info.xs+info.xm,info.ys,info.ys+info.ym,info.zs,info.zs+info.zm);CHKERRQ(ierr); 42 #if !defined(PETSC_USE_COMPLEX) 43 if (da->coordinates) { 44 PetscInt last; 45 const PetscReal *coors; 46 ierr = VecGetArrayRead(da->coordinates,&coors);CHKERRQ(ierr); 47 ierr = VecGetLocalSize(da->coordinates,&last);CHKERRQ(ierr); 48 last = last - 3; 49 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); 50 ierr = VecRestoreArrayRead(da->coordinates,&coors);CHKERRQ(ierr); 51 } 52 #endif 53 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 54 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 55 } else { 56 ierr = DMView_DA_VTK(da,viewer);CHKERRQ(ierr); 57 } 58 } else if (isdraw) { 59 PetscDraw draw; 60 PetscReal ymin = -1.0,ymax = (PetscReal)dd->N; 61 PetscReal xmin = -1.0,xmax = (PetscReal)((dd->M+2)*dd->P),x,y,ycoord,xcoord; 62 PetscInt k,plane,base,*idx; 63 char node[10]; 64 PetscBool isnull; 65 66 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 67 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 68 ierr = PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);CHKERRQ(ierr); 69 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 70 71 /* first processor draw all node lines */ 72 if (!rank) { 73 for (k=0; k<dd->P; k++) { 74 ymin = 0.0; ymax = (PetscReal)(dd->N - 1); 75 for (xmin=(PetscReal)(k*(dd->M+1)); xmin<(PetscReal)(dd->M+(k*(dd->M+1))); xmin++) { 76 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 77 } 78 79 xmin = (PetscReal)(k*(dd->M+1)); xmax = xmin + (PetscReal)(dd->M - 1); 80 for (ymin=0; ymin<(PetscReal)dd->N; ymin++) { 81 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);CHKERRQ(ierr); 82 } 83 } 84 } 85 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 86 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 87 88 for (k=0; k<dd->P; k++) { /*Go through and draw for each plane*/ 89 if ((k >= dd->zs) && (k < dd->ze)) { 90 /* draw my box */ 91 ymin = dd->ys; 92 ymax = dd->ye - 1; 93 xmin = dd->xs/dd->w + (dd->M+1)*k; 94 xmax =(dd->xe-1)/dd->w + (dd->M+1)*k; 95 96 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);CHKERRQ(ierr); 97 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 98 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 99 ierr = PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 100 101 xmin = dd->xs/dd->w; 102 xmax =(dd->xe-1)/dd->w; 103 104 /* put in numbers*/ 105 base = (dd->base+(dd->xe-dd->xs)*(dd->ye-dd->ys)*(k-dd->zs))/dd->w; 106 107 /* Identify which processor owns the box */ 108 sprintf(node,"%d",rank); 109 ierr = PetscDrawString(draw,xmin+(dd->M+1)*k+.2,ymin+.3,PETSC_DRAW_RED,node);CHKERRQ(ierr); 110 111 for (y=ymin; y<=ymax; y++) { 112 for (x=xmin+(dd->M+1)*k; x<=xmax+(dd->M+1)*k; x++) { 113 sprintf(node,"%d",(int)base++); 114 ierr = PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);CHKERRQ(ierr); 115 } 116 } 117 118 } 119 } 120 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 121 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 122 123 for (k=0-dd->s; k<dd->P+dd->s; k++) { 124 /* Go through and draw for each plane */ 125 if ((k >= dd->Zs) && (k < dd->Ze)) { 126 127 /* overlay ghost numbers, useful for error checking */ 128 base = (dd->Xe-dd->Xs)*(dd->Ye-dd->Ys)*(k-dd->Zs); idx = dd->idx; 129 plane=k; 130 /* Keep z wrap around points on the dradrawg */ 131 if (k<0) plane=dd->P+k; 132 if (k>=dd->P) plane=k-dd->P; 133 ymin = dd->Ys; ymax = dd->Ye; 134 xmin = (dd->M+1)*plane*dd->w; 135 xmax = (dd->M+1)*plane*dd->w+dd->M*dd->w; 136 for (y=ymin; y<ymax; y++) { 137 for (x=xmin+dd->Xs; x<xmin+dd->Xe; x+=dd->w) { 138 sprintf(node,"%d",(int)(idx[base]/dd->w)); 139 ycoord = y; 140 /*Keep y wrap around points on drawing */ 141 if (y<0) ycoord = dd->N+y; 142 143 if (y>=dd->N) ycoord = y-dd->N; 144 xcoord = x; /* Keep x wrap points on drawing */ 145 146 if (x<xmin) xcoord = xmax - (xmin-x); 147 if (x>=xmax) xcoord = xmin + (x-xmax); 148 ierr = PetscDrawString(draw,xcoord/dd->w,ycoord,PETSC_DRAW_BLUE,node);CHKERRQ(ierr); 149 base+=dd->w; 150 } 151 } 152 } 153 } 154 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 155 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 156 } else if (isbinary) { 157 ierr = DMView_DA_Binary(da,viewer);CHKERRQ(ierr); 158 #if defined(PETSC_HAVE_MATLAB_ENGINE) 159 } else if (ismatlab) { 160 ierr = DMView_DA_Matlab(da,viewer);CHKERRQ(ierr); 161 #endif 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "DMSetUp_DA_3D" 168 PetscErrorCode DMSetUp_DA_3D(DM da) 169 { 170 DM_DA *dd = (DM_DA*)da->data; 171 const PetscInt M = dd->M; 172 const PetscInt N = dd->N; 173 const PetscInt P = dd->P; 174 PetscInt m = dd->m; 175 PetscInt n = dd->n; 176 PetscInt p = dd->p; 177 const PetscInt dof = dd->w; 178 const PetscInt s = dd->s; 179 DMDABoundaryType bx = dd->bx; 180 DMDABoundaryType by = dd->by; 181 DMDABoundaryType bz = dd->bz; 182 DMDAStencilType stencil_type = dd->stencil_type; 183 PetscInt *lx = dd->lx; 184 PetscInt *ly = dd->ly; 185 PetscInt *lz = dd->lz; 186 MPI_Comm comm; 187 PetscMPIInt rank,size; 188 PetscInt xs = 0,xe,ys = 0,ye,zs = 0,ze,x = 0,y = 0,z = 0; 189 PetscInt Xs,Xe,Ys,Ye,Zs,Ze,IXs,IXe,IYs,IYe,IZs,IZe,start,end,pm; 190 PetscInt left,right,up,down,bottom,top,i,j,k,*idx,*idx_cpy,nn; 191 const PetscInt *idx_full; 192 PetscInt n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n14; 193 PetscInt n15,n16,n17,n18,n19,n20,n21,n22,n23,n24,n25,n26; 194 PetscInt *bases,*ldims,base,x_t,y_t,z_t,s_t,count,s_x,s_y,s_z; 195 PetscInt sn0 = 0,sn1 = 0,sn2 = 0,sn3 = 0,sn5 = 0,sn6 = 0,sn7 = 0; 196 PetscInt sn8 = 0,sn9 = 0,sn11 = 0,sn15 = 0,sn24 = 0,sn25 = 0,sn26 = 0; 197 PetscInt sn17 = 0,sn18 = 0,sn19 = 0,sn20 = 0,sn21 = 0,sn23 = 0; 198 Vec local,global; 199 VecScatter ltog,gtol; 200 IS to,from,ltogis; 201 PetscBool twod; 202 PetscErrorCode ierr; 203 204 205 PetscFunctionBegin; 206 if (stencil_type == DMDA_STENCIL_BOX && (bx == DMDA_BOUNDARY_MIRROR || by == DMDA_BOUNDARY_MIRROR || bz == DMDA_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil"); 207 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 208 #if !defined(PETSC_USE_64BIT_INDICES) 209 if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) P)*((Petsc64bitInt) dof) > (Petsc64bitInt) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof); 210 #endif 211 212 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 213 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 214 215 if (m != PETSC_DECIDE) { 216 if (m < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m); 217 else if (m > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size); 218 } 219 if (n != PETSC_DECIDE) { 220 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n); 221 else if (n > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size); 222 } 223 if (p != PETSC_DECIDE) { 224 if (p < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Z direction: %D",p); 225 else if (p > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Z direction: %D %d",p,size); 226 } 227 if ((m > 0) && (n > 0) && (p > 0) && (m*n*p != size)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m %D * n %D * p %D != size %d",m,n,p,size); 228 229 /* Partition the array among the processors */ 230 if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) { 231 m = size/(n*p); 232 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { 233 n = size/(m*p); 234 } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { 235 p = size/(m*n); 236 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { 237 /* try for squarish distribution */ 238 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N*p))); 239 if (!m) m = 1; 240 while (m > 0) { 241 n = size/(m*p); 242 if (m*n*p == size) break; 243 m--; 244 } 245 if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad p value: p = %D",p); 246 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 247 } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { 248 /* try for squarish distribution */ 249 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n))); 250 if (!m) m = 1; 251 while (m > 0) { 252 p = size/(m*n); 253 if (m*n*p == size) break; 254 m--; 255 } 256 if (!m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad n value: n = %D",n); 257 if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;} 258 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { 259 /* try for squarish distribution */ 260 n = (int)(0.5 + PetscSqrtReal(((PetscReal)N)*((PetscReal)size)/((PetscReal)P*m))); 261 if (!n) n = 1; 262 while (n > 0) { 263 p = size/(m*n); 264 if (m*n*p == size) break; 265 n--; 266 } 267 if (!n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"bad m value: m = %D",n); 268 if (N > P && n < p) {PetscInt _n = n; n = p; p = _n;} 269 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { 270 /* try for squarish distribution */ 271 n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); 272 if (!n) n = 1; 273 while (n > 0) { 274 pm = size/n; 275 if (n*pm == size) break; 276 n--; 277 } 278 if (!n) n = 1; 279 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n))); 280 if (!m) m = 1; 281 while (m > 0) { 282 p = size/(m*n); 283 if (m*n*p == size) break; 284 m--; 285 } 286 if (M > P && m < p) {PetscInt _m = m; m = p; p = _m;} 287 } else if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition"); 288 289 if (m*n*p != size) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_PLIB,"Could not find good partition"); 290 if (M < m) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m); 291 if (N < n) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n); 292 if (P < p) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Partition in z direction is too fine! %D %D",P,p); 293 294 /* 295 Determine locally owned region 296 [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes 297 */ 298 299 if (!lx) { 300 ierr = PetscMalloc(m*sizeof(PetscInt), &dd->lx);CHKERRQ(ierr); 301 lx = dd->lx; 302 for (i=0; i<m; i++) lx[i] = M/m + ((M % m) > (i % m)); 303 } 304 x = lx[rank % m]; 305 xs = 0; 306 for (i=0; i<(rank%m); i++) xs += lx[i]; 307 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); 308 309 if (!ly) { 310 ierr = PetscMalloc(n*sizeof(PetscInt), &dd->ly);CHKERRQ(ierr); 311 ly = dd->ly; 312 for (i=0; i<n; i++) ly[i] = N/n + ((N % n) > (i % n)); 313 } 314 y = ly[(rank % (m*n))/m]; 315 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); 316 317 ys = 0; 318 for (i=0; i<(rank % (m*n))/m; i++) ys += ly[i]; 319 320 if (!lz) { 321 ierr = PetscMalloc(p*sizeof(PetscInt), &dd->lz);CHKERRQ(ierr); 322 lz = dd->lz; 323 for (i=0; i<p; i++) lz[i] = P/p + ((P % p) > (i % p)); 324 } 325 z = lz[rank/(m*n)]; 326 327 /* note this is different than x- and y-, as we will handle as an important special 328 case when p=P=1 and DMDA_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems 329 in a 3D code. Additional code for this case is noted with "2d case" comments */ 330 twod = PETSC_FALSE; 331 if (P == 1) twod = PETSC_TRUE; 332 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); 333 zs = 0; 334 for (i=0; i<(rank/(m*n)); i++) zs += lz[i]; 335 ye = ys + y; 336 xe = xs + x; 337 ze = zs + z; 338 339 /* determine ghost region (Xs) and region scattered into (IXs) */ 340 if (xs-s > 0) { 341 Xs = xs - s; IXs = xs - s; 342 } else { 343 if (bx) Xs = xs - s; 344 else Xs = 0; 345 IXs = 0; 346 } 347 if (xe+s <= M) { 348 Xe = xe + s; IXe = xe + s; 349 } else { 350 if (bx) { 351 Xs = xs - s; Xe = xe + s; 352 } else Xe = M; 353 IXe = M; 354 } 355 356 if (bx == DMDA_BOUNDARY_PERIODIC || bx == DMDA_BOUNDARY_MIRROR) { 357 IXs = xs - s; 358 IXe = xe + s; 359 Xs = xs - s; 360 Xe = xe + s; 361 } 362 363 if (ys-s > 0) { 364 Ys = ys - s; IYs = ys - s; 365 } else { 366 if (by) Ys = ys - s; 367 else Ys = 0; 368 IYs = 0; 369 } 370 if (ye+s <= N) { 371 Ye = ye + s; IYe = ye + s; 372 } else { 373 if (by) Ye = ye + s; 374 else Ye = N; 375 IYe = N; 376 } 377 378 if (by == DMDA_BOUNDARY_PERIODIC || by == DMDA_BOUNDARY_MIRROR) { 379 IYs = ys - s; 380 IYe = ye + s; 381 Ys = ys - s; 382 Ye = ye + s; 383 } 384 385 if (zs-s > 0) { 386 Zs = zs - s; IZs = zs - s; 387 } else { 388 if (bz) Zs = zs - s; 389 else Zs = 0; 390 IZs = 0; 391 } 392 if (ze+s <= P) { 393 Ze = ze + s; IZe = ze + s; 394 } else { 395 if (bz) Ze = ze + s; 396 else Ze = P; 397 IZe = P; 398 } 399 400 if (bz == DMDA_BOUNDARY_PERIODIC || bz == DMDA_BOUNDARY_MIRROR) { 401 IZs = zs - s; 402 IZe = ze + s; 403 Zs = zs - s; 404 Ze = ze + s; 405 } 406 407 /* Resize all X parameters to reflect w */ 408 s_x = s; 409 s_y = s; 410 s_z = s; 411 412 /* determine starting point of each processor */ 413 nn = x*y*z; 414 ierr = PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);CHKERRQ(ierr); 415 ierr = MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);CHKERRQ(ierr); 416 bases[0] = 0; 417 for (i=1; i<=size; i++) bases[i] = ldims[i-1]; 418 for (i=1; i<=size; i++) bases[i] += bases[i-1]; 419 base = bases[rank]*dof; 420 421 /* allocate the base parallel and sequential vectors */ 422 dd->Nlocal = x*y*z*dof; 423 ierr = VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr); 424 dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs)*dof; 425 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);CHKERRQ(ierr); 426 427 /* generate appropriate vector scatters */ 428 /* local to global inserts non-ghost point region into global */ 429 ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr); 430 ierr = ISCreateStride(comm,x*y*z*dof,start,1,&to);CHKERRQ(ierr); 431 432 ierr = PetscMalloc(x*y*z*sizeof(PetscInt),&idx);CHKERRQ(ierr); 433 left = xs - Xs; right = left + x; 434 bottom = ys - Ys; top = bottom + y; 435 down = zs - Zs; up = down + z; 436 count = 0; 437 for (i=down; i<up; i++) { 438 for (j=bottom; j<top; j++) { 439 for (k=left; k<right; k++) { 440 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 441 } 442 } 443 } 444 445 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 446 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 447 ierr = PetscLogObjectParent(da,ltog);CHKERRQ(ierr); 448 ierr = ISDestroy(&from);CHKERRQ(ierr); 449 ierr = ISDestroy(&to);CHKERRQ(ierr); 450 451 /* global to local must include ghost points within the domain, 452 but not ghost points outside the domain that aren't periodic */ 453 if (stencil_type == DMDA_STENCIL_BOX) { 454 count = (IXe-IXs)*(IYe-IYs)*(IZe-IZs); 455 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 456 457 left = IXs - Xs; right = left + (IXe-IXs); 458 bottom = IYs - Ys; top = bottom + (IYe-IYs); 459 down = IZs - Zs; up = down + (IZe-IZs); 460 count = 0; 461 for (i=down; i<up; i++) { 462 for (j=bottom; j<top; j++) { 463 for (k=left; k<right; k++) { 464 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 465 } 466 } 467 } 468 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 469 470 } else { 471 /* This is way ugly! We need to list the funny cross type region */ 472 count = ((ys-IYs) + (IYe-ye))*x*z + ((xs-IXs) + (IXe-xe))*y*z + ((zs-IZs) + (IZe-ze))*x*y + x*y*z; 473 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 474 475 left = xs - Xs; right = left + x; 476 bottom = ys - Ys; top = bottom + y; 477 down = zs - Zs; up = down + z; 478 count = 0; 479 /* the bottom chunck */ 480 for (i=(IZs-Zs); i<down; i++) { 481 for (j=bottom; j<top; j++) { 482 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 483 } 484 } 485 /* the middle piece */ 486 for (i=down; i<up; i++) { 487 /* front */ 488 for (j=(IYs-Ys); j<bottom; j++) { 489 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 490 } 491 /* middle */ 492 for (j=bottom; j<top; j++) { 493 for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 494 } 495 /* back */ 496 for (j=top; j<top+IYe-ye; j++) { 497 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 498 } 499 } 500 /* the top piece */ 501 for (i=up; i<up+IZe-ze; i++) { 502 for (j=bottom; j<top; j++) { 503 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 504 } 505 } 506 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 507 } 508 509 /* determine who lies on each side of use stored in n24 n25 n26 510 n21 n22 n23 511 n18 n19 n20 512 513 n15 n16 n17 514 n12 n14 515 n9 n10 n11 516 517 n6 n7 n8 518 n3 n4 n5 519 n0 n1 n2 520 */ 521 522 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ 523 /* Assume Nodes are Internal to the Cube */ 524 n0 = rank - m*n - m - 1; 525 n1 = rank - m*n - m; 526 n2 = rank - m*n - m + 1; 527 n3 = rank - m*n -1; 528 n4 = rank - m*n; 529 n5 = rank - m*n + 1; 530 n6 = rank - m*n + m - 1; 531 n7 = rank - m*n + m; 532 n8 = rank - m*n + m + 1; 533 534 n9 = rank - m - 1; 535 n10 = rank - m; 536 n11 = rank - m + 1; 537 n12 = rank - 1; 538 n14 = rank + 1; 539 n15 = rank + m - 1; 540 n16 = rank + m; 541 n17 = rank + m + 1; 542 543 n18 = rank + m*n - m - 1; 544 n19 = rank + m*n - m; 545 n20 = rank + m*n - m + 1; 546 n21 = rank + m*n - 1; 547 n22 = rank + m*n; 548 n23 = rank + m*n + 1; 549 n24 = rank + m*n + m - 1; 550 n25 = rank + m*n + m; 551 n26 = rank + m*n + m + 1; 552 553 /* Assume Pieces are on Faces of Cube */ 554 555 if (xs == 0) { /* First assume not corner or edge */ 556 n0 = rank -1 - (m*n); 557 n3 = rank + m -1 - (m*n); 558 n6 = rank + 2*m -1 - (m*n); 559 n9 = rank -1; 560 n12 = rank + m -1; 561 n15 = rank + 2*m -1; 562 n18 = rank -1 + (m*n); 563 n21 = rank + m -1 + (m*n); 564 n24 = rank + 2*m -1 + (m*n); 565 } 566 567 if (xe == M) { /* First assume not corner or edge */ 568 n2 = rank -2*m +1 - (m*n); 569 n5 = rank - m +1 - (m*n); 570 n8 = rank +1 - (m*n); 571 n11 = rank -2*m +1; 572 n14 = rank - m +1; 573 n17 = rank +1; 574 n20 = rank -2*m +1 + (m*n); 575 n23 = rank - m +1 + (m*n); 576 n26 = rank +1 + (m*n); 577 } 578 579 if (ys==0) { /* First assume not corner or edge */ 580 n0 = rank + m * (n-1) -1 - (m*n); 581 n1 = rank + m * (n-1) - (m*n); 582 n2 = rank + m * (n-1) +1 - (m*n); 583 n9 = rank + m * (n-1) -1; 584 n10 = rank + m * (n-1); 585 n11 = rank + m * (n-1) +1; 586 n18 = rank + m * (n-1) -1 + (m*n); 587 n19 = rank + m * (n-1) + (m*n); 588 n20 = rank + m * (n-1) +1 + (m*n); 589 } 590 591 if (ye == N) { /* First assume not corner or edge */ 592 n6 = rank - m * (n-1) -1 - (m*n); 593 n7 = rank - m * (n-1) - (m*n); 594 n8 = rank - m * (n-1) +1 - (m*n); 595 n15 = rank - m * (n-1) -1; 596 n16 = rank - m * (n-1); 597 n17 = rank - m * (n-1) +1; 598 n24 = rank - m * (n-1) -1 + (m*n); 599 n25 = rank - m * (n-1) + (m*n); 600 n26 = rank - m * (n-1) +1 + (m*n); 601 } 602 603 if (zs == 0) { /* First assume not corner or edge */ 604 n0 = size - (m*n) + rank - m - 1; 605 n1 = size - (m*n) + rank - m; 606 n2 = size - (m*n) + rank - m + 1; 607 n3 = size - (m*n) + rank - 1; 608 n4 = size - (m*n) + rank; 609 n5 = size - (m*n) + rank + 1; 610 n6 = size - (m*n) + rank + m - 1; 611 n7 = size - (m*n) + rank + m; 612 n8 = size - (m*n) + rank + m + 1; 613 } 614 615 if (ze == P) { /* First assume not corner or edge */ 616 n18 = (m*n) - (size-rank) - m - 1; 617 n19 = (m*n) - (size-rank) - m; 618 n20 = (m*n) - (size-rank) - m + 1; 619 n21 = (m*n) - (size-rank) - 1; 620 n22 = (m*n) - (size-rank); 621 n23 = (m*n) - (size-rank) + 1; 622 n24 = (m*n) - (size-rank) + m - 1; 623 n25 = (m*n) - (size-rank) + m; 624 n26 = (m*n) - (size-rank) + m + 1; 625 } 626 627 if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */ 628 n0 = size - m*n + rank + m-1 - m; 629 n3 = size - m*n + rank + m-1; 630 n6 = size - m*n + rank + m-1 + m; 631 } 632 633 if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */ 634 n18 = m*n - (size - rank) + m-1 - m; 635 n21 = m*n - (size - rank) + m-1; 636 n24 = m*n - (size - rank) + m-1 + m; 637 } 638 639 if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */ 640 n0 = rank + m*n -1 - m*n; 641 n9 = rank + m*n -1; 642 n18 = rank + m*n -1 + m*n; 643 } 644 645 if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */ 646 n6 = rank - m*(n-1) + m-1 - m*n; 647 n15 = rank - m*(n-1) + m-1; 648 n24 = rank - m*(n-1) + m-1 + m*n; 649 } 650 651 if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */ 652 n2 = size - (m*n-rank) - (m-1) - m; 653 n5 = size - (m*n-rank) - (m-1); 654 n8 = size - (m*n-rank) - (m-1) + m; 655 } 656 657 if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */ 658 n20 = m*n - (size - rank) - (m-1) - m; 659 n23 = m*n - (size - rank) - (m-1); 660 n26 = m*n - (size - rank) - (m-1) + m; 661 } 662 663 if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */ 664 n2 = rank + m*(n-1) - (m-1) - m*n; 665 n11 = rank + m*(n-1) - (m-1); 666 n20 = rank + m*(n-1) - (m-1) + m*n; 667 } 668 669 if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */ 670 n8 = rank - m*n +1 - m*n; 671 n17 = rank - m*n +1; 672 n26 = rank - m*n +1 + m*n; 673 } 674 675 if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */ 676 n0 = size - m + rank -1; 677 n1 = size - m + rank; 678 n2 = size - m + rank +1; 679 } 680 681 if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */ 682 n18 = m*n - (size - rank) + m*(n-1) -1; 683 n19 = m*n - (size - rank) + m*(n-1); 684 n20 = m*n - (size - rank) + m*(n-1) +1; 685 } 686 687 if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */ 688 n6 = size - (m*n-rank) - m * (n-1) -1; 689 n7 = size - (m*n-rank) - m * (n-1); 690 n8 = size - (m*n-rank) - m * (n-1) +1; 691 } 692 693 if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */ 694 n24 = rank - (size-m) -1; 695 n25 = rank - (size-m); 696 n26 = rank - (size-m) +1; 697 } 698 699 /* Check for Corners */ 700 if ((xs==0) && (ys==0) && (zs==0)) n0 = size -1; 701 if ((xs==0) && (ys==0) && (ze==P)) n18 = m*n-1; 702 if ((xs==0) && (ye==N) && (zs==0)) n6 = (size-1)-m*(n-1); 703 if ((xs==0) && (ye==N) && (ze==P)) n24 = m-1; 704 if ((xe==M) && (ys==0) && (zs==0)) n2 = size-m; 705 if ((xe==M) && (ys==0) && (ze==P)) n20 = m*n-m; 706 if ((xe==M) && (ye==N) && (zs==0)) n8 = size-m*n; 707 if ((xe==M) && (ye==N) && (ze==P)) n26 = 0; 708 709 /* Check for when not X,Y, and Z Periodic */ 710 711 /* If not X periodic */ 712 if (bx != DMDA_BOUNDARY_PERIODIC) { 713 if (xs==0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2; 714 if (xe==M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2; 715 } 716 717 /* If not Y periodic */ 718 if (by != DMDA_BOUNDARY_PERIODIC) { 719 if (ys==0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2; 720 if (ye==N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2; 721 } 722 723 /* If not Z periodic */ 724 if (bz != DMDA_BOUNDARY_PERIODIC) { 725 if (zs==0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2; 726 if (ze==P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2; 727 } 728 729 ierr = PetscMalloc(27*sizeof(PetscInt),&dd->neighbors);CHKERRQ(ierr); 730 731 dd->neighbors[0] = n0; 732 dd->neighbors[1] = n1; 733 dd->neighbors[2] = n2; 734 dd->neighbors[3] = n3; 735 dd->neighbors[4] = n4; 736 dd->neighbors[5] = n5; 737 dd->neighbors[6] = n6; 738 dd->neighbors[7] = n7; 739 dd->neighbors[8] = n8; 740 dd->neighbors[9] = n9; 741 dd->neighbors[10] = n10; 742 dd->neighbors[11] = n11; 743 dd->neighbors[12] = n12; 744 dd->neighbors[13] = rank; 745 dd->neighbors[14] = n14; 746 dd->neighbors[15] = n15; 747 dd->neighbors[16] = n16; 748 dd->neighbors[17] = n17; 749 dd->neighbors[18] = n18; 750 dd->neighbors[19] = n19; 751 dd->neighbors[20] = n20; 752 dd->neighbors[21] = n21; 753 dd->neighbors[22] = n22; 754 dd->neighbors[23] = n23; 755 dd->neighbors[24] = n24; 756 dd->neighbors[25] = n25; 757 dd->neighbors[26] = n26; 758 759 /* If star stencil then delete the corner neighbors */ 760 if (stencil_type == DMDA_STENCIL_STAR) { 761 /* save information about corner neighbors */ 762 sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7; 763 sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18; 764 sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25; 765 sn26 = n26; 766 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; 767 } 768 769 ierr = PetscMalloc((Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt),&idx);CHKERRQ(ierr); 770 ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt));CHKERRQ(ierr); 771 772 nn = 0; 773 /* Bottom Level */ 774 for (k=0; k<s_z; k++) { 775 for (i=1; i<=s_y; i++) { 776 if (n0 >= 0) { /* left below */ 777 x_t = lx[n0 % m]; 778 y_t = ly[(n0 % (m*n))/m]; 779 z_t = lz[n0 / (m*n)]; 780 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; 781 if (twod && (s_t < 0)) s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x; /* 2D case */ 782 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 783 } 784 if (n1 >= 0) { /* directly below */ 785 x_t = x; 786 y_t = ly[(n1 % (m*n))/m]; 787 z_t = lz[n1 / (m*n)]; 788 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 789 if (twod && (s_t < 0)) s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 790 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 791 } 792 if (n2 >= 0) { /* right below */ 793 x_t = lx[n2 % m]; 794 y_t = ly[(n2 % (m*n))/m]; 795 z_t = lz[n2 / (m*n)]; 796 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 797 if (twod && (s_t < 0)) s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 798 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 799 } 800 } 801 802 for (i=0; i<y; i++) { 803 if (n3 >= 0) { /* directly left */ 804 x_t = lx[n3 % m]; 805 y_t = y; 806 z_t = lz[n3 / (m*n)]; 807 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 808 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 */ 809 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 810 } 811 812 if (n4 >= 0) { /* middle */ 813 x_t = x; 814 y_t = y; 815 z_t = lz[n4 / (m*n)]; 816 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 817 if (twod && (s_t < 0)) s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 818 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 819 } else if (bz == DMDA_BOUNDARY_MIRROR) { 820 for (j=0; j<x; j++) idx[nn++] = 0; 821 } 822 823 if (n5 >= 0) { /* directly right */ 824 x_t = lx[n5 % m]; 825 y_t = y; 826 z_t = lz[n5 / (m*n)]; 827 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 828 if (twod && (s_t < 0)) s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 829 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 830 } 831 } 832 833 for (i=1; i<=s_y; i++) { 834 if (n6 >= 0) { /* left above */ 835 x_t = lx[n6 % m]; 836 y_t = ly[(n6 % (m*n))/m]; 837 z_t = lz[n6 / (m*n)]; 838 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 839 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 */ 840 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 841 } 842 if (n7 >= 0) { /* directly above */ 843 x_t = x; 844 y_t = ly[(n7 % (m*n))/m]; 845 z_t = lz[n7 / (m*n)]; 846 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 847 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 */ 848 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 849 } 850 if (n8 >= 0) { /* right above */ 851 x_t = lx[n8 % m]; 852 y_t = ly[(n8 % (m*n))/m]; 853 z_t = lz[n8 / (m*n)]; 854 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 855 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 */ 856 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 857 } 858 } 859 } 860 861 /* Middle Level */ 862 for (k=0; k<z; k++) { 863 for (i=1; i<=s_y; i++) { 864 if (n9 >= 0) { /* left below */ 865 x_t = lx[n9 % m]; 866 y_t = ly[(n9 % (m*n))/m]; 867 /* z_t = z; */ 868 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 869 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 870 } 871 if (n10 >= 0) { /* directly below */ 872 x_t = x; 873 y_t = ly[(n10 % (m*n))/m]; 874 /* z_t = z; */ 875 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 876 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 877 } else if (by == DMDA_BOUNDARY_MIRROR) { 878 for (j=0; j<x; j++) idx[nn++] = 0; 879 } 880 if (n11 >= 0) { /* right below */ 881 x_t = lx[n11 % m]; 882 y_t = ly[(n11 % (m*n))/m]; 883 /* z_t = z; */ 884 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 885 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 886 } 887 } 888 889 for (i=0; i<y; i++) { 890 if (n12 >= 0) { /* directly left */ 891 x_t = lx[n12 % m]; 892 y_t = y; 893 /* z_t = z; */ 894 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 895 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 896 } else if (bx == DMDA_BOUNDARY_MIRROR) { 897 for (j=0; j<s_x; j++) idx[nn++] = 0; 898 } 899 900 /* Interior */ 901 s_t = bases[rank] + i*x + k*x*y; 902 for (j=0; j<x; j++) idx[nn++] = s_t++; 903 904 if (n14 >= 0) { /* directly right */ 905 x_t = lx[n14 % m]; 906 y_t = y; 907 /* z_t = z; */ 908 s_t = bases[n14] + i*x_t + k*x_t*y_t; 909 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 910 } else if (bx == DMDA_BOUNDARY_MIRROR) { 911 for (j=0; j<s_x; j++) idx[nn++] = 0; 912 } 913 } 914 915 for (i=1; i<=s_y; i++) { 916 if (n15 >= 0) { /* left above */ 917 x_t = lx[n15 % m]; 918 y_t = ly[(n15 % (m*n))/m]; 919 /* z_t = z; */ 920 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 921 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 922 } 923 if (n16 >= 0) { /* directly above */ 924 x_t = x; 925 y_t = ly[(n16 % (m*n))/m]; 926 /* z_t = z; */ 927 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 928 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 929 } else if (by == DMDA_BOUNDARY_MIRROR) { 930 for (j=0; j<x; j++) idx[nn++] = 0; 931 } 932 if (n17 >= 0) { /* right above */ 933 x_t = lx[n17 % m]; 934 y_t = ly[(n17 % (m*n))/m]; 935 /* z_t = z; */ 936 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 937 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 938 } 939 } 940 } 941 942 /* Upper Level */ 943 for (k=0; k<s_z; k++) { 944 for (i=1; i<=s_y; i++) { 945 if (n18 >= 0) { /* left below */ 946 x_t = lx[n18 % m]; 947 y_t = ly[(n18 % (m*n))/m]; 948 /* z_t = lz[n18 / (m*n)]; */ 949 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 950 if (twod && (s_t >= M*N*P)) s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t; /* 2d case */ 951 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 952 } 953 if (n19 >= 0) { /* directly below */ 954 x_t = x; 955 y_t = ly[(n19 % (m*n))/m]; 956 /* z_t = lz[n19 / (m*n)]; */ 957 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 958 if (twod && (s_t >= M*N*P)) s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 959 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 960 } 961 if (n20 >= 0) { /* right below */ 962 x_t = lx[n20 % m]; 963 y_t = ly[(n20 % (m*n))/m]; 964 /* z_t = lz[n20 / (m*n)]; */ 965 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 966 if (twod && (s_t >= M*N*P)) s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 967 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 968 } 969 } 970 971 for (i=0; i<y; i++) { 972 if (n21 >= 0) { /* directly left */ 973 x_t = lx[n21 % m]; 974 y_t = y; 975 /* z_t = lz[n21 / (m*n)]; */ 976 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 977 if (twod && (s_t >= M*N*P)) s_t = bases[n21] + (i+1)*x_t - s_x; /* 2d case */ 978 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 979 } 980 981 if (n22 >= 0) { /* middle */ 982 x_t = x; 983 y_t = y; 984 /* z_t = lz[n22 / (m*n)]; */ 985 s_t = bases[n22] + i*x_t + k*x_t*y_t; 986 if (twod && (s_t >= M*N*P)) s_t = bases[n22] + i*x_t; /* 2d case */ 987 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 988 } else if (bz == DMDA_BOUNDARY_MIRROR) { 989 for (j=0; j<x; j++) idx[nn++] = 0; 990 } 991 992 if (n23 >= 0) { /* directly right */ 993 x_t = lx[n23 % m]; 994 y_t = y; 995 /* z_t = lz[n23 / (m*n)]; */ 996 s_t = bases[n23] + i*x_t + k*x_t*y_t; 997 if (twod && (s_t >= M*N*P)) s_t = bases[n23] + i*x_t; /* 2d case */ 998 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 999 } 1000 } 1001 1002 for (i=1; i<=s_y; i++) { 1003 if (n24 >= 0) { /* left above */ 1004 x_t = lx[n24 % m]; 1005 y_t = ly[(n24 % (m*n))/m]; 1006 /* z_t = lz[n24 / (m*n)]; */ 1007 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1008 if (twod && (s_t >= M*N*P)) s_t = bases[n24] + i*x_t - s_x; /* 2d case */ 1009 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1010 } 1011 if (n25 >= 0) { /* directly above */ 1012 x_t = x; 1013 y_t = ly[(n25 % (m*n))/m]; 1014 /* z_t = lz[n25 / (m*n)]; */ 1015 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1016 if (twod && (s_t >= M*N*P)) s_t = bases[n25] + (i-1)*x_t; /* 2d case */ 1017 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1018 } 1019 if (n26 >= 0) { /* right above */ 1020 x_t = lx[n26 % m]; 1021 y_t = ly[(n26 % (m*n))/m]; 1022 /* z_t = lz[n26 / (m*n)]; */ 1023 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1024 if (twod && (s_t >= M*N*P)) s_t = bases[n26] + (i-1)*x_t; /* 2d case */ 1025 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1026 } 1027 } 1028 } 1029 1030 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 1031 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1032 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 1033 ierr = ISDestroy(&to);CHKERRQ(ierr); 1034 ierr = ISDestroy(&from);CHKERRQ(ierr); 1035 1036 if (stencil_type == DMDA_STENCIL_STAR) { 1037 n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; 1038 n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; 1039 n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; 1040 n26 = sn26; 1041 } 1042 1043 if (((stencil_type == DMDA_STENCIL_STAR) || 1044 (bx != DMDA_BOUNDARY_PERIODIC && bx) || 1045 (by != DMDA_BOUNDARY_PERIODIC && by) || 1046 (bz != DMDA_BOUNDARY_PERIODIC && bz))) { 1047 /* 1048 Recompute the local to global mappings, this time keeping the 1049 information about the cross corner processor numbers. 1050 */ 1051 nn = 0; 1052 /* Bottom Level */ 1053 for (k=0; k<s_z; k++) { 1054 for (i=1; i<=s_y; i++) { 1055 if (n0 >= 0) { /* left below */ 1056 x_t = lx[n0 % m]; 1057 y_t = ly[(n0 % (m*n))/m]; 1058 z_t = lz[n0 / (m*n)]; 1059 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; 1060 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1061 } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) { 1062 for (j=0; j<s_x; j++) idx[nn++] = -1; 1063 } 1064 if (n1 >= 0) { /* directly below */ 1065 x_t = x; 1066 y_t = ly[(n1 % (m*n))/m]; 1067 z_t = lz[n1 / (m*n)]; 1068 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1069 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1070 } else if (Ys-ys < 0 && Zs-zs < 0) { 1071 for (j=0; j<x; j++) idx[nn++] = -1; 1072 } 1073 if (n2 >= 0) { /* right below */ 1074 x_t = lx[n2 % m]; 1075 y_t = ly[(n2 % (m*n))/m]; 1076 z_t = lz[n2 / (m*n)]; 1077 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1078 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1079 } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) { 1080 for (j=0; j<s_x; j++) idx[nn++] = -1; 1081 } 1082 } 1083 1084 for (i=0; i<y; i++) { 1085 if (n3 >= 0) { /* directly left */ 1086 x_t = lx[n3 % m]; 1087 y_t = y; 1088 z_t = lz[n3 / (m*n)]; 1089 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1090 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1091 } else if (Xs-xs < 0 && Zs-zs < 0) { 1092 for (j=0; j<s_x; j++) idx[nn++] = -1; 1093 } 1094 1095 if (n4 >= 0) { /* middle */ 1096 x_t = x; 1097 y_t = y; 1098 z_t = lz[n4 / (m*n)]; 1099 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1100 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1101 } else if (Zs-zs < 0) { 1102 if (bz == DMDA_BOUNDARY_MIRROR) { 1103 for (j=0; j<x; j++) idx[nn++] = 0; 1104 } else { 1105 for (j=0; j<x; j++) idx[nn++] = -1; 1106 } 1107 } 1108 1109 if (n5 >= 0) { /* directly right */ 1110 x_t = lx[n5 % m]; 1111 y_t = y; 1112 z_t = lz[n5 / (m*n)]; 1113 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1114 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1115 } else if (xe-Xe < 0 && Zs-zs < 0) { 1116 for (j=0; j<s_x; j++) idx[nn++] = -1; 1117 } 1118 } 1119 1120 for (i=1; i<=s_y; i++) { 1121 if (n6 >= 0) { /* left above */ 1122 x_t = lx[n6 % m]; 1123 y_t = ly[(n6 % (m*n))/m]; 1124 z_t = lz[n6 / (m*n)]; 1125 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1126 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1127 } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) { 1128 for (j=0; j<s_x; j++) idx[nn++] = -1; 1129 } 1130 if (n7 >= 0) { /* directly above */ 1131 x_t = x; 1132 y_t = ly[(n7 % (m*n))/m]; 1133 z_t = lz[n7 / (m*n)]; 1134 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1135 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1136 } else if (ye-Ye < 0 && Zs-zs < 0) { 1137 for (j=0; j<x; j++) idx[nn++] = -1; 1138 } 1139 if (n8 >= 0) { /* right above */ 1140 x_t = lx[n8 % m]; 1141 y_t = ly[(n8 % (m*n))/m]; 1142 z_t = lz[n8 / (m*n)]; 1143 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1144 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1145 } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) { 1146 for (j=0; j<s_x; j++) idx[nn++] = -1; 1147 } 1148 } 1149 } 1150 1151 /* Middle Level */ 1152 for (k=0; k<z; k++) { 1153 for (i=1; i<=s_y; i++) { 1154 if (n9 >= 0) { /* left below */ 1155 x_t = lx[n9 % m]; 1156 y_t = ly[(n9 % (m*n))/m]; 1157 /* z_t = z; */ 1158 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1159 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1160 } else if (Xs-xs < 0 && Ys-ys < 0) { 1161 for (j=0; j<s_x; j++) idx[nn++] = -1; 1162 } 1163 if (n10 >= 0) { /* directly below */ 1164 x_t = x; 1165 y_t = ly[(n10 % (m*n))/m]; 1166 /* z_t = z; */ 1167 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1168 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1169 } else if (Ys-ys < 0) { 1170 if (by == DMDA_BOUNDARY_MIRROR) { 1171 for (j=0; j<x; j++) idx[nn++] = -1; 1172 } else { 1173 for (j=0; j<x; j++) idx[nn++] = -1; 1174 } 1175 } 1176 if (n11 >= 0) { /* right below */ 1177 x_t = lx[n11 % m]; 1178 y_t = ly[(n11 % (m*n))/m]; 1179 /* z_t = z; */ 1180 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1181 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1182 } else if (xe-Xe < 0 && Ys-ys < 0) { 1183 for (j=0; j<s_x; j++) idx[nn++] = -1; 1184 } 1185 } 1186 1187 for (i=0; i<y; i++) { 1188 if (n12 >= 0) { /* directly left */ 1189 x_t = lx[n12 % m]; 1190 y_t = y; 1191 /* z_t = z; */ 1192 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 1193 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1194 } else if (Xs-xs < 0) { 1195 if (bx == DMDA_BOUNDARY_MIRROR) { 1196 for (j=0; j<s_x; j++) idx[nn++] = 0; 1197 } else { 1198 for (j=0; j<s_x; j++) idx[nn++] = -1; 1199 } 1200 } 1201 1202 /* Interior */ 1203 s_t = bases[rank] + i*x + k*x*y; 1204 for (j=0; j<x; j++) idx[nn++] = s_t++; 1205 1206 if (n14 >= 0) { /* directly right */ 1207 x_t = lx[n14 % m]; 1208 y_t = y; 1209 /* z_t = z; */ 1210 s_t = bases[n14] + i*x_t + k*x_t*y_t; 1211 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1212 } else if (xe-Xe < 0) { 1213 if (bx == DMDA_BOUNDARY_MIRROR) { 1214 for (j=0; j<s_x; j++) idx[nn++] = 0; 1215 } else { 1216 for (j=0; j<s_x; j++) idx[nn++] = -1; 1217 } 1218 } 1219 } 1220 1221 for (i=1; i<=s_y; i++) { 1222 if (n15 >= 0) { /* left above */ 1223 x_t = lx[n15 % m]; 1224 y_t = ly[(n15 % (m*n))/m]; 1225 /* z_t = z; */ 1226 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 1227 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1228 } else if (Xs-xs < 0 && ye-Ye < 0) { 1229 for (j=0; j<s_x; j++) idx[nn++] = -1; 1230 } 1231 if (n16 >= 0) { /* directly above */ 1232 x_t = x; 1233 y_t = ly[(n16 % (m*n))/m]; 1234 /* z_t = z; */ 1235 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 1236 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1237 } else if (ye-Ye < 0) { 1238 if (by == DMDA_BOUNDARY_MIRROR) { 1239 for (j=0; j<x; j++) idx[nn++] = 0; 1240 } else { 1241 for (j=0; j<x; j++) idx[nn++] = -1; 1242 } 1243 } 1244 if (n17 >= 0) { /* right above */ 1245 x_t = lx[n17 % m]; 1246 y_t = ly[(n17 % (m*n))/m]; 1247 /* z_t = z; */ 1248 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 1249 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1250 } else if (xe-Xe < 0 && ye-Ye < 0) { 1251 for (j=0; j<s_x; j++) idx[nn++] = -1; 1252 } 1253 } 1254 } 1255 1256 /* Upper Level */ 1257 for (k=0; k<s_z; k++) { 1258 for (i=1; i<=s_y; i++) { 1259 if (n18 >= 0) { /* left below */ 1260 x_t = lx[n18 % m]; 1261 y_t = ly[(n18 % (m*n))/m]; 1262 /* z_t = lz[n18 / (m*n)]; */ 1263 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1264 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1265 } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) { 1266 for (j=0; j<s_x; j++) idx[nn++] = -1; 1267 } 1268 if (n19 >= 0) { /* directly below */ 1269 x_t = x; 1270 y_t = ly[(n19 % (m*n))/m]; 1271 /* z_t = lz[n19 / (m*n)]; */ 1272 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1273 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1274 } else if (Ys-ys < 0 && ze-Ze < 0) { 1275 for (j=0; j<x; j++) idx[nn++] = -1; 1276 } 1277 if (n20 >= 0) { /* right below */ 1278 x_t = lx[n20 % m]; 1279 y_t = ly[(n20 % (m*n))/m]; 1280 /* z_t = lz[n20 / (m*n)]; */ 1281 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1282 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1283 } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) { 1284 for (j=0; j<s_x; j++) idx[nn++] = -1; 1285 } 1286 } 1287 1288 for (i=0; i<y; i++) { 1289 if (n21 >= 0) { /* directly left */ 1290 x_t = lx[n21 % m]; 1291 y_t = y; 1292 /* z_t = lz[n21 / (m*n)]; */ 1293 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1294 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1295 } else if (Xs-xs < 0 && ze-Ze < 0) { 1296 for (j=0; j<s_x; j++) idx[nn++] = -1; 1297 } 1298 1299 if (n22 >= 0) { /* middle */ 1300 x_t = x; 1301 y_t = y; 1302 /* z_t = lz[n22 / (m*n)]; */ 1303 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1304 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1305 } else if (ze-Ze < 0) { 1306 if (bz == DMDA_BOUNDARY_MIRROR) { 1307 for (j=0; j<x; j++) idx[nn++] = 0; 1308 } else { 1309 for (j=0; j<x; j++) idx[nn++] = -1; 1310 } 1311 } 1312 1313 if (n23 >= 0) { /* directly right */ 1314 x_t = lx[n23 % m]; 1315 y_t = y; 1316 /* z_t = lz[n23 / (m*n)]; */ 1317 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1318 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1319 } else if (xe-Xe < 0 && ze-Ze < 0) { 1320 for (j=0; j<s_x; j++) idx[nn++] = -1; 1321 } 1322 } 1323 1324 for (i=1; i<=s_y; i++) { 1325 if (n24 >= 0) { /* left above */ 1326 x_t = lx[n24 % m]; 1327 y_t = ly[(n24 % (m*n))/m]; 1328 /* z_t = lz[n24 / (m*n)]; */ 1329 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1330 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1331 } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) { 1332 for (j=0; j<s_x; j++) idx[nn++] = -1; 1333 } 1334 if (n25 >= 0) { /* directly above */ 1335 x_t = x; 1336 y_t = ly[(n25 % (m*n))/m]; 1337 /* z_t = lz[n25 / (m*n)]; */ 1338 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1339 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1340 } else if (ye-Ye < 0 && ze-Ze < 0) { 1341 for (j=0; j<x; j++) idx[nn++] = -1; 1342 } 1343 if (n26 >= 0) { /* right above */ 1344 x_t = lx[n26 % m]; 1345 y_t = ly[(n26 % (m*n))/m]; 1346 /* z_t = lz[n26 / (m*n)]; */ 1347 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1348 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1349 } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) { 1350 for (j=0; j<s_x; j++) idx[nn++] = -1; 1351 } 1352 } 1353 } 1354 } 1355 /* 1356 Set the local to global ordering in the global vector, this allows use 1357 of VecSetValuesLocal(). 1358 */ 1359 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr); 1360 ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr); 1361 ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1362 ierr = ISGetIndices(ltogis, &idx_full);CHKERRQ(ierr); 1363 ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1364 ierr = ISRestoreIndices(ltogis, &idx_full);CHKERRQ(ierr); 1365 ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr); 1366 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1367 ierr = ISDestroy(<ogis);CHKERRQ(ierr); 1368 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 1369 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1370 1371 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1372 dd->m = m; dd->n = n; dd->p = p; 1373 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1374 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; 1375 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; 1376 1377 ierr = VecDestroy(&local);CHKERRQ(ierr); 1378 ierr = VecDestroy(&global);CHKERRQ(ierr); 1379 1380 dd->gtol = gtol; 1381 dd->ltog = ltog; 1382 dd->idx = idx_cpy; 1383 dd->Nl = nn*dof; 1384 dd->base = base; 1385 da->ops->view = DMView_DA_3d; 1386 dd->ltol = NULL; 1387 dd->ao = NULL; 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 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 + -dm_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 = DMViewFromOptions(*da,"-dm_view");CHKERRQ(ierr); 1470 PetscFunctionReturn(0); 1471 } 1472