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