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