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