Actual source code: plexproject.c

      1: #include <petsc/private/dmpleximpl.h>

      3: #include <petsc/private/petscfeimpl.h>

      5: /*@
      6:   DMPlexGetActivePoint - Get the point on which projection is currently working

      8:   Not Collective

     10:   Input Parameter:
     11: . dm - the `DM`

     13:   Output Parameter:
     14: . point - The mesh point involved in the current projection

     16:   Level: developer

     18: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`
     19: @*/
     20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
     21: {
     22:   PetscFunctionBeginHot;
     23:   *point = ((DM_Plex *)dm->data)->activePoint;
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: /*@
     28:   DMPlexSetActivePoint - Set the point on which projection is currently working

     30:   Not Collective

     32:   Input Parameters:
     33: + dm    - the `DM`
     34: - point - The mesh point involved in the current projection

     36:   Level: developer

     38: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetActivePoint()`
     39: @*/
     40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
     41: {
     42:   PetscFunctionBeginHot;
     43:   ((DM_Plex *)dm->data)->activePoint = point;
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: /*
     48:   DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point

     50:   Input Parameters:
     51: + dm     - The output `DM`
     52: . ds     - The output `DS`
     53: . dmIn   - The input `DM`
     54: . dsIn   - The input `DS`
     55: . time   - The time for this evaluation
     56: . fegeom - The FE geometry for this point
     57: . fvgeom - The FV geometry for this point
     58: . isFE   - Flag indicating whether each output field has an FE discretization
     59: . sp     - The output `PetscDualSpace` for each field
     60: . funcs  - The evaluation function for each field
     61: - ctxs   - The user context for each field

     63:   Output Parameter:
     64: . values - The value for each dual basis vector in the output dual space

     66:   Level: developer

     68: .seealso:[](ch_unstructured), `DM`, `DMPLEX`, `PetscDS`, `PetscFEGeom`, `PetscFVCellGeom`, `PetscDualSpace`
     69: */
     70: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, PetscScalar values[])
     71: {
     72:   PetscInt  debug = ((DM_Plex *)dm->data)->printProject;
     73:   PetscInt  coordDim, Nf, *Nc, f, spDim, d, v, tp;
     74:   PetscBool isAffine, isCohesive, transform;

     76:   PetscFunctionBeginHot;
     77:   PetscCall(DMGetCoordinateDim(dmIn, &coordDim));
     78:   PetscCall(DMHasBasisTransform(dmIn, &transform));
     79:   PetscCall(PetscDSGetNumFields(ds, &Nf));
     80:   PetscCall(PetscDSGetComponents(ds, &Nc));
     81:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
     82:   /* Get values for closure */
     83:   isAffine = fegeom->isAffine;
     84:   for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
     85:     void *const ctx = ctxs ? ctxs[f] : NULL;
     86:     PetscBool   cohesive;

     88:     if (!sp[f]) continue;
     89:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
     90:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
     91:     if (funcs[f]) {
     92:       if (isFE[f]) {
     93:         PetscQuadrature  allPoints;
     94:         PetscInt         q, dim, numPoints;
     95:         const PetscReal *points;
     96:         PetscScalar     *pointEval;
     97:         PetscReal       *x;
     98:         DM               rdm;

    100:         PetscCall(PetscDualSpaceGetDM(sp[f], &rdm));
    101:         PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
    102:         PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
    103:         PetscCall(DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    104:         PetscCall(DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x));
    105:         PetscCall(PetscArrayzero(pointEval, numPoints * Nc[f]));
    106:         for (q = 0; q < numPoints; q++, tp++) {
    107:           const PetscReal *v0;

    109:           if (isAffine) {
    110:             const PetscReal *refpoint    = &points[q * dim];
    111:             PetscReal        injpoint[3] = {0., 0., 0.};

    113:             if (dim != fegeom->dim) {
    114:               PetscCheck(isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %" PetscInt_FMT " != %" PetscInt_FMT " dual basis spatial dimension", fegeom->dim, dim);
    115:               /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
    116:               for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
    117:               refpoint = injpoint;
    118:             }
    119:             CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
    120:             v0 = x;
    121:           } else {
    122:             v0 = &fegeom->v[tp * coordDim];
    123:           }
    124:           if (transform) {
    125:             PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx));
    126:             v0 = x;
    127:           }
    128:           if (debug > 3) {
    129:             PetscInt ap;
    130:             PetscCall(DMPlexGetActivePoint(dm, &ap));
    131:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "Project point %" PetscInt_FMT ", analytic: ref (", ap));
    132:             for (PetscInt d = 0; d < dim; ++d) {
    133:               if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    134:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)points[q * dim + d]));
    135:             }
    136:             PetscCall(PetscPrintf(PETSC_COMM_SELF, ") real ("));
    137:             for (PetscInt d = 0; d < dim; ++d) {
    138:               if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    139:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)v0[d]));
    140:             }
    141:             PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
    142:           }
    143:           PetscCall((*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx));
    144:         }
    145:         /* Transform point evaluations pointEval[q,c] */
    146:         PetscCall(PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval));
    147:         PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    148:         PetscCall(DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x));
    149:         PetscCall(DMRestoreWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    150:         v += spDim;
    151:         if (isCohesive && !cohesive) {
    152:           for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    153:         }
    154:       } else {
    155:         for (d = 0; d < spDim; ++d, ++v) PetscCall(PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]));
    156:       }
    157:     } else {
    158:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    159:       if (isCohesive && !cohesive) {
    160:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    161:       }
    162:     }
    163:   }
    164:   PetscFunctionReturn(PETSC_SUCCESS);
    165: }

    167: /*
    168:   DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point

    170:   Input Parameters:
    171: + dm             - The output DM
    172: . ds             - The output DS
    173: . dmIn           - The input DM
    174: . dsIn           - The input DS
    175: . dmAux          - The auxiliary DM, which is always for the input space
    176: . dsAux          - The auxiliary DS, which is always for the input space
    177: . time           - The time for this evaluation
    178: . localU         - The local solution
    179: . localA         - The local auziliary fields
    180: . cgeom          - The FE geometry for this point
    181: . sp             - The output PetscDualSpace for each field
    182: . p              - The point in the output DM
    183: . T              - Input basis and derivatives for each field tabulated on the quadrature points
    184: . TAux           - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
    185: . funcs          - The evaluation function for each field
    186: - ctxs           - The user context for each field

    188:   Output Parameter:
    189: . values         - The value for each dual basis vector in the output dual space

    191:   Level: developer

    193:   Note:
    194:   Not supported for FV

    196: .seealso: `DMProjectPoint_Field_Private()`
    197: */
    198: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
    199: {
    200:   PetscSection       section, sectionAux = NULL;
    201:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
    202:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
    203:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
    204:   const PetscScalar *constants;
    205:   PetscReal         *x;
    206:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
    207:   PetscFEGeom        fegeom, fgeomN[2];
    208:   const PetscInt     dE = cgeom->dimEmbed, *cone, *ornt;
    209:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
    210:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
    211:   DMPolytopeType     qct;

    213:   PetscFunctionBeginHot;
    214:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    215:   PetscCall(PetscDSGetComponents(ds, &Nc));
    216:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    217:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    218:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    219:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
    220:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
    221:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
    222:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
    223:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
    224:   PetscCall(DMHasBasisTransform(dmIn, &transform));
    225:   PetscCall(DMGetLocalSection(dmIn, &section));
    226:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
    227:   // Get cohesive cell hanging off face
    228:   if (isCohesiveIn) {
    229:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
    230:     if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
    231:       DMPolytopeType  ct;
    232:       const PetscInt *support;
    233:       PetscInt        Ns, s;

    235:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
    236:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
    237:       for (s = 0; s < Ns; ++s) {
    238:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
    239:         if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
    240:           inp = support[s];
    241:           break;
    242:         }
    243:       }
    244:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
    245:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
    246:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
    247:       face[0] = 0;
    248:       face[1] = 0;
    249:     }
    250:   }
    251:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
    252:   if (dmAux) {
    253:     PetscInt subp;

    255:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
    256:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    257:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
    258:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    259:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    260:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
    261:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
    262:   }
    263:   /* Get values for closure */
    264:   isAffine        = cgeom->isAffine;
    265:   fegeom.dim      = cgeom->dim;
    266:   fegeom.dimEmbed = cgeom->dimEmbed;
    267:   if (isCohesiveIn) {
    268:     fgeomN[0].dim      = cgeom->dim;
    269:     fgeomN[0].dimEmbed = cgeom->dimEmbed;
    270:     fgeomN[1].dim      = cgeom->dim;
    271:     fgeomN[1].dimEmbed = cgeom->dimEmbed;
    272:   }
    273:   if (isAffine) {
    274:     fegeom.v    = x;
    275:     fegeom.xi   = cgeom->xi;
    276:     fegeom.J    = cgeom->J;
    277:     fegeom.invJ = cgeom->invJ;
    278:     fegeom.detJ = cgeom->detJ;
    279:     if (isCohesiveIn) {
    280:       fgeomN[0].J    = cgeom->suppJ[0];
    281:       fgeomN[0].invJ = cgeom->suppInvJ[0];
    282:       fgeomN[0].detJ = cgeom->suppDetJ[0];
    283:       fgeomN[1].J    = cgeom->suppJ[1];
    284:       fgeomN[1].invJ = cgeom->suppInvJ[1];
    285:       fgeomN[1].detJ = cgeom->suppDetJ[1];
    286:     }
    287:   }
    288:   for (f = 0, v = 0; f < Nf; ++f) {
    289:     PetscQuadrature  allPoints;
    290:     PetscInt         q, dim, numPoints;
    291:     const PetscReal *points;
    292:     PetscScalar     *pointEval;
    293:     PetscBool        cohesive;
    294:     DM               dm;

    296:     if (!sp[f]) continue;
    297:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    298:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    299:     if (!funcs[f]) {
    300:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    301:       if (isCohesive && !cohesive) {
    302:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    303:       }
    304:       continue;
    305:     }
    306:     const PetscInt ***perms;
    307:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
    308:     PetscCall(PetscDualSpaceGetSymmetries(sp[f], &perms, NULL));
    309:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
    310:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
    311:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    312:     for (q = 0; q < numPoints; ++q, ++tp) {
    313:       PetscInt qpt[2];

    315:       if (isCohesiveIn) {
    316:         qpt[0] = perms ? perms[0][ornt[0]][q] : q;
    317:         qpt[1] = perms ? perms[0][DMPolytopeTypeComposeOrientationInv(qct, ornt[1], 0)][q] : q;
    318:       }
    319:       if (isAffine) {
    320:         CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
    321:       } else {
    322:         fegeom.v    = &cgeom->v[tp * dE];
    323:         fegeom.J    = &cgeom->J[tp * dE * dE];
    324:         fegeom.invJ = &cgeom->invJ[tp * dE * dE];
    325:         fegeom.detJ = &cgeom->detJ[tp];
    326:         if (isCohesiveIn) {
    327:           fgeomN[0].J    = &cgeom->suppJ[0][tp * dE * dE];
    328:           fgeomN[0].invJ = &cgeom->suppInvJ[0][tp * dE * dE];
    329:           fgeomN[0].detJ = &cgeom->suppDetJ[0][tp];
    330:           fgeomN[1].J    = &cgeom->suppJ[1][tp * dE * dE];
    331:           fgeomN[1].invJ = &cgeom->suppInvJ[1][tp * dE * dE];
    332:           fgeomN[1].detJ = &cgeom->suppDetJ[1][tp];
    333:         }
    334:       }
    335:       if (coefficients) {
    336:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, fgeomN, coefficients, coefficients_t, u, u_x, u_t));
    337:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
    338:       }
    339:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
    340:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
    341:       (*funcs[f])(dE, NfIn, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, fegeom.v, numConstants, constants, &pointEval[Nc[f] * q]);
    342:     }
    343:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    344:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    345:     v += spDim;
    346:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
    347:     if (isCohesive && !cohesive) {
    348:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    349:     }
    350:   }
    351:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
    352:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
    353:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
    354:   PetscFunctionReturn(PETSC_SUCCESS);
    355: }

    357: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
    358: {
    359:   PetscSection       section, sectionAux = NULL;
    360:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
    361:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
    362:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
    363:   const PetscScalar *constants;
    364:   PetscReal         *x;
    365:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
    366:   PetscFEGeom        fegeom, cgeom, fgeomN[2];
    367:   const PetscInt     dE = fgeom->dimEmbed, *cone, *ornt;
    368:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
    369:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
    370:   DMPolytopeType     qct;

    372:   PetscFunctionBeginHot;
    373:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    374:   PetscCall(PetscDSGetComponents(ds, &Nc));
    375:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    376:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    377:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    378:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
    379:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
    380:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
    381:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
    382:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
    383:   PetscCall(DMHasBasisTransform(dmIn, &transform));
    384:   PetscCall(DMGetLocalSection(dmIn, &section));
    385:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
    386:   // Get cohesive cell hanging off face
    387:   if (isCohesiveIn) {
    388:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
    389:     if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
    390:       DMPolytopeType  ct;
    391:       const PetscInt *support;
    392:       PetscInt        Ns, s;

    394:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
    395:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
    396:       for (s = 0; s < Ns; ++s) {
    397:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
    398:         if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
    399:           inp = support[s];
    400:           break;
    401:         }
    402:       }
    403:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
    404:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
    405:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
    406:       face[0] = 0;
    407:       face[1] = 0;
    408:     }
    409:   }
    410:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
    411:   if (dmAux) {
    412:     PetscInt subp;

    414:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
    415:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    416:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
    417:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    418:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    419:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
    420:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
    421:   }
    422:   /* Get values for closure */
    423:   isAffine        = fgeom->isAffine;
    424:   fegeom.dim      = fgeom->dim;
    425:   fegeom.dimEmbed = fgeom->dimEmbed;
    426:   fegeom.n        = NULL;
    427:   fegeom.J        = NULL;
    428:   fegeom.v        = NULL;
    429:   fegeom.xi       = NULL;
    430:   cgeom.dim       = fgeom->dim;
    431:   cgeom.dimEmbed  = fgeom->dimEmbed;
    432:   if (isCohesiveIn) {
    433:     fgeomN[0].dim      = fgeom->dim;
    434:     fgeomN[0].dimEmbed = fgeom->dimEmbed;
    435:     fgeomN[1].dim      = fgeom->dim;
    436:     fgeomN[1].dimEmbed = fgeom->dimEmbed;
    437:   }
    438:   if (isAffine) {
    439:     fegeom.v    = x;
    440:     fegeom.xi   = fgeom->xi;
    441:     fegeom.J    = fgeom->J;
    442:     fegeom.invJ = fgeom->invJ;
    443:     fegeom.detJ = fgeom->detJ;
    444:     fegeom.n    = fgeom->n;

    446:     cgeom.J    = fgeom->suppJ[0];
    447:     cgeom.invJ = fgeom->suppInvJ[0];
    448:     cgeom.detJ = fgeom->suppDetJ[0];

    450:     if (isCohesiveIn) {
    451:       fgeomN[0].J    = fgeom->suppJ[0];
    452:       fgeomN[0].invJ = fgeom->suppInvJ[0];
    453:       fgeomN[0].detJ = fgeom->suppDetJ[0];
    454:       fgeomN[1].J    = fgeom->suppJ[1];
    455:       fgeomN[1].invJ = fgeom->suppInvJ[1];
    456:       fgeomN[1].detJ = fgeom->suppDetJ[1];
    457:     }
    458:   }
    459:   for (f = 0, v = 0; f < Nf; ++f) {
    460:     PetscQuadrature  allPoints;
    461:     PetscInt         q, dim, numPoints;
    462:     const PetscReal *points;
    463:     PetscScalar     *pointEval;
    464:     PetscBool        cohesive;
    465:     DM               dm;

    467:     if (!sp[f]) continue;
    468:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    469:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    470:     if (!funcs[f]) {
    471:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    472:       if (isCohesive && !cohesive) {
    473:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    474:       }
    475:       continue;
    476:     }
    477:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
    478:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
    479:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
    480:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    481:     for (q = 0; q < numPoints; ++q, ++tp) {
    482:       PetscInt qpt[2];

    484:       if (isCohesiveIn) {
    485:         // These points are not integration quadratures, but dual space quadratures
    486:         // If they had multiple points we should match them from both sides, similar to hybrid residual eval
    487:         qpt[0] = qpt[1] = q;
    488:       }
    489:       if (isAffine) {
    490:         CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
    491:       } else {
    492:         fegeom.v    = &fgeom->v[tp * dE];
    493:         fegeom.J    = &fgeom->J[tp * dE * dE];
    494:         fegeom.invJ = &fgeom->invJ[tp * dE * dE];
    495:         fegeom.detJ = &fgeom->detJ[tp];
    496:         fegeom.n    = &fgeom->n[tp * dE];

    498:         cgeom.J    = &fgeom->suppJ[0][tp * dE * dE];
    499:         cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
    500:         cgeom.detJ = &fgeom->suppDetJ[0][tp];
    501:         if (isCohesiveIn) {
    502:           fgeomN[0].J    = &fgeom->suppJ[0][tp * dE * dE];
    503:           fgeomN[0].invJ = &fgeom->suppInvJ[0][tp * dE * dE];
    504:           fgeomN[0].detJ = &fgeom->suppDetJ[0][tp];
    505:           fgeomN[1].J    = &fgeom->suppJ[1][tp * dE * dE];
    506:           fgeomN[1].invJ = &fgeom->suppInvJ[1][tp * dE * dE];
    507:           fgeomN[1].detJ = &fgeom->suppDetJ[1][tp];
    508:         }
    509:       }
    510:       /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
    511:       if (coefficients) {
    512:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, fgeomN, coefficients, coefficients_t, u, u_x, u_t));
    513:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
    514:       }
    515:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
    516:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
    517:       (*funcs[f])(dE, NfIn, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, fegeom.v, fegeom.n, numConstants, constants, &pointEval[Nc[f] * q]);
    518:     }
    519:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    520:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    521:     v += spDim;
    522:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
    523:     if (isCohesive && !cohesive) {
    524:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    525:     }
    526:   }
    527:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
    528:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
    529:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
    530:   PetscFunctionReturn(PETSC_SUCCESS);
    531: }

    533: static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscFEGeom *fegeom, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, DMBoundaryConditionType type, PetscVoidFn **funcs, void **ctxs, PetscBool fieldActive[], PetscScalar values[])
    534: {
    535:   PetscFVCellGeom fvgeom;
    536:   PetscInt        dim, dimEmbed;

    538:   PetscFunctionBeginHot;
    539:   PetscCall(DMGetDimension(dm, &dim));
    540:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
    541:   if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
    542:   switch (type) {
    543:   case DM_BC_ESSENTIAL:
    544:   case DM_BC_NATURAL:
    545:     PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
    546:     break;
    547:   case DM_BC_ESSENTIAL_FIELD:
    548:   case DM_BC_NATURAL_FIELD:
    549:     PetscCall(DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
    550:     break;
    551:   case DM_BC_ESSENTIAL_BD_FIELD:
    552:     PetscCall(DMProjectPoint_BdField_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
    553:     break;
    554:   default:
    555:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
    556:   }
    557:   PetscFunctionReturn(PETSC_SUCCESS);
    558: }

    560: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, PetscVoidFn **funcs, PetscQuadrature *allPoints)
    561: {
    562:   PetscReal *points;
    563:   PetscInt   f, numPoints;

    565:   PetscFunctionBegin;
    566:   if (!dim) {
    567:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
    568:     PetscFunctionReturn(PETSC_SUCCESS);
    569:   }
    570:   numPoints = 0;
    571:   for (f = 0; f < Nf; ++f) {
    572:     if (funcs[f]) {
    573:       PetscQuadrature fAllPoints;
    574:       PetscInt        fNumPoints;

    576:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
    577:       PetscCall(PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL));
    578:       numPoints += fNumPoints;
    579:     }
    580:   }
    581:   PetscCall(PetscMalloc1(dim * numPoints, &points));
    582:   numPoints = 0;
    583:   for (f = 0; f < Nf; ++f) {
    584:     if (funcs[f]) {
    585:       PetscQuadrature  fAllPoints;
    586:       PetscInt         qdim, fNumPoints, q;
    587:       const PetscReal *fPoints;

    589:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
    590:       PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
    591:       PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %" PetscInt_FMT " for dual basis does not match input dimension %" PetscInt_FMT, qdim, dim);
    592:       for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
    593:       numPoints += fNumPoints;
    594:     }
    595:   }
    596:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
    597:   PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
    598:   PetscFunctionReturn(PETSC_SUCCESS);
    599: }

    601: /*@C
    602:   DMGetFirstLabeledPoint - Find first labeled `point` in `odm` such that the corresponding point in `dm` has the specified `height`. Return `point` and the corresponding `ds`.

    604:   Input Parameters:
    605: + dm     - the `DM`
    606: . odm    - the enclosing `DM`
    607: . label  - label for `DM` domain, or `NULL` for whole domain
    608: . numIds - the number of `ids`
    609: . ids    - An array of the label ids in sequence for the domain
    610: - height - Height of target cells in `DMPLEX` topology

    612:   Output Parameters:
    613: + point - the first labeled point
    614: - ds    - the `PetscDS` corresponding to the first labeled point

    616:   Level: developer

    618: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
    619: @*/
    620: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
    621: {
    622:   DM              plex;
    623:   DMEnclosureType enc;
    624:   PetscInt        ls = -1;

    626:   PetscFunctionBegin;
    627:   if (point) *point = -1;
    628:   if (!label) PetscFunctionReturn(PETSC_SUCCESS);
    629:   PetscCall(DMGetEnclosureRelation(dm, odm, &enc));
    630:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    631:   for (PetscInt i = 0; i < numIds; ++i) {
    632:     IS       labelIS;
    633:     PetscInt num_points, pStart, pEnd;
    634:     PetscCall(DMLabelGetStratumIS(label, ids[i], &labelIS));
    635:     if (!labelIS) continue; /* No points with that id on this process */
    636:     PetscCall(DMPlexGetHeightStratum(plex, height, &pStart, &pEnd));
    637:     PetscCall(ISGetSize(labelIS, &num_points));
    638:     if (num_points) {
    639:       const PetscInt *points;
    640:       PetscCall(ISGetIndices(labelIS, &points));
    641:       for (PetscInt i = 0; i < num_points; i++) {
    642:         PetscInt point;
    643:         PetscCall(DMGetEnclosurePoint(dm, odm, enc, points[i], &point));
    644:         if (pStart <= point && point < pEnd) {
    645:           ls = point;
    646:           if (ds) {
    647:             // If this is a face of a cohesive cell, then prefer that DS
    648:             if (height == 1) {
    649:               const PetscInt *supp;
    650:               PetscInt        suppSize;
    651:               DMPolytopeType  ct;

    653:               PetscCall(DMPlexGetSupport(dm, ls, &supp));
    654:               PetscCall(DMPlexGetSupportSize(dm, ls, &suppSize));
    655:               for (PetscInt s = 0; s < suppSize; ++s) {
    656:                 PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
    657:                 if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
    658:                   ls = supp[s];
    659:                   break;
    660:                 }
    661:               }
    662:             }
    663:             PetscCall(DMGetCellDS(dm, ls, ds, NULL));
    664:           }
    665:           if (ls >= 0) break;
    666:         }
    667:       }
    668:       PetscCall(ISRestoreIndices(labelIS, &points));
    669:     }
    670:     PetscCall(ISDestroy(&labelIS));
    671:     if (ls >= 0) break;
    672:   }
    673:   if (point) *point = ls;
    674:   PetscCall(DMDestroy(&plex));
    675:   PetscFunctionReturn(PETSC_SUCCESS);
    676: }

    678: /*
    679:   This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM

    681:   There are several different scenarios:

    683:   1) Volumetric mesh with volumetric auxiliary data

    685:      Here minHeight=0 since we loop over cells.

    687:   2) Boundary mesh with boundary auxiliary data

    689:      Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.

    691:   3) Volumetric mesh with boundary auxiliary data

    693:      Here minHeight=1 and auxbd=PETSC_TRUE since we loop over faces and use data only supported on those faces. This is common when imposing Dirichlet boundary conditions.

    695:   4) Volumetric input mesh with boundary output mesh

    697:      Here we must get a subspace for the input DS

    699:   The maxHeight is used to support enforcement of constraints in DMForest.

    701:   If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.

    703:   If we are using an input field (DM_BC_ESSENTIAL_FIELD or DM_BC_NATURAL_FIELD), we need to evaluate it at all the quadrature points of the dual basis functionals.
    704:     - We use effectiveHeight to mean the height above our incoming DS. For example, if the DS is for a submesh then the effective height is zero, whereas if the DS
    705:       is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When the effective height is nonzero, we need to extract
    706:       dual spaces for the boundary from our input spaces.
    707:     - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.

    709:   We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration

    711:   If we have a label, we iterate over those points. This will probably break the maxHeight functionality since we do not check the height of those points.
    712: */
    713: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], DMBoundaryConditionType type, PetscVoidFn **funcs, void **ctxs, InsertMode mode, Vec localX)
    714: {
    715:   DM               plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
    716:   DMEnclosureType  encIn, encAux;
    717:   PetscDS          ds = NULL, dsIn = NULL, dsAux = NULL;
    718:   Vec              localA = NULL, tv;
    719:   IS               fieldIS;
    720:   PetscSection     section;
    721:   PetscDualSpace  *sp, *cellsp, *spIn, *cellspIn;
    722:   PetscTabulation *T = NULL, *TAux = NULL;
    723:   PetscInt        *Nc;
    724:   PetscInt         dim, dimEmbed, depth, pStart, pEnd, lStart = PETSC_DETERMINE, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, minHeightIn, minHeightAux, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
    725:   PetscBool       *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, isCohesiveIn = PETSC_FALSE, transform;
    726:   DMField          coordField;
    727:   DMLabel          depthLabel;
    728:   PetscQuadrature  allPoints = NULL;

    730:   PetscFunctionBegin;
    731:   if (localU) PetscCall(VecGetDM(localU, &dmIn));
    732:   else dmIn = dm;
    733:   PetscCall(DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA));
    734:   if (localA) PetscCall(VecGetDM(localA, &dmAux));
    735:   else dmAux = NULL;
    736:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    737:   PetscCall(DMConvert(dmIn, DMPLEX, &plexIn));
    738:   PetscCall(DMGetEnclosureRelation(dmIn, dm, &encIn));
    739:   PetscCall(DMGetEnclosureRelation(dmAux, dm, &encAux));
    740:   PetscCall(DMGetDimension(dm, &dim));
    741:   PetscCall(DMPlexGetVTKCellHeight(plex, &minHeight));
    742:   PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
    743:   PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));
    744:   PetscCall(DMHasBasisTransform(dm, &transform));
    745:   /* Auxiliary information can only be used with interpolation of field functions */
    746:   if (dmAux) {
    747:     PetscCall(DMConvert(dmAux, DMPLEX, &plexAux));
    748:     if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) PetscCheck(localA, PETSC_COMM_SELF, PETSC_ERR_USER, "Missing localA vector");
    749:   }
    750:   if (localU && localU != localX) PetscCall(DMPlexInsertBoundaryValues(plexIn, PETSC_TRUE, localU, time, NULL, NULL, NULL));
    751:   PetscCall(DMGetCoordinateField(dm, &coordField));
    752:   PetscCheck(coordField, PETSC_COMM_SELF, PETSC_ERR_USER, "DM must have a coordinate field");
    753:   /**** No collective calls below this point ****/
    754:   /* Determine height for iteration of all meshes */
    755:   {
    756:     DMPolytopeType ct, ctIn, ctAux;
    757:     PetscInt       p, pStartIn, pStartAux, pEndAux;
    758:     PetscInt       dim = -1, dimIn = -1, dimAux = -1;

    760:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd));
    761:     if (pEnd > pStart) {
    762:       PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL));
    763:       p = lStart < 0 ? pStart : lStart;
    764:       PetscCall(DMPlexGetCellType(plex, p, &ct));
    765:       dim = DMPolytopeTypeGetDim(ct);
    766:       PetscCall(DMPlexGetVTKCellHeight(plexIn, &minHeightIn));
    767:       PetscCall(DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL));
    768:       PetscCall(DMPlexGetCellType(plexIn, pStartIn, &ctIn));
    769:       dimIn = DMPolytopeTypeGetDim(ctIn);
    770:       if (dmAux) {
    771:         PetscCall(DMPlexGetVTKCellHeight(plexAux, &minHeightAux));
    772:         PetscCall(DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux));
    773:         if (pStartAux < pEndAux) {
    774:           PetscCall(DMPlexGetCellType(plexAux, pStartAux, &ctAux));
    775:           dimAux = DMPolytopeTypeGetDim(ctAux);
    776:         }
    777:       } else dimAux = dim;
    778:     } else {
    779:       PetscCall(DMDestroy(&plex));
    780:       PetscCall(DMDestroy(&plexIn));
    781:       if (dmAux) PetscCall(DMDestroy(&plexAux));
    782:       PetscFunctionReturn(PETSC_SUCCESS);
    783:     }
    784:     if (dim < 0) {
    785:       DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;

    787:       /* Fall back to determination based on being a submesh */
    788:       PetscCall(DMPlexGetSubpointMap(plex, &spmap));
    789:       PetscCall(DMPlexGetSubpointMap(plexIn, &spmapIn));
    790:       if (plexAux) PetscCall(DMPlexGetSubpointMap(plexAux, &spmapAux));
    791:       dim    = spmap ? 1 : 0;
    792:       dimIn  = spmapIn ? 1 : 0;
    793:       dimAux = spmapAux ? 1 : 0;
    794:     }
    795:     {
    796:       PetscInt dimProj   = PetscMin(PetscMin(dim, dimIn), dimAux < 0 ? PETSC_INT_MAX : dimAux);
    797:       PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;

    799:       PetscCheck(PetscAbsInt(dimProj - dim) <= 1 && PetscAbsInt(dimProj - dimIn) <= 1 && PetscAbsInt(dimProj - dimAuxEff) <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not currently support differences of more than 1 in dimension");
    800:       if (dimProj < dim) minHeight = 1;
    801:       htInc    = dim - dimProj;
    802:       htIncIn  = dimIn - dimProj;
    803:       htIncAux = dimAuxEff - dimProj;
    804:     }
    805:   }
    806:   PetscCall(DMPlexGetDepth(plex, &depth));
    807:   PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
    808:   PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
    809:   maxHeight = PetscMax(maxHeight, minHeight);
    810:   PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
    811:   PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, NULL, &ds));
    812:   if (!ds) PetscCall(DMGetDS(dm, &ds));
    813:   PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, minHeight, NULL, &dsIn));
    814:   if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
    815:   if (!dsIn) {
    816:     if (encIn == DM_ENC_SUPERMESH) {
    817:       PetscInt p = pStart, pIn;

    819:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? p : lStart, &pIn));
    820:       // If the input mesh is higher dimensional than the output mesh, get a cell from the output mesh
    821:       if (htIncIn) PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &p, NULL));
    822:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? p : lStart, &pIn));
    823:       PetscCall(DMGetCellDS(dmIn, pIn, &dsIn, NULL));
    824:     } else PetscCall(DMGetDS(dmIn, &dsIn));
    825:   }
    826:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    827:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    828:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    829:   if (isCohesiveIn) --htIncIn; // Should be rearranged
    830:   PetscCall(DMGetNumFields(dm, &NfTot));
    831:   PetscCall(DMFindRegionNum(dm, ds, &regionNum));
    832:   PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
    833:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    834:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
    835:   PetscCall(DMGetLocalSection(dm, &section));
    836:   if (dmAux) {
    837:     if (encAux == DM_ENC_SUPERMESH) {
    838:       PetscInt p = pStart, pAux;

    840:       // If the auxiliary mesh is higher dimensional than the output mesh, get a cell from the output mesh
    841:       if (htIncAux) PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &p, NULL));
    842:       PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, lStart < 0 ? p : lStart, &pAux));
    843:       PetscCall(DMGetCellDS(dmAux, pAux, &dsAux, NULL));
    844:     } else PetscCall(DMGetDS(dmAux, &dsAux));
    845:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    846:   }
    847:   PetscCall(PetscDSGetComponents(ds, &Nc));
    848:   PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
    849:   if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
    850:   else {
    851:     cellsp   = sp;
    852:     cellspIn = spIn;
    853:   }
    854:   /* Get cell dual spaces */
    855:   for (f = 0; f < Nf; ++f) {
    856:     PetscDiscType disctype;

    858:     PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
    859:     if (disctype == PETSC_DISC_FE) {
    860:       PetscFE fe;

    862:       isFE[f] = PETSC_TRUE;
    863:       hasFE   = PETSC_TRUE;
    864:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
    865:       PetscCall(PetscFEGetDualSpace(fe, &cellsp[f]));
    866:     } else if (disctype == PETSC_DISC_FV) {
    867:       PetscFV fv;

    869:       isFE[f] = PETSC_FALSE;
    870:       hasFV   = PETSC_TRUE;
    871:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fv));
    872:       PetscCall(PetscFVGetDualSpace(fv, &cellsp[f]));
    873:     } else {
    874:       isFE[f]   = PETSC_FALSE;
    875:       cellsp[f] = NULL;
    876:     }
    877:   }
    878:   for (f = 0; f < NfIn; ++f) {
    879:     PetscDiscType disctype;

    881:     PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
    882:     if (disctype == PETSC_DISC_FE) {
    883:       PetscFE fe;

    885:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
    886:       PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
    887:     } else if (disctype == PETSC_DISC_FV) {
    888:       PetscFV fv;

    890:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv));
    891:       PetscCall(PetscFVGetDualSpace(fv, &cellspIn[f]));
    892:     } else {
    893:       cellspIn[f] = NULL;
    894:     }
    895:   }
    896:   for (f = 0; f < Nf; ++f) {
    897:     if (!htInc) {
    898:       sp[f] = cellsp[f];
    899:     } else PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]));
    900:   }
    901:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
    902:     PetscFE          fem, subfem;
    903:     PetscDiscType    disctype;
    904:     const PetscReal *points;
    905:     PetscInt         numPoints, k;

    907:     PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
    908:     PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
    909:     PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
    910:     PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
    911:     for (f = 0; f < NfIn; ++f) {
    912:       if (!htIncIn) {
    913:         spIn[f] = cellspIn[f];
    914:       } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));

    916:       PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
    917:       if (disctype != PETSC_DISC_FE) continue;
    918:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem));
    919:       if (!htIncIn) {
    920:         subfem = fem;
    921:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncIn, &subfem));
    922:       PetscCall(PetscDSGetJetDegree(dsIn, f, &k));
    923:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, k, &T[f]));
    924:     }
    925:     for (f = 0; f < NfAux; ++f) {
    926:       PetscCall(PetscDSGetDiscType_Internal(dsAux, f, &disctype));
    927:       if (disctype != PETSC_DISC_FE) continue;
    928:       PetscCall(PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem));
    929:       if (!htIncAux) {
    930:         subfem = fem;
    931:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncAux, &subfem));
    932:       PetscCall(PetscDSGetJetDegree(dsAux, f, &k));
    933:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, k, &TAux[f]));
    934:     }
    935:   }
    936:   /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
    937:   for (h = minHeight; h <= maxHeight; h++) {
    938:     PetscInt     hEff     = h - minHeight + htInc;
    939:     PetscInt     hEffIn   = h - minHeight + htIncIn;
    940:     PetscInt     hEffAux  = h - minHeight + htIncAux;
    941:     PetscDS      dsEff    = ds;
    942:     PetscDS      dsEffIn  = dsIn;
    943:     PetscDS      dsEffAux = dsAux;
    944:     PetscScalar *values;
    945:     PetscBool   *fieldActive;
    946:     PetscInt     maxDegree;
    947:     PetscInt     p, spDim, totDim, numValues;
    948:     IS           heightIS;

    950:     if (h > minHeight) {
    951:       for (f = 0; f < Nf; ++f) PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]));
    952:     }
    953:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd));
    954:     PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL));
    955:     PetscCall(DMLabelGetStratumIS(depthLabel, depth - h, &heightIS));
    956:     if (pEnd <= pStart) {
    957:       PetscCall(ISDestroy(&heightIS));
    958:       continue;
    959:     }
    960:     /* Compute totDim, the number of dofs in the closure of a point at this height */
    961:     totDim = 0;
    962:     for (f = 0; f < Nf; ++f) {
    963:       PetscBool cohesive;

    965:       if (!sp[f]) continue;
    966:       PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    967:       PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    968:       totDim += spDim;
    969:       if (isCohesive && !cohesive) totDim += spDim;
    970:     }
    971:     p = lStart < 0 ? pStart : lStart;
    972:     PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
    973:     PetscCheck(numValues == totDim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The output section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT " in [%" PetscInt_FMT ", %" PetscInt_FMT "]", p, numValues, totDim, h, minHeight, maxHeight);
    974:     if (!totDim) {
    975:       PetscCall(ISDestroy(&heightIS));
    976:       continue;
    977:     }
    978:     if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
    979:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
    980:     if (localU) {
    981:       PetscInt totDimIn, pIn, numValuesIn;

    983:       totDimIn = 0;
    984:       for (f = 0; f < NfIn; ++f) {
    985:         PetscBool cohesive;

    987:         if (!spIn[f]) continue;
    988:         PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
    989:         PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
    990:         totDimIn += spDim;
    991:         if (isCohesiveIn && !cohesive) totDimIn += spDim;
    992:       }
    993:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
    994:       PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
    995:       // TODO We could check that pIn is a cohesive cell for this check
    996:       PetscCheck(isCohesiveIn || (numValuesIn == totDimIn), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The input section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT, pIn, numValuesIn, totDimIn, htIncIn);
    997:       if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
    998:     }
    999:     if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
    1000:     /* Loop over points at this height */
    1001:     PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
    1002:     PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
    1003:     {
    1004:       const PetscInt *fields;

    1006:       PetscCall(ISGetIndices(fieldIS, &fields));
    1007:       for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
    1008:       for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
    1009:       PetscCall(ISRestoreIndices(fieldIS, &fields));
    1010:     }
    1011:     if (label) {
    1012:       PetscInt i;

    1014:       for (i = 0; i < numIds; ++i) {
    1015:         IS              pointIS, isectIS;
    1016:         const PetscInt *points;
    1017:         PetscInt        n;
    1018:         PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
    1019:         PetscQuadrature quad = NULL;

    1021:         PetscCall(DMLabelGetStratumIS(label, ids[i], &pointIS));
    1022:         if (!pointIS) continue; /* No points with that id on this process */
    1023:         PetscCall(ISIntersect(pointIS, heightIS, &isectIS));
    1024:         PetscCall(ISDestroy(&pointIS));
    1025:         if (!isectIS) continue;
    1026:         PetscCall(ISGetLocalSize(isectIS, &n));
    1027:         PetscCall(ISGetIndices(isectIS, &points));
    1028:         PetscCall(DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree));
    1029:         if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad));
    1030:         if (!quad) {
    1031:           if (!h && allPoints) {
    1032:             quad      = allPoints;
    1033:             allPoints = NULL;
    1034:           } else {
    1035:             PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad));
    1036:           }
    1037:         }
    1038:         PetscFEGeomMode geommode = htInc && h == minHeight ? PETSC_FEGEOM_BOUNDARY : PETSC_FEGEOM_BASIC;

    1040:         if (n) {
    1041:           PetscInt depth, dep;

    1043:           PetscCall(DMPlexGetDepth(dm, &depth));
    1044:           PetscCall(DMPlexGetPointDepth(dm, points[0], &dep));
    1045:           if (dep < depth && h == minHeight) geommode = PETSC_FEGEOM_BOUNDARY;
    1046:         }
    1047:         PetscCall(DMFieldCreateFEGeom(coordField, isectIS, quad, geommode, &fegeom));
    1048:         for (p = 0; p < n; ++p) {
    1049:           const PetscInt point = points[p];

    1051:           PetscCall(PetscArrayzero(values, numValues));
    1052:           PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
    1053:           PetscCall(DMPlexSetActivePoint(dm, point));
    1054:           PetscCall(DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsEffIn, plexAux, encAux, dsEffAux, chunkgeom, htInc, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values));
    1055:           if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
    1056:           PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
    1057:         }
    1058:         PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
    1059:         PetscCall(PetscFEGeomDestroy(&fegeom));
    1060:         PetscCall(PetscQuadratureDestroy(&quad));
    1061:         PetscCall(ISRestoreIndices(isectIS, &points));
    1062:         PetscCall(ISDestroy(&isectIS));
    1063:       }
    1064:     } else {
    1065:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
    1066:       PetscQuadrature quad = NULL;
    1067:       IS              pointIS;

    1069:       PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS));
    1070:       PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
    1071:       if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
    1072:       if (!quad) {
    1073:         if (!h && allPoints) {
    1074:           quad      = allPoints;
    1075:           allPoints = NULL;
    1076:         } else {
    1077:           PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad));
    1078:         }
    1079:       }
    1080:       PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_FEGEOM_BOUNDARY : PETSC_FEGEOM_BASIC, &fegeom));
    1081:       for (p = pStart; p < pEnd; ++p) {
    1082:         PetscCall(PetscArrayzero(values, numValues));
    1083:         PetscCall(PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom));
    1084:         PetscCall(DMPlexSetActivePoint(dm, p));
    1085:         PetscCall(DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsEffIn, plexAux, encAux, dsEffAux, chunkgeom, htInc, time, localU, localA, hasFE, hasFV, isFE, sp, p, T, TAux, type, funcs, ctxs, fieldActive, values));
    1086:         if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values));
    1087:         PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode));
    1088:       }
    1089:       PetscCall(PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom));
    1090:       PetscCall(PetscFEGeomDestroy(&fegeom));
    1091:       PetscCall(PetscQuadratureDestroy(&quad));
    1092:       PetscCall(ISDestroy(&pointIS));
    1093:     }
    1094:     PetscCall(ISDestroy(&heightIS));
    1095:     PetscCall(DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values));
    1096:     PetscCall(DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive));
    1097:   }
    1098:   /* Cleanup */
    1099:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
    1100:     for (f = 0; f < NfIn; ++f) PetscCall(PetscTabulationDestroy(&T[f]));
    1101:     for (f = 0; f < NfAux; ++f) PetscCall(PetscTabulationDestroy(&TAux[f]));
    1102:     PetscCall(PetscFree2(T, TAux));
    1103:   }
    1104:   PetscCall(PetscQuadratureDestroy(&allPoints));
    1105:   PetscCall(PetscFree3(isFE, sp, spIn));
    1106:   if (maxHeight > 0) PetscCall(PetscFree2(cellsp, cellspIn));
    1107:   PetscCall(DMDestroy(&plex));
    1108:   PetscCall(DMDestroy(&plexIn));
    1109:   if (dmAux) PetscCall(DMDestroy(&plexAux));
    1110:   PetscFunctionReturn(PETSC_SUCCESS);
    1111: }

    1113: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
    1114: {
    1115:   PetscFunctionBegin;
    1116:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (PetscVoidFn **)funcs, ctxs, mode, localX));
    1117:   PetscFunctionReturn(PETSC_SUCCESS);
    1118: }

    1120: PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
    1121: {
    1122:   PetscFunctionBegin;
    1123:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (PetscVoidFn **)funcs, ctxs, mode, localX));
    1124:   PetscFunctionReturn(PETSC_SUCCESS);
    1125: }

    1127: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
    1128: {
    1129:   PetscFunctionBegin;
    1130:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (PetscVoidFn **)funcs, NULL, mode, localX));
    1131:   PetscFunctionReturn(PETSC_SUCCESS);
    1132: }

    1134: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
    1135: {
    1136:   PetscFunctionBegin;
    1137:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (PetscVoidFn **)funcs, NULL, mode, localX));
    1138:   PetscFunctionReturn(PETSC_SUCCESS);
    1139: }

    1141: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
    1142: {
    1143:   PetscFunctionBegin;
    1144:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (PetscVoidFn **)funcs, NULL, mode, localX));
    1145:   PetscFunctionReturn(PETSC_SUCCESS);
    1146: }