Actual source code: stagstencil.c

      1: /* Functions concerning getting and setting Vec and Mat values with DMStagStencil */
      2: #include <petsc/private/dmstagimpl.h>

      4: /* Strings corresponding to the types defined in $PETSC_DIR/include/petscdmstag.h */
      5: const char *const DMStagStencilTypes[] = {"NONE", "STAR", "BOX", "DMStagStencilType", "DM_STAG_STENCIL_", NULL};

      7: /* Strings corresponding the positions in $PETSC_DIR/include/petscdmstag.h */
      8: const char *const DMStagStencilLocations[] = {"NONE", "BACK_DOWN_LEFT", "BACK_DOWN", "BACK_DOWN_RIGHT", "BACK_LEFT", "BACK", "BACK_RIGHT", "BACK_UP_LEFT", "BACK_UP", "BACK_UP_RIGHT", "DOWN_LEFT", "DOWN", "DOWN_RIGHT", "LEFT", "ELEMENT", "RIGHT", "UP_LEFT", "UP", "UP_RIGHT", "FRONT_DOWN_LEFT", "FRONT_DOWN", "FRONT_DOWN_RIGHT", "FRONT_LEFT", "FRONT", "FRONT_RIGHT", "FRONT_UP_LEFT", "FRONT_UP", "FRONT_UP_RIGHT", "DMStagStencilLocation", "", NULL};

     10: /*@C
     11:   DMStagCreateISFromStencils - Create an `IS`, using global numberings, for a subset of DOF in a `DMSTAG` object

     13:   Collective

     15:   Input Parameters:
     16: + dm        - the `DMSTAG` object
     17: . n_stencil - the number of stencils provided
     18: - stencils  - an array of `DMStagStencil` objects (`i`, `j`, and `k` are ignored)

     20:   Output Parameter:
     21: . is - the global `IS`

     23:   Note:
     24:   Redundant entries in the stencils argument are ignored

     26:   Level: advanced

     28: .seealso: [](ch_stag), `DMSTAG`, `IS`, `DMStagStencil`, `DMCreateGlobalVector`
     29: @*/
     30: PetscErrorCode DMStagCreateISFromStencils(DM dm, PetscInt n_stencil, DMStagStencil stencils[], IS *is)
     31: {
     32:   PetscInt              *stencil_active;
     33:   DMStagStencil         *stencils_ordered_unique;
     34:   PetscInt              *idx, *idxLocal;
     35:   const PetscInt        *ltogidx;
     36:   PetscInt               n_stencil_unique, dim, count, nidx, nc_max;
     37:   ISLocalToGlobalMapping ltog;
     38:   PetscInt               start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM];

     40:   PetscFunctionBegin;
     41:   PetscCall(DMGetDimension(dm, &dim));
     42:   PetscCheck(dim >= 1 && dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);

     44:   /* To assert that the resulting IS has unique, sorted, entries, we perform
     45:      a bucket sort, taking advantage of the fact that DMStagStencilLocation
     46:      enum values are integers starting with 1, in canonical order */
     47:   nc_max           = 1; // maximum number of components to represent these stencils
     48:   n_stencil_unique = 0;
     49:   for (PetscInt p = 0; p < n_stencil; ++p) nc_max = PetscMax(nc_max, stencils[p].c + 1);
     50:   PetscCall(PetscCalloc1(DMSTAG_NUMBER_LOCATIONS * nc_max, &stencil_active));
     51:   for (PetscInt p = 0; p < n_stencil; ++p) {
     52:     DMStagStencilLocation loc_canonical;
     53:     PetscInt              slot;

     55:     PetscCall(DMStagStencilLocationCanonicalize(stencils[p].loc, &loc_canonical));
     56:     slot = nc_max * ((PetscInt)loc_canonical) + stencils[p].c;
     57:     if (stencil_active[slot] == 0) {
     58:       stencil_active[slot] = 1;
     59:       ++n_stencil_unique;
     60:     }
     61:   }
     62:   PetscCall(PetscMalloc1(n_stencil_unique, &stencils_ordered_unique));
     63:   {
     64:     PetscInt p = 0;

     66:     for (PetscInt i = 1; i < DMSTAG_NUMBER_LOCATIONS; ++i) {
     67:       for (PetscInt c = 0; c < nc_max; ++c) {
     68:         if (stencil_active[nc_max * i + c] != 0) {
     69:           stencils_ordered_unique[p].loc = (DMStagStencilLocation)i;
     70:           stencils_ordered_unique[p].c   = c;
     71:           ++p;
     72:         }
     73:       }
     74:     }
     75:   }
     76:   PetscCall(PetscFree(stencil_active));

     78:   PetscCall(PetscMalloc1(n_stencil_unique, &idxLocal));
     79:   PetscCall(DMGetLocalToGlobalMapping(dm, &ltog));
     80:   PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &ltogidx));
     81:   PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &extraPoint[0], &extraPoint[1], &extraPoint[2]));
     82:   for (PetscInt d = dim; d < DMSTAG_MAX_DIM; ++d) {
     83:     start[d]      = 0;
     84:     n[d]          = 1; /* To allow for a single loop nest below */
     85:     extraPoint[d] = 0;
     86:   }
     87:   nidx = n_stencil_unique;
     88:   for (PetscInt d = 0; d < dim; ++d) nidx *= (n[d] + 1); /* Overestimate (always assumes extraPoint) */
     89:   PetscCall(PetscMalloc1(nidx, &idx));
     90:   count = 0;
     91:   /* Note that unused loop variables are not accessed, for lower dimensions */
     92:   for (PetscInt k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
     93:     for (PetscInt j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
     94:       for (PetscInt i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
     95:         for (PetscInt p = 0; p < n_stencil_unique; ++p) {
     96:           stencils_ordered_unique[p].i = i;
     97:           stencils_ordered_unique[p].j = j;
     98:           stencils_ordered_unique[p].k = k;
     99:         }
    100:         PetscCall(DMStagStencilToIndexLocal(dm, dim, n_stencil_unique, stencils_ordered_unique, idxLocal));
    101:         for (PetscInt p = 0; p < n_stencil_unique; ++p) {
    102:           const PetscInt gidx = ltogidx[idxLocal[p]];
    103:           if (gidx >= 0) {
    104:             idx[count] = gidx;
    105:             ++count;
    106:           }
    107:         }
    108:       }
    109:     }
    110:   }

    112:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &ltogidx));
    113:   PetscCall(PetscFree(stencils_ordered_unique));
    114:   PetscCall(PetscFree(idxLocal));

    116:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), count, idx, PETSC_OWN_POINTER, is));
    117:   PetscCall(ISSetInfo(*is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
    118:   PetscCall(ISSetInfo(*is, IS_UNIQUE, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: /*@C
    123:   DMStagGetLocationDOF - Get number of DOF associated with a given point in a `DMSTAG` grid

    125:   Not Collective

    127:   Input Parameters:
    128: + dm  - the `DMSTAG` object
    129: - loc - grid point (see `DMStagStencilLocation`)

    131:   Output Parameter:
    132: . dof - the number of DOF (components) living at `loc` in `dm`

    134:   Level: intermediate

    136: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMDAGetDof()`
    137: @*/
    138: PetscErrorCode DMStagGetLocationDOF(DM dm, DMStagStencilLocation loc, PetscInt *dof)
    139: {
    140:   const DM_Stag *const stag = (DM_Stag *)dm->data;
    141:   PetscInt             dim;

    143:   PetscFunctionBegin;
    145:   PetscCall(DMGetDimension(dm, &dim));
    146:   switch (dim) {
    147:   case 1:
    148:     switch (loc) {
    149:     case DMSTAG_LEFT:
    150:     case DMSTAG_RIGHT:
    151:       *dof = stag->dof[0];
    152:       break;
    153:     case DMSTAG_ELEMENT:
    154:       *dof = stag->dof[1];
    155:       break;
    156:     default:
    157:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
    158:     }
    159:     break;
    160:   case 2:
    161:     switch (loc) {
    162:     case DMSTAG_DOWN_LEFT:
    163:     case DMSTAG_DOWN_RIGHT:
    164:     case DMSTAG_UP_LEFT:
    165:     case DMSTAG_UP_RIGHT:
    166:       *dof = stag->dof[0];
    167:       break;
    168:     case DMSTAG_LEFT:
    169:     case DMSTAG_RIGHT:
    170:     case DMSTAG_UP:
    171:     case DMSTAG_DOWN:
    172:       *dof = stag->dof[1];
    173:       break;
    174:     case DMSTAG_ELEMENT:
    175:       *dof = stag->dof[2];
    176:       break;
    177:     default:
    178:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
    179:     }
    180:     break;
    181:   case 3:
    182:     switch (loc) {
    183:     case DMSTAG_BACK_DOWN_LEFT:
    184:     case DMSTAG_BACK_DOWN_RIGHT:
    185:     case DMSTAG_BACK_UP_LEFT:
    186:     case DMSTAG_BACK_UP_RIGHT:
    187:     case DMSTAG_FRONT_DOWN_LEFT:
    188:     case DMSTAG_FRONT_DOWN_RIGHT:
    189:     case DMSTAG_FRONT_UP_LEFT:
    190:     case DMSTAG_FRONT_UP_RIGHT:
    191:       *dof = stag->dof[0];
    192:       break;
    193:     case DMSTAG_BACK_DOWN:
    194:     case DMSTAG_BACK_LEFT:
    195:     case DMSTAG_BACK_RIGHT:
    196:     case DMSTAG_BACK_UP:
    197:     case DMSTAG_DOWN_LEFT:
    198:     case DMSTAG_DOWN_RIGHT:
    199:     case DMSTAG_UP_LEFT:
    200:     case DMSTAG_UP_RIGHT:
    201:     case DMSTAG_FRONT_DOWN:
    202:     case DMSTAG_FRONT_LEFT:
    203:     case DMSTAG_FRONT_RIGHT:
    204:     case DMSTAG_FRONT_UP:
    205:       *dof = stag->dof[1];
    206:       break;
    207:     case DMSTAG_LEFT:
    208:     case DMSTAG_RIGHT:
    209:     case DMSTAG_DOWN:
    210:     case DMSTAG_UP:
    211:     case DMSTAG_BACK:
    212:     case DMSTAG_FRONT:
    213:       *dof = stag->dof[2];
    214:       break;
    215:     case DMSTAG_ELEMENT:
    216:       *dof = stag->dof[3];
    217:       break;
    218:     default:
    219:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
    220:     }
    221:     break;
    222:   default:
    223:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
    224:   }
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: /*
    229: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved
    230: */
    231: PETSC_INTERN PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc, DMStagStencilLocation *locCanonical)
    232: {
    233:   PetscFunctionBegin;
    234:   switch (loc) {
    235:   case DMSTAG_ELEMENT:
    236:     *locCanonical = DMSTAG_ELEMENT;
    237:     break;
    238:   case DMSTAG_LEFT:
    239:   case DMSTAG_RIGHT:
    240:     *locCanonical = DMSTAG_LEFT;
    241:     break;
    242:   case DMSTAG_DOWN:
    243:   case DMSTAG_UP:
    244:     *locCanonical = DMSTAG_DOWN;
    245:     break;
    246:   case DMSTAG_BACK:
    247:   case DMSTAG_FRONT:
    248:     *locCanonical = DMSTAG_BACK;
    249:     break;
    250:   case DMSTAG_DOWN_LEFT:
    251:   case DMSTAG_DOWN_RIGHT:
    252:   case DMSTAG_UP_LEFT:
    253:   case DMSTAG_UP_RIGHT:
    254:     *locCanonical = DMSTAG_DOWN_LEFT;
    255:     break;
    256:   case DMSTAG_BACK_LEFT:
    257:   case DMSTAG_BACK_RIGHT:
    258:   case DMSTAG_FRONT_LEFT:
    259:   case DMSTAG_FRONT_RIGHT:
    260:     *locCanonical = DMSTAG_BACK_LEFT;
    261:     break;
    262:   case DMSTAG_BACK_DOWN:
    263:   case DMSTAG_BACK_UP:
    264:   case DMSTAG_FRONT_DOWN:
    265:   case DMSTAG_FRONT_UP:
    266:     *locCanonical = DMSTAG_BACK_DOWN;
    267:     break;
    268:   case DMSTAG_BACK_DOWN_LEFT:
    269:   case DMSTAG_BACK_DOWN_RIGHT:
    270:   case DMSTAG_BACK_UP_LEFT:
    271:   case DMSTAG_BACK_UP_RIGHT:
    272:   case DMSTAG_FRONT_DOWN_LEFT:
    273:   case DMSTAG_FRONT_DOWN_RIGHT:
    274:   case DMSTAG_FRONT_UP_LEFT:
    275:   case DMSTAG_FRONT_UP_RIGHT:
    276:     *locCanonical = DMSTAG_BACK_DOWN_LEFT;
    277:     break;
    278:   default:
    279:     *locCanonical = DMSTAG_NULL_LOCATION;
    280:     break;
    281:   }
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: /*@C
    286:   DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing

    288:   Not Collective

    290:   Input Parameters:
    291: + dm     - the `DMSTAG` object
    292: . mat    - the matrix
    293: . nRow   - number of rows
    294: . posRow - grid locations (including components) of rows
    295: . nCol   - number of columns
    296: - posCol - grid locations (including components) of columns

    298:   Output Parameter:
    299: . val - logically two-dimensional array of values

    301:   Level: advanced

    303: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
    304: @*/
    305: PetscErrorCode DMStagMatGetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, PetscScalar *val)
    306: {
    307:   PetscInt  dim;
    308:   PetscInt *ir, *ic;

    310:   PetscFunctionBegin;
    313:   PetscCall(DMGetDimension(dm, &dim));
    314:   PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
    315:   PetscCall(DMStagStencilToIndexLocal(dm, dim, nRow, posRow, ir));
    316:   PetscCall(DMStagStencilToIndexLocal(dm, dim, nCol, posCol, ic));
    317:   PetscCall(MatGetValuesLocal(mat, nRow, ir, nCol, ic, val));
    318:   PetscCall(PetscFree2(ir, ic));
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: /*@C
    323:   DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing

    325:   Not Collective

    327:   Input Parameters:
    328: + dm         - the `DMSTAG` object
    329: . mat        - the matrix
    330: . nRow       - number of rows
    331: . posRow     - grid locations (including components) of rows
    332: . nCol       - number of columns
    333: . posCol     - grid locations (including components) of columns
    334: . val        - logically two-dimensional array of values
    335: - insertMode - `INSERT_VALUES` or `ADD_VALUES`

    337:   Notes:
    338:   See notes for `MatSetValuesStencil()`

    340:   Level: intermediate

    342: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatGetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
    343: @*/
    344: PetscErrorCode DMStagMatSetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, const PetscScalar *val, InsertMode insertMode)
    345: {
    346:   PetscInt *ir, *ic;

    348:   PetscFunctionBegin;
    351:   PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
    352:   PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nRow, posRow, ir));
    353:   PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nCol, posCol, ic));
    354:   PetscCall(MatSetValuesLocal(mat, nRow, ir, nCol, ic, val, insertMode));
    355:   PetscCall(PetscFree2(ir, ic));
    356:   PetscFunctionReturn(PETSC_SUCCESS);
    357: }

    359: /*@C
    360:   DMStagStencilToIndexLocal - Convert an array of `DMStagStenci`l objects to an array of indices into a local vector.

    362:   Not Collective

    364:   Input Parameters:
    365: + dm  - the `DMSTAG` object
    366: . dim - the dimension of the `DMSTAG` object
    367: . n   - the number of `DMStagStencil` objects
    368: - pos - an array of `n` `DMStagStencil` objects

    370:   Output Parameter:
    371: . ix - output array of `n` indices

    373:   Notes:
    374:   The `DMStagStencil` objects in `pos` use global element indices.

    376:   The `.c` fields in `pos` must always be set (even if to `0`).

    378:   Developer Notes:
    379:   This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.

    381:   Level: developer

    383: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMGetLocalVector`, `DMCreateLocalVector`
    384: @*/
    385: PetscErrorCode DMStagStencilToIndexLocal(DM dm, PetscInt dim, PetscInt n, const DMStagStencil *pos, PetscInt *ix)
    386: {
    387:   const DM_Stag *const stag = (DM_Stag *)dm->data;
    388:   const PetscInt       epe  = stag->entriesPerElement;

    390:   PetscFunctionBeginHot;
    391:   if (dim == 1) {
    392:     for (PetscInt idx = 0; idx < n; ++idx) {
    393:       const PetscInt eLocal = pos[idx].i - stag->startGhost[0];

    395:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
    396:     }
    397:   } else if (dim == 2) {
    398:     const PetscInt epr = stag->nGhost[0];

    400:     for (PetscInt idx = 0; idx < n; ++idx) {
    401:       const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
    402:       const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
    403:       const PetscInt eLocal  = eLocalx + epr * eLocaly;

    405:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
    406:     }
    407:   } else if (dim == 3) {
    408:     const PetscInt epr = stag->nGhost[0];
    409:     const PetscInt epl = stag->nGhost[0] * stag->nGhost[1];

    411:     for (PetscInt idx = 0; idx < n; ++idx) {
    412:       const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
    413:       const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
    414:       const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
    415:       const PetscInt eLocal  = epl * eLocalz + epr * eLocaly + eLocalx;

    417:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
    418:     }
    419:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    420:   PetscFunctionReturn(PETSC_SUCCESS);
    421: }

    423: /*@C
    424:   DMStagVecGetValuesStencil - get vector values using grid indexing

    426:   Not Collective

    428:   Input Parameters:
    429: + dm  - the `DMSTAG` object
    430: . vec - the vector object
    431: . n   - the number of values to obtain
    432: - pos - locations to obtain values from (as an array of `DMStagStencil` values)

    434:   Output Parameter:
    435: . val - value at the point

    437:   Notes:
    438:   Accepts stencils which refer to global element numbers, but
    439:   only allows access to entries in the local representation (including ghosts).

    441:   This approach is not as efficient as getting values directly with `DMStagVecGetArray()`,
    442:   which is recommended for matrix-free operators.

    444:   Level: advanced

    446: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMStagVecGetArray()`
    447: @*/
    448: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, PetscScalar *val)
    449: {
    450:   DM_Stag *const     stag = (DM_Stag *)dm->data;
    451:   PetscInt           nLocal, idx;
    452:   PetscInt          *ix;
    453:   PetscScalar const *arr;

    455:   PetscFunctionBegin;
    458:   PetscCall(VecGetLocalSize(vec, &nLocal));
    459:   PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector should be a local vector. Local size %" PetscInt_FMT " does not match expected %" PetscInt_FMT, nLocal, stag->entriesGhost);
    460:   PetscCall(PetscMalloc1(n, &ix));
    461:   PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
    462:   PetscCall(VecGetArrayRead(vec, &arr));
    463:   for (idx = 0; idx < n; ++idx) val[idx] = arr[ix[idx]];
    464:   PetscCall(VecRestoreArrayRead(vec, &arr));
    465:   PetscCall(PetscFree(ix));
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: /*@C
    470:   DMStagVecSetValuesStencil - Set `Vec` values using global grid indexing

    472:   Not Collective

    474:   Input Parameters:
    475: + dm         - the `DMSTAG` object
    476: . vec        - the `Vec`
    477: . n          - the number of values to set
    478: . pos        - the locations to set values, as an array of `DMStagStencil` structs
    479: . val        - the values to set
    480: - insertMode - `INSERT_VALUES` or `ADD_VALUES`

    482:   Notes:
    483:   The vector is expected to be a global vector compatible with the DM (usually obtained by `DMGetGlobalVector()` or `DMCreateGlobalVector()`).

    485:   This approach is not as efficient as setting values directly with `DMStagVecGetArray()`, which is recommended for matrix-free operators.
    486:   For assembling systems, where overhead may be less important than convenience, this routine could be helpful in assembling a righthand side and a matrix (using `DMStagMatSetValuesStencil()`).

    488:   Level: advanced

    490: .seealso: [](ch_stag), `DMSTAG`, `Vec`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMCreateGlobalVector()`, `DMGetLocalVector()`, `DMStagVecGetArray()`
    491: @*/
    492: PetscErrorCode DMStagVecSetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, const PetscScalar *val, InsertMode insertMode)
    493: {
    494:   DM_Stag *const stag = (DM_Stag *)dm->data;
    495:   PetscInt       nLocal;
    496:   PetscInt      *ix;

    498:   PetscFunctionBegin;
    501:   PetscCall(VecGetLocalSize(vec, &nLocal));
    502:   PetscCheck(nLocal == stag->entries, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Provided vec has a different number of local entries (%" PetscInt_FMT ") than expected (%" PetscInt_FMT "). It should be a global vector", nLocal, stag->entries);
    503:   PetscCall(PetscMalloc1(n, &ix));
    504:   PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
    505:   PetscCall(VecSetValuesLocal(vec, n, ix, val, insertMode));
    506:   PetscCall(PetscFree(ix));
    507:   PetscFunctionReturn(PETSC_SUCCESS);
    508: }