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