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, void (**funcs)(void), 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, void (**funcs)(void), 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, void (**funcs)(void), 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, (void (**)(void))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, (void (**)(void))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, (void (**)(void))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, (void (**)(void))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, (void (**)(void))funcs, NULL, mode, localX));
1145:   PetscFunctionReturn(PETSC_SUCCESS);
1146: }