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