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