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