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