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