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