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