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