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