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