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