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