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