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