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