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