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