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