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