Actual source code: stagda.c

      1: /* Routines to convert between a (subset of) DMStag and DMDA */

      3: #include <petscdmda.h>
      4: #include <petsc/private/dmstagimpl.h>
      5: #include <petscdmdatypes.h>

      7: static PetscErrorCode DMStagCreateCompatibleDMDA(DM dm, DMStagStencilLocation loc, PetscInt c, DM *dmda)
      8: {
      9:   DM_Stag *const  stag = (DM_Stag *)dm->data;
     10:   PetscInt        dim, i, j, stencilWidth, dof, N[DMSTAG_MAX_DIM];
     11:   DMDAStencilType stencilType;
     12:   PetscInt       *l[DMSTAG_MAX_DIM];

     14:   PetscFunctionBegin;
     16:   PetscCall(DMGetDimension(dm, &dim));

     18:   /* Create grid decomposition (to be adjusted later) */
     19:   for (i = 0; i < dim; ++i) {
     20:     PetscCall(PetscMalloc1(stag->nRanks[i], &l[i]));
     21:     for (j = 0; j < stag->nRanks[i]; ++j) l[i][j] = stag->l[i][j];
     22:     N[i] = stag->N[i];
     23:   }

     25:   /* dof */
     26:   dof = c < 0 ? -c : 1;

     28:   /* Determine/adjust sizes */
     29:   switch (loc) {
     30:   case DMSTAG_ELEMENT:
     31:     break;
     32:   case DMSTAG_LEFT:
     33:   case DMSTAG_RIGHT:
     34:     PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     35:     l[0][stag->nRanks[0] - 1] += 1; /* extra vertex in direction 0 on last rank in dimension 0 */
     36:     N[0] += 1;
     37:     break;
     38:   case DMSTAG_UP:
     39:   case DMSTAG_DOWN:
     40:     PetscCheck(dim >= 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     41:     l[1][stag->nRanks[1] - 1] += 1; /* extra vertex in direction 1 on last rank in dimension 1 */
     42:     N[1] += 1;
     43:     break;
     44:   case DMSTAG_BACK:
     45:   case DMSTAG_FRONT:
     46:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     47:     l[2][stag->nRanks[2] - 1] += 1; /* extra vertex in direction 2 on last rank in dimension 2 */
     48:     N[2] += 1;
     49:     break;
     50:   case DMSTAG_DOWN_LEFT:
     51:   case DMSTAG_DOWN_RIGHT:
     52:   case DMSTAG_UP_LEFT:
     53:   case DMSTAG_UP_RIGHT:
     54:     PetscCheck(dim >= 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     55:     for (i = 0; i < 2; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1 */
     56:       l[i][stag->nRanks[i] - 1] += 1;
     57:       N[i] += 1;
     58:     }
     59:     break;
     60:   case DMSTAG_BACK_LEFT:
     61:   case DMSTAG_BACK_RIGHT:
     62:   case DMSTAG_FRONT_LEFT:
     63:   case DMSTAG_FRONT_RIGHT:
     64:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     65:     for (i = 0; i < 3; i += 2) { /* extra vertex in direction i on last rank in dimension i = 0,2 */
     66:       l[i][stag->nRanks[i] - 1] += 1;
     67:       N[i] += 1;
     68:     }
     69:     break;
     70:   case DMSTAG_BACK_DOWN:
     71:   case DMSTAG_BACK_UP:
     72:   case DMSTAG_FRONT_DOWN:
     73:   case DMSTAG_FRONT_UP:
     74:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     75:     for (i = 1; i < 3; ++i) { /* extra vertex in direction i on last rank in dimension i = 1,2 */
     76:       l[i][stag->nRanks[i] - 1] += 1;
     77:       N[i] += 1;
     78:     }
     79:     break;
     80:   case DMSTAG_BACK_DOWN_LEFT:
     81:   case DMSTAG_BACK_DOWN_RIGHT:
     82:   case DMSTAG_BACK_UP_LEFT:
     83:   case DMSTAG_BACK_UP_RIGHT:
     84:   case DMSTAG_FRONT_DOWN_LEFT:
     85:   case DMSTAG_FRONT_DOWN_RIGHT:
     86:   case DMSTAG_FRONT_UP_LEFT:
     87:   case DMSTAG_FRONT_UP_RIGHT:
     88:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
     89:     for (i = 0; i < 3; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1,2 */
     90:       l[i][stag->nRanks[i] - 1] += 1;
     91:       N[i] += 1;
     92:     }
     93:     break;
     94:   default:
     95:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
     96:   }

     98:   /* Use the same stencil type */
     99:   switch (stag->stencilType) {
    100:   case DMSTAG_STENCIL_STAR:
    101:     stencilType  = DMDA_STENCIL_STAR;
    102:     stencilWidth = stag->stencilWidth;
    103:     break;
    104:   case DMSTAG_STENCIL_BOX:
    105:     stencilType  = DMDA_STENCIL_BOX;
    106:     stencilWidth = stag->stencilWidth;
    107:     break;
    108:   default:
    109:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported Stencil Type %d", stag->stencilType);
    110:   }

    112:   /* Create DMDA, using same boundary type */
    113:   switch (dim) {
    114:   case 1:
    115:     PetscCall(DMDACreate1d(PetscObjectComm((PetscObject)dm), stag->boundaryType[0], N[0], dof, stencilWidth, l[0], dmda));
    116:     break;
    117:   case 2:
    118:     PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)dm), stag->boundaryType[0], stag->boundaryType[1], stencilType, N[0], N[1], stag->nRanks[0], stag->nRanks[1], dof, stencilWidth, l[0], l[1], dmda));
    119:     break;
    120:   case 3:
    121:     PetscCall(DMDACreate3d(PetscObjectComm((PetscObject)dm), stag->boundaryType[0], stag->boundaryType[1], stag->boundaryType[2], stencilType, N[0], N[1], N[2], stag->nRanks[0], stag->nRanks[1], stag->nRanks[2], dof, stencilWidth, l[0], l[1], l[2], dmda));
    122:     break;
    123:   default:
    124:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim %" PetscInt_FMT, dim);
    125:   }
    126:   for (i = 0; i < dim; ++i) PetscCall(PetscFree(l[i]));
    127:   PetscFunctionReturn(PETSC_SUCCESS);
    128: }

    130: /*
    131: Helper function to get the number of extra points in a DMDA representation for a given canonical location.
    132: */
    133: static PetscErrorCode DMStagDMDAGetExtraPoints(DM dm, DMStagStencilLocation locCanonical, PetscInt *extraPoint)
    134: {
    135:   PetscInt dim, d, nExtra[DMSTAG_MAX_DIM];

    137:   PetscFunctionBegin;
    139:   PetscCall(DMGetDimension(dm, &dim));
    140:   PetscCheck(dim <= DMSTAG_MAX_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimensions", dim);
    141:   PetscCall(DMStagGetCorners(dm, NULL, NULL, NULL, NULL, NULL, NULL, &nExtra[0], &nExtra[1], &nExtra[2]));
    142:   for (d = 0; d < dim; ++d) extraPoint[d] = 0;
    143:   switch (locCanonical) {
    144:   case DMSTAG_ELEMENT:
    145:     break; /* no extra points */
    146:   case DMSTAG_LEFT:
    147:     extraPoint[0] = nExtra[0];
    148:     break; /* only extra point in x */
    149:   case DMSTAG_DOWN:
    150:     extraPoint[1] = nExtra[1];
    151:     break; /* only extra point in y */
    152:   case DMSTAG_BACK:
    153:     extraPoint[2] = nExtra[2];
    154:     break; /* only extra point in z */
    155:   case DMSTAG_DOWN_LEFT:
    156:     extraPoint[0] = nExtra[0];
    157:     extraPoint[1] = nExtra[1];
    158:     break; /* extra point in both x and y  */
    159:   case DMSTAG_BACK_LEFT:
    160:     extraPoint[0] = nExtra[0];
    161:     extraPoint[2] = nExtra[2];
    162:     break; /* extra point in both x and z  */
    163:   case DMSTAG_BACK_DOWN:
    164:     extraPoint[1] = nExtra[1];
    165:     extraPoint[2] = nExtra[2];
    166:     break; /* extra point in both y and z  */
    167:   case DMSTAG_BACK_DOWN_LEFT:
    168:     extraPoint[0] = nExtra[0];
    169:     extraPoint[1] = nExtra[1];
    170:     extraPoint[2] = nExtra[2];
    171:     break; /* extra points in x,y,z */
    172:   default:
    173:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location (perhaps not canonical) %s", DMStagStencilLocations[locCanonical]);
    174:   }
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: /*
    179: Function much like DMStagMigrateVec(), but which accepts an additional position argument to disambiguate which
    180: type of DMDA to migrate to.
    181: */

    183: static PetscErrorCode DMStagMigrateVecDMDA(DM dm, Vec vec, DMStagStencilLocation loc, PetscInt c, DM dmTo, Vec vecTo)
    184: {
    185:   PetscInt i, j, k, d, dim, dof, dofToMax, start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM];
    186:   Vec      vecLocal;

    188:   PetscFunctionBegin;
    193:   PetscCall(DMGetDimension(dm, &dim));
    194:   PetscCall(DMDAGetDof(dmTo, &dofToMax));
    195:   PetscCheck(-c <= dofToMax, PetscObjectComm((PetscObject)dmTo), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative component value. Must be >= -%" PetscInt_FMT, dofToMax);
    196:   PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], NULL, NULL, NULL));
    197:   PetscCall(DMStagDMDAGetExtraPoints(dm, loc, extraPoint));
    198:   PetscCall(DMStagGetLocationDOF(dm, loc, &dof));
    199:   PetscCall(DMGetLocalVector(dm, &vecLocal));
    200:   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
    201:   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
    202:   if (dim == 1) {
    203:     PetscScalar **arrTo;
    204:     PetscCall(DMDAVecGetArrayDOF(dmTo, vecTo, &arrTo));
    205:     if (c < 0) {
    206:       const PetscInt dofTo = -c;
    207:       for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    208:         for (d = 0; d < PetscMin(dof, dofTo); ++d) {
    209:           DMStagStencil pos;
    210:           pos.i   = i;
    211:           pos.loc = loc;
    212:           pos.c   = d;
    213:           PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[i][d]));
    214:         }
    215:         for (; d < dofTo; ++d) arrTo[i][d] = 0.0; /* Pad extra dof with zeros */
    216:       }
    217:     } else {
    218:       for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    219:         DMStagStencil pos;
    220:         pos.i   = i;
    221:         pos.loc = loc;
    222:         pos.c   = c;
    223:         PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[i][0]));
    224:       }
    225:     }
    226:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
    227:   } else if (dim == 2) {
    228:     PetscScalar ***arrTo;
    229:     PetscCall(DMDAVecGetArrayDOF(dmTo, vecTo, &arrTo));
    230:     if (c < 0) {
    231:       const PetscInt dofTo = -c;
    232:       for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
    233:         for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    234:           for (d = 0; d < PetscMin(dof, dofTo); ++d) {
    235:             DMStagStencil pos;
    236:             pos.i   = i;
    237:             pos.j   = j;
    238:             pos.loc = loc;
    239:             pos.c   = d;
    240:             PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[j][i][d]));
    241:           }
    242:           for (; d < dofTo; ++d) arrTo[j][i][d] = 0.0; /* Pad extra dof with zeros */
    243:         }
    244:       }
    245:     } else {
    246:       for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
    247:         for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    248:           DMStagStencil pos;
    249:           pos.i   = i;
    250:           pos.j   = j;
    251:           pos.loc = loc;
    252:           pos.c   = c;
    253:           PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[j][i][0]));
    254:         }
    255:       }
    256:     }
    257:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
    258:   } else if (dim == 3) {
    259:     PetscScalar ****arrTo;
    260:     PetscCall(DMDAVecGetArrayDOF(dmTo, vecTo, &arrTo));
    261:     if (c < 0) {
    262:       const PetscInt dofTo = -c;
    263:       for (k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
    264:         for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
    265:           for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    266:             for (d = 0; d < PetscMin(dof, dofTo); ++d) {
    267:               DMStagStencil pos;
    268:               pos.i   = i;
    269:               pos.j   = j;
    270:               pos.k   = k;
    271:               pos.loc = loc;
    272:               pos.c   = d;
    273:               PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[k][j][i][d]));
    274:             }
    275:             for (; d < dofTo; ++d) arrTo[k][j][i][d] = 0.0; /* Pad extra dof with zeros */
    276:           }
    277:         }
    278:       }
    279:     } else {
    280:       for (k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
    281:         for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
    282:           for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
    283:             DMStagStencil pos;
    284:             pos.i   = i;
    285:             pos.j   = j;
    286:             pos.k   = k;
    287:             pos.loc = loc;
    288:             pos.c   = c;
    289:             PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[k][j][i][0]));
    290:           }
    291:         }
    292:       }
    293:     }
    294:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
    295:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
    296:   PetscCall(DMRestoreLocalVector(dm, &vecLocal));
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: /* Transfer coordinates from a DMStag to a DMDA, specifying which location */
    301: static PetscErrorCode DMStagTransferCoordinatesToDMDA(DM dmstag, DMStagStencilLocation loc, DM dmda)
    302: {
    303:   PetscInt  dim, start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM], d;
    304:   DM        dmstagCoord, dmdaCoord;
    305:   DMType    dmstagCoordType;
    306:   Vec       stagCoord, daCoord;
    307:   PetscBool daCoordIsStag, daCoordIsProduct;

    309:   PetscFunctionBegin;
    312:   PetscCall(DMGetDimension(dmstag, &dim));
    313:   PetscCall(DMGetCoordinateDM(dmstag, &dmstagCoord));
    314:   PetscCall(DMGetCoordinatesLocal(dmstag, &stagCoord)); /* Note local */
    315:   PetscCall(DMGetCoordinateDM(dmda, &dmdaCoord));
    316:   daCoord = NULL;
    317:   PetscCall(DMGetCoordinates(dmda, &daCoord));
    318:   if (!daCoord) {
    319:     PetscCall(DMCreateGlobalVector(dmdaCoord, &daCoord));
    320:     PetscCall(DMSetCoordinates(dmda, daCoord));
    321:     PetscCall(VecDestroy(&daCoord));
    322:     PetscCall(DMGetCoordinates(dmda, &daCoord));
    323:   }
    324:   PetscCall(DMGetType(dmstagCoord, &dmstagCoordType));
    325:   PetscCall(PetscStrcmp(dmstagCoordType, DMSTAG, &daCoordIsStag));
    326:   PetscCall(PetscStrcmp(dmstagCoordType, DMPRODUCT, &daCoordIsProduct));
    327:   PetscCall(DMStagGetCorners(dmstag, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], NULL, NULL, NULL));
    328:   PetscCall(DMStagDMDAGetExtraPoints(dmstag, loc, extraPoint));
    329:   if (dim == 1) {
    330:     PetscInt      ex;
    331:     PetscScalar **cArrDa;
    332:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
    333:     if (daCoordIsStag) {
    334:       PetscInt      slot;
    335:       PetscScalar **cArrStag;
    336:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
    337:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
    338:       for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) cArrDa[ex][0] = cArrStag[ex][slot];
    339:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
    340:     } else if (daCoordIsProduct) {
    341:       PetscScalar **cArrX;
    342:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, NULL, NULL));
    343:       for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) cArrDa[ex][0] = cArrX[ex][0];
    344:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, NULL, NULL));
    345:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
    346:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
    347:   } else if (dim == 2) {
    348:     PetscInt       ex, ey;
    349:     PetscScalar ***cArrDa;
    350:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
    351:     if (daCoordIsStag) {
    352:       PetscInt       slot;
    353:       PetscScalar ***cArrStag;
    354:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
    355:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
    356:       for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
    357:         for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
    358:           for (d = 0; d < 2; ++d) cArrDa[ey][ex][d] = cArrStag[ey][ex][slot + d];
    359:         }
    360:       }
    361:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
    362:     } else if (daCoordIsProduct) {
    363:       PetscScalar **cArrX, **cArrY;
    364:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, NULL));
    365:       for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
    366:         for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
    367:           cArrDa[ey][ex][0] = cArrX[ex][0];
    368:           cArrDa[ey][ex][1] = cArrY[ey][0];
    369:         }
    370:       }
    371:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, NULL));
    372:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
    373:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
    374:   } else if (dim == 3) {
    375:     PetscInt        ex, ey, ez;
    376:     PetscScalar ****cArrDa;
    377:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
    378:     if (daCoordIsStag) {
    379:       PetscInt        slot;
    380:       PetscScalar ****cArrStag;
    381:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
    382:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
    383:       for (ez = start[2]; ez < start[2] + n[2] + extraPoint[2]; ++ez) {
    384:         for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
    385:           for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
    386:             for (d = 0; d < 3; ++d) cArrDa[ez][ey][ex][d] = cArrStag[ez][ey][ex][slot + d];
    387:           }
    388:         }
    389:       }
    390:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
    391:     } else if (daCoordIsProduct) {
    392:       PetscScalar **cArrX, **cArrY, **cArrZ;
    393:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, &cArrZ));
    394:       for (ez = start[2]; ez < start[2] + n[2] + extraPoint[2]; ++ez) {
    395:         for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
    396:           for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
    397:             cArrDa[ez][ey][ex][0] = cArrX[ex][0];
    398:             cArrDa[ez][ey][ex][1] = cArrY[ey][0];
    399:             cArrDa[ez][ey][ex][2] = cArrZ[ez][0];
    400:           }
    401:         }
    402:       }
    403:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, &cArrZ));
    404:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
    405:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
    406:   } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
    407:   PetscFunctionReturn(PETSC_SUCCESS);
    408: }

    410: /*@
    411:   DMStagVecSplitToDMDA - create a `DMDA` and `Vec` from a subgrid of a `DMSTAG` and its `Vec`

    413:   Collective

    415:   Input Parameters:
    416: + dm  - the `DMSTAG` object
    417: . vec - `Vec` object associated with `dm`
    418: . loc - which subgrid to extract (see `DMStagStencilLocation`)
    419: - c   - which component to extract (see note below)

    421:   Output Parameters:
    422: + pda    - the `DMDA`
    423: - pdavec - the new `Vec`

    425:   Level: advanced

    427:   Notes:
    428:   If a `c` value of `-k` is provided, the first `k` DOF for that position are extracted,
    429:   padding with zero values if needed. If a non-negative value is provided, a single
    430:   DOF is extracted.

    432:   The caller is responsible for destroying the created `DMDA` and `Vec`.

    434: .seealso: [](ch_stag), `DMSTAG`, `DMDA`, `DMStagStencilLocation`, `DM`, `Vec`, `DMStagMigrateVec()`, `DMStagCreateCompatibleDMStag()`
    435: @*/
    436: PetscErrorCode DMStagVecSplitToDMDA(DM dm, Vec vec, DMStagStencilLocation loc, PetscInt c, DM *pda, Vec *pdavec)
    437: {
    438:   PetscInt              dim, locdof;
    439:   DM                    da, coordDM;
    440:   Vec                   davec;
    441:   DMStagStencilLocation locCanonical;

    443:   PetscFunctionBegin;
    446:   PetscCall(DMGetDimension(dm, &dim));
    447:   PetscCall(DMStagGetLocationDOF(dm, loc, &locdof));
    448:   PetscCheck(c < locdof, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Location %s has %" PetscInt_FMT " dof, but component %" PetscInt_FMT " requested", DMStagStencilLocations[loc], locdof, c);
    449:   PetscCall(DMStagStencilLocationCanonicalize(loc, &locCanonical));
    450:   PetscCall(DMStagCreateCompatibleDMDA(dm, locCanonical, c, pda));
    451:   da = *pda;
    452:   PetscCall(DMSetUp(*pda));
    453:   if (dm->coordinates[0].dm != NULL) {
    454:     PetscCall(DMGetCoordinateDM(dm, &coordDM));
    455:     PetscCall(DMStagTransferCoordinatesToDMDA(dm, locCanonical, da));
    456:   }
    457:   PetscCall(DMCreateGlobalVector(da, pdavec));
    458:   davec = *pdavec;
    459:   PetscCall(DMStagMigrateVecDMDA(dm, vec, locCanonical, c, da, davec));
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }