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