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