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