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