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