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  coordDim, Nf, *Nc, f, spDim, d, v, tp;
     73:   PetscBool isAffine, isCohesive, transform;

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

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

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

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

    112:             if (dim != fegeom->dim) {
    113:               if (isCohesive) {
    114:                 /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
    115:                 for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
    116:                 refpoint = injpoint;
    117:               } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %" PetscInt_FMT " != %" PetscInt_FMT " dual basis spatial dimension", fegeom->dim, dim);
    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:           PetscCall((*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx));
    129:         }
    130:         /* Transform point evaluations pointEval[q,c] */
    131:         PetscCall(PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval));
    132:         PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    133:         PetscCall(DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x));
    134:         PetscCall(DMRestoreWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    135:         v += spDim;
    136:         if (isCohesive && !cohesive) {
    137:           for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    138:         }
    139:       } else {
    140:         for (d = 0; d < spDim; ++d, ++v) PetscCall(PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]));
    141:       }
    142:     } else {
    143:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    144:       if (isCohesive && !cohesive) {
    145:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    146:       }
    147:     }
    148:   }
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

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

    155:   Input Parameters:
    156: + dm             - The output DM
    157: . ds             - The output DS
    158: . dmIn           - The input DM
    159: . dsIn           - The input DS
    160: . dmAux          - The auxiliary DM, which is always for the input space
    161: . dsAux          - The auxiliary DS, which is always for the input space
    162: . time           - The time for this evaluation
    163: . localU         - The local solution
    164: . localA         - The local auziliary fields
    165: . cgeom          - The FE geometry for this point
    166: . sp             - The output PetscDualSpace for each field
    167: . p              - The point in the output DM
    168: . T              - Input basis and derivatives for each field tabulated on the quadrature points
    169: . TAux           - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
    170: . funcs          - The evaluation function for each field
    171: - ctxs           - The user context for each field

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

    176:   Level: developer

    178:   Note:
    179:   Not supported for FV

    181: .seealso: `DMProjectPoint_Field_Private()`
    182: */
    183: 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[])
    184: {
    185:   PetscSection       section, sectionAux = NULL;
    186:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
    187:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
    188:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
    189:   const PetscScalar *constants;
    190:   PetscReal         *x;
    191:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
    192:   PetscFEGeom        fegeom;
    193:   const PetscInt     dE = cgeom->dimEmbed, *cone, *ornt;
    194:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
    195:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
    196:   DMPolytopeType     qct;

    198:   PetscFunctionBeginHot;
    199:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    200:   PetscCall(PetscDSGetComponents(ds, &Nc));
    201:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    202:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    203:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    204:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
    205:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
    206:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
    207:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
    208:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
    209:   PetscCall(DMHasBasisTransform(dmIn, &transform));
    210:   PetscCall(DMGetLocalSection(dmIn, &section));
    211:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
    212:   // Get cohesive cell hanging off face
    213:   if (isCohesiveIn) {
    214:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
    215:     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)) {
    216:       DMPolytopeType  ct;
    217:       const PetscInt *support;
    218:       PetscInt        Ns, s;

    220:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
    221:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
    222:       for (s = 0; s < Ns; ++s) {
    223:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
    224:         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)) {
    225:           inp = support[s];
    226:           break;
    227:         }
    228:       }
    229:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
    230:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
    231:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
    232:       face[0] = 0;
    233:       face[1] = 0;
    234:     }
    235:   }
    236:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
    237:   if (dmAux) {
    238:     PetscInt subp;

    240:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
    241:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    242:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
    243:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    244:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    245:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
    246:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
    247:   }
    248:   /* Get values for closure */
    249:   isAffine        = cgeom->isAffine;
    250:   fegeom.dim      = cgeom->dim;
    251:   fegeom.dimEmbed = cgeom->dimEmbed;
    252:   if (isAffine) {
    253:     fegeom.v    = x;
    254:     fegeom.xi   = cgeom->xi;
    255:     fegeom.J    = cgeom->J;
    256:     fegeom.invJ = cgeom->invJ;
    257:     fegeom.detJ = cgeom->detJ;
    258:   }
    259:   for (f = 0, v = 0; f < Nf; ++f) {
    260:     PetscQuadrature  allPoints;
    261:     PetscInt         q, dim, numPoints;
    262:     const PetscReal *points;
    263:     PetscScalar     *pointEval;
    264:     PetscBool        cohesive;
    265:     DM               dm;

    267:     if (!sp[f]) continue;
    268:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    269:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    270:     if (!funcs[f]) {
    271:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    272:       if (isCohesive && !cohesive) {
    273:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    274:       }
    275:       continue;
    276:     }
    277:     const PetscInt ***perms;
    278:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
    279:     PetscCall(PetscDualSpaceGetSymmetries(sp[f], &perms, NULL));
    280:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
    281:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
    282:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    283:     for (q = 0; q < numPoints; ++q, ++tp) {
    284:       PetscInt qpt[2];

    286:       if (isCohesiveIn) {
    287:         qpt[0] = perms ? perms[0][ornt[0]][q] : q;
    288:         qpt[1] = perms ? perms[0][DMPolytopeTypeComposeOrientationInv(qct, ornt[1], 0)][q] : q;
    289:       }
    290:       if (isAffine) {
    291:         CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
    292:       } else {
    293:         fegeom.v    = &cgeom->v[tp * dE];
    294:         fegeom.J    = &cgeom->J[tp * dE * dE];
    295:         fegeom.invJ = &cgeom->invJ[tp * dE * dE];
    296:         fegeom.detJ = &cgeom->detJ[tp];
    297:       }
    298:       if (coefficients) {
    299:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
    300:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
    301:       }
    302:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
    303:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
    304:       (*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]);
    305:     }
    306:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    307:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    308:     v += spDim;
    309:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
    310:     if (isCohesive && !cohesive) {
    311:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    312:     }
    313:   }
    314:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
    315:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
    316:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
    317:   PetscFunctionReturn(PETSC_SUCCESS);
    318: }

    320: 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[])
    321: {
    322:   PetscSection       section, sectionAux = NULL;
    323:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
    324:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
    325:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
    326:   const PetscScalar *constants;
    327:   PetscReal         *x;
    328:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
    329:   PetscFEGeom        fegeom, cgeom;
    330:   const PetscInt     dE = fgeom->dimEmbed, *cone, *ornt;
    331:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
    332:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
    333:   DMPolytopeType     qct;

    335:   PetscFunctionBeginHot;
    336:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    337:   PetscCall(PetscDSGetComponents(ds, &Nc));
    338:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    339:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    340:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    341:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
    342:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
    343:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
    344:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
    345:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
    346:   PetscCall(DMHasBasisTransform(dmIn, &transform));
    347:   PetscCall(DMGetLocalSection(dmIn, &section));
    348:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
    349:   // Get cohesive cell hanging off face
    350:   if (isCohesiveIn) {
    351:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
    352:     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)) {
    353:       DMPolytopeType  ct;
    354:       const PetscInt *support;
    355:       PetscInt        Ns, s;

    357:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
    358:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
    359:       for (s = 0; s < Ns; ++s) {
    360:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
    361:         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)) {
    362:           inp = support[s];
    363:           break;
    364:         }
    365:       }
    366:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
    367:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
    368:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
    369:       face[0] = 0;
    370:       face[1] = 0;
    371:     }
    372:   }
    373:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
    374:   if (dmAux) {
    375:     PetscInt subp;

    377:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
    378:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    379:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
    380:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    381:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    382:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
    383:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
    384:   }
    385:   /* Get values for closure */
    386:   isAffine       = fgeom->isAffine;
    387:   fegeom.n       = NULL;
    388:   fegeom.J       = NULL;
    389:   fegeom.v       = NULL;
    390:   fegeom.xi      = NULL;
    391:   cgeom.dim      = fgeom->dim;
    392:   cgeom.dimEmbed = fgeom->dimEmbed;
    393:   if (isAffine) {
    394:     fegeom.v    = x;
    395:     fegeom.xi   = fgeom->xi;
    396:     fegeom.J    = fgeom->J;
    397:     fegeom.invJ = fgeom->invJ;
    398:     fegeom.detJ = fgeom->detJ;
    399:     fegeom.n    = fgeom->n;

    401:     cgeom.J    = fgeom->suppJ[0];
    402:     cgeom.invJ = fgeom->suppInvJ[0];
    403:     cgeom.detJ = fgeom->suppDetJ[0];
    404:   }
    405:   for (f = 0, v = 0; f < Nf; ++f) {
    406:     PetscQuadrature  allPoints;
    407:     PetscInt         q, dim, numPoints;
    408:     const PetscReal *points;
    409:     PetscScalar     *pointEval;
    410:     PetscBool        cohesive;
    411:     DM               dm;

    413:     if (!sp[f]) continue;
    414:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    415:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    416:     if (!funcs[f]) {
    417:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    418:       if (isCohesive && !cohesive) {
    419:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
    420:       }
    421:       continue;
    422:     }
    423:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
    424:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
    425:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
    426:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    427:     for (q = 0; q < numPoints; ++q, ++tp) {
    428:       PetscInt qpt[2];

    430:       if (isCohesiveIn) {
    431:         // These points are not integration quadratures, but dual space quadratures
    432:         // If they had multiple points we should match them from both sides, similar to hybrid residual eval
    433:         qpt[0] = qpt[1] = q;
    434:       }
    435:       if (isAffine) {
    436:         CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
    437:       } else {
    438:         fegeom.v    = &fgeom->v[tp * dE];
    439:         fegeom.J    = &fgeom->J[tp * dE * dE];
    440:         fegeom.invJ = &fgeom->invJ[tp * dE * dE];
    441:         fegeom.detJ = &fgeom->detJ[tp];
    442:         fegeom.n    = &fgeom->n[tp * dE];

    444:         cgeom.J    = &fgeom->suppJ[0][tp * dE * dE];
    445:         cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
    446:         cgeom.detJ = &fgeom->suppDetJ[0][tp];
    447:       }
    448:       /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
    449:       if (coefficients) {
    450:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
    451:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
    452:       }
    453:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
    454:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
    455:       (*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]);
    456:     }
    457:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
    458:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
    459:     v += spDim;
    460:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
    461:     if (isCohesive && !cohesive) {
    462:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
    463:     }
    464:   }
    465:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
    466:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
    467:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
    468:   PetscFunctionReturn(PETSC_SUCCESS);
    469: }

    471: 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, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
    472: {
    473:   PetscFVCellGeom fvgeom;
    474:   PetscInt        dim, dimEmbed;

    476:   PetscFunctionBeginHot;
    477:   PetscCall(DMGetDimension(dm, &dim));
    478:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
    479:   if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
    480:   switch (type) {
    481:   case DM_BC_ESSENTIAL:
    482:   case DM_BC_NATURAL:
    483:     PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
    484:     break;
    485:   case DM_BC_ESSENTIAL_FIELD:
    486:   case DM_BC_NATURAL_FIELD:
    487:     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));
    488:     break;
    489:   case DM_BC_ESSENTIAL_BD_FIELD:
    490:     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));
    491:     break;
    492:   default:
    493:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
    494:   }
    495:   PetscFunctionReturn(PETSC_SUCCESS);
    496: }

    498: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
    499: {
    500:   PetscReal *points;
    501:   PetscInt   f, numPoints;

    503:   PetscFunctionBegin;
    504:   if (!dim) {
    505:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
    506:     PetscFunctionReturn(PETSC_SUCCESS);
    507:   }
    508:   numPoints = 0;
    509:   for (f = 0; f < Nf; ++f) {
    510:     if (funcs[f]) {
    511:       PetscQuadrature fAllPoints;
    512:       PetscInt        fNumPoints;

    514:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
    515:       PetscCall(PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL));
    516:       numPoints += fNumPoints;
    517:     }
    518:   }
    519:   PetscCall(PetscMalloc1(dim * numPoints, &points));
    520:   numPoints = 0;
    521:   for (f = 0; f < Nf; ++f) {
    522:     if (funcs[f]) {
    523:       PetscQuadrature  fAllPoints;
    524:       PetscInt         qdim, fNumPoints, q;
    525:       const PetscReal *fPoints;

    527:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
    528:       PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
    529:       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);
    530:       for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
    531:       numPoints += fNumPoints;
    532:     }
    533:   }
    534:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
    535:   PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

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

    542:   Input Parameters:
    543: + dm     - the `DM`
    544: . odm    - the enclosing `DM`
    545: . label  - label for `DM` domain, or `NULL` for whole domain
    546: . numIds - the number of `ids`
    547: . ids    - An array of the label ids in sequence for the domain
    548: - height - Height of target cells in `DMPLEX` topology

    550:   Output Parameters:
    551: + point - the first labeled point
    552: - ds    - the `PetscDS` corresponding to the first labeled point

    554:   Level: developer

    556: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
    557: @*/
    558: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
    559: {
    560:   DM              plex;
    561:   DMEnclosureType enc;
    562:   PetscInt        ls = -1;

    564:   PetscFunctionBegin;
    565:   if (point) *point = -1;
    566:   if (!label) PetscFunctionReturn(PETSC_SUCCESS);
    567:   PetscCall(DMGetEnclosureRelation(dm, odm, &enc));
    568:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    569:   for (PetscInt i = 0; i < numIds; ++i) {
    570:     IS       labelIS;
    571:     PetscInt num_points, pStart, pEnd;
    572:     PetscCall(DMLabelGetStratumIS(label, ids[i], &labelIS));
    573:     if (!labelIS) continue; /* No points with that id on this process */
    574:     PetscCall(DMPlexGetHeightStratum(plex, height, &pStart, &pEnd));
    575:     PetscCall(ISGetSize(labelIS, &num_points));
    576:     if (num_points) {
    577:       const PetscInt *points;
    578:       PetscCall(ISGetIndices(labelIS, &points));
    579:       for (PetscInt i = 0; i < num_points; i++) {
    580:         PetscInt point;
    581:         PetscCall(DMGetEnclosurePoint(dm, odm, enc, points[i], &point));
    582:         if (pStart <= point && point < pEnd) {
    583:           ls = point;
    584:           if (ds) {
    585:             // If this is a face of a cohesive cell, then prefer that DS
    586:             if (height == 1) {
    587:               const PetscInt *supp;
    588:               PetscInt        suppSize;
    589:               DMPolytopeType  ct;

    591:               PetscCall(DMPlexGetSupport(dm, ls, &supp));
    592:               PetscCall(DMPlexGetSupportSize(dm, ls, &suppSize));
    593:               for (PetscInt s = 0; s < suppSize; ++s) {
    594:                 PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
    595:                 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)) {
    596:                   ls = supp[s];
    597:                   break;
    598:                 }
    599:               }
    600:             }
    601:             PetscCall(DMGetCellDS(dm, ls, ds, NULL));
    602:           }
    603:           if (ls >= 0) break;
    604:         }
    605:       }
    606:       PetscCall(ISRestoreIndices(labelIS, &points));
    607:     }
    608:     PetscCall(ISDestroy(&labelIS));
    609:     if (ls >= 0) break;
    610:   }
    611:   if (point) *point = ls;
    612:   PetscCall(DMDestroy(&plex));
    613:   PetscFunctionReturn(PETSC_SUCCESS);
    614: }

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

    619:   There are several different scenarios:

    621:   1) Volumetric mesh with volumetric auxiliary data

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

    625:   2) Boundary mesh with boundary auxiliary data

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

    629:   3) Volumetric mesh with boundary auxiliary data

    631:      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.

    633:   4) Volumetric input mesh with boundary output mesh

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

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

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

    641:   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.
    642:     - 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
    643:       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
    644:       dual spaces for the boundary from our input spaces.
    645:     - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.

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

    649:   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.
    650: */
    651: 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, void (**funcs)(void), void **ctxs, InsertMode mode, Vec localX)
    652: {
    653:   DM               plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
    654:   DMEnclosureType  encIn, encAux;
    655:   PetscDS          ds = NULL, dsIn = NULL, dsAux = NULL;
    656:   Vec              localA = NULL, tv;
    657:   IS               fieldIS;
    658:   PetscSection     section;
    659:   PetscDualSpace  *sp, *cellsp, *spIn, *cellspIn;
    660:   PetscTabulation *T = NULL, *TAux = NULL;
    661:   PetscInt        *Nc;
    662:   PetscInt         dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, minHeightIn, minHeightAux, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
    663:   PetscBool       *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, isCohesiveIn = PETSC_FALSE, transform;
    664:   DMField          coordField;
    665:   DMLabel          depthLabel;
    666:   PetscQuadrature  allPoints = NULL;

    668:   PetscFunctionBegin;
    669:   if (localU) PetscCall(VecGetDM(localU, &dmIn));
    670:   else dmIn = dm;
    671:   PetscCall(DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA));
    672:   if (localA) PetscCall(VecGetDM(localA, &dmAux));
    673:   else dmAux = NULL;
    674:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    675:   PetscCall(DMConvert(dmIn, DMPLEX, &plexIn));
    676:   PetscCall(DMGetEnclosureRelation(dmIn, dm, &encIn));
    677:   PetscCall(DMGetEnclosureRelation(dmAux, dm, &encAux));
    678:   PetscCall(DMGetDimension(dm, &dim));
    679:   PetscCall(DMPlexGetVTKCellHeight(plex, &minHeight));
    680:   PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
    681:   PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));
    682:   PetscCall(DMHasBasisTransform(dm, &transform));
    683:   /* Auxiliary information can only be used with interpolation of field functions */
    684:   if (dmAux) {
    685:     PetscCall(DMConvert(dmAux, DMPLEX, &plexAux));
    686:     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");
    687:   }
    688:   if (localU && localU != localX) PetscCall(DMPlexInsertBoundaryValues(plexIn, PETSC_TRUE, localU, time, NULL, NULL, NULL));
    689:   PetscCall(DMGetCoordinateField(dm, &coordField));
    690:   PetscCheck(coordField, PETSC_COMM_SELF, PETSC_ERR_USER, "DM must have a coordinate field");
    691:   /**** No collective calls below this point ****/
    692:   /* Determine height for iteration of all meshes */
    693:   {
    694:     DMPolytopeType ct, ctIn, ctAux;
    695:     PetscInt       lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
    696:     PetscInt       dim = -1, dimIn = -1, dimAux = -1;

    698:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd));
    699:     if (pEnd > pStart) {
    700:       PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL));
    701:       p = lStart < 0 ? pStart : lStart;
    702:       PetscCall(DMPlexGetCellType(plex, p, &ct));
    703:       dim = DMPolytopeTypeGetDim(ct);
    704:       PetscCall(DMPlexGetVTKCellHeight(plexIn, &minHeightIn));
    705:       PetscCall(DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL));
    706:       PetscCall(DMPlexGetCellType(plexIn, pStartIn, &ctIn));
    707:       dimIn = DMPolytopeTypeGetDim(ctIn);
    708:       if (dmAux) {
    709:         PetscCall(DMPlexGetVTKCellHeight(plexAux, &minHeightAux));
    710:         PetscCall(DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux));
    711:         if (pStartAux < pEndAux) {
    712:           PetscCall(DMPlexGetCellType(plexAux, pStartAux, &ctAux));
    713:           dimAux = DMPolytopeTypeGetDim(ctAux);
    714:         }
    715:       } else dimAux = dim;
    716:     } else {
    717:       PetscCall(DMDestroy(&plex));
    718:       PetscCall(DMDestroy(&plexIn));
    719:       if (dmAux) PetscCall(DMDestroy(&plexAux));
    720:       PetscFunctionReturn(PETSC_SUCCESS);
    721:     }
    722:     if (dim < 0) {
    723:       DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;

    725:       /* Fall back to determination based on being a submesh */
    726:       PetscCall(DMPlexGetSubpointMap(plex, &spmap));
    727:       PetscCall(DMPlexGetSubpointMap(plexIn, &spmapIn));
    728:       if (plexAux) PetscCall(DMPlexGetSubpointMap(plexAux, &spmapAux));
    729:       dim    = spmap ? 1 : 0;
    730:       dimIn  = spmapIn ? 1 : 0;
    731:       dimAux = spmapAux ? 1 : 0;
    732:     }
    733:     {
    734:       PetscInt dimProj   = PetscMin(PetscMin(dim, dimIn), dimAux < 0 ? PETSC_INT_MAX : dimAux);
    735:       PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;

    737:       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");
    738:       if (dimProj < dim) minHeight = 1;
    739:       htInc    = dim - dimProj;
    740:       htIncIn  = dimIn - dimProj;
    741:       htIncAux = dimAuxEff - dimProj;
    742:     }
    743:   }
    744:   PetscCall(DMPlexGetDepth(plex, &depth));
    745:   PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
    746:   PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
    747:   maxHeight = PetscMax(maxHeight, minHeight);
    748:   PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
    749:   PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, NULL, &ds));
    750:   if (!ds) PetscCall(DMGetDS(dm, &ds));
    751:   PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, minHeight, NULL, &dsIn));
    752:   if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
    753:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    754:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
    755:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
    756:   if (isCohesiveIn) --htIncIn; // Should be rearranged
    757:   PetscCall(DMGetNumFields(dm, &NfTot));
    758:   PetscCall(DMFindRegionNum(dm, ds, &regionNum));
    759:   PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
    760:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    761:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
    762:   PetscCall(DMGetLocalSection(dm, &section));
    763:   if (dmAux) {
    764:     PetscCall(DMGetDS(dmAux, &dsAux));
    765:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    766:   }
    767:   PetscCall(PetscDSGetComponents(ds, &Nc));
    768:   PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
    769:   if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
    770:   else {
    771:     cellsp   = sp;
    772:     cellspIn = spIn;
    773:   }
    774:   /* Get cell dual spaces */
    775:   for (f = 0; f < Nf; ++f) {
    776:     PetscDiscType disctype;

    778:     PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
    779:     if (disctype == PETSC_DISC_FE) {
    780:       PetscFE fe;

    782:       isFE[f] = PETSC_TRUE;
    783:       hasFE   = PETSC_TRUE;
    784:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
    785:       PetscCall(PetscFEGetDualSpace(fe, &cellsp[f]));
    786:     } else if (disctype == PETSC_DISC_FV) {
    787:       PetscFV fv;

    789:       isFE[f] = PETSC_FALSE;
    790:       hasFV   = PETSC_TRUE;
    791:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fv));
    792:       PetscCall(PetscFVGetDualSpace(fv, &cellsp[f]));
    793:     } else {
    794:       isFE[f]   = PETSC_FALSE;
    795:       cellsp[f] = NULL;
    796:     }
    797:   }
    798:   for (f = 0; f < NfIn; ++f) {
    799:     PetscDiscType disctype;

    801:     PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
    802:     if (disctype == PETSC_DISC_FE) {
    803:       PetscFE fe;

    805:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
    806:       PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
    807:     } else if (disctype == PETSC_DISC_FV) {
    808:       PetscFV fv;

    810:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv));
    811:       PetscCall(PetscFVGetDualSpace(fv, &cellspIn[f]));
    812:     } else {
    813:       cellspIn[f] = NULL;
    814:     }
    815:   }
    816:   for (f = 0; f < Nf; ++f) {
    817:     if (!htInc) {
    818:       sp[f] = cellsp[f];
    819:     } else PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]));
    820:   }
    821:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
    822:     PetscFE          fem, subfem;
    823:     PetscDiscType    disctype;
    824:     const PetscReal *points;
    825:     PetscInt         numPoints;

    827:     PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
    828:     PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
    829:     PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
    830:     PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
    831:     for (f = 0; f < NfIn; ++f) {
    832:       if (!htIncIn) {
    833:         spIn[f] = cellspIn[f];
    834:       } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));

    836:       PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
    837:       if (disctype != PETSC_DISC_FE) continue;
    838:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem));
    839:       if (!htIncIn) {
    840:         subfem = fem;
    841:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncIn, &subfem));
    842:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]));
    843:     }
    844:     for (f = 0; f < NfAux; ++f) {
    845:       PetscCall(PetscDSGetDiscType_Internal(dsAux, f, &disctype));
    846:       if (disctype != PETSC_DISC_FE) continue;
    847:       PetscCall(PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem));
    848:       if (!htIncAux) {
    849:         subfem = fem;
    850:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncAux, &subfem));
    851:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]));
    852:     }
    853:   }
    854:   /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
    855:   for (h = minHeight; h <= maxHeight; h++) {
    856:     PetscInt     hEff     = h - minHeight + htInc;
    857:     PetscInt     hEffIn   = h - minHeight + htIncIn;
    858:     PetscInt     hEffAux  = h - minHeight + htIncAux;
    859:     PetscDS      dsEff    = ds;
    860:     PetscDS      dsEffIn  = dsIn;
    861:     PetscDS      dsEffAux = dsAux;
    862:     PetscScalar *values;
    863:     PetscBool   *fieldActive;
    864:     PetscInt     maxDegree;
    865:     PetscInt     pStart, pEnd, p, lStart, spDim, totDim, numValues;
    866:     IS           heightIS;

    868:     if (h > minHeight) {
    869:       for (f = 0; f < Nf; ++f) PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]));
    870:     }
    871:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd));
    872:     PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL));
    873:     PetscCall(DMLabelGetStratumIS(depthLabel, depth - h, &heightIS));
    874:     if (pEnd <= pStart) {
    875:       PetscCall(ISDestroy(&heightIS));
    876:       continue;
    877:     }
    878:     /* Compute totDim, the number of dofs in the closure of a point at this height */
    879:     totDim = 0;
    880:     for (f = 0; f < Nf; ++f) {
    881:       PetscBool cohesive;

    883:       if (!sp[f]) continue;
    884:       PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
    885:       PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
    886:       totDim += spDim;
    887:       if (isCohesive && !cohesive) totDim += spDim;
    888:     }
    889:     p = lStart < 0 ? pStart : lStart;
    890:     PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
    891:     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);
    892:     if (!totDim) {
    893:       PetscCall(ISDestroy(&heightIS));
    894:       continue;
    895:     }
    896:     if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
    897:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
    898:     if (localU) {
    899:       PetscInt totDimIn, pIn, numValuesIn;

    901:       totDimIn = 0;
    902:       for (f = 0; f < NfIn; ++f) {
    903:         PetscBool cohesive;

    905:         if (!spIn[f]) continue;
    906:         PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
    907:         PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
    908:         totDimIn += spDim;
    909:         if (isCohesiveIn && !cohesive) totDimIn += spDim;
    910:       }
    911:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
    912:       PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
    913:       // TODO We could check that pIn is a cohesive cell for this check
    914:       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);
    915:       if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
    916:     }
    917:     if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
    918:     /* Loop over points at this height */
    919:     PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
    920:     PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
    921:     {
    922:       const PetscInt *fields;

    924:       PetscCall(ISGetIndices(fieldIS, &fields));
    925:       for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
    926:       for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
    927:       PetscCall(ISRestoreIndices(fieldIS, &fields));
    928:     }
    929:     if (label) {
    930:       PetscInt i;

    932:       for (i = 0; i < numIds; ++i) {
    933:         IS              pointIS, isectIS;
    934:         const PetscInt *points;
    935:         PetscInt        n;
    936:         PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
    937:         PetscQuadrature quad = NULL;

    939:         PetscCall(DMLabelGetStratumIS(label, ids[i], &pointIS));
    940:         if (!pointIS) continue; /* No points with that id on this process */
    941:         PetscCall(ISIntersect(pointIS, heightIS, &isectIS));
    942:         PetscCall(ISDestroy(&pointIS));
    943:         if (!isectIS) continue;
    944:         PetscCall(ISGetLocalSize(isectIS, &n));
    945:         PetscCall(ISGetIndices(isectIS, &points));
    946:         PetscCall(DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree));
    947:         if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad));
    948:         if (!quad) {
    949:           if (!h && allPoints) {
    950:             quad      = allPoints;
    951:             allPoints = NULL;
    952:           } else {
    953:             PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad));
    954:           }
    955:         }
    956:         PetscBool computeFaceGeom = htInc && h == minHeight ? PETSC_TRUE : PETSC_FALSE;

    958:         if (n) {
    959:           PetscInt depth, dep;

    961:           PetscCall(DMPlexGetDepth(dm, &depth));
    962:           PetscCall(DMPlexGetPointDepth(dm, points[0], &dep));
    963:           if (dep < depth && h == minHeight) computeFaceGeom = PETSC_TRUE;
    964:         }
    965:         PetscCall(DMFieldCreateFEGeom(coordField, isectIS, quad, computeFaceGeom, &fegeom));
    966:         for (p = 0; p < n; ++p) {
    967:           const PetscInt point = points[p];

    969:           PetscCall(PetscArrayzero(values, numValues));
    970:           PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
    971:           PetscCall(DMPlexSetActivePoint(dm, point));
    972:           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));
    973:           if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
    974:           PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
    975:         }
    976:         PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
    977:         PetscCall(PetscFEGeomDestroy(&fegeom));
    978:         PetscCall(PetscQuadratureDestroy(&quad));
    979:         PetscCall(ISRestoreIndices(isectIS, &points));
    980:         PetscCall(ISDestroy(&isectIS));
    981:       }
    982:     } else {
    983:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
    984:       PetscQuadrature quad = NULL;
    985:       IS              pointIS;

    987:       PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS));
    988:       PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
    989:       if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
    990:       if (!quad) {
    991:         if (!h && allPoints) {
    992:           quad      = allPoints;
    993:           allPoints = NULL;
    994:         } else {
    995:           PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad));
    996:         }
    997:       }
    998:       PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom));
    999:       for (p = pStart; p < pEnd; ++p) {
    1000:         PetscCall(PetscArrayzero(values, numValues));
    1001:         PetscCall(PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom));
    1002:         PetscCall(DMPlexSetActivePoint(dm, p));
    1003:         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));
    1004:         if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values));
    1005:         PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode));
    1006:       }
    1007:       PetscCall(PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom));
    1008:       PetscCall(PetscFEGeomDestroy(&fegeom));
    1009:       PetscCall(PetscQuadratureDestroy(&quad));
    1010:       PetscCall(ISDestroy(&pointIS));
    1011:     }
    1012:     PetscCall(ISDestroy(&heightIS));
    1013:     PetscCall(DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values));
    1014:     PetscCall(DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive));
    1015:   }
    1016:   /* Cleanup */
    1017:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
    1018:     for (f = 0; f < NfIn; ++f) PetscCall(PetscTabulationDestroy(&T[f]));
    1019:     for (f = 0; f < NfAux; ++f) PetscCall(PetscTabulationDestroy(&TAux[f]));
    1020:     PetscCall(PetscFree2(T, TAux));
    1021:   }
    1022:   PetscCall(PetscQuadratureDestroy(&allPoints));
    1023:   PetscCall(PetscFree3(isFE, sp, spIn));
    1024:   if (maxHeight > 0) PetscCall(PetscFree2(cellsp, cellspIn));
    1025:   PetscCall(DMDestroy(&plex));
    1026:   PetscCall(DMDestroy(&plexIn));
    1027:   if (dmAux) PetscCall(DMDestroy(&plexAux));
    1028:   PetscFunctionReturn(PETSC_SUCCESS);
    1029: }

    1031: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
    1032: {
    1033:   PetscFunctionBegin;
    1034:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
    1035:   PetscFunctionReturn(PETSC_SUCCESS);
    1036: }

    1038: 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)
    1039: {
    1040:   PetscFunctionBegin;
    1041:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
    1042:   PetscFunctionReturn(PETSC_SUCCESS);
    1043: }

    1045: 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)
    1046: {
    1047:   PetscFunctionBegin;
    1048:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
    1049:   PetscFunctionReturn(PETSC_SUCCESS);
    1050: }

    1052: 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)
    1053: {
    1054:   PetscFunctionBegin;
    1055:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
    1056:   PetscFunctionReturn(PETSC_SUCCESS);
    1057: }

    1059: 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)
    1060: {
    1061:   PetscFunctionBegin;
    1062:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX));
    1063:   PetscFunctionReturn(PETSC_SUCCESS);
    1064: }