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 appropriate vector scatters */ 437 /* local to global inserts non-ghost point region into global */ 438 ierr = PetscMalloc1((IXe-IXs)*(IYe-IYs)*(IZe-IZs),&idx);CHKERRQ(ierr); 439 left = xs - Xs; right = left + x; 440 bottom = ys - Ys; top = bottom + y; 441 down = zs - Zs; up = down + z; 442 count = 0; 443 for (i=down; i<up; i++) { 444 for (j=bottom; j<top; j++) { 445 for (k=left; k<right; k++) { 446 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 447 } 448 } 449 } 450 451 /* global to local must include ghost points within the domain, 452 but not ghost points outside the domain that aren't periodic */ 453 if (stencil_type == DMDA_STENCIL_BOX) { 454 left = IXs - Xs; right = left + (IXe-IXs); 455 bottom = IYs - Ys; top = bottom + (IYe-IYs); 456 down = IZs - Zs; up = down + (IZe-IZs); 457 count = 0; 458 for (i=down; i<up; i++) { 459 for (j=bottom; j<top; j++) { 460 for (k=left; k<right; k++) { 461 idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 462 } 463 } 464 } 465 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 466 } else { 467 /* This is way ugly! We need to list the funny cross type region */ 468 left = xs - Xs; right = left + x; 469 bottom = ys - Ys; top = bottom + y; 470 down = zs - Zs; up = down + z; 471 count = 0; 472 /* the bottom chunck */ 473 for (i=(IZs-Zs); i<down; i++) { 474 for (j=bottom; j<top; j++) { 475 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 476 } 477 } 478 /* the middle piece */ 479 for (i=down; i<up; i++) { 480 /* front */ 481 for (j=(IYs-Ys); j<bottom; j++) { 482 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 483 } 484 /* middle */ 485 for (j=bottom; j<top; j++) { 486 for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 487 } 488 /* back */ 489 for (j=top; j<top+IYe-ye; j++) { 490 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 491 } 492 } 493 /* the top piece */ 494 for (i=up; i<up+IZe-ze; i++) { 495 for (j=bottom; j<top; j++) { 496 for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k; 497 } 498 } 499 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 500 } 501 502 /* determine who lies on each side of use stored in n24 n25 n26 503 n21 n22 n23 504 n18 n19 n20 505 506 n15 n16 n17 507 n12 n14 508 n9 n10 n11 509 510 n6 n7 n8 511 n3 n4 n5 512 n0 n1 n2 513 */ 514 515 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ 516 /* Assume Nodes are Internal to the Cube */ 517 n0 = rank - m*n - m - 1; 518 n1 = rank - m*n - m; 519 n2 = rank - m*n - m + 1; 520 n3 = rank - m*n -1; 521 n4 = rank - m*n; 522 n5 = rank - m*n + 1; 523 n6 = rank - m*n + m - 1; 524 n7 = rank - m*n + m; 525 n8 = rank - m*n + m + 1; 526 527 n9 = rank - m - 1; 528 n10 = rank - m; 529 n11 = rank - m + 1; 530 n12 = rank - 1; 531 n14 = rank + 1; 532 n15 = rank + m - 1; 533 n16 = rank + m; 534 n17 = rank + m + 1; 535 536 n18 = rank + m*n - m - 1; 537 n19 = rank + m*n - m; 538 n20 = rank + m*n - m + 1; 539 n21 = rank + m*n - 1; 540 n22 = rank + m*n; 541 n23 = rank + m*n + 1; 542 n24 = rank + m*n + m - 1; 543 n25 = rank + m*n + m; 544 n26 = rank + m*n + m + 1; 545 546 /* Assume Pieces are on Faces of Cube */ 547 548 if (xs == 0) { /* First assume not corner or edge */ 549 n0 = rank -1 - (m*n); 550 n3 = rank + m -1 - (m*n); 551 n6 = rank + 2*m -1 - (m*n); 552 n9 = rank -1; 553 n12 = rank + m -1; 554 n15 = rank + 2*m -1; 555 n18 = rank -1 + (m*n); 556 n21 = rank + m -1 + (m*n); 557 n24 = rank + 2*m -1 + (m*n); 558 } 559 560 if (xe == M) { /* First assume not corner or edge */ 561 n2 = rank -2*m +1 - (m*n); 562 n5 = rank - m +1 - (m*n); 563 n8 = rank +1 - (m*n); 564 n11 = rank -2*m +1; 565 n14 = rank - m +1; 566 n17 = rank +1; 567 n20 = rank -2*m +1 + (m*n); 568 n23 = rank - m +1 + (m*n); 569 n26 = rank +1 + (m*n); 570 } 571 572 if (ys==0) { /* First assume not corner or edge */ 573 n0 = rank + m * (n-1) -1 - (m*n); 574 n1 = rank + m * (n-1) - (m*n); 575 n2 = rank + m * (n-1) +1 - (m*n); 576 n9 = rank + m * (n-1) -1; 577 n10 = rank + m * (n-1); 578 n11 = rank + m * (n-1) +1; 579 n18 = rank + m * (n-1) -1 + (m*n); 580 n19 = rank + m * (n-1) + (m*n); 581 n20 = rank + m * (n-1) +1 + (m*n); 582 } 583 584 if (ye == N) { /* First assume not corner or edge */ 585 n6 = rank - m * (n-1) -1 - (m*n); 586 n7 = rank - m * (n-1) - (m*n); 587 n8 = rank - m * (n-1) +1 - (m*n); 588 n15 = rank - m * (n-1) -1; 589 n16 = rank - m * (n-1); 590 n17 = rank - m * (n-1) +1; 591 n24 = rank - m * (n-1) -1 + (m*n); 592 n25 = rank - m * (n-1) + (m*n); 593 n26 = rank - m * (n-1) +1 + (m*n); 594 } 595 596 if (zs == 0) { /* First assume not corner or edge */ 597 n0 = size - (m*n) + rank - m - 1; 598 n1 = size - (m*n) + rank - m; 599 n2 = size - (m*n) + rank - m + 1; 600 n3 = size - (m*n) + rank - 1; 601 n4 = size - (m*n) + rank; 602 n5 = size - (m*n) + rank + 1; 603 n6 = size - (m*n) + rank + m - 1; 604 n7 = size - (m*n) + rank + m; 605 n8 = size - (m*n) + rank + m + 1; 606 } 607 608 if (ze == P) { /* First assume not corner or edge */ 609 n18 = (m*n) - (size-rank) - m - 1; 610 n19 = (m*n) - (size-rank) - m; 611 n20 = (m*n) - (size-rank) - m + 1; 612 n21 = (m*n) - (size-rank) - 1; 613 n22 = (m*n) - (size-rank); 614 n23 = (m*n) - (size-rank) + 1; 615 n24 = (m*n) - (size-rank) + m - 1; 616 n25 = (m*n) - (size-rank) + m; 617 n26 = (m*n) - (size-rank) + m + 1; 618 } 619 620 if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */ 621 n0 = size - m*n + rank + m-1 - m; 622 n3 = size - m*n + rank + m-1; 623 n6 = size - m*n + rank + m-1 + m; 624 } 625 626 if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */ 627 n18 = m*n - (size - rank) + m-1 - m; 628 n21 = m*n - (size - rank) + m-1; 629 n24 = m*n - (size - rank) + m-1 + m; 630 } 631 632 if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */ 633 n0 = rank + m*n -1 - m*n; 634 n9 = rank + m*n -1; 635 n18 = rank + m*n -1 + m*n; 636 } 637 638 if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */ 639 n6 = rank - m*(n-1) + m-1 - m*n; 640 n15 = rank - m*(n-1) + m-1; 641 n24 = rank - m*(n-1) + m-1 + m*n; 642 } 643 644 if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */ 645 n2 = size - (m*n-rank) - (m-1) - m; 646 n5 = size - (m*n-rank) - (m-1); 647 n8 = size - (m*n-rank) - (m-1) + m; 648 } 649 650 if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */ 651 n20 = m*n - (size - rank) - (m-1) - m; 652 n23 = m*n - (size - rank) - (m-1); 653 n26 = m*n - (size - rank) - (m-1) + m; 654 } 655 656 if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */ 657 n2 = rank + m*(n-1) - (m-1) - m*n; 658 n11 = rank + m*(n-1) - (m-1); 659 n20 = rank + m*(n-1) - (m-1) + m*n; 660 } 661 662 if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */ 663 n8 = rank - m*n +1 - m*n; 664 n17 = rank - m*n +1; 665 n26 = rank - m*n +1 + m*n; 666 } 667 668 if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */ 669 n0 = size - m + rank -1; 670 n1 = size - m + rank; 671 n2 = size - m + rank +1; 672 } 673 674 if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */ 675 n18 = m*n - (size - rank) + m*(n-1) -1; 676 n19 = m*n - (size - rank) + m*(n-1); 677 n20 = m*n - (size - rank) + m*(n-1) +1; 678 } 679 680 if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */ 681 n6 = size - (m*n-rank) - m * (n-1) -1; 682 n7 = size - (m*n-rank) - m * (n-1); 683 n8 = size - (m*n-rank) - m * (n-1) +1; 684 } 685 686 if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */ 687 n24 = rank - (size-m) -1; 688 n25 = rank - (size-m); 689 n26 = rank - (size-m) +1; 690 } 691 692 /* Check for Corners */ 693 if ((xs==0) && (ys==0) && (zs==0)) n0 = size -1; 694 if ((xs==0) && (ys==0) && (ze==P)) n18 = m*n-1; 695 if ((xs==0) && (ye==N) && (zs==0)) n6 = (size-1)-m*(n-1); 696 if ((xs==0) && (ye==N) && (ze==P)) n24 = m-1; 697 if ((xe==M) && (ys==0) && (zs==0)) n2 = size-m; 698 if ((xe==M) && (ys==0) && (ze==P)) n20 = m*n-m; 699 if ((xe==M) && (ye==N) && (zs==0)) n8 = size-m*n; 700 if ((xe==M) && (ye==N) && (ze==P)) n26 = 0; 701 702 /* Check for when not X,Y, and Z Periodic */ 703 704 /* If not X periodic */ 705 if (bx != DM_BOUNDARY_PERIODIC) { 706 if (xs==0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2; 707 if (xe==M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2; 708 } 709 710 /* If not Y periodic */ 711 if (by != DM_BOUNDARY_PERIODIC) { 712 if (ys==0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2; 713 if (ye==N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2; 714 } 715 716 /* If not Z periodic */ 717 if (bz != DM_BOUNDARY_PERIODIC) { 718 if (zs==0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2; 719 if (ze==P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2; 720 } 721 722 ierr = PetscMalloc1(27,&dd->neighbors);CHKERRQ(ierr); 723 724 dd->neighbors[0] = n0; 725 dd->neighbors[1] = n1; 726 dd->neighbors[2] = n2; 727 dd->neighbors[3] = n3; 728 dd->neighbors[4] = n4; 729 dd->neighbors[5] = n5; 730 dd->neighbors[6] = n6; 731 dd->neighbors[7] = n7; 732 dd->neighbors[8] = n8; 733 dd->neighbors[9] = n9; 734 dd->neighbors[10] = n10; 735 dd->neighbors[11] = n11; 736 dd->neighbors[12] = n12; 737 dd->neighbors[13] = rank; 738 dd->neighbors[14] = n14; 739 dd->neighbors[15] = n15; 740 dd->neighbors[16] = n16; 741 dd->neighbors[17] = n17; 742 dd->neighbors[18] = n18; 743 dd->neighbors[19] = n19; 744 dd->neighbors[20] = n20; 745 dd->neighbors[21] = n21; 746 dd->neighbors[22] = n22; 747 dd->neighbors[23] = n23; 748 dd->neighbors[24] = n24; 749 dd->neighbors[25] = n25; 750 dd->neighbors[26] = n26; 751 752 /* If star stencil then delete the corner neighbors */ 753 if (stencil_type == DMDA_STENCIL_STAR) { 754 /* save information about corner neighbors */ 755 sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7; 756 sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18; 757 sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25; 758 sn26 = n26; 759 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; 760 } 761 762 ierr = PetscMalloc1((Xe-Xs)*(Ye-Ys)*(Ze-Zs),&idx);CHKERRQ(ierr); 763 764 nn = 0; 765 /* Bottom Level */ 766 for (k=0; k<s_z; k++) { 767 for (i=1; i<=s_y; i++) { 768 if (n0 >= 0) { /* left below */ 769 x_t = lx[n0 % m]; 770 y_t = ly[(n0 % (m*n))/m]; 771 z_t = lz[n0 / (m*n)]; 772 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; 773 if (twod && (s_t < 0)) s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x; /* 2D case */ 774 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 775 } 776 if (n1 >= 0) { /* directly below */ 777 x_t = x; 778 y_t = ly[(n1 % (m*n))/m]; 779 z_t = lz[n1 / (m*n)]; 780 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 781 if (twod && (s_t < 0)) s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 782 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 783 } 784 if (n2 >= 0) { /* right below */ 785 x_t = lx[n2 % m]; 786 y_t = ly[(n2 % (m*n))/m]; 787 z_t = lz[n2 / (m*n)]; 788 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 789 if (twod && (s_t < 0)) s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t; /* 2D case */ 790 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 791 } 792 } 793 794 for (i=0; i<y; i++) { 795 if (n3 >= 0) { /* directly left */ 796 x_t = lx[n3 % m]; 797 y_t = y; 798 z_t = lz[n3 / (m*n)]; 799 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 800 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 */ 801 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 802 } 803 804 if (n4 >= 0) { /* middle */ 805 x_t = x; 806 y_t = y; 807 z_t = lz[n4 / (m*n)]; 808 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 809 if (twod && (s_t < 0)) s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 810 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 811 } else if (bz == DM_BOUNDARY_MIRROR) { 812 for (j=0; j<x; j++) idx[nn++] = 0; 813 } 814 815 if (n5 >= 0) { /* directly right */ 816 x_t = lx[n5 % m]; 817 y_t = y; 818 z_t = lz[n5 / (m*n)]; 819 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 820 if (twod && (s_t < 0)) s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t; /* 2D case */ 821 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 822 } 823 } 824 825 for (i=1; i<=s_y; i++) { 826 if (n6 >= 0) { /* left above */ 827 x_t = lx[n6 % m]; 828 y_t = ly[(n6 % (m*n))/m]; 829 z_t = lz[n6 / (m*n)]; 830 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 831 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 */ 832 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 833 } 834 if (n7 >= 0) { /* directly above */ 835 x_t = x; 836 y_t = ly[(n7 % (m*n))/m]; 837 z_t = lz[n7 / (m*n)]; 838 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 839 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 */ 840 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 841 } 842 if (n8 >= 0) { /* right above */ 843 x_t = lx[n8 % m]; 844 y_t = ly[(n8 % (m*n))/m]; 845 z_t = lz[n8 / (m*n)]; 846 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 847 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 */ 848 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 849 } 850 } 851 } 852 853 /* Middle Level */ 854 for (k=0; k<z; k++) { 855 for (i=1; i<=s_y; i++) { 856 if (n9 >= 0) { /* left below */ 857 x_t = lx[n9 % m]; 858 y_t = ly[(n9 % (m*n))/m]; 859 /* z_t = z; */ 860 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 861 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 862 } 863 if (n10 >= 0) { /* directly below */ 864 x_t = x; 865 y_t = ly[(n10 % (m*n))/m]; 866 /* z_t = z; */ 867 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 868 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 869 } else if (by == DM_BOUNDARY_MIRROR) { 870 for (j=0; j<x; j++) idx[nn++] = 0; 871 } 872 if (n11 >= 0) { /* right below */ 873 x_t = lx[n11 % m]; 874 y_t = ly[(n11 % (m*n))/m]; 875 /* z_t = z; */ 876 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 877 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 878 } 879 } 880 881 for (i=0; i<y; i++) { 882 if (n12 >= 0) { /* directly left */ 883 x_t = lx[n12 % m]; 884 y_t = y; 885 /* z_t = z; */ 886 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 887 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 888 } else if (bx == DM_BOUNDARY_MIRROR) { 889 for (j=0; j<s_x; j++) idx[nn++] = 0; 890 } 891 892 /* Interior */ 893 s_t = bases[rank] + i*x + k*x*y; 894 for (j=0; j<x; j++) idx[nn++] = s_t++; 895 896 if (n14 >= 0) { /* directly right */ 897 x_t = lx[n14 % m]; 898 y_t = y; 899 /* z_t = z; */ 900 s_t = bases[n14] + i*x_t + k*x_t*y_t; 901 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 902 } else if (bx == DM_BOUNDARY_MIRROR) { 903 for (j=0; j<s_x; j++) idx[nn++] = 0; 904 } 905 } 906 907 for (i=1; i<=s_y; i++) { 908 if (n15 >= 0) { /* left above */ 909 x_t = lx[n15 % m]; 910 y_t = ly[(n15 % (m*n))/m]; 911 /* z_t = z; */ 912 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 913 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 914 } 915 if (n16 >= 0) { /* directly above */ 916 x_t = x; 917 y_t = ly[(n16 % (m*n))/m]; 918 /* z_t = z; */ 919 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 920 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 921 } else if (by == DM_BOUNDARY_MIRROR) { 922 for (j=0; j<x; j++) idx[nn++] = 0; 923 } 924 if (n17 >= 0) { /* right above */ 925 x_t = lx[n17 % m]; 926 y_t = ly[(n17 % (m*n))/m]; 927 /* z_t = z; */ 928 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 929 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 930 } 931 } 932 } 933 934 /* Upper Level */ 935 for (k=0; k<s_z; k++) { 936 for (i=1; i<=s_y; i++) { 937 if (n18 >= 0) { /* left below */ 938 x_t = lx[n18 % m]; 939 y_t = ly[(n18 % (m*n))/m]; 940 /* z_t = lz[n18 / (m*n)]; */ 941 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 942 if (twod && (s_t >= M*N*P)) s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t; /* 2d case */ 943 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 944 } 945 if (n19 >= 0) { /* directly below */ 946 x_t = x; 947 y_t = ly[(n19 % (m*n))/m]; 948 /* z_t = lz[n19 / (m*n)]; */ 949 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 950 if (twod && (s_t >= M*N*P)) s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 951 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 952 } 953 if (n20 >= 0) { /* right below */ 954 x_t = lx[n20 % m]; 955 y_t = ly[(n20 % (m*n))/m]; 956 /* z_t = lz[n20 / (m*n)]; */ 957 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 958 if (twod && (s_t >= M*N*P)) s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t; /* 2d case */ 959 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 960 } 961 } 962 963 for (i=0; i<y; i++) { 964 if (n21 >= 0) { /* directly left */ 965 x_t = lx[n21 % m]; 966 y_t = y; 967 /* z_t = lz[n21 / (m*n)]; */ 968 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 969 if (twod && (s_t >= M*N*P)) s_t = bases[n21] + (i+1)*x_t - s_x; /* 2d case */ 970 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 971 } 972 973 if (n22 >= 0) { /* middle */ 974 x_t = x; 975 y_t = y; 976 /* z_t = lz[n22 / (m*n)]; */ 977 s_t = bases[n22] + i*x_t + k*x_t*y_t; 978 if (twod && (s_t >= M*N*P)) s_t = bases[n22] + i*x_t; /* 2d case */ 979 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 980 } else if (bz == DM_BOUNDARY_MIRROR) { 981 for (j=0; j<x; j++) idx[nn++] = 0; 982 } 983 984 if (n23 >= 0) { /* directly right */ 985 x_t = lx[n23 % m]; 986 y_t = y; 987 /* z_t = lz[n23 / (m*n)]; */ 988 s_t = bases[n23] + i*x_t + k*x_t*y_t; 989 if (twod && (s_t >= M*N*P)) s_t = bases[n23] + i*x_t; /* 2d case */ 990 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 991 } 992 } 993 994 for (i=1; i<=s_y; i++) { 995 if (n24 >= 0) { /* left above */ 996 x_t = lx[n24 % m]; 997 y_t = ly[(n24 % (m*n))/m]; 998 /* z_t = lz[n24 / (m*n)]; */ 999 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1000 if (twod && (s_t >= M*N*P)) s_t = bases[n24] + i*x_t - s_x; /* 2d case */ 1001 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1002 } 1003 if (n25 >= 0) { /* directly above */ 1004 x_t = x; 1005 y_t = ly[(n25 % (m*n))/m]; 1006 /* z_t = lz[n25 / (m*n)]; */ 1007 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1008 if (twod && (s_t >= M*N*P)) s_t = bases[n25] + (i-1)*x_t; /* 2d case */ 1009 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1010 } 1011 if (n26 >= 0) { /* right above */ 1012 x_t = lx[n26 % m]; 1013 y_t = ly[(n26 % (m*n))/m]; 1014 /* z_t = lz[n26 / (m*n)]; */ 1015 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1016 if (twod && (s_t >= M*N*P)) s_t = bases[n26] + (i-1)*x_t; /* 2d case */ 1017 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1018 } 1019 } 1020 } 1021 1022 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 1023 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1024 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);CHKERRQ(ierr); 1025 ierr = ISDestroy(&to);CHKERRQ(ierr); 1026 ierr = ISDestroy(&from);CHKERRQ(ierr); 1027 1028 if (stencil_type == DMDA_STENCIL_STAR) { 1029 n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; 1030 n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; 1031 n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; 1032 n26 = sn26; 1033 } 1034 1035 if (((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz))) { 1036 /* 1037 Recompute the local to global mappings, this time keeping the 1038 information about the cross corner processor numbers. 1039 */ 1040 nn = 0; 1041 /* Bottom Level */ 1042 for (k=0; k<s_z; k++) { 1043 for (i=1; i<=s_y; i++) { 1044 if (n0 >= 0) { /* left below */ 1045 x_t = lx[n0 % m]; 1046 y_t = ly[(n0 % (m*n))/m]; 1047 z_t = lz[n0 / (m*n)]; 1048 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; 1049 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1050 } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) { 1051 for (j=0; j<s_x; j++) idx[nn++] = -1; 1052 } 1053 if (n1 >= 0) { /* directly below */ 1054 x_t = x; 1055 y_t = ly[(n1 % (m*n))/m]; 1056 z_t = lz[n1 / (m*n)]; 1057 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1058 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1059 } else if (Ys-ys < 0 && Zs-zs < 0) { 1060 for (j=0; j<x; j++) idx[nn++] = -1; 1061 } 1062 if (n2 >= 0) { /* right below */ 1063 x_t = lx[n2 % m]; 1064 y_t = ly[(n2 % (m*n))/m]; 1065 z_t = lz[n2 / (m*n)]; 1066 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1067 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1068 } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) { 1069 for (j=0; j<s_x; j++) idx[nn++] = -1; 1070 } 1071 } 1072 1073 for (i=0; i<y; i++) { 1074 if (n3 >= 0) { /* directly left */ 1075 x_t = lx[n3 % m]; 1076 y_t = y; 1077 z_t = lz[n3 / (m*n)]; 1078 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1079 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1080 } else if (Xs-xs < 0 && Zs-zs < 0) { 1081 for (j=0; j<s_x; j++) idx[nn++] = -1; 1082 } 1083 1084 if (n4 >= 0) { /* middle */ 1085 x_t = x; 1086 y_t = y; 1087 z_t = lz[n4 / (m*n)]; 1088 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1089 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1090 } else if (Zs-zs < 0) { 1091 if (bz == DM_BOUNDARY_MIRROR) { 1092 for (j=0; j<x; j++) idx[nn++] = 0; 1093 } else { 1094 for (j=0; j<x; j++) idx[nn++] = -1; 1095 } 1096 } 1097 1098 if (n5 >= 0) { /* directly right */ 1099 x_t = lx[n5 % m]; 1100 y_t = y; 1101 z_t = lz[n5 / (m*n)]; 1102 s_t = bases[n5] + i*x_t + 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 (xe-Xe < 0 && Zs-zs < 0) { 1105 for (j=0; j<s_x; j++) idx[nn++] = -1; 1106 } 1107 } 1108 1109 for (i=1; i<=s_y; i++) { 1110 if (n6 >= 0) { /* left above */ 1111 x_t = lx[n6 % m]; 1112 y_t = ly[(n6 % (m*n))/m]; 1113 z_t = lz[n6 / (m*n)]; 1114 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1115 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1116 } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) { 1117 for (j=0; j<s_x; j++) idx[nn++] = -1; 1118 } 1119 if (n7 >= 0) { /* directly above */ 1120 x_t = x; 1121 y_t = ly[(n7 % (m*n))/m]; 1122 z_t = lz[n7 / (m*n)]; 1123 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1124 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1125 } else if (ye-Ye < 0 && Zs-zs < 0) { 1126 for (j=0; j<x; j++) idx[nn++] = -1; 1127 } 1128 if (n8 >= 0) { /* right above */ 1129 x_t = lx[n8 % m]; 1130 y_t = ly[(n8 % (m*n))/m]; 1131 z_t = lz[n8 / (m*n)]; 1132 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1133 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1134 } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) { 1135 for (j=0; j<s_x; j++) idx[nn++] = -1; 1136 } 1137 } 1138 } 1139 1140 /* Middle Level */ 1141 for (k=0; k<z; k++) { 1142 for (i=1; i<=s_y; i++) { 1143 if (n9 >= 0) { /* left below */ 1144 x_t = lx[n9 % m]; 1145 y_t = ly[(n9 % (m*n))/m]; 1146 /* z_t = z; */ 1147 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1148 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1149 } else if (Xs-xs < 0 && Ys-ys < 0) { 1150 for (j=0; j<s_x; j++) idx[nn++] = -1; 1151 } 1152 if (n10 >= 0) { /* directly below */ 1153 x_t = x; 1154 y_t = ly[(n10 % (m*n))/m]; 1155 /* z_t = z; */ 1156 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1157 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1158 } else if (Ys-ys < 0) { 1159 if (by == DM_BOUNDARY_MIRROR) { 1160 for (j=0; j<x; j++) idx[nn++] = -1; 1161 } else { 1162 for (j=0; j<x; j++) idx[nn++] = -1; 1163 } 1164 } 1165 if (n11 >= 0) { /* right below */ 1166 x_t = lx[n11 % m]; 1167 y_t = ly[(n11 % (m*n))/m]; 1168 /* z_t = z; */ 1169 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1170 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1171 } else if (xe-Xe < 0 && Ys-ys < 0) { 1172 for (j=0; j<s_x; j++) idx[nn++] = -1; 1173 } 1174 } 1175 1176 for (i=0; i<y; i++) { 1177 if (n12 >= 0) { /* directly left */ 1178 x_t = lx[n12 % m]; 1179 y_t = y; 1180 /* z_t = z; */ 1181 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 1182 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1183 } else if (Xs-xs < 0) { 1184 if (bx == DM_BOUNDARY_MIRROR) { 1185 for (j=0; j<s_x; j++) idx[nn++] = 0; 1186 } else { 1187 for (j=0; j<s_x; j++) idx[nn++] = -1; 1188 } 1189 } 1190 1191 /* Interior */ 1192 s_t = bases[rank] + i*x + k*x*y; 1193 for (j=0; j<x; j++) idx[nn++] = s_t++; 1194 1195 if (n14 >= 0) { /* directly right */ 1196 x_t = lx[n14 % m]; 1197 y_t = y; 1198 /* z_t = z; */ 1199 s_t = bases[n14] + i*x_t + k*x_t*y_t; 1200 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1201 } else if (xe-Xe < 0) { 1202 if (bx == DM_BOUNDARY_MIRROR) { 1203 for (j=0; j<s_x; j++) idx[nn++] = 0; 1204 } else { 1205 for (j=0; j<s_x; j++) idx[nn++] = -1; 1206 } 1207 } 1208 } 1209 1210 for (i=1; i<=s_y; i++) { 1211 if (n15 >= 0) { /* left above */ 1212 x_t = lx[n15 % m]; 1213 y_t = ly[(n15 % (m*n))/m]; 1214 /* z_t = z; */ 1215 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 1216 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1217 } else if (Xs-xs < 0 && ye-Ye < 0) { 1218 for (j=0; j<s_x; j++) idx[nn++] = -1; 1219 } 1220 if (n16 >= 0) { /* directly above */ 1221 x_t = x; 1222 y_t = ly[(n16 % (m*n))/m]; 1223 /* z_t = z; */ 1224 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 1225 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1226 } else if (ye-Ye < 0) { 1227 if (by == DM_BOUNDARY_MIRROR) { 1228 for (j=0; j<x; j++) idx[nn++] = 0; 1229 } else { 1230 for (j=0; j<x; j++) idx[nn++] = -1; 1231 } 1232 } 1233 if (n17 >= 0) { /* right above */ 1234 x_t = lx[n17 % m]; 1235 y_t = ly[(n17 % (m*n))/m]; 1236 /* z_t = z; */ 1237 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 1238 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1239 } else if (xe-Xe < 0 && ye-Ye < 0) { 1240 for (j=0; j<s_x; j++) idx[nn++] = -1; 1241 } 1242 } 1243 } 1244 1245 /* Upper Level */ 1246 for (k=0; k<s_z; k++) { 1247 for (i=1; i<=s_y; i++) { 1248 if (n18 >= 0) { /* left below */ 1249 x_t = lx[n18 % m]; 1250 y_t = ly[(n18 % (m*n))/m]; 1251 /* z_t = lz[n18 / (m*n)]; */ 1252 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1253 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1254 } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) { 1255 for (j=0; j<s_x; j++) idx[nn++] = -1; 1256 } 1257 if (n19 >= 0) { /* directly below */ 1258 x_t = x; 1259 y_t = ly[(n19 % (m*n))/m]; 1260 /* z_t = lz[n19 / (m*n)]; */ 1261 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1262 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1263 } else if (Ys-ys < 0 && ze-Ze < 0) { 1264 for (j=0; j<x; j++) idx[nn++] = -1; 1265 } 1266 if (n20 >= 0) { /* right below */ 1267 x_t = lx[n20 % m]; 1268 y_t = ly[(n20 % (m*n))/m]; 1269 /* z_t = lz[n20 / (m*n)]; */ 1270 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1271 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1272 } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) { 1273 for (j=0; j<s_x; j++) idx[nn++] = -1; 1274 } 1275 } 1276 1277 for (i=0; i<y; i++) { 1278 if (n21 >= 0) { /* directly left */ 1279 x_t = lx[n21 % m]; 1280 y_t = y; 1281 /* z_t = lz[n21 / (m*n)]; */ 1282 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1283 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1284 } else if (Xs-xs < 0 && ze-Ze < 0) { 1285 for (j=0; j<s_x; j++) idx[nn++] = -1; 1286 } 1287 1288 if (n22 >= 0) { /* middle */ 1289 x_t = x; 1290 y_t = y; 1291 /* z_t = lz[n22 / (m*n)]; */ 1292 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1293 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1294 } else if (ze-Ze < 0) { 1295 if (bz == DM_BOUNDARY_MIRROR) { 1296 for (j=0; j<x; j++) idx[nn++] = 0; 1297 } else { 1298 for (j=0; j<x; j++) idx[nn++] = -1; 1299 } 1300 } 1301 1302 if (n23 >= 0) { /* directly right */ 1303 x_t = lx[n23 % m]; 1304 y_t = y; 1305 /* z_t = lz[n23 / (m*n)]; */ 1306 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1307 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1308 } else if (xe-Xe < 0 && ze-Ze < 0) { 1309 for (j=0; j<s_x; j++) idx[nn++] = -1; 1310 } 1311 } 1312 1313 for (i=1; i<=s_y; i++) { 1314 if (n24 >= 0) { /* left above */ 1315 x_t = lx[n24 % m]; 1316 y_t = ly[(n24 % (m*n))/m]; 1317 /* z_t = lz[n24 / (m*n)]; */ 1318 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1319 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1320 } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) { 1321 for (j=0; j<s_x; j++) idx[nn++] = -1; 1322 } 1323 if (n25 >= 0) { /* directly above */ 1324 x_t = x; 1325 y_t = ly[(n25 % (m*n))/m]; 1326 /* z_t = lz[n25 / (m*n)]; */ 1327 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1328 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 1329 } else if (ye-Ye < 0 && ze-Ze < 0) { 1330 for (j=0; j<x; j++) idx[nn++] = -1; 1331 } 1332 if (n26 >= 0) { /* right above */ 1333 x_t = lx[n26 % m]; 1334 y_t = ly[(n26 % (m*n))/m]; 1335 /* z_t = lz[n26 / (m*n)]; */ 1336 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1337 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 1338 } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) { 1339 for (j=0; j<s_x; j++) idx[nn++] = -1; 1340 } 1341 } 1342 } 1343 } 1344 /* 1345 Set the local to global ordering in the global vector, this allows use 1346 of VecSetValuesLocal(). 1347 */ 1348 ierr = ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr); 1349 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);CHKERRQ(ierr); 1350 1351 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1352 dd->m = m; dd->n = n; dd->p = p; 1353 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1354 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; 1355 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; 1356 1357 ierr = VecDestroy(&local);CHKERRQ(ierr); 1358 ierr = VecDestroy(&global);CHKERRQ(ierr); 1359 1360 dd->gtol = gtol; 1361 dd->base = base; 1362 da->ops->view = DMView_DA_3d; 1363 dd->ltol = NULL; 1364 dd->ao = NULL; 1365 PetscFunctionReturn(0); 1366 } 1367 1368 1369 /*@C 1370 DMDACreate3d - Creates an object that will manage the communication of three-dimensional 1371 regular array data that is distributed across some processors. 1372 1373 Collective on MPI_Comm 1374 1375 Input Parameters: 1376 + comm - MPI communicator 1377 . bx,by,bz - type of ghost nodes the array have. 1378 Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC. 1379 . stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX) 1380 . M,N,P - global dimension in each direction of the array 1381 . m,n,p - corresponding number of processors in each dimension 1382 (or PETSC_DECIDE to have calculated) 1383 . dof - number of degrees of freedom per node 1384 . s - stencil width 1385 - lx, ly, lz - arrays containing the number of nodes in each cell along 1386 the x, y, and z coordinates, or NULL. If non-null, these 1387 must be of length as m,n,p and the corresponding 1388 m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of 1389 the ly[] must N, sum of the lz[] must be P 1390 1391 Output Parameter: 1392 . da - the resulting distributed array object 1393 1394 Options Database Key: 1395 + -dm_view - Calls DMView() at the conclusion of DMDACreate3d() 1396 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 1397 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 1398 . -da_grid_z <nz> - number of grid points in z direction, if P < 0 1399 . -da_processors_x <MX> - number of processors in x direction 1400 . -da_processors_y <MY> - number of processors in y direction 1401 . -da_processors_z <MZ> - number of processors in z direction 1402 . -da_refine_x <rx> - refinement ratio in x direction 1403 . -da_refine_y <ry> - refinement ratio in y direction 1404 . -da_refine_z <rz>- refinement ratio in z directio 1405 - -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0 1406 1407 Level: beginner 1408 1409 Notes: 1410 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 1411 standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 1412 the standard 27-pt stencil. 1413 1414 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 1415 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 1416 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 1417 1418 You must call DMSetUp() after this call before using this DM. 1419 1420 If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call 1421 but before DMSetUp(). 1422 1423 .keywords: distributed array, create, three-dimensional 1424 1425 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 1426 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(), 1427 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 1428 1429 @*/ 1430 PetscErrorCode DMDACreate3d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz,DMDAStencilType stencil_type,PetscInt M, 1431 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) 1432 { 1433 PetscErrorCode ierr; 1434 1435 PetscFunctionBegin; 1436 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 1437 ierr = DMSetDimension(*da, 3);CHKERRQ(ierr); 1438 ierr = DMDASetSizes(*da, M, N, P);CHKERRQ(ierr); 1439 ierr = DMDASetNumProcs(*da, m, n, p);CHKERRQ(ierr); 1440 ierr = DMDASetBoundaryType(*da, bx, by, bz);CHKERRQ(ierr); 1441 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 1442 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 1443 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 1444 ierr = DMDASetOwnershipRanges(*da, lx, ly, lz);CHKERRQ(ierr); 1445 PetscFunctionReturn(0); 1446 } 1447