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