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