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