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