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