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