/* Code for manipulating distributed regular 3d arrays in parallel. File created by Peter Mell 7/14/95 */ #include /*I "petscdmda.h" I*/ #include static PetscErrorCode DMView_DA_3d(DM da, PetscViewer viewer) { PetscMPIInt rank; PetscBool isascii, isdraw, isglvis, isbinary; DM_DA *dd = (DM_DA *)da->data; Vec coordinates; #if defined(PETSC_HAVE_MATLAB) PetscBool ismatlab; #endif PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); #if defined(PETSC_HAVE_MATLAB) PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab)); #endif if (isascii) { PetscViewerFormat format; PetscCall(PetscViewerASCIIPushSynchronized(viewer)); PetscCall(PetscViewerGetFormat(viewer, &format)); if (format == PETSC_VIEWER_LOAD_BALANCE) { PetscInt i, nmax = 0, nmin = PETSC_INT_MAX, navg = 0, *nz, nzlocal; DMDALocalInfo info; PetscMPIInt size; PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size)); PetscCall(DMDAGetLocalInfo(da, &info)); nzlocal = info.xm * info.ym * info.zm; PetscCall(PetscMalloc1(size, &nz)); PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da))); for (i = 0; i < size; i++) { nmax = PetscMax(nmax, nz[i]); nmin = PetscMin(nmin, nz[i]); navg += nz[i]; } PetscCall(PetscFree(nz)); navg = navg / size; PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax)); PetscFunctionReturn(PETSC_SUCCESS); } if (format != PETSC_VIEWER_ASCII_GLVIS) { DMDALocalInfo info; PetscCall(DMDAGetLocalInfo(da, &info)); 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)); 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, info.xs + info.xm, info.ys, info.ys + info.ym, info.zs, info.zs + info.zm)); PetscCall(DMGetCoordinates(da, &coordinates)); #if !defined(PETSC_USE_COMPLEX) if (coordinates) { PetscInt last; const PetscReal *coors; PetscCall(VecGetArrayRead(coordinates, &coors)); PetscCall(VecGetLocalSize(coordinates, &last)); last = last - 3; 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])); PetscCall(VecRestoreArrayRead(coordinates, &coors)); } #endif PetscCall(PetscViewerFlush(viewer)); PetscCall(PetscViewerASCIIPopSynchronized(viewer)); } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer)); } else if (isdraw) { PetscDraw draw; PetscReal ymin = -1.0, ymax = (PetscReal)dd->N; PetscReal xmin = -1.0, xmax = (PetscReal)((dd->M + 2) * dd->P), x, y, ycoord, xcoord; PetscInt k, plane, base; const PetscInt *idx; char node[10]; PetscBool isnull; PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); PetscCall(PetscDrawIsNull(draw, &isnull)); if (isnull) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(PetscDrawCheckResizedWindow(draw)); PetscCall(PetscDrawClear(draw)); PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax)); PetscDrawCollectiveBegin(draw); /* first processor draw all node lines */ if (rank == 0) { for (k = 0; k < dd->P; k++) { ymin = 0.0; ymax = (PetscReal)(dd->N - 1); 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)); xmin = (PetscReal)(k * (dd->M + 1)); xmax = xmin + (PetscReal)(dd->M - 1); for (ymin = 0; ymin < (PetscReal)dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK)); } } PetscDrawCollectiveEnd(draw); PetscCall(PetscDrawFlush(draw)); PetscCall(PetscDrawPause(draw)); PetscDrawCollectiveBegin(draw); /*Go through and draw for each plane*/ for (k = 0; k < dd->P; k++) { if ((k >= dd->zs) && (k < dd->ze)) { /* draw my box */ ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs / dd->w + (dd->M + 1) * k; xmax = (dd->xe - 1) / dd->w + (dd->M + 1) * k; PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED)); PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED)); PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED)); PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED)); xmin = dd->xs / dd->w; xmax = (dd->xe - 1) / dd->w; /* identify which processor owns the box */ PetscCall(PetscSNPrintf(node, sizeof(node), "%d", rank)); PetscCall(PetscDrawString(draw, xmin + (dd->M + 1) * k + .2, ymin + .3, PETSC_DRAW_RED, node)); /* put in numbers*/ base = (dd->base + (dd->xe - dd->xs) * (dd->ye - dd->ys) * (k - dd->zs)) / dd->w; for (y = ymin; y <= ymax; y++) { for (x = xmin + (dd->M + 1) * k; x <= xmax + (dd->M + 1) * k; x++) { PetscCall(PetscSNPrintf(node, sizeof(node), "%" PetscInt_FMT, base++)); PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node)); } } } } PetscDrawCollectiveEnd(draw); PetscCall(PetscDrawFlush(draw)); PetscCall(PetscDrawPause(draw)); PetscDrawCollectiveBegin(draw); for (k = 0 - dd->s; k < dd->P + dd->s; k++) { /* Go through and draw for each plane */ if ((k >= dd->Zs) && (k < dd->Ze)) { /* overlay ghost numbers, useful for error checking */ base = (dd->Xe - dd->Xs) * (dd->Ye - dd->Ys) * (k - dd->Zs) / dd->w; PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx)); plane = k; /* Keep z wrap around points on the drawing */ if (k < 0) plane = dd->P + k; if (k >= dd->P) plane = k - dd->P; ymin = dd->Ys; ymax = dd->Ye; xmin = (dd->M + 1) * plane * dd->w; xmax = (dd->M + 1) * plane * dd->w + dd->M * dd->w; for (y = ymin; y < ymax; y++) { for (x = xmin + dd->Xs; x < xmin + dd->Xe; x += dd->w) { PetscCall(PetscSNPrintf(node, PETSC_STATIC_ARRAY_LENGTH(node), "%" PetscInt_FMT, idx[base])); ycoord = y; /*Keep y wrap around points on drawing */ if (y < 0) ycoord = dd->N + y; if (y >= dd->N) ycoord = y - dd->N; xcoord = x; /* Keep x wrap points on drawing */ if (x < xmin) xcoord = xmax - (xmin - x); if (x >= xmax) xcoord = xmin + (x - xmax); PetscCall(PetscDrawString(draw, xcoord / dd->w, ycoord, PETSC_DRAW_BLUE, node)); base++; } } PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx)); } } PetscDrawCollectiveEnd(draw); PetscCall(PetscDrawFlush(draw)); PetscCall(PetscDrawPause(draw)); PetscCall(PetscDrawSave(draw)); } else if (isglvis) { PetscCall(DMView_DA_GLVis(da, viewer)); } else if (isbinary) { PetscCall(DMView_DA_Binary(da, viewer)); #if defined(PETSC_HAVE_MATLAB) } else if (ismatlab) { PetscCall(DMView_DA_Matlab(da, viewer)); #endif } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMSetUp_DA_3D(DM da) { DM_DA *dd = (DM_DA *)da->data; const PetscInt M = dd->M; const PetscInt N = dd->N; const PetscInt P = dd->P; PetscMPIInt m, n, p; const PetscInt dof = dd->w; const PetscInt s = dd->s; DMBoundaryType bx = dd->bx; DMBoundaryType by = dd->by; DMBoundaryType bz = dd->bz; DMDAStencilType stencil_type = dd->stencil_type; PetscInt *lx = dd->lx; PetscInt *ly = dd->ly; PetscInt *lz = dd->lz; MPI_Comm comm; PetscMPIInt rank, size; PetscInt xs = 0, xe, ys = 0, ye, zs = 0, ze, x = 0, y = 0, z = 0; PetscInt Xs, Xe, Ys, Ye, Zs, Ze, IXs, IXe, IYs, IYe, IZs, IZe, pm; PetscInt left, right, up, down, bottom, top, i, j, k, *idx, nn; PetscMPIInt n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n14; PetscMPIInt n15, n16, n17, n18, n19, n20, n21, n22, n23, n24, n25, n26; PetscInt *bases, *ldims, base, x_t, y_t, z_t, s_t, count, s_x, s_y, s_z; PetscMPIInt sn0 = 0, sn1 = 0, sn2 = 0, sn3 = 0, sn5 = 0, sn6 = 0, sn7 = 0; PetscMPIInt sn8 = 0, sn9 = 0, sn11 = 0, sn15 = 0, sn24 = 0, sn25 = 0, sn26 = 0; PetscMPIInt sn17 = 0, sn18 = 0, sn19 = 0, sn20 = 0, sn21 = 0, sn23 = 0; Vec local, global; VecScatter gtol; IS to, from; PetscBool twod; PetscFunctionBegin; 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"); PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); #if !defined(PETSC_USE_64BIT_INDICES) 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); #endif PetscCall(PetscMPIIntCast(dd->m, &m)); PetscCall(PetscMPIIntCast(dd->n, &n)); PetscCall(PetscMPIIntCast(dd->p, &p)); PetscCallMPI(MPI_Comm_size(comm, &size)); PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (m != PETSC_DECIDE) { PetscCheck(m >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %d", m); PetscCheck(m <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %d %d", m, size); } if (n != PETSC_DECIDE) { PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %d", n); PetscCheck(n <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %d %d", n, size); } if (p != PETSC_DECIDE) { PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Z direction: %d", p); PetscCheck(p <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Z direction: %d %d", p, size); } 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); /* Partition the array among the processors */ if (m == PETSC_DECIDE && n != PETSC_DECIDE && p != PETSC_DECIDE) { m = size / (n * p); } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { n = size / (m * p); } else if (m != PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { p = size / (m * n); } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p != PETSC_DECIDE) { /* try for squarish distribution */ m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N * p))); if (!m) m = 1; while (m > 0) { n = size / (m * p); if (m * n * p == size) break; m--; } PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad p value: p = %d", p); if (M > N && m < n) { PetscMPIInt _m = m; m = n; n = _m; } } else if (m == PETSC_DECIDE && n != PETSC_DECIDE && p == PETSC_DECIDE) { /* try for squarish distribution */ m = (int)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n))); if (!m) m = 1; while (m > 0) { p = size / (m * n); if (m * n * p == size) break; m--; } PetscCheck(m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad n value: n = %d", n); if (M > P && m < p) { PetscMPIInt _m = m; m = p; p = _m; } } else if (m != PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { /* try for squarish distribution */ n = (int)(0.5 + PetscSqrtReal(((PetscReal)N) * ((PetscReal)size) / ((PetscReal)P * m))); if (!n) n = 1; while (n > 0) { p = size / (m * n); if (m * n * p == size) break; n--; } PetscCheck(n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad m value: m = %d", n); if (N > P && n < p) { PetscMPIInt _n = n; n = p; p = _n; } } else if (m == PETSC_DECIDE && n == PETSC_DECIDE && p == PETSC_DECIDE) { /* try for squarish distribution */ n = (PetscMPIInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.))); if (!n) n = 1; while (n > 0) { pm = size / n; if (n * pm == size) break; n--; } if (!n) n = 1; m = (PetscMPIInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n))); if (!m) m = 1; while (m > 0) { p = size / (m * n); if (m * n * p == size) break; m--; } if (M > P && m < p) { PetscMPIInt _m = m; m = p; p = _m; } } else PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition"); PetscCheck(m * n * p == size, PetscObjectComm((PetscObject)da), PETSC_ERR_PLIB, "Could not find good partition"); PetscCheck(M >= m, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %d", M, m); PetscCheck(N >= n, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %d", N, n); PetscCheck(P >= p, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Partition in z direction is too fine! %" PetscInt_FMT " %d", P, p); /* Determine locally owned region [x, y, or z]s is the first local node number, [x, y, z] is the number of local nodes */ if (!lx) { PetscCall(PetscMalloc1(m, &dd->lx)); lx = dd->lx; for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > (i % m)); } x = lx[rank % m]; xs = 0; for (i = 0; i < (rank % m); i++) xs += lx[i]; 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); if (!ly) { PetscCall(PetscMalloc1(n, &dd->ly)); ly = dd->ly; for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > (i % n)); } y = ly[(rank % (m * n)) / m]; 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); ys = 0; for (i = 0; i < (rank % (m * n)) / m; i++) ys += ly[i]; if (!lz) { PetscCall(PetscMalloc1(p, &dd->lz)); lz = dd->lz; for (i = 0; i < p; i++) lz[i] = P / p + ((P % p) > (i % p)); } z = lz[rank / (m * n)]; 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); 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); 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); /* note this is different than x- and y-, as we will handle as an important special case when p=P=1 and DM_BOUNDARY_PERIODIC and s > z. This is to deal with 2D problems in a 3D code. Additional code for this case is noted with "2d case" comments */ twod = PETSC_FALSE; if (P == 1) twod = PETSC_TRUE; 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); zs = 0; for (i = 0; i < (rank / (m * n)); i++) zs += lz[i]; ye = ys + y; xe = xs + x; ze = zs + z; /* determine ghost region (Xs) and region scattered into (IXs) */ if (xs - s > 0) { Xs = xs - s; IXs = xs - s; } else { if (bx) Xs = xs - s; else Xs = 0; IXs = 0; } if (xe + s <= M) { Xe = xe + s; IXe = xe + s; } else { if (bx) { Xs = xs - s; Xe = xe + s; } else Xe = M; IXe = M; } if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { IXs = xs - s; IXe = xe + s; Xs = xs - s; Xe = xe + s; } if (ys - s > 0) { Ys = ys - s; IYs = ys - s; } else { if (by) Ys = ys - s; else Ys = 0; IYs = 0; } if (ye + s <= N) { Ye = ye + s; IYe = ye + s; } else { if (by) Ye = ye + s; else Ye = N; IYe = N; } if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) { IYs = ys - s; IYe = ye + s; Ys = ys - s; Ye = ye + s; } if (zs - s > 0) { Zs = zs - s; IZs = zs - s; } else { if (bz) Zs = zs - s; else Zs = 0; IZs = 0; } if (ze + s <= P) { Ze = ze + s; IZe = ze + s; } else { if (bz) Ze = ze + s; else Ze = P; IZe = P; } if (bz == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_MIRROR) { IZs = zs - s; IZe = ze + s; Zs = zs - s; Ze = ze + s; } /* Resize all X parameters to reflect w */ s_x = s; s_y = s; s_z = s; /* determine starting point of each processor */ nn = x * y * z; PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims)); PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm)); bases[0] = 0; for (i = 1; i <= size; i++) bases[i] = ldims[i - 1]; for (i = 1; i <= size; i++) bases[i] += bases[i - 1]; base = bases[rank] * dof; /* allocate the base parallel and sequential vectors */ dd->Nlocal = x * y * z * dof; PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global)); dd->nlocal = (Xe - Xs) * (Ye - Ys) * (Ze - Zs) * dof; PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local)); /* generate global to local vector scatter and local to global mapping*/ /* global to local must include ghost points within the domain, but not ghost points outside the domain that aren't periodic */ PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs) * (IZe - IZs), &idx)); if (stencil_type == DMDA_STENCIL_BOX) { left = IXs - Xs; right = left + (IXe - IXs); bottom = IYs - Ys; top = bottom + (IYe - IYs); down = IZs - Zs; up = down + (IZe - IZs); count = 0; for (i = down; i < up; i++) { for (j = bottom; j < top; j++) { for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } } PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); } else { /* This is way ugly! We need to list the funny cross type region */ left = xs - Xs; right = left + x; bottom = ys - Ys; top = bottom + y; down = zs - Zs; up = down + z; count = 0; /* the bottom chunk */ for (i = (IZs - Zs); i < down; i++) { for (j = bottom; j < top; j++) { for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } } /* the middle piece */ for (i = down; i < up; i++) { /* front */ for (j = (IYs - Ys); j < bottom; j++) { for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } /* middle */ for (j = bottom; j < top; j++) { for (k = IXs - Xs; k < IXe - Xs; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } /* back */ for (j = top; j < top + IYe - ye; j++) { for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } } /* the top piece */ for (i = up; i < up + IZe - ze; i++) { for (j = bottom; j < top; j++) { for (k = left; k < right; k++) idx[count++] = (i * (Ye - Ys) + j) * (Xe - Xs) + k; } } PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); } /* determine who lies on each side of use stored in n24 n25 n26 n21 n22 n23 n18 n19 n20 n15 n16 n17 n12 n14 n9 n10 n11 n6 n7 n8 n3 n4 n5 n0 n1 n2 */ /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */ /* Assume Nodes are Internal to the Cube */ n0 = rank - m * n - m - 1; n1 = rank - m * n - m; n2 = rank - m * n - m + 1; n3 = rank - m * n - 1; n4 = rank - m * n; n5 = rank - m * n + 1; n6 = rank - m * n + m - 1; n7 = rank - m * n + m; n8 = rank - m * n + m + 1; n9 = rank - m - 1; n10 = rank - m; n11 = rank - m + 1; n12 = rank - 1; n14 = rank + 1; n15 = rank + m - 1; n16 = rank + m; n17 = rank + m + 1; n18 = rank + m * n - m - 1; n19 = rank + m * n - m; n20 = rank + m * n - m + 1; n21 = rank + m * n - 1; n22 = rank + m * n; n23 = rank + m * n + 1; n24 = rank + m * n + m - 1; n25 = rank + m * n + m; n26 = rank + m * n + m + 1; /* Assume Pieces are on Faces of Cube */ if (xs == 0) { /* First assume not corner or edge */ n0 = rank - 1 - (m * n); n3 = rank + m - 1 - (m * n); n6 = rank + 2 * m - 1 - (m * n); n9 = rank - 1; n12 = rank + m - 1; n15 = rank + 2 * m - 1; n18 = rank - 1 + (m * n); n21 = rank + m - 1 + (m * n); n24 = rank + 2 * m - 1 + (m * n); } if (xe == M) { /* First assume not corner or edge */ n2 = rank - 2 * m + 1 - (m * n); n5 = rank - m + 1 - (m * n); n8 = rank + 1 - (m * n); n11 = rank - 2 * m + 1; n14 = rank - m + 1; n17 = rank + 1; n20 = rank - 2 * m + 1 + (m * n); n23 = rank - m + 1 + (m * n); n26 = rank + 1 + (m * n); } if (ys == 0) { /* First assume not corner or edge */ n0 = rank + m * (n - 1) - 1 - (m * n); n1 = rank + m * (n - 1) - (m * n); n2 = rank + m * (n - 1) + 1 - (m * n); n9 = rank + m * (n - 1) - 1; n10 = rank + m * (n - 1); n11 = rank + m * (n - 1) + 1; n18 = rank + m * (n - 1) - 1 + (m * n); n19 = rank + m * (n - 1) + (m * n); n20 = rank + m * (n - 1) + 1 + (m * n); } if (ye == N) { /* First assume not corner or edge */ n6 = rank - m * (n - 1) - 1 - (m * n); n7 = rank - m * (n - 1) - (m * n); n8 = rank - m * (n - 1) + 1 - (m * n); n15 = rank - m * (n - 1) - 1; n16 = rank - m * (n - 1); n17 = rank - m * (n - 1) + 1; n24 = rank - m * (n - 1) - 1 + (m * n); n25 = rank - m * (n - 1) + (m * n); n26 = rank - m * (n - 1) + 1 + (m * n); } if (zs == 0) { /* First assume not corner or edge */ n0 = size - (m * n) + rank - m - 1; n1 = size - (m * n) + rank - m; n2 = size - (m * n) + rank - m + 1; n3 = size - (m * n) + rank - 1; n4 = size - (m * n) + rank; n5 = size - (m * n) + rank + 1; n6 = size - (m * n) + rank + m - 1; n7 = size - (m * n) + rank + m; n8 = size - (m * n) + rank + m + 1; } if (ze == P) { /* First assume not corner or edge */ n18 = (m * n) - (size - rank) - m - 1; n19 = (m * n) - (size - rank) - m; n20 = (m * n) - (size - rank) - m + 1; n21 = (m * n) - (size - rank) - 1; n22 = (m * n) - (size - rank); n23 = (m * n) - (size - rank) + 1; n24 = (m * n) - (size - rank) + m - 1; n25 = (m * n) - (size - rank) + m; n26 = (m * n) - (size - rank) + m + 1; } if ((xs == 0) && (zs == 0)) { /* Assume an edge, not corner */ n0 = size - m * n + rank + m - 1 - m; n3 = size - m * n + rank + m - 1; n6 = size - m * n + rank + m - 1 + m; } if ((xs == 0) && (ze == P)) { /* Assume an edge, not corner */ n18 = m * n - (size - rank) + m - 1 - m; n21 = m * n - (size - rank) + m - 1; n24 = m * n - (size - rank) + m - 1 + m; } if ((xs == 0) && (ys == 0)) { /* Assume an edge, not corner */ n0 = rank + m * n - 1 - m * n; n9 = rank + m * n - 1; n18 = rank + m * n - 1 + m * n; } if ((xs == 0) && (ye == N)) { /* Assume an edge, not corner */ n6 = rank - m * (n - 1) + m - 1 - m * n; n15 = rank - m * (n - 1) + m - 1; n24 = rank - m * (n - 1) + m - 1 + m * n; } if ((xe == M) && (zs == 0)) { /* Assume an edge, not corner */ n2 = size - (m * n - rank) - (m - 1) - m; n5 = size - (m * n - rank) - (m - 1); n8 = size - (m * n - rank) - (m - 1) + m; } if ((xe == M) && (ze == P)) { /* Assume an edge, not corner */ n20 = m * n - (size - rank) - (m - 1) - m; n23 = m * n - (size - rank) - (m - 1); n26 = m * n - (size - rank) - (m - 1) + m; } if ((xe == M) && (ys == 0)) { /* Assume an edge, not corner */ n2 = rank + m * (n - 1) - (m - 1) - m * n; n11 = rank + m * (n - 1) - (m - 1); n20 = rank + m * (n - 1) - (m - 1) + m * n; } if ((xe == M) && (ye == N)) { /* Assume an edge, not corner */ n8 = rank - m * n + 1 - m * n; n17 = rank - m * n + 1; n26 = rank - m * n + 1 + m * n; } if ((ys == 0) && (zs == 0)) { /* Assume an edge, not corner */ n0 = size - m + rank - 1; n1 = size - m + rank; n2 = size - m + rank + 1; } if ((ys == 0) && (ze == P)) { /* Assume an edge, not corner */ n18 = m * n - (size - rank) + m * (n - 1) - 1; n19 = m * n - (size - rank) + m * (n - 1); n20 = m * n - (size - rank) + m * (n - 1) + 1; } if ((ye == N) && (zs == 0)) { /* Assume an edge, not corner */ n6 = size - (m * n - rank) - m * (n - 1) - 1; n7 = size - (m * n - rank) - m * (n - 1); n8 = size - (m * n - rank) - m * (n - 1) + 1; } if ((ye == N) && (ze == P)) { /* Assume an edge, not corner */ n24 = rank - (size - m) - 1; n25 = rank - (size - m); n26 = rank - (size - m) + 1; } /* Check for Corners */ if ((xs == 0) && (ys == 0) && (zs == 0)) n0 = size - 1; if ((xs == 0) && (ys == 0) && (ze == P)) n18 = m * n - 1; if ((xs == 0) && (ye == N) && (zs == 0)) n6 = (size - 1) - m * (n - 1); if ((xs == 0) && (ye == N) && (ze == P)) n24 = m - 1; if ((xe == M) && (ys == 0) && (zs == 0)) n2 = size - m; if ((xe == M) && (ys == 0) && (ze == P)) n20 = m * n - m; if ((xe == M) && (ye == N) && (zs == 0)) n8 = size - m * n; if ((xe == M) && (ye == N) && (ze == P)) n26 = 0; /* Check for when not X,Y, and Z Periodic */ /* If not X periodic */ if (bx != DM_BOUNDARY_PERIODIC) { if (xs == 0) n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2; if (xe == M) n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2; } /* If not Y periodic */ if (by != DM_BOUNDARY_PERIODIC) { if (ys == 0) n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2; if (ye == N) n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2; } /* If not Z periodic */ if (bz != DM_BOUNDARY_PERIODIC) { if (zs == 0) n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2; if (ze == P) n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2; } PetscCall(PetscMalloc1(27, &dd->neighbors)); dd->neighbors[0] = n0; dd->neighbors[1] = n1; dd->neighbors[2] = n2; dd->neighbors[3] = n3; dd->neighbors[4] = n4; dd->neighbors[5] = n5; dd->neighbors[6] = n6; dd->neighbors[7] = n7; dd->neighbors[8] = n8; dd->neighbors[9] = n9; dd->neighbors[10] = n10; dd->neighbors[11] = n11; dd->neighbors[12] = n12; dd->neighbors[13] = rank; dd->neighbors[14] = n14; dd->neighbors[15] = n15; dd->neighbors[16] = n16; dd->neighbors[17] = n17; dd->neighbors[18] = n18; dd->neighbors[19] = n19; dd->neighbors[20] = n20; dd->neighbors[21] = n21; dd->neighbors[22] = n22; dd->neighbors[23] = n23; dd->neighbors[24] = n24; dd->neighbors[25] = n25; dd->neighbors[26] = n26; /* If star stencil then delete the corner neighbors */ if (stencil_type == DMDA_STENCIL_STAR) { /* save information about corner neighbors */ sn0 = n0; sn1 = n1; sn2 = n2; sn3 = n3; sn5 = n5; sn6 = n6; sn7 = n7; sn8 = n8; sn9 = n9; sn11 = n11; sn15 = n15; sn17 = n17; sn18 = n18; sn19 = n19; sn20 = n20; sn21 = n21; sn23 = n23; sn24 = n24; sn25 = n25; sn26 = n26; n0 = n1 = n2 = n3 = n5 = n6 = n7 = n8 = n9 = n11 = n15 = n17 = n18 = n19 = n20 = n21 = n23 = n24 = n25 = n26 = -1; } PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys) * (Ze - Zs), &idx)); nn = 0; /* Bottom Level */ for (k = 0; k < s_z; k++) { for (i = 1; i <= s_y; i++) { if (n0 >= 0) { /* left below */ x_t = lx[n0 % m]; y_t = ly[(n0 % (m * n)) / m]; z_t = lz[n0 / (m * n)]; 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; if (twod && (s_t < 0)) s_t = bases[n0] + x_t * y_t * z_t - (s_y - i) * x_t - s_x; /* 2D case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n1 >= 0) { /* directly below */ x_t = x; y_t = ly[(n1 % (m * n)) / m]; z_t = lz[n1 / (m * n)]; s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; if (twod && (s_t < 0)) s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } if (n2 >= 0) { /* right below */ x_t = lx[n2 % m]; y_t = ly[(n2 % (m * n)) / m]; z_t = lz[n2 / (m * n)]; s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; if (twod && (s_t < 0)) s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t; /* 2D case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } for (i = 0; i < y; i++) { if (n3 >= 0) { /* directly left */ x_t = lx[n3 % m]; y_t = y; z_t = lz[n3 / (m * n)]; s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 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 */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n4 >= 0) { /* middle */ x_t = x; y_t = y; z_t = lz[n4 / (m * n)]; s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; if (twod && (s_t < 0)) s_t = bases[n4] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (bz == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k - 1) * x * y; } if (n5 >= 0) { /* directly right */ x_t = lx[n5 % m]; y_t = y; z_t = lz[n5 / (m * n)]; s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; if (twod && (s_t < 0)) s_t = bases[n5] + i * x_t + x_t * y_t * z_t - x_t * y_t; /* 2D case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } for (i = 1; i <= s_y; i++) { if (n6 >= 0) { /* left above */ x_t = lx[n6 % m]; y_t = ly[(n6 % (m * n)) / m]; z_t = lz[n6 / (m * n)]; s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 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 */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n7 >= 0) { /* directly above */ x_t = x; y_t = ly[(n7 % (m * n)) / m]; z_t = lz[n7 / (m * n)]; s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 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 */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } if (n8 >= 0) { /* right above */ x_t = lx[n8 % m]; y_t = ly[(n8 % (m * n)) / m]; z_t = lz[n8 / (m * n)]; s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; 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 */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } } /* Middle Level */ for (k = 0; k < z; k++) { for (i = 1; i <= s_y; i++) { if (n9 >= 0) { /* left below */ x_t = lx[n9 % m]; y_t = ly[(n9 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n10 >= 0) { /* directly below */ x_t = x; y_t = ly[(n10 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (by == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i) * x; } if (n11 >= 0) { /* right below */ x_t = lx[n11 % m]; y_t = ly[(n11 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } for (i = 0; i < y; i++) { if (n12 >= 0) { /* directly left */ x_t = lx[n12 % m]; y_t = y; /* z_t = z; */ s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (bx == DM_BOUNDARY_MIRROR) { for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x; } /* Interior */ s_t = bases[rank] + i * x + k * x * y; for (j = 0; j < x; j++) idx[nn++] = s_t++; if (n14 >= 0) { /* directly right */ x_t = lx[n14 % m]; y_t = y; /* z_t = z; */ s_t = bases[n14] + i * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (bx == DM_BOUNDARY_MIRROR) { for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x; } } for (i = 1; i <= s_y; i++) { if (n15 >= 0) { /* left above */ x_t = lx[n15 % m]; y_t = ly[(n15 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n16 >= 0) { /* directly above */ x_t = x; y_t = ly[(n16 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (by == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i) * x; } if (n17 >= 0) { /* right above */ x_t = lx[n17 % m]; y_t = ly[(n17 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } } /* Upper Level */ for (k = 0; k < s_z; k++) { for (i = 1; i <= s_y; i++) { if (n18 >= 0) { /* left below */ x_t = lx[n18 % m]; y_t = ly[(n18 % (m * n)) / m]; /* z_t = lz[n18 / (m*n)]; */ s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n18] - (s_y - i) * x_t - s_x + x_t * y_t; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n19 >= 0) { /* directly below */ x_t = x; y_t = ly[(n19 % (m * n)) / m]; /* z_t = lz[n19 / (m*n)]; */ s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n19] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } if (n20 >= 0) { /* right below */ x_t = lx[n20 % m]; y_t = ly[(n20 % (m * n)) / m]; /* z_t = lz[n20 / (m*n)]; */ s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n20] - (s_y + 1 - i) * x_t + x_t * y_t; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } for (i = 0; i < y; i++) { if (n21 >= 0) { /* directly left */ x_t = lx[n21 % m]; y_t = y; /* z_t = lz[n21 / (m*n)]; */ s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n21] + (i + 1) * x_t - s_x; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n22 >= 0) { /* middle */ x_t = x; y_t = y; /* z_t = lz[n22 / (m*n)]; */ s_t = bases[n22] + i * x_t + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n22] + i * x_t; /* 2d case */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (bz == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 1) * x * y + i * x; } if (n23 >= 0) { /* directly right */ x_t = lx[n23 % m]; y_t = y; /* z_t = lz[n23 / (m*n)]; */ s_t = bases[n23] + i * x_t + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n23] + i * x_t; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } for (i = 1; i <= s_y; i++) { if (n24 >= 0) { /* left above */ x_t = lx[n24 % m]; y_t = ly[(n24 % (m * n)) / m]; /* z_t = lz[n24 / (m*n)]; */ s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n24] + i * x_t - s_x; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } if (n25 >= 0) { /* directly above */ x_t = x; y_t = ly[(n25 % (m * n)) / m]; /* z_t = lz[n25 / (m*n)]; */ s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n25] + (i - 1) * x_t; /* 2d case */ for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } if (n26 >= 0) { /* right above */ x_t = lx[n26 % m]; y_t = ly[(n26 % (m * n)) / m]; /* z_t = lz[n26 / (m*n)]; */ s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t; if (twod && (s_t >= M * N * P)) s_t = bases[n26] + (i - 1) * x_t; /* 2d case */ for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } } } PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from)); PetscCall(VecScatterCreate(global, from, local, to, >ol)); PetscCall(ISDestroy(&to)); PetscCall(ISDestroy(&from)); if (stencil_type == DMDA_STENCIL_STAR) { n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7; n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18; n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25; n26 = sn26; } if ((stencil_type == DMDA_STENCIL_STAR) || (bx != DM_BOUNDARY_PERIODIC && bx) || (by != DM_BOUNDARY_PERIODIC && by) || (bz != DM_BOUNDARY_PERIODIC && bz)) { /* Recompute the local to global mappings, this time keeping the information about the cross corner processor numbers. */ nn = 0; /* Bottom Level */ for (k = 0; k < s_z; k++) { for (i = 1; i <= s_y; i++) { if (n0 >= 0) { /* left below */ x_t = lx[n0 % m]; y_t = ly[(n0 % (m * n)) / m]; z_t = lz[n0 / (m * n)]; 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; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && Ys - ys < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n1 >= 0) { /* directly below */ x_t = x; y_t = ly[(n1 % (m * n)) / m]; z_t = lz[n1 / (m * n)]; s_t = bases[n1] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (Ys - ys < 0 && Zs - zs < 0) { for (j = 0; j < x; j++) idx[nn++] = -1; } if (n2 >= 0) { /* right below */ x_t = lx[n2 % m]; y_t = ly[(n2 % (m * n)) / m]; z_t = lz[n2 / (m * n)]; s_t = bases[n2] + x_t * y_t * z_t - (s_y + 1 - i) * x_t - (s_z - k - 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && Ys - ys < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } for (i = 0; i < y; i++) { if (n3 >= 0) { /* directly left */ x_t = lx[n3 % m]; y_t = y; z_t = lz[n3 / (m * n)]; s_t = bases[n3] + (i + 1) * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n4 >= 0) { /* middle */ x_t = x; y_t = y; z_t = lz[n4 / (m * n)]; s_t = bases[n4] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (Zs - zs < 0) { if (bz == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + x * i + (s_z - k) * x * y; } else { for (j = 0; j < x; j++) idx[nn++] = -1; } } if (n5 >= 0) { /* directly right */ x_t = lx[n5 % m]; y_t = y; z_t = lz[n5 / (m * n)]; s_t = bases[n5] + i * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } for (i = 1; i <= s_y; i++) { if (n6 >= 0) { /* left above */ x_t = lx[n6 % m]; y_t = ly[(n6 % (m * n)) / m]; z_t = lz[n6 / (m * n)]; s_t = bases[n6] + i * x_t - s_x + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && ye - Ye < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n7 >= 0) { /* directly above */ x_t = x; y_t = ly[(n7 % (m * n)) / m]; z_t = lz[n7 / (m * n)]; s_t = bases[n7] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (ye - Ye < 0 && Zs - zs < 0) { for (j = 0; j < x; j++) idx[nn++] = -1; } if (n8 >= 0) { /* right above */ x_t = lx[n8 % m]; y_t = ly[(n8 % (m * n)) / m]; z_t = lz[n8 / (m * n)]; s_t = bases[n8] + (i - 1) * x_t + x_t * y_t * z_t - (s_z - k) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && ye - Ye < 0 && Zs - zs < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } } /* Middle Level */ for (k = 0; k < z; k++) { for (i = 1; i <= s_y; i++) { if (n9 >= 0) { /* left below */ x_t = lx[n9 % m]; y_t = ly[(n9 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n9] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && Ys - ys < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n10 >= 0) { /* directly below */ x_t = x; y_t = ly[(n10 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n10] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (Ys - ys < 0) { if (by == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (s_y - i + 1) * x; } else { for (j = 0; j < x; j++) idx[nn++] = -1; } } if (n11 >= 0) { /* right below */ x_t = lx[n11 % m]; y_t = ly[(n11 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n11] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && Ys - ys < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } for (i = 0; i < y; i++) { if (n12 >= 0) { /* directly left */ x_t = lx[n12 % m]; y_t = y; /* z_t = z; */ s_t = bases[n12] + (i + 1) * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0) { if (bx == DM_BOUNDARY_MIRROR) { for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + s_x - j - 1 + k * x * y + i * x + 1; } else { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } /* Interior */ s_t = bases[rank] + i * x + k * x * y; for (j = 0; j < x; j++) idx[nn++] = s_t++; if (n14 >= 0) { /* directly right */ x_t = lx[n14 % m]; y_t = y; /* z_t = z; */ s_t = bases[n14] + i * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0) { if (bx == DM_BOUNDARY_MIRROR) { for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x - j - 1 + k * x * y + i * x - 1; } else { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } } for (i = 1; i <= s_y; i++) { if (n15 >= 0) { /* left above */ x_t = lx[n15 % m]; y_t = ly[(n15 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n15] + i * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && ye - Ye < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n16 >= 0) { /* directly above */ x_t = x; y_t = ly[(n16 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n16] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (ye - Ye < 0) { if (by == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + k * x * y + (y - i - 1) * x; } else { for (j = 0; j < x; j++) idx[nn++] = -1; } } if (n17 >= 0) { /* right above */ x_t = lx[n17 % m]; y_t = ly[(n17 % (m * n)) / m]; /* z_t = z; */ s_t = bases[n17] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && ye - Ye < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } } /* Upper Level */ for (k = 0; k < s_z; k++) { for (i = 1; i <= s_y; i++) { if (n18 >= 0) { /* left below */ x_t = lx[n18 % m]; y_t = ly[(n18 % (m * n)) / m]; /* z_t = lz[n18 / (m*n)]; */ s_t = bases[n18] - (s_y - i) * x_t - s_x + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && Ys - ys < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n19 >= 0) { /* directly below */ x_t = x; y_t = ly[(n19 % (m * n)) / m]; /* z_t = lz[n19 / (m*n)]; */ s_t = bases[n19] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (Ys - ys < 0 && ze - Ze < 0) { for (j = 0; j < x; j++) idx[nn++] = -1; } if (n20 >= 0) { /* right below */ x_t = lx[n20 % m]; y_t = ly[(n20 % (m * n)) / m]; /* z_t = lz[n20 / (m*n)]; */ s_t = bases[n20] - (s_y + 1 - i) * x_t + (k + 1) * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && Ys - ys < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } for (i = 0; i < y; i++) { if (n21 >= 0) { /* directly left */ x_t = lx[n21 % m]; y_t = y; /* z_t = lz[n21 / (m*n)]; */ s_t = bases[n21] + (i + 1) * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n22 >= 0) { /* middle */ x_t = x; y_t = y; /* z_t = lz[n22 / (m*n)]; */ s_t = bases[n22] + i * x_t + k * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (ze - Ze < 0) { if (bz == DM_BOUNDARY_MIRROR) { for (j = 0; j < x; j++) idx[nn++] = bases[rank] + j + (z - k - 2) * x * y + i * x; } else { for (j = 0; j < x; j++) idx[nn++] = -1; } } if (n23 >= 0) { /* directly right */ x_t = lx[n23 % m]; y_t = y; /* z_t = lz[n23 / (m*n)]; */ s_t = bases[n23] + i * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } for (i = 1; i <= s_y; i++) { if (n24 >= 0) { /* left above */ x_t = lx[n24 % m]; y_t = ly[(n24 % (m * n)) / m]; /* z_t = lz[n24 / (m*n)]; */ s_t = bases[n24] + i * x_t - s_x + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (Xs - xs < 0 && ye - Ye < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } if (n25 >= 0) { /* directly above */ x_t = x; y_t = ly[(n25 % (m * n)) / m]; /* z_t = lz[n25 / (m*n)]; */ s_t = bases[n25] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < x_t; j++) idx[nn++] = s_t++; } else if (ye - Ye < 0 && ze - Ze < 0) { for (j = 0; j < x; j++) idx[nn++] = -1; } if (n26 >= 0) { /* right above */ x_t = lx[n26 % m]; y_t = ly[(n26 % (m * n)) / m]; /* z_t = lz[n26 / (m*n)]; */ s_t = bases[n26] + (i - 1) * x_t + k * x_t * y_t; for (j = 0; j < s_x; j++) idx[nn++] = s_t++; } else if (xe - Xe < 0 && ye - Ye < 0 && ze - Ze < 0) { for (j = 0; j < s_x; j++) idx[nn++] = -1; } } } } /* Set the local to global ordering in the global vector, this allows use of VecSetValuesLocal(). */ PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap)); PetscCall(PetscFree2(bases, ldims)); dd->m = m; dd->n = n; dd->p = p; /* note PETSc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ dd->xs = xs * dof; dd->xe = xe * dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze; dd->Xs = Xs * dof; dd->Xe = Xe * dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze; PetscCall(VecDestroy(&local)); PetscCall(VecDestroy(&global)); dd->gtol = gtol; dd->base = base; da->ops->view = DMView_DA_3d; dd->ltol = NULL; dd->ao = NULL; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMDACreate3d - Creates an object that will manage the communication of three-dimensional regular array data that is distributed across one or more MPI processes. Collective Input Parameters: + comm - MPI communicator . bx - type of x ghost nodes the array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. . by - type of y ghost nodes the array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. . bz - type of z ghost nodes the array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. . stencil_type - Type of stencil (`DMDA_STENCIL_STAR` or `DMDA_STENCIL_BOX`) . M - global dimension in x direction of the array . N - global dimension in y direction of the array . P - global dimension in z direction of the array . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated) . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated) . p - corresponding number of processors in z dimension (or `PETSC_DECIDE` to have calculated) . dof - number of degrees of freedom per node . s - stencil width . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`. . ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`. - lz - arrays containing the number of nodes in each cell along the z coordinates, or `NULL`. Output Parameter: . da - the resulting distributed array object Options Database Keys: + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate3d()` . -da_grid_x - number of grid points in x direction . -da_grid_y - number of grid points in y direction . -da_grid_z - number of grid points in z direction . -da_processors_x - number of processors in x direction . -da_processors_y - number of processors in y direction . -da_processors_z - number of processors in z direction . -da_bd_x - boundary type in x direction . -da_bd_y - boundary type in y direction . -da_bd_z - boundary type in x direction . -da_bd_all - boundary type in all directions . -da_refine_x - refinement ratio in x direction . -da_refine_y - refinement ratio in y direction . -da_refine_z - refinement ratio in z directio - -da_refine - refine the `DMDA` n times before creating it Level: beginner Notes: If `lx`, `ly`, or `lz` are non-null, these must be of length as `m`, `n`, `p` and the corresponding `m`, `n`, or `p` cannot be `PETSC_DECIDE`. Sum of the `lx` entries must be `M`, sum of the `ly` must `N`, sum of the `lz` must be `P`. The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the standard 7-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes the standard 27-pt stencil. The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects; The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()` and `DMCreateLocalVector()` and calls to `VecDuplicate()` if more are needed. You must call `DMSetUp()` after this call before using this `DM`. To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call but before `DMSetUp()`. .seealso: [](sec_struct), `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, `DMStagCreate3d()`, `DMBoundaryType` @*/ 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) { PetscFunctionBegin; PetscCall(DMDACreate(comm, da)); PetscCall(DMSetDimension(*da, 3)); PetscCall(DMDASetSizes(*da, M, N, P)); PetscCall(DMDASetNumProcs(*da, m, n, p)); PetscCall(DMDASetBoundaryType(*da, bx, by, bz)); PetscCall(DMDASetDof(*da, dof)); PetscCall(DMDASetStencilType(*da, stencil_type)); PetscCall(DMDASetStencilWidth(*da, s)); PetscCall(DMDASetOwnershipRanges(*da, lx, ly, lz)); PetscFunctionReturn(PETSC_SUCCESS); }