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 = PetscViewerASCIIPushSynchronized(viewer);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 = PetscViewerASCIIPopSynchronized(viewer);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)/dd->w; 130 ierr = ISLocalToGlobalMappingGetBlockIndices(da->ltogmap,&idx);CHKERRQ(ierr); 131 plane=k; 132 /* Keep z wrap around points on the drawing */ 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++; 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,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 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 = PetscMalloc1((IXe-IXs)*(IYe-IYs)*(IZe-IZs),&idx);CHKERRQ(ierr); 432 left = xs - Xs; right = left + x; 433 bottom = ys - Ys; top = bottom + y; 434 down = zs - Zs; up = down + z; 435 count = 0; 436 for (i=down; i<up; i++) { 437 for (j=bottom; j<top; j++) { 438 for (k=left; k<right; k++) { 439 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 440 } 441 } 442 } 443 444 /* global to local must include ghost points within the domain, 445 but not ghost points outside the domain that aren't periodic */ 446 if (stencil_type == DMDA_STENCIL_BOX) { 447 left = IXs - Xs; right = left + (IXe-IXs); 448 bottom = IYs - Ys; top = bottom + (IYe-IYs); 449 down = IZs - Zs; up = down + (IZe-IZs); 450 count = 0; 451 for (i=down; i<up; i++) { 452 for (j=bottom; j<top; j++) { 453 for (k=left; k<right; k++) { 454 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 455 } 456 } 457 } 458 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 459 } else { 460 /* This is way ugly! We need to list the funny cross type region */ 461 left = xs - Xs; right = left + x; 462 bottom = ys - Ys; top = bottom + y; 463 down = zs - Zs; up = down + z; 464 count = 0; 465 /* the bottom chunck */ 466 for (i=(IZs-Zs); i<down; i++) { 467 for (j=bottom; j<top; j++) { 468 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 469 } 470 } 471 /* the middle piece */ 472 for (i=down; i<up; i++) { 473 /* front */ 474 for (j=(IYs-Ys); j<bottom; j++) { 475 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 476 } 477 /* middle */ 478 for (j=bottom; j<top; j++) { 479 for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 480 } 481 /* back */ 482 for (j=top; j<top+IYe-ye; j++) { 483 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 484 } 485 } 486 /* the top piece */ 487 for (i=up; i<up+IZe-ze; i++) { 488 for (j=bottom; j<top; j++) { 489 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 490 } 491 } 492 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 493 } 494 495 /* determine who lies on each side of use stored in n24 n25 n26 496 n21 n22 n23 497 n18 n19 n20 498 499 n15 n16 n17 500 n12 n14 501 n9 n10 n11 502 503 n6 n7 n8 504 n3 n4 n5 505 n0 n1 n2 506 */ 507 508 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ 509 /* Assume Nodes are Internal to the Cube */ 510 n0 = rank - m*n - m - 1; 511 n1 = rank - m*n - m; 512 n2 = rank - m*n - m + 1; 513 n3 = rank - m*n -1; 514 n4 = rank - m*n; 515 n5 = rank - m*n + 1; 516 n6 = rank - m*n + m - 1; 517 n7 = rank - m*n + m; 518 n8 = rank - m*n + m + 1; 519 520 n9 = rank - m - 1; 521 n10 = rank - m; 522 n11 = rank - m + 1; 523 n12 = rank - 1; 524 n14 = rank + 1; 525 n15 = rank + m - 1; 526 n16 = rank + m; 527 n17 = rank + m + 1; 528 529 n18 = rank + m*n - m - 1; 530 n19 = rank + m*n - m; 531 n20 = rank + m*n - m + 1; 532 n21 = rank + m*n - 1; 533 n22 = rank + m*n; 534 n23 = rank + m*n + 1; 535 n24 = rank + m*n + m - 1; 536 n25 = rank + m*n + m; 537 n26 = rank + m*n + m + 1; 538 539 /* Assume Pieces are on Faces of Cube */ 540 541 if (xs == 0) { /* First assume not corner or edge */ 542 n0 = rank -1 - (m*n); 543 n3 = rank + m -1 - (m*n); 544 n6 = rank + 2*m -1 - (m*n); 545 n9 = rank -1; 546 n12 = rank + m -1; 547 n15 = rank + 2*m -1; 548 n18 = rank -1 + (m*n); 549 n21 = rank + m -1 + (m*n); 550 n24 = rank + 2*m -1 + (m*n); 551 } 552 553 if (xe == M) { /* First assume not corner or edge */ 554 n2 = rank -2*m +1 - (m*n); 555 n5 = rank - m +1 - (m*n); 556 n8 = rank +1 - (m*n); 557 n11 = rank -2*m +1; 558 n14 = rank - m +1; 559 n17 = rank +1; 560 n20 = rank -2*m +1 + (m*n); 561 n23 = rank - m +1 + (m*n); 562 n26 = rank +1 + (m*n); 563 } 564 565 if (ys==0) { /* First assume not corner or edge */ 566 n0 = rank + m * (n-1) -1 - (m*n); 567 n1 = rank + m * (n-1) - (m*n); 568 n2 = rank + m * (n-1) +1 - (m*n); 569 n9 = rank + m * (n-1) -1; 570 n10 = rank + m * (n-1); 571 n11 = rank + m * (n-1) +1; 572 n18 = rank + m * (n-1) -1 + (m*n); 573 n19 = rank + m * (n-1) + (m*n); 574 n20 = rank + m * (n-1) +1 + (m*n); 575 } 576 577 if (ye == N) { /* First assume not corner or edge */ 578 n6 = rank - m * (n-1) -1 - (m*n); 579 n7 = rank - m * (n-1) - (m*n); 580 n8 = rank - m * (n-1) +1 - (m*n); 581 n15 = rank - m * (n-1) -1; 582 n16 = rank - m * (n-1); 583 n17 = rank - m * (n-1) +1; 584 n24 = rank - m * (n-1) -1 + (m*n); 585 n25 = rank - m * (n-1) + (m*n); 586 n26 = rank - m * (n-1) +1 + (m*n); 587 } 588 589 if (zs == 0) { /* First assume not corner or edge */ 590 n0 = size - (m*n) + rank - m - 1; 591 n1 = size - (m*n) + rank - m; 592 n2 = size - (m*n) + rank - m + 1; 593 n3 = size - (m*n) + rank - 1; 594 n4 = size - (m*n) + rank; 595 n5 = size - (m*n) + rank + 1; 596 n6 = size - (m*n) + rank + m - 1; 597 n7 = size - (m*n) + rank + m; 598 n8 = size - (m*n) + rank + m + 1; 599 } 600 601 if (ze == P) { /* First assume not corner or edge */ 602 n18 = (m*n) - (size-rank) - m - 1; 603 n19 = (m*n) - (size-rank) - m; 604 n20 = (m*n) - (size-rank) - m + 1; 605 n21 = (m*n) - (size-rank) - 1; 606 n22 = (m*n) - (size-rank); 607 n23 = (m*n) - (size-rank) + 1; 608 n24 = (m*n) - (size-rank) + m - 1; 609 n25 = (m*n) - (size-rank) + m; 610 n26 = (m*n) - (size-rank) + m + 1; 611 } 612 613 if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */ 614 n0 = size - m*n + rank + m-1 - m; 615 n3 = size - m*n + rank + m-1; 616 n6 = size - m*n + rank + m-1 + m; 617 } 618 619 if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */ 620 n18 = m*n - (size - rank) + m-1 - m; 621 n21 = m*n - (size - rank) + m-1; 622 n24 = m*n - (size - rank) + m-1 + m; 623 } 624 625 if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */ 626 n0 = rank + m*n -1 - m*n; 627 n9 = rank + m*n -1; 628 n18 = rank + m*n -1 + m*n; 629 } 630 631 if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */ 632 n6 = rank - m*(n-1) + m-1 - m*n; 633 n15 = rank - m*(n-1) + m-1; 634 n24 = rank - m*(n-1) + m-1 + m*n; 635 } 636 637 if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */ 638 n2 = size - (m*n-rank) - (m-1) - m; 639 n5 = size - (m*n-rank) - (m-1); 640 n8 = size - (m*n-rank) - (m-1) + m; 641 } 642 643 if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */ 644 n20 = m*n - (size - rank) - (m-1) - m; 645 n23 = m*n - (size - rank) - (m-1); 646 n26 = m*n - (size - rank) - (m-1) + m; 647 } 648 649 if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */ 650 n2 = rank + m*(n-1) - (m-1) - m*n; 651 n11 = rank + m*(n-1) - (m-1); 652 n20 = rank + m*(n-1) - (m-1) + m*n; 653 } 654 655 if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */ 656 n8 = rank - m*n +1 - m*n; 657 n17 = rank - m*n +1; 658 n26 = rank - m*n +1 + m*n; 659 } 660 661 if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */ 662 n0 = size - m + rank -1; 663 n1 = size - m + rank; 664 n2 = size - m + rank +1; 665 } 666 667 if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */ 668 n18 = m*n - (size - rank) + m*(n-1) -1; 669 n19 = m*n - (size - rank) + m*(n-1); 670 n20 = m*n - (size - rank) + m*(n-1) +1; 671 } 672 673 if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */ 674 n6 = size - (m*n-rank) - m * (n-1) -1; 675 n7 = size - (m*n-rank) - m * (n-1); 676 n8 = size - (m*n-rank) - m * (n-1) +1; 677 } 678 679 if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */ 680 n24 = rank - (size-m) -1; 681 n25 = rank - (size-m); 682 n26 = rank - (size-m) +1; 683 } 684 685 /* Check for Corners */ 686 if ((xs==0) && (ys==0) && (zs==0)) n0 = size -1; 687 if ((xs==0) && (ys==0) && (ze==P)) n18 = m*n-1; 688 if ((xs==0) && (ye==N) && (zs==0)) n6 = (size-1)-m*(n-1); 689 if ((xs==0) && (ye==N) && (ze==P)) n24 = m-1; 690 if ((xe==M) && (ys==0) && (zs==0)) n2 = size-m; 691 if ((xe==M) && (ys==0) && (ze==P)) n20 = m*n-m; 692 if ((xe==M) && (ye==N) && (zs==0)) n8 = size-m*n; 693 if ((xe==M) && (ye==N) && (ze==P)) n26 = 0; 694 695 /* Check for when not X,Y, and Z Periodic */ 696 697 /* If not X periodic */ 698 if (bx != DM_BOUNDARY_PERIODIC) { 699 if (xs==0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2; 700 if (xe==M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2; 701 } 702 703 /* If not Y periodic */ 704 if (by != DM_BOUNDARY_PERIODIC) { 705 if (ys==0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2; 706 if (ye==N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2; 707 } 708 709 /* If not Z periodic */ 710 if (bz != DM_BOUNDARY_PERIODIC) { 711 if (zs==0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2; 712 if (ze==P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2; 713 } 714 715 ierr = PetscMalloc1(27,&dd->neighbors);CHKERRQ(ierr); 716 717 dd->neighbors[0] = n0; 718 dd->neighbors[1] = n1; 719 dd->neighbors[2] = n2; 720 dd->neighbors[3] = n3; 721 dd->neighbors[4] = n4; 722 dd->neighbors[5] = n5; 723 dd->neighbors[6] = n6; 724 dd->neighbors[7] = n7; 725 dd->neighbors[8] = n8; 726 dd->neighbors[9] = n9; 727 dd->neighbors[10] = n10; 728 dd->neighbors[11] = n11; 729 dd->neighbors[12] = n12; 730 dd->neighbors[13] = rank; 731 dd->neighbors[14] = n14; 732 dd->neighbors[15] = n15; 733 dd->neighbors[16] = n16; 734 dd->neighbors[17] = n17; 735 dd->neighbors[18] = n18; 736 dd->neighbors[19] = n19; 737 dd->neighbors[20] = n20; 738 dd->neighbors[21] = n21; 739 dd->neighbors[22] = n22; 740 dd->neighbors[23] = n23; 741 dd->neighbors[24] = n24; 742 dd->neighbors[25] = n25; 743 dd->neighbors[26] = n26; 744 745 /* If star stencil then delete the corner neighbors */ 746 if (stencil_type == DMDA_STENCIL_STAR) { 747 /* save information about corner neighbors */ 748 sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7; 749 sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18; 750 sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25; 751 sn26 = n26; 752 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; 753 } 754 755 ierr = PetscMalloc1((Xe-Xs)*(Ye-Ys)*(Ze-Zs),&idx);CHKERRQ(ierr); 756 757 nn = 0; 758 /* Bottom Level */ 759 for (k=0; k<s_z; k++) { 760 for (i=1; i<=s_y; i++) { 761 if (n0 >= 0) { /* left below */ 762 x_t = lx[n0 % m]; 763 y_t = ly[(n0 % (m*n))/m]; 764 z_t = lz[n0 / (m*n)]; 765 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; 766 if (twod && (s_t < 0)) s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x; /* 2D case */ 767 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 768 } 769 if (n1 >= 0) { /* directly below */ 770 x_t = x; 771 y_t = ly[(n1 % (m*n))/m]; 772 z_t = lz[n1 / (m*n)]; 773 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 774 if (twod && (s_t < 0)) s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 775 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 776 } 777 if (n2 >= 0) { /* right below */ 778 x_t = lx[n2 % m]; 779 y_t = ly[(n2 % (m*n))/m]; 780 z_t = lz[n2 / (m*n)]; 781 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 782 if (twod && (s_t < 0)) s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 783 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 784 } 785 } 786 787 for (i=0; i<y; i++) { 788 if (n3 >= 0) { /* directly left */ 789 x_t = lx[n3 % m]; 790 y_t = y; 791 z_t = lz[n3 / (m*n)]; 792 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 793 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 */ 794 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 795 } 796 797 if (n4 >= 0) { /* middle */ 798 x_t = x; 799 y_t = y; 800 z_t = lz[n4 / (m*n)]; 801 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 802 if (twod && (s_t < 0)) s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 803 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 804 } else if (bz == DM_BOUNDARY_MIRROR) { 805 for (j=0; j<x; j++) idx[nn++] = 0; 806 } 807 808 if (n5 >= 0) { /* directly right */ 809 x_t = lx[n5 % m]; 810 y_t = y; 811 z_t = lz[n5 / (m*n)]; 812 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 813 if (twod && (s_t < 0)) s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 814 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 815 } 816 } 817 818 for (i=1; i<=s_y; i++) { 819 if (n6 >= 0) { /* left above */ 820 x_t = lx[n6 % m]; 821 y_t = ly[(n6 % (m*n))/m]; 822 z_t = lz[n6 / (m*n)]; 823 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 824 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 */ 825 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 826 } 827 if (n7 >= 0) { /* directly above */ 828 x_t = x; 829 y_t = ly[(n7 % (m*n))/m]; 830 z_t = lz[n7 / (m*n)]; 831 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 832 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 */ 833 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 834 } 835 if (n8 >= 0) { /* right above */ 836 x_t = lx[n8 % m]; 837 y_t = ly[(n8 % (m*n))/m]; 838 z_t = lz[n8 / (m*n)]; 839 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 840 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 */ 841 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 842 } 843 } 844 } 845 846 /* Middle Level */ 847 for (k=0; k<z; k++) { 848 for (i=1; i<=s_y; i++) { 849 if (n9 >= 0) { /* left below */ 850 x_t = lx[n9 % m]; 851 y_t = ly[(n9 % (m*n))/m]; 852 /* z_t = z; */ 853 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 854 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 855 } 856 if (n10 >= 0) { /* directly below */ 857 x_t = x; 858 y_t = ly[(n10 % (m*n))/m]; 859 /* z_t = z; */ 860 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 861 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 862 } else if (by == DM_BOUNDARY_MIRROR) { 863 for (j=0; j<x; j++) idx[nn++] = 0; 864 } 865 if (n11 >= 0) { /* right below */ 866 x_t = lx[n11 % m]; 867 y_t = ly[(n11 % (m*n))/m]; 868 /* z_t = z; */ 869 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 870 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 871 } 872 } 873 874 for (i=0; i<y; i++) { 875 if (n12 >= 0) { /* directly left */ 876 x_t = lx[n12 % m]; 877 y_t = y; 878 /* z_t = z; */ 879 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 880 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 881 } else if (bx == DM_BOUNDARY_MIRROR) { 882 for (j=0; j<s_x; j++) idx[nn++] = 0; 883 } 884 885 /* Interior */ 886 s_t = bases[rank] + i*x + k*x*y; 887 for (j=0; j<x; j++) idx[nn++] = s_t++; 888 889 if (n14 >= 0) { /* directly right */ 890 x_t = lx[n14 % m]; 891 y_t = y; 892 /* z_t = z; */ 893 s_t = bases[n14] + i*x_t + k*x_t*y_t; 894 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 895 } else if (bx == DM_BOUNDARY_MIRROR) { 896 for (j=0; j<s_x; j++) idx[nn++] = 0; 897 } 898 } 899 900 for (i=1; i<=s_y; i++) { 901 if (n15 >= 0) { /* left above */ 902 x_t = lx[n15 % m]; 903 y_t = ly[(n15 % (m*n))/m]; 904 /* z_t = z; */ 905 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 906 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 907 } 908 if (n16 >= 0) { /* directly above */ 909 x_t = x; 910 y_t = ly[(n16 % (m*n))/m]; 911 /* z_t = z; */ 912 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 913 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 914 } else if (by == DM_BOUNDARY_MIRROR) { 915 for (j=0; j<x; j++) idx[nn++] = 0; 916 } 917 if (n17 >= 0) { /* right above */ 918 x_t = lx[n17 % m]; 919 y_t = ly[(n17 % (m*n))/m]; 920 /* z_t = z; */ 921 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 922 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 923 } 924 } 925 } 926 927 /* Upper Level */ 928 for (k=0; k<s_z; k++) { 929 for (i=1; i<=s_y; i++) { 930 if (n18 >= 0) { /* left below */ 931 x_t = lx[n18 % m]; 932 y_t = ly[(n18 % (m*n))/m]; 933 /* z_t = lz[n18 / (m*n)]; */ 934 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 935 if (twod && (s_t >= M*N*P)) s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t; /* 2d case */ 936 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 937 } 938 if (n19 >= 0) { /* directly below */ 939 x_t = x; 940 y_t = ly[(n19 % (m*n))/m]; 941 /* z_t = lz[n19 / (m*n)]; */ 942 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 943 if (twod && (s_t >= M*N*P)) s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 944 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 945 } 946 if (n20 >= 0) { /* right below */ 947 x_t = lx[n20 % m]; 948 y_t = ly[(n20 % (m*n))/m]; 949 /* z_t = lz[n20 / (m*n)]; */ 950 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 951 if (twod && (s_t >= M*N*P)) s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 952 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 953 } 954 } 955 956 for (i=0; i<y; i++) { 957 if (n21 >= 0) { /* directly left */ 958 x_t = lx[n21 % m]; 959 y_t = y; 960 /* z_t = lz[n21 / (m*n)]; */ 961 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 962 if (twod && (s_t >= M*N*P)) s_t = bases[n21] + (i+1)*x_t - s_x; /* 2d case */ 963 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 964 } 965 966 if (n22 >= 0) { /* middle */ 967 x_t = x; 968 y_t = y; 969 /* z_t = lz[n22 / (m*n)]; */ 970 s_t = bases[n22] + i*x_t + k*x_t*y_t; 971 if (twod && (s_t >= M*N*P)) s_t = bases[n22] + i*x_t; /* 2d case */ 972 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 973 } else if (bz == DM_BOUNDARY_MIRROR) { 974 for (j=0; j<x; j++) idx[nn++] = 0; 975 } 976 977 if (n23 >= 0) { /* directly right */ 978 x_t = lx[n23 % m]; 979 y_t = y; 980 /* z_t = lz[n23 / (m*n)]; */ 981 s_t = bases[n23] + i*x_t + k*x_t*y_t; 982 if (twod && (s_t >= M*N*P)) s_t = bases[n23] + i*x_t; /* 2d case */ 983 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 984 } 985 } 986 987 for (i=1; i<=s_y; i++) { 988 if (n24 >= 0) { /* left above */ 989 x_t = lx[n24 % m]; 990 y_t = ly[(n24 % (m*n))/m]; 991 /* z_t = lz[n24 / (m*n)]; */ 992 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 993 if (twod && (s_t >= M*N*P)) s_t = bases[n24] + i*x_t - s_x; /* 2d case */ 994 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 995 } 996 if (n25 >= 0) { /* directly above */ 997 x_t = x; 998 y_t = ly[(n25 % (m*n))/m]; 999 /* z_t = lz[n25 / (m*n)]; */ 1000 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1001 if (twod && (s_t >= M*N*P)) s_t = bases[n25] + (i-1)*x_t; /* 2d case */ 1002 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1003 } 1004 if (n26 >= 0) { /* right above */ 1005 x_t = lx[n26 % m]; 1006 y_t = ly[(n26 % (m*n))/m]; 1007 /* z_t = lz[n26 / (m*n)]; */ 1008 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1009 if (twod && (s_t >= M*N*P)) s_t = bases[n26] + (i-1)*x_t; /* 2d case */ 1010 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1011 } 1012 } 1013 } 1014 1015 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 1016 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1017 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);CHKERRQ(ierr); 1018 ierr = ISDestroy(&to);CHKERRQ(ierr); 1019 ierr = ISDestroy(&from);CHKERRQ(ierr); 1020 1021 if (stencil_type == DMDA_STENCIL_STAR) { 1022 n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; 1023 n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; 1024 n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; 1025 n26 = sn26; 1026 } 1027 1028 if (((stencil_type == DMDA_STENCIL_STAR) || 1029 (bx != DM_BOUNDARY_PERIODIC && bx) || 1030 (by != DM_BOUNDARY_PERIODIC && by) || 1031 (bz != DM_BOUNDARY_PERIODIC && bz))) { 1032 /* 1033 Recompute the local to global mappings, this time keeping the 1034 information about the cross corner processor numbers. 1035 */ 1036 nn = 0; 1037 /* Bottom Level */ 1038 for (k=0; k<s_z; k++) { 1039 for (i=1; i<=s_y; i++) { 1040 if (n0 >= 0) { /* left below */ 1041 x_t = lx[n0 % m]; 1042 y_t = ly[(n0 % (m*n))/m]; 1043 z_t = lz[n0 / (m*n)]; 1044 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; 1045 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1046 } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) { 1047 for (j=0; j<s_x; j++) idx[nn++] = -1; 1048 } 1049 if (n1 >= 0) { /* directly below */ 1050 x_t = x; 1051 y_t = ly[(n1 % (m*n))/m]; 1052 z_t = lz[n1 / (m*n)]; 1053 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1054 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1055 } else if (Ys-ys < 0 && Zs-zs < 0) { 1056 for (j=0; j<x; j++) idx[nn++] = -1; 1057 } 1058 if (n2 >= 0) { /* right below */ 1059 x_t = lx[n2 % m]; 1060 y_t = ly[(n2 % (m*n))/m]; 1061 z_t = lz[n2 / (m*n)]; 1062 s_t = bases[n2] + 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<s_x; j++) idx[nn++] = s_t++; 1064 } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) { 1065 for (j=0; j<s_x; j++) idx[nn++] = -1; 1066 } 1067 } 1068 1069 for (i=0; i<y; i++) { 1070 if (n3 >= 0) { /* directly left */ 1071 x_t = lx[n3 % m]; 1072 y_t = y; 1073 z_t = lz[n3 / (m*n)]; 1074 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1075 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1076 } else if (Xs-xs < 0 && Zs-zs < 0) { 1077 for (j=0; j<s_x; j++) idx[nn++] = -1; 1078 } 1079 1080 if (n4 >= 0) { /* middle */ 1081 x_t = x; 1082 y_t = y; 1083 z_t = lz[n4 / (m*n)]; 1084 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1085 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1086 } else if (Zs-zs < 0) { 1087 if (bz == DM_BOUNDARY_MIRROR) { 1088 for (j=0; j<x; j++) idx[nn++] = 0; 1089 } else { 1090 for (j=0; j<x; j++) idx[nn++] = -1; 1091 } 1092 } 1093 1094 if (n5 >= 0) { /* directly right */ 1095 x_t = lx[n5 % m]; 1096 y_t = y; 1097 z_t = lz[n5 / (m*n)]; 1098 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1099 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1100 } else if (xe-Xe < 0 && Zs-zs < 0) { 1101 for (j=0; j<s_x; j++) idx[nn++] = -1; 1102 } 1103 } 1104 1105 for (i=1; i<=s_y; i++) { 1106 if (n6 >= 0) { /* left above */ 1107 x_t = lx[n6 % m]; 1108 y_t = ly[(n6 % (m*n))/m]; 1109 z_t = lz[n6 / (m*n)]; 1110 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1111 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1112 } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) { 1113 for (j=0; j<s_x; j++) idx[nn++] = -1; 1114 } 1115 if (n7 >= 0) { /* directly above */ 1116 x_t = x; 1117 y_t = ly[(n7 % (m*n))/m]; 1118 z_t = lz[n7 / (m*n)]; 1119 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1120 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1121 } else if (ye-Ye < 0 && Zs-zs < 0) { 1122 for (j=0; j<x; j++) idx[nn++] = -1; 1123 } 1124 if (n8 >= 0) { /* right above */ 1125 x_t = lx[n8 % m]; 1126 y_t = ly[(n8 % (m*n))/m]; 1127 z_t = lz[n8 / (m*n)]; 1128 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1129 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1130 } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) { 1131 for (j=0; j<s_x; j++) idx[nn++] = -1; 1132 } 1133 } 1134 } 1135 1136 /* Middle Level */ 1137 for (k=0; k<z; k++) { 1138 for (i=1; i<=s_y; i++) { 1139 if (n9 >= 0) { /* left below */ 1140 x_t = lx[n9 % m]; 1141 y_t = ly[(n9 % (m*n))/m]; 1142 /* z_t = z; */ 1143 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1144 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1145 } else if (Xs-xs < 0 && Ys-ys < 0) { 1146 for (j=0; j<s_x; j++) idx[nn++] = -1; 1147 } 1148 if (n10 >= 0) { /* directly below */ 1149 x_t = x; 1150 y_t = ly[(n10 % (m*n))/m]; 1151 /* z_t = z; */ 1152 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1153 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1154 } else if (Ys-ys < 0) { 1155 if (by == DM_BOUNDARY_MIRROR) { 1156 for (j=0; j<x; j++) idx[nn++] = -1; 1157 } else { 1158 for (j=0; j<x; j++) idx[nn++] = -1; 1159 } 1160 } 1161 if (n11 >= 0) { /* right below */ 1162 x_t = lx[n11 % m]; 1163 y_t = ly[(n11 % (m*n))/m]; 1164 /* z_t = z; */ 1165 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1166 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1167 } else if (xe-Xe < 0 && Ys-ys < 0) { 1168 for (j=0; j<s_x; j++) idx[nn++] = -1; 1169 } 1170 } 1171 1172 for (i=0; i<y; i++) { 1173 if (n12 >= 0) { /* directly left */ 1174 x_t = lx[n12 % m]; 1175 y_t = y; 1176 /* z_t = z; */ 1177 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 1178 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1179 } else if (Xs-xs < 0) { 1180 if (bx == DM_BOUNDARY_MIRROR) { 1181 for (j=0; j<s_x; j++) idx[nn++] = 0; 1182 } else { 1183 for (j=0; j<s_x; j++) idx[nn++] = -1; 1184 } 1185 } 1186 1187 /* Interior */ 1188 s_t = bases[rank] + i*x + k*x*y; 1189 for (j=0; j<x; j++) idx[nn++] = s_t++; 1190 1191 if (n14 >= 0) { /* directly right */ 1192 x_t = lx[n14 % m]; 1193 y_t = y; 1194 /* z_t = z; */ 1195 s_t = bases[n14] + i*x_t + k*x_t*y_t; 1196 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1197 } else if (xe-Xe < 0) { 1198 if (bx == DM_BOUNDARY_MIRROR) { 1199 for (j=0; j<s_x; j++) idx[nn++] = 0; 1200 } else { 1201 for (j=0; j<s_x; j++) idx[nn++] = -1; 1202 } 1203 } 1204 } 1205 1206 for (i=1; i<=s_y; i++) { 1207 if (n15 >= 0) { /* left above */ 1208 x_t = lx[n15 % m]; 1209 y_t = ly[(n15 % (m*n))/m]; 1210 /* z_t = z; */ 1211 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 1212 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1213 } else if (Xs-xs < 0 && ye-Ye < 0) { 1214 for (j=0; j<s_x; j++) idx[nn++] = -1; 1215 } 1216 if (n16 >= 0) { /* directly above */ 1217 x_t = x; 1218 y_t = ly[(n16 % (m*n))/m]; 1219 /* z_t = z; */ 1220 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 1221 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1222 } else if (ye-Ye < 0) { 1223 if (by == DM_BOUNDARY_MIRROR) { 1224 for (j=0; j<x; j++) idx[nn++] = 0; 1225 } else { 1226 for (j=0; j<x; j++) idx[nn++] = -1; 1227 } 1228 } 1229 if (n17 >= 0) { /* right above */ 1230 x_t = lx[n17 % m]; 1231 y_t = ly[(n17 % (m*n))/m]; 1232 /* z_t = z; */ 1233 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 1234 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1235 } else if (xe-Xe < 0 && ye-Ye < 0) { 1236 for (j=0; j<s_x; j++) idx[nn++] = -1; 1237 } 1238 } 1239 } 1240 1241 /* Upper Level */ 1242 for (k=0; k<s_z; k++) { 1243 for (i=1; i<=s_y; i++) { 1244 if (n18 >= 0) { /* left below */ 1245 x_t = lx[n18 % m]; 1246 y_t = ly[(n18 % (m*n))/m]; 1247 /* z_t = lz[n18 / (m*n)]; */ 1248 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1249 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1250 } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) { 1251 for (j=0; j<s_x; j++) idx[nn++] = -1; 1252 } 1253 if (n19 >= 0) { /* directly below */ 1254 x_t = x; 1255 y_t = ly[(n19 % (m*n))/m]; 1256 /* z_t = lz[n19 / (m*n)]; */ 1257 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1258 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1259 } else if (Ys-ys < 0 && ze-Ze < 0) { 1260 for (j=0; j<x; j++) idx[nn++] = -1; 1261 } 1262 if (n20 >= 0) { /* right below */ 1263 x_t = lx[n20 % m]; 1264 y_t = ly[(n20 % (m*n))/m]; 1265 /* z_t = lz[n20 / (m*n)]; */ 1266 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1267 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1268 } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) { 1269 for (j=0; j<s_x; j++) idx[nn++] = -1; 1270 } 1271 } 1272 1273 for (i=0; i<y; i++) { 1274 if (n21 >= 0) { /* directly left */ 1275 x_t = lx[n21 % m]; 1276 y_t = y; 1277 /* z_t = lz[n21 / (m*n)]; */ 1278 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1279 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1280 } else if (Xs-xs < 0 && ze-Ze < 0) { 1281 for (j=0; j<s_x; j++) idx[nn++] = -1; 1282 } 1283 1284 if (n22 >= 0) { /* middle */ 1285 x_t = x; 1286 y_t = y; 1287 /* z_t = lz[n22 / (m*n)]; */ 1288 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1289 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1290 } else if (ze-Ze < 0) { 1291 if (bz == DM_BOUNDARY_MIRROR) { 1292 for (j=0; j<x; j++) idx[nn++] = 0; 1293 } else { 1294 for (j=0; j<x; j++) idx[nn++] = -1; 1295 } 1296 } 1297 1298 if (n23 >= 0) { /* directly right */ 1299 x_t = lx[n23 % m]; 1300 y_t = y; 1301 /* z_t = lz[n23 / (m*n)]; */ 1302 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1303 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1304 } else if (xe-Xe < 0 && ze-Ze < 0) { 1305 for (j=0; j<s_x; j++) idx[nn++] = -1; 1306 } 1307 } 1308 1309 for (i=1; i<=s_y; i++) { 1310 if (n24 >= 0) { /* left above */ 1311 x_t = lx[n24 % m]; 1312 y_t = ly[(n24 % (m*n))/m]; 1313 /* z_t = lz[n24 / (m*n)]; */ 1314 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1315 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1316 } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) { 1317 for (j=0; j<s_x; j++) idx[nn++] = -1; 1318 } 1319 if (n25 >= 0) { /* directly above */ 1320 x_t = x; 1321 y_t = ly[(n25 % (m*n))/m]; 1322 /* z_t = lz[n25 / (m*n)]; */ 1323 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1324 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1325 } else if (ye-Ye < 0 && ze-Ze < 0) { 1326 for (j=0; j<x; j++) idx[nn++] = -1; 1327 } 1328 if (n26 >= 0) { /* right above */ 1329 x_t = lx[n26 % m]; 1330 y_t = ly[(n26 % (m*n))/m]; 1331 /* z_t = lz[n26 / (m*n)]; */ 1332 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1333 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1334 } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) { 1335 for (j=0; j<s_x; j++) idx[nn++] = -1; 1336 } 1337 } 1338 } 1339 } 1340 /* 1341 Set the local to global ordering in the global vector, this allows use 1342 of VecSetValuesLocal(). 1343 */ 1344 ierr = ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr); 1345 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);CHKERRQ(ierr); 1346 1347 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1348 dd->m = m; dd->n = n; dd->p = p; 1349 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1350 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; 1351 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; 1352 1353 ierr = VecDestroy(&local);CHKERRQ(ierr); 1354 ierr = VecDestroy(&global);CHKERRQ(ierr); 1355 1356 dd->gtol = gtol; 1357 dd->base = base; 1358 da->ops->view = DMView_DA_3d; 1359 dd->ltol = NULL; 1360 dd->ao = NULL; 1361 PetscFunctionReturn(0); 1362 } 1363 1364 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "DMDACreate3d" 1367 /*@C 1368 DMDACreate3d - Creates an object that will manage the communication of three-dimensional 1369 regular array data that is distributed across some processors. 1370 1371 Collective on MPI_Comm 1372 1373 Input Parameters: 1374 + comm - MPI communicator 1375 . bx,by,bz - type of ghost nodes the array have. 1376 Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC. 1377 . stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX) 1378 . 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 1379 from the command line with -da_grid_x <M> -da_grid_y <N> -da_grid_z <P>) 1380 . m,n,p - corresponding number of processors in each dimension 1381 (or PETSC_DECIDE to have calculated) 1382 . dof - number of degrees of freedom per node 1383 . s - stencil width 1384 - lx, ly, lz - arrays containing the number of nodes in each cell along 1385 the x, y, and z coordinates, or NULL. If non-null, these 1386 must be of length as m,n,p and the corresponding 1387 m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of 1388 the ly[] must N, sum of the lz[] must be P 1389 1390 Output Parameter: 1391 . da - the resulting distributed array object 1392 1393 Options Database Key: 1394 + -dm_view - Calls DMView() at the conclusion of DMDACreate3d() 1395 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 1396 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 1397 . -da_grid_z <nz> - number of grid points in z direction, if P < 0 1398 . -da_processors_x <MX> - number of processors in x direction 1399 . -da_processors_y <MY> - number of processors in y direction 1400 . -da_processors_z <MZ> - number of processors in z direction 1401 . -da_refine_x <rx> - refinement ratio in x direction 1402 . -da_refine_y <ry> - refinement ratio in y direction 1403 . -da_refine_z <rz>- refinement ratio in z directio 1404 - -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0 1405 1406 Level: beginner 1407 1408 Notes: 1409 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 1410 standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 1411 the standard 27-pt stencil. 1412 1413 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 1414 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 1415 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 1416 1417 .keywords: distributed array, create, three-dimensional 1418 1419 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 1420 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(), 1421 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 1422 1423 @*/ 1424 PetscErrorCode DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M, 1425 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) 1426 { 1427 PetscErrorCode ierr; 1428 1429 PetscFunctionBegin; 1430 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 1431 ierr = DMSetDimension(*da, 3);CHKERRQ(ierr); 1432 ierr = DMDASetSizes(*da, M, N, P);CHKERRQ(ierr); 1433 ierr = DMDASetNumProcs(*da, m, n, p);CHKERRQ(ierr); 1434 ierr = DMDASetBoundaryType(*da, bx, by, bz);CHKERRQ(ierr); 1435 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 1436 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 1437 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 1438 ierr = DMDASetOwnershipRanges(*da, lx, ly, lz);CHKERRQ(ierr); 1439 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 1440 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 1441 ierr = DMSetUp(*da);CHKERRQ(ierr); 1442 PetscFunctionReturn(0); 1443 } 1444