1 /* 2 Code for manipulating distributed regular 3d arrays in parallel. 3 File created by Peter Mell 7/14/95 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 #include <petscdraw.h> 9 static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer) 10 { 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) 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) 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, 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) 184 } else if (ismatlab) { 185 PetscCall(DMView_DA_Matlab(da, viewer)); 186 #endif 187 } 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 PetscErrorCode DMSetUp_DA_3D(DM da) 192 { 193 DM_DA *dd = (DM_DA *)da->data; 194 const PetscInt M = dd->M; 195 const PetscInt N = dd->N; 196 const PetscInt P = dd->P; 197 PetscInt m = dd->m; 198 PetscInt n = dd->n; 199 PetscInt p = dd->p; 200 const PetscInt dof = dd->w; 201 const PetscInt s = dd->s; 202 DMBoundaryType bx = dd->bx; 203 DMBoundaryType by = dd->by; 204 DMBoundaryType bz = dd->bz; 205 DMDAStencilType stencil_type = dd->stencil_type; 206 PetscInt *lx = dd->lx; 207 PetscInt *ly = dd->ly; 208 PetscInt *lz = dd->lz; 209 MPI_Comm comm; 210 PetscMPIInt rank, size; 211 PetscInt xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0; 212 PetscInt Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm; 213 PetscInt left, right, up, down, bottom, top, i, j, k, *idx, nn; 214 PetscInt n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14; 215 PetscInt n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26; 216 PetscInt *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z; 217 PetscInt sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0; 218 PetscInt sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0; 219 PetscInt sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0; 220 Vec local, global; 221 VecScatter gtol; 222 IS to, from; 223 PetscBool twod; 224 225 PetscFunctionBegin; 226 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"); 227 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 228 #if !defined(PETSC_USE_64BIT_INDICES) 229 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); 230 #endif 231 232 PetscCallMPI(MPI_Comm_size(comm, &size)); 233 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 234 235 if (m != PETSC_DECIDE) { 236 PetscCheck(m >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m); 237 PetscCheck(m <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %d", m, size); 238 } 239 if (n != PETSC_DECIDE) { 240 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %" PetscInt_FMT, n); 241 PetscCheck(n <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %d", n, size); 242 } 243 if (p != PETSC_DECIDE) { 244 PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Z direction: %" PetscInt_FMT, p); 245 PetscCheck(p <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Z direction: %" PetscInt_FMT " %d", p, size); 246 } 247 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); 248 249 /* Partition the array among the processors */ 250 if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) { 251 m = size / (n * p); 252 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { 253 n = size / (m * p); 254 } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { 255 p = size / (m * n); 256 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { 257 /* try for squarish distribution */ 258 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N * p))); 259 if (!m) m = 1; 260 while (m > 0) { 261 n = size / (m * p); 262 if (m * n * p == size) break; 263 m--; 264 } 265 PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad p value: p = %" PetscInt_FMT, p); 266 if (M > N && m < n) { 267 PetscInt _m = m; 268 m = n; 269 n = _m; 270 } 271 } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { 272 /* try for squarish distribution */ 273 m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n))); 274 if (!m) m = 1; 275 while (m > 0) { 276 p = size / (m * n); 277 if (m * n * p == size) break; 278 m--; 279 } 280 PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad n value: n = %" PetscInt_FMT, n); 281 if (M > P && m < p) { 282 PetscInt _m = m; 283 m = p; 284 p = _m; 285 } 286 } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { 287 /* try for squarish distribution */ 288 n = (int)(0.5 + PetscSqrtReal(((PetscReal)N) * ((PetscReal)size) / ((PetscReal)P * m))); 289 if (!n) n = 1; 290 while (n > 0) { 291 p = size / (m * n); 292 if (m * n * p == size) break; 293 n--; 294 } 295 PetscCheck(n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad m value: m = %" PetscInt_FMT, n); 296 if (N > P && n < p) { 297 PetscInt _n = n; 298 n = p; 299 p = _n; 300 } 301 } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { 302 /* try for squarish distribution */ 303 n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.))); 304 if (!n) n = 1; 305 while (n > 0) { 306 pm = size / n; 307 if (n * pm == size) break; 308 n--; 309 } 310 if (!n) n = 1; 311 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n))); 312 if (!m) m = 1; 313 while (m > 0) { 314 p = size / (m * n); 315 if (m * n * p == size) break; 316 m--; 317 } 318 if (M > P && m < p) { 319 PetscInt _m = m; 320 m = p; 321 p = _m; 322 } 323 } else PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition"); 324 325 PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_PLIB, "Could not find good partition"); 326 PetscCheck(M >= m, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m); 327 PetscCheck(N >= n, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n); 328 PetscCheck(P >= p, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in z direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, P, p); 329 330 /* 331 Determine locally owned region 332 [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes 333 */ 334 335 if (!lx) { 336 PetscCall(PetscMalloc1(m, &dd->lx)); 337 lx = dd->lx; 338 for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m)); 339 } 340 x = lx[rank % m]; 341 xs = 0; 342 for (i = 0; i < (rank % m); i++) xs += lx[i]; 343 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); 344 345 if (!ly) { 346 PetscCall(PetscMalloc1(n, &dd->ly)); 347 ly = dd->ly; 348 for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n)); 349 } 350 y = ly[(rank % (m * n)) / m]; 351 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); 352 353 ys = 0; 354 for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i]; 355 356 if (!lz) { 357 PetscCall(PetscMalloc1(p, &dd->lz)); 358 lz = dd->lz; 359 for (i = 0; i < p; i++) lz[i] = P / p + ((P % p) > (i % p)); 360 } 361 z = lz[rank / (m * n)]; 362 363 /* note this is different than x- and y-, as we will handle as an important special 364 case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems 365 in a 3D code. Additional code for this case is noted with "2d case" comments */ 366 twod = PETSC_FALSE; 367 if (P == 1) twod = PETSC_TRUE; 368 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); 369 zs = 0; 370 for (i = 0; i < (rank / (m * n)); i++) zs += lz[i]; 371 ye = ys + y; 372 xe = xs + x; 373 ze = zs + z; 374 375 /* determine ghost region (Xs) and region scattered into (IXs) */ 376 if (xs - s > 0) { 377 Xs = xs - s; 378 IXs = xs - s; 379 } else { 380 if (bx) Xs = xs - s; 381 else Xs = 0; 382 IXs = 0; 383 } 384 if (xe + s <= M) { 385 Xe = xe + s; 386 IXe = xe + s; 387 } else { 388 if (bx) { 389 Xs = xs - s; 390 Xe = xe + s; 391 } else Xe = M; 392 IXe = M; 393 } 394 395 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { 396 IXs = xs - s; 397 IXe = xe + s; 398 Xs = xs - s; 399 Xe = xe + s; 400 } 401 402 if (ys - s > 0) { 403 Ys = ys - s; 404 IYs = ys - s; 405 } else { 406 if (by) Ys = ys - s; 407 else Ys = 0; 408 IYs = 0; 409 } 410 if (ye + s <= N) { 411 Ye = ye + s; 412 IYe = ye + s; 413 } else { 414 if (by) Ye = ye + s; 415 else Ye = N; 416 IYe = N; 417 } 418 419 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) { 420 IYs = ys - s; 421 IYe = ye + s; 422 Ys = ys - s; 423 Ye = ye + s; 424 } 425 426 if (zs - s > 0) { 427 Zs = zs - s; 428 IZs = zs - s; 429 } else { 430 if (bz) Zs = zs - s; 431 else Zs = 0; 432 IZs = 0; 433 } 434 if (ze + s <= P) { 435 Ze = ze + s; 436 IZe = ze + s; 437 } else { 438 if (bz) Ze = ze + s; 439 else Ze = P; 440 IZe = P; 441 } 442 443 if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) { 444 IZs = zs - s; 445 IZe = ze + s; 446 Zs = zs - s; 447 Ze = ze + s; 448 } 449 450 /* Resize all X parameters to reflect w */ 451 s_x = s; 452 s_y = s; 453 s_z = s; 454 455 /* determine starting point of each processor */ 456 nn = x * y * z; 457 PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims)); 458 PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm)); 459 bases[0] = 0; 460 for (i = 1; i <= size; i++) bases[i] = ldims[i - 1]; 461 for (i = 1; i <= size; i++) bases[i] += bases[i - 1]; 462 base = bases[rank] * dof; 463 464 /* allocate the base parallel and sequential vectors */ 465 dd->Nlocal = x * y * z * dof; 466 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global)); 467 dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof; 468 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local)); 469 470 /* generate global to local vector scatter and local to global mapping*/ 471 472 /* global to local must include ghost points within the domain, 473 but not ghost points outside the domain that aren't periodic */ 474 PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx)); 475 if (stencil_type == DMDA_STENCIL_BOX) { 476 left = IXs - Xs; 477 right = left + (IXe - IXs); 478 bottom = IYs - Ys; 479 top = bottom + (IYe - IYs); 480 down = IZs - Zs; 481 up = down + (IZe - IZs); 482 count = 0; 483 for (i = down; i < up; i++) { 484 for (j = bottom; j < top; j++) { 485 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 486 } 487 } 488 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 489 } else { 490 /* This is way ugly! We need to list the funny cross type region */ 491 left = xs - Xs; 492 right = left + x; 493 bottom = ys - Ys; 494 top = bottom + y; 495 down = zs - Zs; 496 up = down + z; 497 count = 0; 498 /* the bottom chunk */ 499 for (i = (IZs - Zs); i < down; i++) { 500 for (j = bottom; j < top; j++) { 501 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 502 } 503 } 504 /* the middle piece */ 505 for (i = down; i < up; i++) { 506 /* front */ 507 for (j = (IYs - Ys); j < bottom; j++) { 508 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 509 } 510 /* middle */ 511 for (j = bottom; j < top; j++) { 512 for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 513 } 514 /* back */ 515 for (j = top; j < top + IYe - ye; j++) { 516 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 517 } 518 } 519 /* the top piece */ 520 for (i = up; i < up + IZe - ze; i++) { 521 for (j = bottom; j < top; j++) { 522 for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; 523 } 524 } 525 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 526 } 527 528 /* determine who lies on each side of use stored in n24 n25 n26 529 n21 n22 n23 530 n18 n19 n20 531 532 n15 n16 n17 533 n12 n14 534 n9 n10 n11 535 536 n6 n7 n8 537 n3 n4 n5 538 n0 n1 n2 539 */ 540 541 /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ 542 /* Assume Nodes are Internal to the Cube */ 543 n0 = rank - m * n - m - 1; 544 n1 = rank - m * n - m; 545 n2 = rank - m * n - m + 1; 546 n3 = rank - m * n - 1; 547 n4 = rank - m * n; 548 n5 = rank - m * n + 1; 549 n6 = rank - m * n + m - 1; 550 n7 = rank - m * n + m; 551 n8 = rank - m * n + m + 1; 552 553 n9 = rank - m - 1; 554 n10 = rank - m; 555 n11 = rank - m + 1; 556 n12 = rank - 1; 557 n14 = rank + 1; 558 n15 = rank + m - 1; 559 n16 = rank + m; 560 n17 = rank + m + 1; 561 562 n18 = rank + m * n - m - 1; 563 n19 = rank + m * n - m; 564 n20 = rank + m * n - m + 1; 565 n21 = rank + m * n - 1; 566 n22 = rank + m * n; 567 n23 = rank + m * n + 1; 568 n24 = rank + m * n + m - 1; 569 n25 = rank + m * n + m; 570 n26 = rank + m * n + m + 1; 571 572 /* Assume Pieces are on Faces of Cube */ 573 574 if (xs == 0) { /* First assume not corner or edge */ 575 n0 = rank - 1 - (m * n); 576 n3 = rank + m - 1 - (m * n); 577 n6 = rank + 2 * m - 1 - (m * n); 578 n9 = rank - 1; 579 n12 = rank + m - 1; 580 n15 = rank + 2 * m - 1; 581 n18 = rank - 1 + (m * n); 582 n21 = rank + m - 1 + (m * n); 583 n24 = rank + 2 * m - 1 + (m * n); 584 } 585 586 if (xe == M) { /* First assume not corner or edge */ 587 n2 = rank - 2 * m + 1 - (m * n); 588 n5 = rank - m + 1 - (m * n); 589 n8 = rank + 1 - (m * n); 590 n11 = rank - 2 * m + 1; 591 n14 = rank - m + 1; 592 n17 = rank + 1; 593 n20 = rank - 2 * m + 1 + (m * n); 594 n23 = rank - m + 1 + (m * n); 595 n26 = rank + 1 + (m * n); 596 } 597 598 if (ys == 0) { /* First assume not corner or edge */ 599 n0 = rank + m * (n - 1) - 1 - (m * n); 600 n1 = rank + m * (n - 1) - (m * n); 601 n2 = rank + m * (n - 1) + 1 - (m * n); 602 n9 = rank + m * (n - 1) - 1; 603 n10 = rank + m * (n - 1); 604 n11 = rank + m * (n - 1) + 1; 605 n18 = rank + m * (n - 1) - 1 + (m * n); 606 n19 = rank + m * (n - 1) + (m * n); 607 n20 = rank + m * (n - 1) + 1 + (m * n); 608 } 609 610 if (ye == N) { /* First assume not corner or edge */ 611 n6 = rank - m * (n - 1) - 1 - (m * n); 612 n7 = rank - m * (n - 1) - (m * n); 613 n8 = rank - m * (n - 1) + 1 - (m * n); 614 n15 = rank - m * (n - 1) - 1; 615 n16 = rank - m * (n - 1); 616 n17 = rank - m * (n - 1) + 1; 617 n24 = rank - m * (n - 1) - 1 + (m * n); 618 n25 = rank - m * (n - 1) + (m * n); 619 n26 = rank - m * (n - 1) + 1 + (m * n); 620 } 621 622 if (zs == 0) { /* First assume not corner or edge */ 623 n0 = size - (m * n) + rank - m - 1; 624 n1 = size - (m * n) + rank - m; 625 n2 = size - (m * n) + rank - m + 1; 626 n3 = size - (m * n) + rank - 1; 627 n4 = size - (m * n) + rank; 628 n5 = size - (m * n) + rank + 1; 629 n6 = size - (m * n) + rank + m - 1; 630 n7 = size - (m * n) + rank + m; 631 n8 = size - (m * n) + rank + m + 1; 632 } 633 634 if (ze == P) { /* First assume not corner or edge */ 635 n18 = (m * n) - (size - rank) - m - 1; 636 n19 = (m * n) - (size - rank) - m; 637 n20 = (m * n) - (size - rank) - m + 1; 638 n21 = (m * n) - (size - rank) - 1; 639 n22 = (m * n) - (size - rank); 640 n23 = (m * n) - (size - rank) + 1; 641 n24 = (m * n) - (size - rank) + m - 1; 642 n25 = (m * n) - (size - rank) + m; 643 n26 = (m * n) - (size - rank) + m + 1; 644 } 645 646 if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */ 647 n0 = size - m * n + rank + m - 1 - m; 648 n3 = size - m * n + rank + m - 1; 649 n6 = size - m * n + rank + m - 1 + m; 650 } 651 652 if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */ 653 n18 = m * n - (size - rank) + m - 1 - m; 654 n21 = m * n - (size - rank) + m - 1; 655 n24 = m * n - (size - rank) + m - 1 + m; 656 } 657 658 if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */ 659 n0 = rank + m * n - 1 - m * n; 660 n9 = rank + m * n - 1; 661 n18 = rank + m * n - 1 + m * n; 662 } 663 664 if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */ 665 n6 = rank - m * (n - 1) + m - 1 - m * n; 666 n15 = rank - m * (n - 1) + m - 1; 667 n24 = rank - m * (n - 1) + m - 1 + m * n; 668 } 669 670 if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */ 671 n2 = size - (m * n - rank) - (m - 1) - m; 672 n5 = size - (m * n - rank) - (m - 1); 673 n8 = size - (m * n - rank) - (m - 1) + m; 674 } 675 676 if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */ 677 n20 = m * n - (size - rank) - (m - 1) - m; 678 n23 = m * n - (size - rank) - (m - 1); 679 n26 = m * n - (size - rank) - (m - 1) + m; 680 } 681 682 if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */ 683 n2 = rank + m * (n - 1) - (m - 1) - m * n; 684 n11 = rank + m * (n - 1) - (m - 1); 685 n20 = rank + m * (n - 1) - (m - 1) + m * n; 686 } 687 688 if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */ 689 n8 = rank - m * n + 1 - m * n; 690 n17 = rank - m * n + 1; 691 n26 = rank - m * n + 1 + m * n; 692 } 693 694 if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */ 695 n0 = size - m + rank - 1; 696 n1 = size - m + rank; 697 n2 = size - m + rank + 1; 698 } 699 700 if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */ 701 n18 = m * n - (size - rank) + m * (n - 1) - 1; 702 n19 = m * n - (size - rank) + m * (n - 1); 703 n20 = m * n - (size - rank) + m * (n - 1) + 1; 704 } 705 706 if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */ 707 n6 = size - (m * n - rank) - m * (n - 1) - 1; 708 n7 = size - (m * n - rank) - m * (n - 1); 709 n8 = size - (m * n - rank) - m * (n - 1) + 1; 710 } 711 712 if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */ 713 n24 = rank - (size - m) - 1; 714 n25 = rank - (size - m); 715 n26 = rank - (size - m) + 1; 716 } 717 718 /* Check for Corners */ 719 if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1; 720 if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1; 721 if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1); 722 if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1; 723 if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m; 724 if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m; 725 if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n; 726 if ((xe == M) && (ye == N) && (ze == P)) n26 = 0; 727 728 /* Check for when not X,Y, and Z Periodic */ 729 730 /* If not X periodic */ 731 if (bx != DM_BOUNDARY_PERIODIC) { 732 if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2; 733 if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2; 734 } 735 736 /* If not Y periodic */ 737 if (by != DM_BOUNDARY_PERIODIC) { 738 if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2; 739 if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2; 740 } 741 742 /* If not Z periodic */ 743 if (bz != DM_BOUNDARY_PERIODIC) { 744 if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2; 745 if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2; 746 } 747 748 PetscCall(PetscMalloc1(27, &dd->neighbors)); 749 750 dd->neighbors[0] = n0; 751 dd->neighbors[1] = n1; 752 dd->neighbors[2] = n2; 753 dd->neighbors[3] = n3; 754 dd->neighbors[4] = n4; 755 dd->neighbors[5] = n5; 756 dd->neighbors[6] = n6; 757 dd->neighbors[7] = n7; 758 dd->neighbors[8] = n8; 759 dd->neighbors[9] = n9; 760 dd->neighbors[10] = n10; 761 dd->neighbors[11] = n11; 762 dd->neighbors[12] = n12; 763 dd->neighbors[13] = rank; 764 dd->neighbors[14] = n14; 765 dd->neighbors[15] = n15; 766 dd->neighbors[16] = n16; 767 dd->neighbors[17] = n17; 768 dd->neighbors[18] = n18; 769 dd->neighbors[19] = n19; 770 dd->neighbors[20] = n20; 771 dd->neighbors[21] = n21; 772 dd->neighbors[22] = n22; 773 dd->neighbors[23] = n23; 774 dd->neighbors[24] = n24; 775 dd->neighbors[25] = n25; 776 dd->neighbors[26] = n26; 777 778 /* If star stencil then delete the corner neighbors */ 779 if (stencil_type == DMDA_STENCIL_STAR) { 780 /* save information about corner neighbors */ 781 sn0 = n0; 782 sn1 = n1; 783 sn2 = n2; 784 sn3 = n3; 785 sn5 = n5; 786 sn6 = n6; 787 sn7 = n7; 788 sn8 = n8; 789 sn9 = n9; 790 sn11 = n11; 791 sn15 = n15; 792 sn17 = n17; 793 sn18 = n18; 794 sn19 = n19; 795 sn20 = n20; 796 sn21 = n21; 797 sn23 = n23; 798 sn24 = n24; 799 sn25 = n25; 800 sn26 = n26; 801 n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; 802 } 803 804 PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx)); 805 806 nn = 0; 807 /* Bottom Level */ 808 for (k = 0; k < s_z; k++) { 809 for (i = 1; i <= s_y; i++) { 810 if (n0 >= 0) { /* left below */ 811 x_t = lx[n0 % m]; 812 y_t = ly[(n0 % (m * n)) / m]; 813 z_t = lz[n0 / (m * n)]; 814 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; 815 if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */ 816 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 817 } 818 if (n1 >= 0) { /* directly below */ 819 x_t = x; 820 y_t = ly[(n1 % (m * n)) / m]; 821 z_t = lz[n1 / (m * n)]; 822 s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; 823 if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */ 824 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 825 } 826 if (n2 >= 0) { /* right below */ 827 x_t = lx[n2 % m]; 828 y_t = ly[(n2 % (m * n)) / m]; 829 z_t = lz[n2 / (m * n)]; 830 s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; 831 if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */ 832 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 833 } 834 } 835 836 for (i = 0; i < y; i++) { 837 if (n3 >= 0) { /* directly left */ 838 x_t = lx[n3 % m]; 839 y_t = y; 840 z_t = lz[n3 / (m * n)]; 841 s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 842 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 */ 843 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 844 } 845 846 if (n4 >= 0) { /* middle */ 847 x_t = x; 848 y_t = y; 849 z_t = lz[n4 / (m * n)]; 850 s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 851 if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */ 852 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 853 } else if (bz == DM_BOUNDARY_MIRROR) { 854 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y; 855 } 856 857 if (n5 >= 0) { /* directly right */ 858 x_t = lx[n5 % m]; 859 y_t = y; 860 z_t = lz[n5 / (m * n)]; 861 s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 862 if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */ 863 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 864 } 865 } 866 867 for (i = 1; i <= s_y; i++) { 868 if (n6 >= 0) { /* left above */ 869 x_t = lx[n6 % m]; 870 y_t = ly[(n6 % (m * n)) / m]; 871 z_t = lz[n6 / (m * n)]; 872 s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 873 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 */ 874 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 875 } 876 if (n7 >= 0) { /* directly above */ 877 x_t = x; 878 y_t = ly[(n7 % (m * n)) / m]; 879 z_t = lz[n7 / (m * n)]; 880 s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 881 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 */ 882 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 883 } 884 if (n8 >= 0) { /* right above */ 885 x_t = lx[n8 % m]; 886 y_t = ly[(n8 % (m * n)) / m]; 887 z_t = lz[n8 / (m * n)]; 888 s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 889 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 */ 890 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 891 } 892 } 893 } 894 895 /* Middle Level */ 896 for (k = 0; k < z; k++) { 897 for (i = 1; i <= s_y; i++) { 898 if (n9 >= 0) { /* left below */ 899 x_t = lx[n9 % m]; 900 y_t = ly[(n9 % (m * n)) / m]; 901 /* z_t = z; */ 902 s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; 903 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 904 } 905 if (n10 >= 0) { /* directly below */ 906 x_t = x; 907 y_t = ly[(n10 % (m * n)) / m]; 908 /* z_t = z; */ 909 s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; 910 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 911 } else if (by == DM_BOUNDARY_MIRROR) { 912 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x; 913 } 914 if (n11 >= 0) { /* right below */ 915 x_t = lx[n11 % m]; 916 y_t = ly[(n11 % (m * n)) / m]; 917 /* z_t = z; */ 918 s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; 919 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 920 } 921 } 922 923 for (i = 0; i < y; i++) { 924 if (n12 >= 0) { /* directly left */ 925 x_t = lx[n12 % m]; 926 y_t = y; 927 /* z_t = z; */ 928 s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t; 929 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 930 } else if (bx == DM_BOUNDARY_MIRROR) { 931 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x; 932 } 933 934 /* Interior */ 935 s_t = bases[rank] + i * x + k * x * y; 936 for (j = 0; j < x; j++) idx[nn++] = s_t++; 937 938 if (n14 >= 0) { /* directly right */ 939 x_t = lx[n14 % m]; 940 y_t = y; 941 /* z_t = z; */ 942 s_t = bases[n14] + i * x_t + k * x_t * y_t; 943 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 944 } else if (bx == DM_BOUNDARY_MIRROR) { 945 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x; 946 } 947 } 948 949 for (i = 1; i <= s_y; i++) { 950 if (n15 >= 0) { /* left above */ 951 x_t = lx[n15 % m]; 952 y_t = ly[(n15 % (m * n)) / m]; 953 /* z_t = z; */ 954 s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t; 955 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 956 } 957 if (n16 >= 0) { /* directly above */ 958 x_t = x; 959 y_t = ly[(n16 % (m * n)) / m]; 960 /* z_t = z; */ 961 s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t; 962 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 963 } else if (by == DM_BOUNDARY_MIRROR) { 964 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x; 965 } 966 if (n17 >= 0) { /* right above */ 967 x_t = lx[n17 % m]; 968 y_t = ly[(n17 % (m * n)) / m]; 969 /* z_t = z; */ 970 s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t; 971 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 972 } 973 } 974 } 975 976 /* Upper Level */ 977 for (k = 0; k < s_z; k++) { 978 for (i = 1; i <= s_y; i++) { 979 if (n18 >= 0) { /* left below */ 980 x_t = lx[n18 % m]; 981 y_t = ly[(n18 % (m * n)) / m]; 982 /* z_t = lz[n18 / (m*n)]; */ 983 s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; 984 if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */ 985 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 986 } 987 if (n19 >= 0) { /* directly below */ 988 x_t = x; 989 y_t = ly[(n19 % (m * n)) / m]; 990 /* z_t = lz[n19 / (m*n)]; */ 991 s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; 992 if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */ 993 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 994 } 995 if (n20 >= 0) { /* right below */ 996 x_t = lx[n20 % m]; 997 y_t = ly[(n20 % (m * n)) / m]; 998 /* z_t = lz[n20 / (m*n)]; */ 999 s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; 1000 if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */ 1001 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 1002 } 1003 } 1004 1005 for (i = 0; i < y; i++) { 1006 if (n21 >= 0) { /* directly left */ 1007 x_t = lx[n21 % m]; 1008 y_t = y; 1009 /* z_t = lz[n21 / (m*n)]; */ 1010 s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t; 1011 if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */ 1012 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 1013 } 1014 1015 if (n22 >= 0) { /* middle */ 1016 x_t = x; 1017 y_t = y; 1018 /* z_t = lz[n22 / (m*n)]; */ 1019 s_t = bases[n22] + i * x_t + k * x_t * y_t; 1020 if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */ 1021 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 1022 } else if (bz == DM_BOUNDARY_MIRROR) { 1023 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x; 1024 } 1025 1026 if (n23 >= 0) { /* directly right */ 1027 x_t = lx[n23 % m]; 1028 y_t = y; 1029 /* z_t = lz[n23 / (m*n)]; */ 1030 s_t = bases[n23] + i * x_t + k * x_t * y_t; 1031 if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */ 1032 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 1033 } 1034 } 1035 1036 for (i = 1; i <= s_y; i++) { 1037 if (n24 >= 0) { /* left above */ 1038 x_t = lx[n24 % m]; 1039 y_t = ly[(n24 % (m * n)) / m]; 1040 /* z_t = lz[n24 / (m*n)]; */ 1041 s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t; 1042 if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */ 1043 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 1044 } 1045 if (n25 >= 0) { /* directly above */ 1046 x_t = x; 1047 y_t = ly[(n25 % (m * n)) / m]; 1048 /* z_t = lz[n25 / (m*n)]; */ 1049 s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t; 1050 if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */ 1051 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 1052 } 1053 if (n26 >= 0) { /* right above */ 1054 x_t = lx[n26 % m]; 1055 y_t = ly[(n26 % (m * n)) / m]; 1056 /* z_t = lz[n26 / (m*n)]; */ 1057 s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t; 1058 if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */ 1059 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 1060 } 1061 } 1062 } 1063 1064 PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from)); 1065 PetscCall(VecScatterCreate(global, from, local, to, >ol)); 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 1407 PetscCall(PetscFree2(bases, ldims)); 1408 dd->m = m; 1409 dd->n = n; 1410 dd->p = p; 1411 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1412 dd->xs = xs * dof; 1413 dd->xe = xe * dof; 1414 dd->ys = ys; 1415 dd->ye = ye; 1416 dd->zs = zs; 1417 dd->ze = ze; 1418 dd->Xs = Xs * dof; 1419 dd->Xe = Xe * dof; 1420 dd->Ys = Ys; 1421 dd->Ye = Ye; 1422 dd->Zs = Zs; 1423 dd->Ze = Ze; 1424 1425 PetscCall(VecDestroy(&local)); 1426 PetscCall(VecDestroy(&global)); 1427 1428 dd->gtol = gtol; 1429 dd->base = base; 1430 da->ops->view = DMView_DA_3d; 1431 dd->ltol = NULL; 1432 dd->ao = NULL; 1433 PetscFunctionReturn(PETSC_SUCCESS); 1434 } 1435 1436 /*@C 1437 DMDACreate3d - Creates an object that will manage the communication of three-dimensional 1438 regular array data that is distributed across some processors. 1439 1440 Collective 1441 1442 Input Parameters: 1443 + comm - MPI communicator 1444 . bx - type of x ghost nodes the array have. 1445 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 1446 . by - type of y ghost nodes the array have. 1447 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 1448 . bz - type of z ghost nodes the array have. 1449 Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 1450 . stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`) 1451 . M - global dimension in x direction of the array 1452 . N - global dimension in y direction of the array 1453 . P - global dimension in z direction of the array 1454 . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated) 1455 . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated) 1456 . p - corresponding number of processors in z dimension (or `PETSC_DECIDE` to have calculated) 1457 . dof - number of degrees of freedom per node 1458 . s - stencil width 1459 . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`. 1460 . ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`. 1461 - lz - arrays containing the number of nodes in each cell along the z coordinates, or `NULL`. 1462 1463 Output Parameter: 1464 . da - the resulting distributed array object 1465 1466 Options Database Keys: 1467 + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()` 1468 . -da_grid_x <nx> - number of grid points in x direction 1469 . -da_grid_y <ny> - number of grid points in y direction 1470 . -da_grid_z <nz> - number of grid points in z direction 1471 . -da_processors_x <MX> - number of processors in x direction 1472 . -da_processors_y <MY> - number of processors in y direction 1473 . -da_processors_z <MZ> - number of processors in z direction 1474 . -da_refine_x <rx> - refinement ratio in x direction 1475 . -da_refine_y <ry> - refinement ratio in y direction 1476 . -da_refine_z <rz> - refinement ratio in z directio 1477 - -da_refine <n> - refine the DMDA n times before creating it 1478 1479 Level: beginner 1480 1481 Notes: 1482 If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the 1483 corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`, 1484 sum of the `ly` must `N`, sum of the `lz` must be `P`. 1485 1486 The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the 1487 standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes 1488 the standard 27-pt stencil. 1489 1490 The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects; 1491 The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()` 1492 and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed. 1493 1494 You must call `DMSetUp()` after this call before using this `DM`. 1495 1496 If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call 1497 but before `DMSetUp()`. 1498 1499 .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, 1500 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, 1501 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 1502 `DMStagCreate3d()` 1503 @*/ 1504 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) 1505 { 1506 PetscFunctionBegin; 1507 PetscCall(DMDACreate(comm, da)); 1508 PetscCall(DMSetDimension(*da, 3)); 1509 PetscCall(DMDASetSizes(*da, M, N, P)); 1510 PetscCall(DMDASetNumProcs(*da, m, n, p)); 1511 PetscCall(DMDASetBoundaryType(*da, bx, by, bz)); 1512 PetscCall(DMDASetDof(*da, dof)); 1513 PetscCall(DMDASetStencilType(*da, stencil_type)); 1514 PetscCall(DMDASetStencilWidth(*da, s)); 1515 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz)); 1516 PetscFunctionReturn(PETSC_SUCCESS); 1517 } 1518