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