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