Actual source code: stag1d.c

      1: /*
      2:    Functions specific to the 1-dimensional implementation of DMStag
      3: */
      4: #include <petsc/private/dmstagimpl.h>

      6: /*@
      7:   DMStagCreate1d - Create an object to manage data living on the elements and vertices of a parallelized regular 1D grid.

      9:   Collective

     11:   Input Parameters:
     12: + comm         - MPI communicator
     13: . bndx         - boundary type: `DM_BOUNDARY_NONE`, `DM_BOUNDARY_PERIODIC`, or `DM_BOUNDARY_GHOSTED`
     14: . M            - global number of elements
     15: . dof0         - number of degrees of freedom per vertex/0-cell
     16: . dof1         - number of degrees of freedom per element/1-cell
     17: . stencilType  - ghost/halo region type: `DMSTAG_STENCIL_BOX` or `DMSTAG_STENCIL_NONE`
     18: . stencilWidth - width, in elements, of halo/ghost region
     19: - lx           - array of local sizes, of length equal to the comm size, summing to `M` or `NULL`

     21:   Output Parameter:
     22: . dm - the new `DMSTAG` object

     24:   Options Database Keys:
     25: + -dm_view                                      - calls `DMViewFromOptions()` at the conclusion of `DMSetUp()`
     26: . -stag_grid_x nx                               - number of elements in the x direction
     27: . -stag_ghost_stencil_width                     - width of ghost region, in elements
     28: - -stag_boundary_type_x (none|ghosted|periodic) - `DMBoundaryType` value

     30:   Level: beginner

     32:   Notes:
     33:   You must call `DMSetUp()` after this call before using the `DM`.
     34:   If you wish to use the options database (see the keys above) to change values in the `DMSTAG`, you must call
     35:   `DMSetFromOptions()` after this function but before `DMSetUp()`.

     37: .seealso: [](ch_stag), `DMSTAG`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDestroy()`, `DMView()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMLocalToGlobalBegin()`, `DMDACreate1d()`
     38: @*/
     39: PetscErrorCode DMStagCreate1d(MPI_Comm comm, DMBoundaryType bndx, PetscInt M, PetscInt dof0, PetscInt dof1, DMStagStencilType stencilType, PetscInt stencilWidth, const PetscInt lx[], DM *dm)
     40: {
     41:   PetscMPIInt size;

     43:   PetscFunctionBegin;
     44:   PetscCallMPI(MPI_Comm_size(comm, &size));
     45:   PetscCall(DMCreate(comm, dm));
     46:   PetscCall(DMSetDimension(*dm, 1));
     47:   PetscCall(DMStagInitialize(bndx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, M, 0, 0, size, 0, 0, dof0, dof1, 0, 0, stencilType, stencilWidth, lx, NULL, NULL, *dm));
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: PETSC_INTERN PetscErrorCode DMStagRestrictSimple_1d(DM dmf, Vec xf_local, DM dmc, Vec xc_local)
     52: {
     53:   PetscInt            Mf, Mc, factorx, dof[2];
     54:   PetscInt            xc, mc, nExtraxc, i, d;
     55:   PetscInt            ileftf, ielemf, ileftc, ielemc;
     56:   const PetscScalar **arrf;
     57:   PetscScalar       **arrc;

     59:   PetscFunctionBegin;
     60:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, NULL, NULL));
     61:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, NULL, NULL));
     62:   factorx = Mf / Mc;
     63:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));

     65:   PetscCall(DMStagGetCorners(dmc, &xc, NULL, NULL, &mc, NULL, NULL, &nExtraxc, NULL, NULL));
     66:   PetscCall(VecZeroEntries(xc_local));
     67:   PetscCall(DMStagVecGetArray(dmf, xf_local, &arrf));
     68:   PetscCall(DMStagVecGetArray(dmc, xc_local, &arrc));
     69:   PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_LEFT, 0, &ileftf));
     70:   PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_ELEMENT, 0, &ielemf));
     71:   PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_LEFT, 0, &ileftc));
     72:   PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_ELEMENT, 0, &ielemc));

     74:   for (d = 0; d < dof[0]; ++d)
     75:     for (i = xc; i < xc + mc + nExtraxc; ++i) arrc[i][ileftc + d] = arrf[factorx * i][ileftf + d];

     77:   for (d = 0; d < dof[1]; ++d)
     78:     for (i = xc; i < xc + mc; ++i) {
     79:       if (factorx % 2 == 0) arrc[i][ielemc + d] = 0.5 * (arrf[factorx * i + factorx / 2 - 1][ielemf + d] + arrf[factorx * i + factorx / 2][ielemf + d]);
     80:       else arrc[i][ielemc + d] = arrf[factorx * i + factorx / 2][ielemf + d];
     81:     }

     83:   PetscCall(DMStagVecRestoreArray(dmf, xf_local, &arrf));
     84:   PetscCall(DMStagVecRestoreArray(dmc, xc_local, &arrc));
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm, PetscReal xmin, PetscReal xmax)
     89: {
     90:   DM_Stag      *stagCoord;
     91:   DM            dmCoord;
     92:   Vec           coordLocal;
     93:   PetscReal     h, min;
     94:   PetscScalar **arr;
     95:   PetscInt      start_ghost, n_ghost, s;
     96:   PetscInt      ileft, ielement;

     98:   PetscFunctionBegin;
     99:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
    100:   stagCoord = (DM_Stag *)dmCoord->data;
    101:   for (s = 0; s < 2; ++s) {
    102:     PetscCheck(stagCoord->dof[s] == 0 || stagCoord->dof[s] == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Coordinate DM in 1 dimensions must have 0 or 1 dof on each stratum, but stratum %" PetscInt_FMT " has %" PetscInt_FMT " dof", s,
    103:                stagCoord->dof[s]);
    104:   }
    105:   PetscCall(DMCreateLocalVector(dmCoord, &coordLocal));

    107:   PetscCall(DMStagVecGetArray(dmCoord, coordLocal, &arr));
    108:   if (stagCoord->dof[0]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_LEFT, 0, &ileft));
    109:   if (stagCoord->dof[1]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_ELEMENT, 0, &ielement));
    110:   PetscCall(DMStagGetGhostCorners(dmCoord, &start_ghost, NULL, NULL, &n_ghost, NULL, NULL));

    112:   min = xmin;
    113:   h   = (xmax - xmin) / stagCoord->N[0];

    115:   for (PetscInt ind = start_ghost; ind < start_ghost + n_ghost; ++ind) {
    116:     if (stagCoord->dof[0]) {
    117:       const PetscReal off = 0.0;
    118:       arr[ind][ileft]     = min + ((PetscReal)ind + off) * h;
    119:     }
    120:     if (stagCoord->dof[1]) {
    121:       const PetscReal off = 0.5;
    122:       arr[ind][ielement]  = min + ((PetscReal)ind + off) * h;
    123:     }
    124:   }
    125:   PetscCall(DMStagVecRestoreArray(dmCoord, coordLocal, &arr));
    126:   PetscCall(DMSetCoordinatesLocal(dm, coordLocal));
    127:   PetscCall(VecDestroy(&coordLocal));
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: /* Helper functions used in DMSetUp_Stag() */
    132: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);

    134: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
    135: {
    136:   DM_Stag *const stag = (DM_Stag *)dm->data;
    137:   PetscMPIInt    size, rank;
    138:   MPI_Comm       comm;
    139:   PetscInt       j;

    141:   PetscFunctionBegin;
    142:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    143:   PetscCallMPI(MPI_Comm_size(comm, &size));
    144:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    146:   /* Check Global size */
    147:   PetscCheck(stag->N[0] >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Global grid size of %" PetscInt_FMT " < 1 specified", stag->N[0]);

    149:   /* Local sizes */
    150:   PetscCheck(stag->N[0] >= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "More ranks (%d) than elements (%" PetscInt_FMT ") specified", size, stag->N[0]);
    151:   if (!stag->l[0]) {
    152:     /* Divide equally, giving an extra elements to higher ranks */
    153:     PetscCall(PetscMalloc1(stag->nRanks[0], &stag->l[0]));
    154:     for (j = 0; j < stag->nRanks[0]; ++j) stag->l[0][j] = stag->N[0] / stag->nRanks[0] + (stag->N[0] % stag->nRanks[0] > j ? 1 : 0);
    155:   }
    156:   {
    157:     PetscInt Nchk = 0;
    158:     for (j = 0; j < size; ++j) Nchk += stag->l[0][j];
    159:     PetscCheck(Nchk == stag->N[0], comm, PETSC_ERR_ARG_OUTOFRANGE, "Sum of specified local sizes (%" PetscInt_FMT ") is not equal to global size (%" PetscInt_FMT ")", Nchk, stag->N[0]);
    160:   }
    161:   stag->n[0] = stag->l[0][rank];

    163:   /* Rank (trivial in 1d) */
    164:   stag->rank[0]      = rank;
    165:   stag->firstRank[0] = (PetscBool)(rank == 0);
    166:   stag->lastRank[0]  = (PetscBool)(rank == size - 1);

    168:   /* Local (unghosted) numbers of entries */
    169:   stag->entriesPerElement = stag->dof[0] + stag->dof[1];
    170:   switch (stag->boundaryType[0]) {
    171:   case DM_BOUNDARY_NONE:
    172:   case DM_BOUNDARY_GHOSTED:
    173:     stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
    174:     break;
    175:   case DM_BOUNDARY_PERIODIC:
    176:     stag->entries = stag->n[0] * stag->entriesPerElement;
    177:     break;
    178:   default:
    179:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
    180:   }

    182:   /* Starting element */
    183:   stag->start[0] = 0;
    184:   for (j = 0; j < stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];

    186:   /* Local/ghosted size and starting element */
    187:   switch (stag->boundaryType[0]) {
    188:   case DM_BOUNDARY_NONE:
    189:     switch (stag->stencilType) {
    190:     case DMSTAG_STENCIL_NONE: /* Only dummy cells on the right */
    191:       stag->startGhost[0] = stag->start[0];
    192:       stag->nGhost[0]     = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
    193:       break;
    194:     case DMSTAG_STENCIL_STAR:
    195:     case DMSTAG_STENCIL_BOX:
    196:       stag->startGhost[0] = stag->firstRank[0] ? stag->start[0] : stag->start[0] - stag->stencilWidth;
    197:       stag->nGhost[0]     = stag->n[0];
    198:       stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
    199:       stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
    200:       break;
    201:     default:
    202:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
    203:     }
    204:     break;
    205:   case DM_BOUNDARY_GHOSTED:
    206:     switch (stag->stencilType) {
    207:     case DMSTAG_STENCIL_NONE:
    208:       stag->startGhost[0] = stag->start[0];
    209:       stag->nGhost[0]     = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
    210:       break;
    211:     case DMSTAG_STENCIL_STAR:
    212:     case DMSTAG_STENCIL_BOX:
    213:       stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
    214:       stag->nGhost[0]     = stag->n[0] + 2 * stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
    215:       break;
    216:     default:
    217:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
    218:     }
    219:     break;
    220:   case DM_BOUNDARY_PERIODIC:
    221:     switch (stag->stencilType) {
    222:     case DMSTAG_STENCIL_NONE:
    223:       stag->startGhost[0] = stag->start[0];
    224:       stag->nGhost[0]     = stag->n[0];
    225:       break;
    226:     case DMSTAG_STENCIL_STAR:
    227:     case DMSTAG_STENCIL_BOX:
    228:       stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
    229:       stag->nGhost[0]     = stag->n[0] + 2 * stag->stencilWidth;
    230:       break;
    231:     default:
    232:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
    233:     }
    234:     break;
    235:   default:
    236:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
    237:   }

    239:   /* Total size of ghosted/local representation */
    240:   stag->entriesGhost = stag->nGhost[0] * stag->entriesPerElement;

    242:   /* Define neighbors */
    243:   PetscCall(PetscMalloc1(3, &stag->neighbors));
    244:   if (stag->firstRank[0]) {
    245:     switch (stag->boundaryType[0]) {
    246:     case DM_BOUNDARY_GHOSTED:
    247:     case DM_BOUNDARY_NONE:
    248:       stag->neighbors[0] = -1;
    249:       break;
    250:     case DM_BOUNDARY_PERIODIC:
    251:       stag->neighbors[0] = stag->nRanks[0] - 1;
    252:       break;
    253:     default:
    254:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
    255:     }
    256:   } else {
    257:     stag->neighbors[0] = stag->rank[0] - 1;
    258:   }
    259:   stag->neighbors[1] = stag->rank[0];
    260:   if (stag->lastRank[0]) {
    261:     switch (stag->boundaryType[0]) {
    262:     case DM_BOUNDARY_GHOSTED:
    263:     case DM_BOUNDARY_NONE:
    264:       stag->neighbors[2] = -1;
    265:       break;
    266:     case DM_BOUNDARY_PERIODIC:
    267:       stag->neighbors[2] = 0;
    268:       break;
    269:     default:
    270:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
    271:     }
    272:   } else {
    273:     stag->neighbors[2] = stag->rank[0] + 1;
    274:   }

    276:   PetscCheck(stag->n[0] >= stag->stencilWidth, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMStag 1d setup does not support local sizes (%" PetscInt_FMT ") smaller than the elementwise stencil width (%" PetscInt_FMT ")", stag->n[0], stag->stencilWidth);

    278:   /* Create global->local VecScatter and ISLocalToGlobalMapping */
    279:   {
    280:     PetscInt *idxLocal, *idxGlobal, *idxGlobalAll;
    281:     PetscInt  i, iLocal, d, entriesToTransferTotal, ghostOffsetStart, ghostOffsetEnd, nNonDummyGhost;
    282:     IS        isLocal, isGlobal;

    284:     /* The offset on the right (may not be equal to the stencil width, as we
    285:        always have at least one ghost element, to account for the boundary
    286:        point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
    287:     ghostOffsetStart = stag->start[0] - stag->startGhost[0];
    288:     ghostOffsetEnd   = stag->startGhost[0] + stag->nGhost[0] - (stag->start[0] + stag->n[0]);
    289:     nNonDummyGhost   = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);

    291:     /* Compute the number of non-dummy entries in the local representation
    292:        This is equal to the number of non-dummy elements in the local (ghosted) representation,
    293:        plus some extra entries on the right boundary on the last rank*/
    294:     switch (stag->boundaryType[0]) {
    295:     case DM_BOUNDARY_GHOSTED:
    296:     case DM_BOUNDARY_NONE:
    297:       entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
    298:       break;
    299:     case DM_BOUNDARY_PERIODIC:
    300:       entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
    301:       break;
    302:     default:
    303:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
    304:     }

    306:     PetscCall(PetscMalloc1(entriesToTransferTotal, &idxLocal));
    307:     PetscCall(PetscMalloc1(entriesToTransferTotal, &idxGlobal));
    308:     PetscCall(PetscMalloc1(stag->entriesGhost, &idxGlobalAll));
    309:     if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
    310:       PetscInt count = 0, countAll = 0;
    311:       /* Left ghost points and native points */
    312:       for (i = stag->startGhost[0], iLocal = 0; iLocal < nNonDummyGhost; ++i, ++iLocal) {
    313:         for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
    314:           idxLocal[count]        = iLocal * stag->entriesPerElement + d;
    315:           idxGlobal[count]       = i * stag->entriesPerElement + d;
    316:           idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
    317:         }
    318:       }
    319:       /* Ghost points on the right
    320:          Special case for last (partial dummy) element on the last rank */
    321:       if (stag->lastRank[0]) {
    322:         i      = stag->N[0];
    323:         iLocal = (stag->nGhost[0] - ghostOffsetEnd);
    324:         /* Only vertex (0-cell) dofs in global representation */
    325:         for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) {
    326:           idxGlobal[count]       = i * stag->entriesPerElement + d;
    327:           idxLocal[count]        = iLocal * stag->entriesPerElement + d;
    328:           idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
    329:         }
    330:         for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
    331:           idxGlobalAll[countAll] = -1;
    332:         }
    333:       }
    334:     } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
    335:       PetscInt       count = 0, iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
    336:       const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
    337:       const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
    338:       /* Ghost points on the left */
    339:       if (stag->firstRank[0]) {
    340:         for (i = stag->N[0] - stag->stencilWidth; iLocal < stag->stencilWidth; ++i, ++iLocal) {
    341:           for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
    342:             idxGlobal[count]    = i * stag->entriesPerElement + d;
    343:             idxLocal[count]     = iLocal * stag->entriesPerElement + d;
    344:             idxGlobalAll[count] = idxGlobal[count];
    345:           }
    346:         }
    347:       }
    348:       /* Native points */
    349:       for (i = iMin; i < iMax; ++i, ++iLocal) {
    350:         for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
    351:           idxGlobal[count]    = i * stag->entriesPerElement + d;
    352:           idxLocal[count]     = iLocal * stag->entriesPerElement + d;
    353:           idxGlobalAll[count] = idxGlobal[count];
    354:         }
    355:       }
    356:       /* Ghost points on the right */
    357:       if (stag->lastRank[0]) {
    358:         for (i = 0; iLocal < stag->nGhost[0]; ++i, ++iLocal) {
    359:           for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
    360:             idxGlobal[count]    = i * stag->entriesPerElement + d;
    361:             idxLocal[count]     = iLocal * stag->entriesPerElement + d;
    362:             idxGlobalAll[count] = idxGlobal[count];
    363:           }
    364:         }
    365:       }
    366:     } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
    367:       PetscInt count = 0, countAll = 0;
    368:       /* Dummy elements on the left, on the first rank */
    369:       if (stag->firstRank[0]) {
    370:         for (iLocal = 0; iLocal < ghostOffsetStart; ++iLocal) {
    371:           /* Complete elements full of dummy entries */
    372:           for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
    373:         }
    374:         i = 0; /* nonDummy entries start with global entry 0 */
    375:       } else {
    376:         /* nonDummy entries start as usual */
    377:         i      = stag->startGhost[0];
    378:         iLocal = 0;
    379:       }

    381:       /* non-Dummy entries */
    382:       {
    383:         PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
    384:         for (; iLocal < iLocalNonDummyMax; ++i, ++iLocal) {
    385:           for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
    386:             idxLocal[count]        = iLocal * stag->entriesPerElement + d;
    387:             idxGlobal[count]       = i * stag->entriesPerElement + d;
    388:             idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
    389:           }
    390:         }
    391:       }

    393:       /* (partial) dummy elements on the right, on the last rank */
    394:       if (stag->lastRank[0]) {
    395:         /* First one is partial dummy */
    396:         i      = stag->N[0];
    397:         iLocal = (stag->nGhost[0] - ghostOffsetEnd);
    398:         for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) { /* Only vertex (0-cell) dofs in global representation */
    399:           idxLocal[count]        = iLocal * stag->entriesPerElement + d;
    400:           idxGlobal[count]       = i * stag->entriesPerElement + d;
    401:           idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
    402:         }
    403:         for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
    404:           idxGlobalAll[countAll] = -1;
    405:         }
    406:         for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
    407:           /* Additional dummy elements */
    408:           for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
    409:         }
    410:       }
    411:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);

    413:     /* Create Local IS (transferring pointer ownership) */
    414:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxLocal, PETSC_OWN_POINTER, &isLocal));

    416:     /* Create Global IS (transferring pointer ownership) */
    417:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxGlobal, PETSC_OWN_POINTER, &isGlobal));

    419:     /* Create stag->gtol, which doesn't include dummy entries */
    420:     {
    421:       Vec local, global;
    422:       PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
    423:       PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, PetscMax(stag->entriesPerElement, 1), stag->entriesGhost, NULL, &local));
    424:       PetscCall(VecScatterCreate(global, isGlobal, local, isLocal, &stag->gtol));
    425:       PetscCall(VecDestroy(&global));
    426:       PetscCall(VecDestroy(&local));
    427:     }

    429:     /* In special cases, create a dedicated injective local-to-global map */
    430:     if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) PetscCall(DMStagPopulateLocalToGlobalInjective(dm));

    432:     /* Destroy ISs */
    433:     PetscCall(ISDestroy(&isLocal));
    434:     PetscCall(ISDestroy(&isGlobal));

    436:     /* Create local-to-global map (transferring pointer ownership) */
    437:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stag->entriesGhost, idxGlobalAll, PETSC_OWN_POINTER, &dm->ltogmap));
    438:   }

    440:   /* Precompute location offsets */
    441:   PetscCall(DMStagComputeLocationOffsets_1d(dm));

    443:   /* View from Options */
    444:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
    449: {
    450:   DM_Stag *const stag = (DM_Stag *)dm->data;
    451:   const PetscInt epe  = stag->entriesPerElement;

    453:   PetscFunctionBegin;
    454:   PetscCall(PetscMalloc1(DMSTAG_NUMBER_LOCATIONS, &stag->locationOffsets));
    455:   stag->locationOffsets[DMSTAG_LEFT]    = 0;
    456:   stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
    457:   stag->locationOffsets[DMSTAG_RIGHT]   = stag->locationOffsets[DMSTAG_LEFT] + epe;
    458:   PetscFunctionReturn(PETSC_SUCCESS);
    459: }

    461: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
    462: {
    463:   DM_Stag *const stag = (DM_Stag *)dm->data;
    464:   PetscInt      *idxLocal, *idxGlobal;
    465:   PetscInt       i, iLocal, d, count;
    466:   IS             isLocal, isGlobal;

    468:   PetscFunctionBegin;
    469:   PetscCall(PetscMalloc1(stag->entries, &idxLocal));
    470:   PetscCall(PetscMalloc1(stag->entries, &idxGlobal));
    471:   count  = 0;
    472:   iLocal = stag->start[0] - stag->startGhost[0];
    473:   for (i = stag->start[0]; i < stag->start[0] + stag->n[0]; ++i, ++iLocal) {
    474:     for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
    475:       idxGlobal[count] = i * stag->entriesPerElement + d;
    476:       idxLocal[count]  = iLocal * stag->entriesPerElement + d;
    477:     }
    478:   }
    479:   if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
    480:     i      = stag->start[0] + stag->n[0];
    481:     iLocal = stag->start[0] - stag->startGhost[0] + stag->n[0];
    482:     for (d = 0; d < stag->dof[0]; ++d, ++count) {
    483:       idxGlobal[count] = i * stag->entriesPerElement + d;
    484:       idxLocal[count]  = iLocal * stag->entriesPerElement + d;
    485:     }
    486:   }
    487:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxLocal, PETSC_OWN_POINTER, &isLocal));
    488:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
    489:   {
    490:     Vec local, global;
    491:     PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
    492:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, PetscMax(stag->entriesPerElement, 1), stag->entriesGhost, NULL, &local));
    493:     PetscCall(VecScatterCreate(local, isLocal, global, isGlobal, &stag->ltog_injective));
    494:     PetscCall(VecDestroy(&global));
    495:     PetscCall(VecDestroy(&local));
    496:   }
    497:   PetscCall(ISDestroy(&isLocal));
    498:   PetscCall(ISDestroy(&isGlobal));
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToLocal1d_Internal(DM dm)
    503: {
    504:   DM_Stag *const stag = (DM_Stag *)dm->data;
    505:   PetscInt      *idxRemap;
    506:   PetscInt       i, leftGhostEntries;

    508:   PetscFunctionBegin;
    509:   PetscCall(VecScatterCopy(stag->gtol, &stag->ltol));
    510:   PetscCall(PetscMalloc1(stag->entries, &idxRemap));

    512:   leftGhostEntries = (stag->start[0] - stag->startGhost[0]) * stag->entriesPerElement;
    513:   for (i = 0; i < stag->entries; ++i) idxRemap[i] = leftGhostEntries + i;

    515:   PetscCall(VecScatterRemap(stag->ltol, idxRemap, NULL));
    516:   PetscCall(PetscFree(idxRemap));
    517:   PetscFunctionReturn(PETSC_SUCCESS);
    518: }

    520: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_1D_AIJ_Assemble(DM dm, Mat A)
    521: {
    522:   DMStagStencilType stencil_type;
    523:   PetscInt          dof[2], start, n, n_extra, stencil_width, N, epe;
    524:   DMBoundaryType    boundary_type_x;

    526:   PetscFunctionBegin;
    527:   PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL));
    528:   PetscCall(DMStagGetStencilType(dm, &stencil_type));
    529:   PetscCall(DMStagGetStencilWidth(dm, &stencil_width));
    530:   PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL));
    531:   PetscCall(DMStagGetGlobalSizes(dm, &N, NULL, NULL));
    532:   PetscCall(DMStagGetEntriesPerElement(dm, &epe));
    533:   PetscCall(DMStagGetBoundaryTypes(dm, &boundary_type_x, NULL, NULL));
    534:   if (stencil_type == DMSTAG_STENCIL_NONE) {
    535:     /* Couple all DOF at each location to each other */
    536:     DMStagStencil *row_vertex, *row_element;

    538:     PetscCall(PetscMalloc1(dof[0], &row_vertex));
    539:     for (PetscInt c = 0; c < dof[0]; ++c) {
    540:       row_vertex[c].loc = DMSTAG_LEFT;
    541:       row_vertex[c].c   = c;
    542:     }

    544:     PetscCall(PetscMalloc1(dof[1], &row_element));
    545:     for (PetscInt c = 0; c < dof[1]; ++c) {
    546:       row_element[c].loc = DMSTAG_ELEMENT;
    547:       row_element[c].c   = c;
    548:     }

    550:     for (PetscInt e = start; e < start + n + n_extra; ++e) {
    551:       {
    552:         for (PetscInt c = 0; c < dof[0]; ++c) row_vertex[c].i = e;
    553:         PetscCall(DMStagMatSetValuesStencil(dm, A, dof[0], row_vertex, dof[0], row_vertex, NULL, INSERT_VALUES));
    554:       }
    555:       if (e < N) {
    556:         for (PetscInt c = 0; c < dof[1]; ++c) row_element[c].i = e;
    557:         PetscCall(DMStagMatSetValuesStencil(dm, A, dof[1], row_element, dof[1], row_element, NULL, INSERT_VALUES));
    558:       }
    559:     }
    560:     PetscCall(PetscFree(row_vertex));
    561:     PetscCall(PetscFree(row_element));
    562:   } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
    563:     DMStagStencil *col, *row;

    565:     PetscCall(PetscMalloc1(epe, &row));
    566:     {
    567:       PetscInt nrows = 0;
    568:       for (PetscInt c = 0; c < dof[0]; ++c) {
    569:         row[nrows].c   = c;
    570:         row[nrows].loc = DMSTAG_LEFT;
    571:         ++nrows;
    572:       }
    573:       for (PetscInt c = 0; c < dof[1]; ++c) {
    574:         row[nrows].c   = c;
    575:         row[nrows].loc = DMSTAG_ELEMENT;
    576:         ++nrows;
    577:       }
    578:     }
    579:     PetscCall(PetscMalloc1(epe, &col));
    580:     {
    581:       PetscInt ncols = 0;
    582:       for (PetscInt c = 0; c < dof[0]; ++c) {
    583:         col[ncols].c   = c;
    584:         col[ncols].loc = DMSTAG_LEFT;
    585:         ++ncols;
    586:       }
    587:       for (PetscInt c = 0; c < dof[1]; ++c) {
    588:         col[ncols].c   = c;
    589:         col[ncols].loc = DMSTAG_ELEMENT;
    590:         ++ncols;
    591:       }
    592:     }
    593:     for (PetscInt e = start; e < start + n + n_extra; ++e) {
    594:       for (PetscInt i = 0; i < epe; ++i) row[i].i = e;
    595:       for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
    596:         const PetscInt e_offset = e + offset;

    598:         /* Only set values corresponding to elements which can have non-dummy entries,
    599:            meaning those that map to unknowns in the global representation. In the periodic
    600:            case, this is the entire stencil, but in all other cases, only includes a single
    601:            "extra" element which is partially outside the physical domain (those points in the
    602:            global representation */
    603:         if (boundary_type_x == DM_BOUNDARY_PERIODIC || (e_offset < N + 1 && e_offset >= 0)) {
    604:           for (PetscInt i = 0; i < epe; ++i) col[i].i = e_offset;
    605:           PetscCall(DMStagMatSetValuesStencil(dm, A, epe, row, epe, col, NULL, INSERT_VALUES));
    606:         }
    607:       }
    608:     }
    609:     PetscCall(PetscFree(row));
    610:     PetscCall(PetscFree(col));
    611:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported stencil type %s", DMStagStencilTypes[stencil_type]);
    612:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    613:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    614:   PetscFunctionReturn(PETSC_SUCCESS);
    615: }