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 defined(foo) 747 if (s_t < 0) {s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x;} /* 2D case */ 748 #endif 749 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 750 } 751 if (n1 >= 0) { /* directly below */ 752 x_t = x; 753 y_t = ly[(n1 % (m*n))/m]; 754 z_t = lz[n1 / (m*n)]; 755 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 756 #if defined(foo) 757 if (s_t < 0) {s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t;} /* 2D case */ 758 #endif 759 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 760 } 761 if (n2 >= 0) { /* right below */ 762 x_t = lx[n2 % m]; 763 y_t = ly[(n2 % (m*n))/m]; 764 z_t = lz[n2 / (m*n)]; 765 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 766 #if defined(foo) 767 if (s_t < 0) {s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t;} /* 2D case */ 768 #endif 769 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 770 } 771 } 772 773 for (i=0; i<y; i++) { 774 if (n3 >= 0) { /* directly left */ 775 x_t = lx[n3 % m]; 776 y_t = y; 777 z_t = lz[n3 / (m*n)]; 778 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 779 #if defined(foo) 780 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 */ 781 #endif 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 defined(foo) 791 if (s_t < 0) {s_t = bases[n4] + i*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 792 #endif 793 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 794 } 795 796 if (n5 >= 0) { /* directly right */ 797 x_t = lx[n5 % m]; 798 y_t = y; 799 z_t = lz[n5 / (m*n)]; 800 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 801 #if defined(foo) 802 if (s_t < 0) {s_t = bases[n5] + i*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 803 #endif 804 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 805 } 806 } 807 808 for (i=1; i<=s_y; i++) { 809 if (n6 >= 0) { /* left above */ 810 x_t = lx[n6 % m]; 811 y_t = ly[(n6 % (m*n))/m]; 812 z_t = lz[n6 / (m*n)]; 813 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 814 #if defined(foo) 815 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 */ 816 #endif 817 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 818 } 819 if (n7 >= 0) { /* directly above */ 820 x_t = x; 821 y_t = ly[(n7 % (m*n))/m]; 822 z_t = lz[n7 / (m*n)]; 823 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 824 #if defined(foo) 825 if (s_t < 0) {s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 826 #endif 827 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 828 } 829 if (n8 >= 0) { /* right above */ 830 x_t = lx[n8 % m]; 831 y_t = ly[(n8 % (m*n))/m]; 832 z_t = lz[n8 / (m*n)]; 833 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 834 #if defined(foo) 835 if (s_t < 0) {s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - x_t*y_t;} /* 2D case */ 836 #endif 837 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 838 } 839 } 840 } 841 842 /* Middle Level */ 843 for (k=0; k<z; k++) { 844 for (i=1; i<=s_y; i++) { 845 if (n9 >= 0) { /* left below */ 846 x_t = lx[n9 % m]; 847 y_t = ly[(n9 % (m*n))/m]; 848 /* z_t = z; */ 849 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 850 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 851 } 852 if (n10 >= 0) { /* directly below */ 853 x_t = x; 854 y_t = ly[(n10 % (m*n))/m]; 855 /* z_t = z; */ 856 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 857 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 858 } 859 if (n11 >= 0) { /* right below */ 860 x_t = lx[n11 % m]; 861 y_t = ly[(n11 % (m*n))/m]; 862 /* z_t = z; */ 863 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 864 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 865 } 866 } 867 868 for (i=0; i<y; i++) { 869 if (n12 >= 0) { /* directly left */ 870 x_t = lx[n12 % m]; 871 y_t = y; 872 /* z_t = z; */ 873 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 874 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 875 } 876 877 /* Interior */ 878 s_t = bases[rank] + i*x + k*x*y; 879 for (j=0; j<x; j++) { idx[nn++] = s_t++;} 880 881 if (n14 >= 0) { /* directly right */ 882 x_t = lx[n14 % m]; 883 y_t = y; 884 /* z_t = z; */ 885 s_t = bases[n14] + i*x_t + k*x_t*y_t; 886 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 887 } 888 } 889 890 for (i=1; i<=s_y; i++) { 891 if (n15 >= 0) { /* left above */ 892 x_t = lx[n15 % m]; 893 y_t = ly[(n15 % (m*n))/m]; 894 /* z_t = z; */ 895 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 896 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 897 } 898 if (n16 >= 0) { /* directly above */ 899 x_t = x; 900 y_t = ly[(n16 % (m*n))/m]; 901 /* z_t = z; */ 902 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 903 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 904 } 905 if (n17 >= 0) { /* right above */ 906 x_t = lx[n17 % m]; 907 y_t = ly[(n17 % (m*n))/m]; 908 /* z_t = z; */ 909 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 910 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 911 } 912 } 913 } 914 915 /* Upper Level */ 916 for (k=0; k<s_z; k++) { 917 for (i=1; i<=s_y; i++) { 918 if (n18 >= 0) { /* left below */ 919 x_t = lx[n18 % m]; 920 y_t = ly[(n18 % (m*n))/m]; 921 /* z_t = lz[n18 / (m*n)]; */ 922 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 923 #if defined(foo) 924 if (s_t >= x*y*z) {s_t = bases[n18] - (s_y-i)*x_t -s_x + x_t*y_t;} /* 2d case */ 925 #endif 926 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 927 } 928 if (n19 >= 0) { /* directly below */ 929 x_t = x; 930 y_t = ly[(n19 % (m*n))/m]; 931 /* z_t = lz[n19 / (m*n)]; */ 932 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 933 #if defined(foo) 934 if (s_t >= x*y*z) {s_t = bases[n19] - (s_y+1-i)*x_t + x_t*y_t;} /* 2d case */ 935 #endif 936 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 937 } 938 if (n20 >= 0) { /* right below */ 939 x_t = lx[n20 % m]; 940 y_t = ly[(n20 % (m*n))/m]; 941 /* z_t = lz[n20 / (m*n)]; */ 942 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 943 #if defined(foo) 944 if (s_t >= x*y*z) {s_t = bases[n20] - (s_y+1-i)*x_t + x_t*y_t;} /* 2d case */ 945 #endif 946 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 947 } 948 } 949 950 for (i=0; i<y; i++) { 951 if (n21 >= 0) { /* directly left */ 952 x_t = lx[n21 % m]; 953 y_t = y; 954 /* z_t = lz[n21 / (m*n)]; */ 955 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 956 #if defined(foo) 957 if (s_t >= x*y*z) {s_t = bases[n21] + (i+1)*x_t - s_x;} /* 2d case */ 958 #endif 959 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 960 } 961 962 if (n22 >= 0) { /* middle */ 963 x_t = x; 964 y_t = y; 965 /* z_t = lz[n22 / (m*n)]; */ 966 s_t = bases[n22] + i*x_t + k*x_t*y_t; 967 #if defined(foo) 968 if (s_t >= x*y*z) {s_t = bases[n22] + i*x_t;} /* 2d case */ 969 #endif 970 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 971 } 972 973 if (n23 >= 0) { /* directly right */ 974 x_t = lx[n23 % m]; 975 y_t = y; 976 /* z_t = lz[n23 / (m*n)]; */ 977 s_t = bases[n23] + i*x_t + k*x_t*y_t; 978 #if defined(foo) 979 if (s_t >= x*y*z) {s_t = bases[n23] + i*x_t;} /* 2d case */ 980 #endif 981 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 982 } 983 } 984 985 for (i=1; i<=s_y; i++) { 986 if (n24 >= 0) { /* left above */ 987 x_t = lx[n24 % m]; 988 y_t = ly[(n24 % (m*n))/m]; 989 /* z_t = lz[n24 / (m*n)]; */ 990 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 991 #if defined(foo) 992 if (s_t >= x*y*z) {s_t = bases[n24] + i*x_t - s_x;} /* 2d case */ 993 #endif 994 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 995 } 996 if (n25 >= 0) { /* directly above */ 997 x_t = x; 998 y_t = ly[(n25 % (m*n))/m]; 999 /* z_t = lz[n25 / (m*n)]; */ 1000 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1001 #if defined(foo) 1002 if (s_t >= x*y*z) {s_t = bases[n25] + (i-1)*x_t;} /* 2d case */ 1003 #endif 1004 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1005 } 1006 if (n26 >= 0) { /* right above */ 1007 x_t = lx[n26 % m]; 1008 y_t = ly[(n26 % (m*n))/m]; 1009 /* z_t = lz[n26 / (m*n)]; */ 1010 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1011 #if defined(foo) 1012 if (s_t >= x*y*z) {s_t = bases[n26] + (i-1)*x_t;} /* 2d case */ 1013 #endif 1014 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1015 } 1016 } 1017 } 1018 1019 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 1020 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1021 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 1022 ierr = ISDestroy(&to);CHKERRQ(ierr); 1023 ierr = ISDestroy(&from);CHKERRQ(ierr); 1024 1025 if (stencil_type == DMDA_STENCIL_STAR) { 1026 n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; 1027 n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; 1028 n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; 1029 n26 = sn26; 1030 } 1031 1032 if ((stencil_type == DMDA_STENCIL_STAR) || 1033 (bx != DMDA_BOUNDARY_PERIODIC && bx) || 1034 (by != DMDA_BOUNDARY_PERIODIC && by) || 1035 (bz != DMDA_BOUNDARY_PERIODIC && bz)) { 1036 /* 1037 Recompute the local to global mappings, this time keeping the 1038 information about the cross corner processor numbers. 1039 */ 1040 nn = 0; 1041 /* Bottom Level */ 1042 for (k=0; k<s_z; k++) { 1043 for (i=1; i<=s_y; i++) { 1044 if (n0 >= 0) { /* left below */ 1045 x_t = lx[n0 % m]; 1046 y_t = ly[(n0 % (m*n))/m]; 1047 z_t = lz[n0 / (m*n)]; 1048 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; 1049 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1050 } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) { 1051 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1052 } 1053 if (n1 >= 0) { /* directly below */ 1054 x_t = x; 1055 y_t = ly[(n1 % (m*n))/m]; 1056 z_t = lz[n1 / (m*n)]; 1057 s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1058 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1059 } else if (Ys-ys < 0 && Zs-zs < 0) { 1060 for (j=0; j<x; j++) { idx[nn++] = -1;} 1061 } 1062 if (n2 >= 0) { /* right below */ 1063 x_t = lx[n2 % m]; 1064 y_t = ly[(n2 % (m*n))/m]; 1065 z_t = lz[n2 / (m*n)]; 1066 s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t; 1067 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1068 } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) { 1069 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1070 } 1071 } 1072 1073 for (i=0; i<y; i++) { 1074 if (n3 >= 0) { /* directly left */ 1075 x_t = lx[n3 % m]; 1076 y_t = y; 1077 z_t = lz[n3 / (m*n)]; 1078 s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1079 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1080 } else if (Xs-xs < 0 && Zs-zs < 0) { 1081 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1082 } 1083 1084 if (n4 >= 0) { /* middle */ 1085 x_t = x; 1086 y_t = y; 1087 z_t = lz[n4 / (m*n)]; 1088 s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1089 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1090 } else if (Zs-zs < 0) { 1091 for (j=0; j<x; j++) { idx[nn++] = -1;} 1092 } 1093 1094 if (n5 >= 0) { /* directly right */ 1095 x_t = lx[n5 % m]; 1096 y_t = y; 1097 z_t = lz[n5 / (m*n)]; 1098 s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1099 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1100 } else if (xe-Xe < 0 && Zs-zs < 0) { 1101 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1102 } 1103 } 1104 1105 for (i=1; i<=s_y; i++) { 1106 if (n6 >= 0) { /* left above */ 1107 x_t = lx[n6 % m]; 1108 y_t = ly[(n6 % (m*n))/m]; 1109 z_t = lz[n6 / (m*n)]; 1110 s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1111 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1112 } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) { 1113 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1114 } 1115 if (n7 >= 0) { /* directly above */ 1116 x_t = x; 1117 y_t = ly[(n7 % (m*n))/m]; 1118 z_t = lz[n7 / (m*n)]; 1119 s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1120 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1121 } else if (ye-Ye < 0 && Zs-zs < 0) { 1122 for (j=0; j<x; j++) { idx[nn++] = -1;} 1123 } 1124 if (n8 >= 0) { /* right above */ 1125 x_t = lx[n8 % m]; 1126 y_t = ly[(n8 % (m*n))/m]; 1127 z_t = lz[n8 / (m*n)]; 1128 s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t; 1129 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1130 } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) { 1131 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1132 } 1133 } 1134 } 1135 1136 /* Middle Level */ 1137 for (k=0; k<z; k++) { 1138 for (i=1; i<=s_y; i++) { 1139 if (n9 >= 0) { /* left below */ 1140 x_t = lx[n9 % m]; 1141 y_t = ly[(n9 % (m*n))/m]; 1142 /* z_t = z; */ 1143 s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1144 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1145 } else if (Xs-xs < 0 && Ys-ys < 0) { 1146 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1147 } 1148 if (n10 >= 0) { /* directly below */ 1149 x_t = x; 1150 y_t = ly[(n10 % (m*n))/m]; 1151 /* z_t = z; */ 1152 s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1153 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1154 } else if (Ys-ys < 0) { 1155 for (j=0; j<x; j++) { idx[nn++] = -1;} 1156 } 1157 if (n11 >= 0) { /* right below */ 1158 x_t = lx[n11 % m]; 1159 y_t = ly[(n11 % (m*n))/m]; 1160 /* z_t = z; */ 1161 s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1162 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1163 } else if (xe-Xe < 0 && Ys-ys < 0) { 1164 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1165 } 1166 } 1167 1168 for (i=0; i<y; i++) { 1169 if (n12 >= 0) { /* directly left */ 1170 x_t = lx[n12 % m]; 1171 y_t = y; 1172 /* z_t = z; */ 1173 s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t; 1174 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1175 } else if (Xs-xs < 0) { 1176 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1177 } 1178 1179 /* Interior */ 1180 s_t = bases[rank] + i*x + k*x*y; 1181 for (j=0; j<x; j++) { idx[nn++] = s_t++;} 1182 1183 if (n14 >= 0) { /* directly right */ 1184 x_t = lx[n14 % m]; 1185 y_t = y; 1186 /* z_t = z; */ 1187 s_t = bases[n14] + i*x_t + k*x_t*y_t; 1188 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1189 } else if (xe-Xe < 0) { 1190 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1191 } 1192 } 1193 1194 for (i=1; i<=s_y; i++) { 1195 if (n15 >= 0) { /* left above */ 1196 x_t = lx[n15 % m]; 1197 y_t = ly[(n15 % (m*n))/m]; 1198 /* z_t = z; */ 1199 s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t; 1200 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1201 } else if (Xs-xs < 0 && ye-Ye < 0) { 1202 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1203 } 1204 if (n16 >= 0) { /* directly above */ 1205 x_t = x; 1206 y_t = ly[(n16 % (m*n))/m]; 1207 /* z_t = z; */ 1208 s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t; 1209 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1210 } else if (ye-Ye < 0) { 1211 for (j=0; j<x; j++) { idx[nn++] = -1;} 1212 } 1213 if (n17 >= 0) { /* right above */ 1214 x_t = lx[n17 % m]; 1215 y_t = ly[(n17 % (m*n))/m]; 1216 /* z_t = z; */ 1217 s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t; 1218 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1219 } else if (xe-Xe < 0 && ye-Ye < 0) { 1220 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1221 } 1222 } 1223 } 1224 1225 /* Upper Level */ 1226 for (k=0; k<s_z; k++) { 1227 for (i=1; i<=s_y; i++) { 1228 if (n18 >= 0) { /* left below */ 1229 x_t = lx[n18 % m]; 1230 y_t = ly[(n18 % (m*n))/m]; 1231 /* z_t = lz[n18 / (m*n)]; */ 1232 s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t; 1233 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1234 } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) { 1235 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1236 } 1237 if (n19 >= 0) { /* directly below */ 1238 x_t = x; 1239 y_t = ly[(n19 % (m*n))/m]; 1240 /* z_t = lz[n19 / (m*n)]; */ 1241 s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1242 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1243 } else if (Ys-ys < 0 && ze-Ze < 0) { 1244 for (j=0; j<x; j++) { idx[nn++] = -1;} 1245 } 1246 if (n20 >= 0) { /* right below */ 1247 x_t = lx[n20 % m]; 1248 y_t = ly[(n20 % (m*n))/m]; 1249 /* z_t = lz[n20 / (m*n)]; */ 1250 s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t; 1251 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1252 } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) { 1253 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1254 } 1255 } 1256 1257 for (i=0; i<y; i++) { 1258 if (n21 >= 0) { /* directly left */ 1259 x_t = lx[n21 % m]; 1260 y_t = y; 1261 /* z_t = lz[n21 / (m*n)]; */ 1262 s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t; 1263 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1264 } else if (Xs-xs < 0 && ze-Ze < 0) { 1265 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1266 } 1267 1268 if (n22 >= 0) { /* middle */ 1269 x_t = x; 1270 y_t = y; 1271 /* z_t = lz[n22 / (m*n)]; */ 1272 s_t = bases[n22] + i*x_t + k*x_t*y_t; 1273 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1274 } else if (ze-Ze < 0) { 1275 for (j=0; j<x; j++) { idx[nn++] = -1;} 1276 } 1277 1278 if (n23 >= 0) { /* directly right */ 1279 x_t = lx[n23 % m]; 1280 y_t = y; 1281 /* z_t = lz[n23 / (m*n)]; */ 1282 s_t = bases[n23] + i*x_t + k*x_t*y_t; 1283 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1284 } else if (xe-Xe < 0 && ze-Ze < 0) { 1285 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1286 } 1287 } 1288 1289 for (i=1; i<=s_y; i++) { 1290 if (n24 >= 0) { /* left above */ 1291 x_t = lx[n24 % m]; 1292 y_t = ly[(n24 % (m*n))/m]; 1293 /* z_t = lz[n24 / (m*n)]; */ 1294 s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t; 1295 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1296 } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) { 1297 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1298 } 1299 if (n25 >= 0) { /* directly above */ 1300 x_t = x; 1301 y_t = ly[(n25 % (m*n))/m]; 1302 /* z_t = lz[n25 / (m*n)]; */ 1303 s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t; 1304 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1305 } else if (ye-Ye < 0 && ze-Ze < 0) { 1306 for (j=0; j<x; j++) { idx[nn++] = -1;} 1307 } 1308 if (n26 >= 0) { /* right above */ 1309 x_t = lx[n26 % m]; 1310 y_t = ly[(n26 % (m*n))/m]; 1311 /* z_t = lz[n26 / (m*n)]; */ 1312 s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t; 1313 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1314 } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) { 1315 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1316 } 1317 } 1318 } 1319 } 1320 /* 1321 Set the local to global ordering in the global vector, this allows use 1322 of VecSetValuesLocal(). 1323 */ 1324 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr); 1325 ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr); 1326 ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1327 ierr = ISGetIndices(ltogis, &idx_full); 1328 ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1329 ierr = ISRestoreIndices(ltogis, &idx_full); 1330 ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr); 1331 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1332 ierr = ISDestroy(<ogis);CHKERRQ(ierr); 1333 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 1334 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1335 1336 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1337 dd->m = m; dd->n = n; dd->p = p; 1338 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1339 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; 1340 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; 1341 1342 ierr = VecDestroy(&local);CHKERRQ(ierr); 1343 ierr = VecDestroy(&global);CHKERRQ(ierr); 1344 1345 dd->gtol = gtol; 1346 dd->ltog = ltog; 1347 dd->idx = idx_cpy; 1348 dd->Nl = nn*dof; 1349 dd->base = base; 1350 da->ops->view = DMView_DA_3d; 1351 dd->ltol = PETSC_NULL; 1352 dd->ao = PETSC_NULL; 1353 1354 PetscFunctionReturn(0); 1355 } 1356 1357 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "DMDACreate3d" 1360 /*@C 1361 DMDACreate3d - Creates an object that will manage the communication of three-dimensional 1362 regular array data that is distributed across some processors. 1363 1364 Collective on MPI_Comm 1365 1366 Input Parameters: 1367 + comm - MPI communicator 1368 . bx,by,bz - type of ghost nodes the array have. 1369 Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC. 1370 . stencil_type - Type of stencil (DMDA_STENCIL_STAR or DMDA_STENCIL_BOX) 1371 . 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 1372 from the command line with -da_grid_x <M> -da_grid_y <N> -da_grid_z <P>) 1373 . m,n,p - corresponding number of processors in each dimension 1374 (or PETSC_DECIDE to have calculated) 1375 . dof - number of degrees of freedom per node 1376 . lx, ly, lz - arrays containing the number of nodes in each cell along 1377 the x, y, and z coordinates, or PETSC_NULL. If non-null, these 1378 must be of length as m,n,p and the corresponding 1379 m,n, or p cannot be PETSC_DECIDE. Sum of the lx[] entries must be M, sum of 1380 the ly[] must N, sum of the lz[] must be P 1381 - s - stencil width 1382 1383 Output Parameter: 1384 . da - the resulting distributed array object 1385 1386 Options Database Key: 1387 + -da_view - Calls DMView() at the conclusion of DMDACreate3d() 1388 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 1389 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 1390 . -da_grid_z <nz> - number of grid points in z direction, if P < 0 1391 . -da_processors_x <MX> - number of processors in x direction 1392 . -da_processors_y <MY> - number of processors in y direction 1393 . -da_processors_z <MZ> - number of processors in z direction 1394 . -da_refine_x <rx> - refinement ratio in x direction 1395 . -da_refine_y <ry> - refinement ratio in y direction 1396 . -da_refine_z <rz>- refinement ratio in z directio 1397 - -da_refine <n> - refine the DMDA n times before creating it, , if M, N, or P < 0 1398 1399 Level: beginner 1400 1401 Notes: 1402 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 1403 standard 7-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 1404 the standard 27-pt stencil. 1405 1406 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 1407 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 1408 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 1409 1410 .keywords: distributed array, create, three-dimensional 1411 1412 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate2d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 1413 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(), 1414 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 1415 1416 @*/ 1417 PetscErrorCode DMDACreate3d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz,DMDAStencilType stencil_type,PetscInt M, 1418 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) 1419 { 1420 PetscErrorCode ierr; 1421 1422 PetscFunctionBegin; 1423 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 1424 ierr = DMDASetDim(*da, 3);CHKERRQ(ierr); 1425 ierr = DMDASetSizes(*da, M, N, P);CHKERRQ(ierr); 1426 ierr = DMDASetNumProcs(*da, m, n, p);CHKERRQ(ierr); 1427 ierr = DMDASetBoundaryType(*da, bx, by, bz);CHKERRQ(ierr); 1428 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 1429 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 1430 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 1431 ierr = DMDASetOwnershipRanges(*da, lx, ly, lz);CHKERRQ(ierr); 1432 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 1433 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 1434 ierr = DMSetUp(*da);CHKERRQ(ierr); 1435 ierr = DMView_DA_Private(*da);CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438