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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`
 19: @*/
 20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
 21: {
 23:   *point = ((DM_Plex *)dm->data)->activePoint;
 24:   return 0;
 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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexGetActivePoint()`
 39: @*/
 40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
 41: {
 43:   ((DM_Plex *)dm->data)->activePoint = point;
 44:   return 0;
 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:[](chapter_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;

 76:   DMGetCoordinateDim(dmIn, &coordDim);
 77:   DMHasBasisTransform(dmIn, &transform);
 78:   PetscDSGetNumFields(ds, &Nf);
 79:   PetscDSGetComponents(ds, &Nc);
 80:   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:     PetscDSGetCohesive(ds, f, &cohesive);
 89:     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:         PetscDualSpaceGetDM(sp[f], &rdm);
100:         PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
101:         PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
102:         DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
103:         DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x);
104:         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:             DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx);
126:             v0 = x;
127:           }
128:           (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx);
129:         }
130:         /* Transform point evaluations pointEval[q,c] */
131:         PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
132:         PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
133:         DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x);
134:         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) 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:   return 0;
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;
192:   PetscFEGeom        fegeom;
193:   const PetscInt     dE = cgeom->dimEmbed;
194:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
195:   PetscBool          isAffine, isCohesive, transform;

198:   PetscDSGetNumFields(ds, &Nf);
199:   PetscDSGetComponents(ds, &Nc);
200:   PetscDSIsCohesive(ds, &isCohesive);
201:   PetscDSGetNumFields(dsIn, &NfIn);
202:   PetscDSGetComponentOffsets(dsIn, &uOff);
203:   PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
204:   PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
205:   PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
206:   PetscDSGetConstants(dsIn, &numConstants, &constants);
207:   DMHasBasisTransform(dmIn, &transform);
208:   DMGetLocalSection(dmIn, &section);
209:   DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
210:   DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
211:   if (dmAux) {
212:     PetscInt subp;

214:     DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
215:     PetscDSGetNumFields(dsAux, &NfAux);
216:     DMGetLocalSection(dmAux, &sectionAux);
217:     PetscDSGetComponentOffsets(dsAux, &aOff);
218:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
219:     PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
220:     DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
221:   }
222:   /* Get values for closure */
223:   isAffine        = cgeom->isAffine;
224:   fegeom.dim      = cgeom->dim;
225:   fegeom.dimEmbed = cgeom->dimEmbed;
226:   if (isAffine) {
227:     fegeom.v    = x;
228:     fegeom.xi   = cgeom->xi;
229:     fegeom.J    = cgeom->J;
230:     fegeom.invJ = cgeom->invJ;
231:     fegeom.detJ = cgeom->detJ;
232:   }
233:   for (f = 0, v = 0; f < Nf; ++f) {
234:     PetscQuadrature  allPoints;
235:     PetscInt         q, dim, numPoints;
236:     const PetscReal *points;
237:     PetscScalar     *pointEval;
238:     PetscBool        cohesive;
239:     DM               dm;

241:     if (!sp[f]) continue;
242:     PetscDSGetCohesive(ds, f, &cohesive);
243:     PetscDualSpaceGetDimension(sp[f], &spDim);
244:     if (!funcs[f]) {
245:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
246:       if (isCohesive && !cohesive) {
247:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
248:       }
249:       continue;
250:     }
251:     PetscDualSpaceGetDM(sp[f], &dm);
252:     PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
253:     PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
254:     DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
255:     for (q = 0; q < numPoints; ++q, ++tp) {
256:       if (isAffine) {
257:         CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
258:       } else {
259:         fegeom.v    = &cgeom->v[tp * dE];
260:         fegeom.J    = &cgeom->J[tp * dE * dE];
261:         fegeom.invJ = &cgeom->invJ[tp * dE * dE];
262:         fegeom.detJ = &cgeom->detJ[tp];
263:       }
264:       PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
265:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
266:       if (transform) DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);
267:       (*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]);
268:     }
269:     PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
270:     DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
271:     v += spDim;
272:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
273:     if (isCohesive && !cohesive) {
274:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
275:     }
276:   }
277:   DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
278:   if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
279:   return 0;
280: }

282: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, 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[])
283: {
284:   PetscSection       section, sectionAux = NULL;
285:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
286:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
287:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
288:   const PetscScalar *constants;
289:   PetscReal         *x;
290:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
291:   PetscFEGeom        fegeom, cgeom;
292:   const PetscInt     dE = fgeom->dimEmbed;
293:   PetscInt           numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
294:   PetscBool          isAffine;

298:   PetscDSGetNumFields(ds, &Nf);
299:   PetscDSGetComponents(ds, &Nc);
300:   PetscDSGetComponentOffsets(ds, &uOff);
301:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
302:   PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
303:   PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
304:   PetscDSGetConstants(ds, &numConstants, &constants);
305:   DMGetLocalSection(dm, &section);
306:   DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
307:   if (dmAux) {
308:     PetscInt subp;

310:     DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
311:     PetscDSGetNumFields(dsAux, &NfAux);
312:     DMGetLocalSection(dmAux, &sectionAux);
313:     PetscDSGetComponentOffsets(dsAux, &aOff);
314:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
315:     PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
316:     DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
317:   }
318:   /* Get values for closure */
319:   isAffine       = fgeom->isAffine;
320:   fegeom.n       = NULL;
321:   fegeom.J       = NULL;
322:   fegeom.v       = NULL;
323:   fegeom.xi      = NULL;
324:   cgeom.dim      = fgeom->dim;
325:   cgeom.dimEmbed = fgeom->dimEmbed;
326:   if (isAffine) {
327:     fegeom.v    = x;
328:     fegeom.xi   = fgeom->xi;
329:     fegeom.J    = fgeom->J;
330:     fegeom.invJ = fgeom->invJ;
331:     fegeom.detJ = fgeom->detJ;
332:     fegeom.n    = fgeom->n;

334:     cgeom.J    = fgeom->suppJ[0];
335:     cgeom.invJ = fgeom->suppInvJ[0];
336:     cgeom.detJ = fgeom->suppDetJ[0];
337:   }
338:   for (f = 0, v = 0; f < Nf; ++f) {
339:     PetscQuadrature  allPoints;
340:     PetscInt         q, dim, numPoints;
341:     const PetscReal *points;
342:     PetscScalar     *pointEval;
343:     DM               dm;

345:     if (!sp[f]) continue;
346:     PetscDualSpaceGetDimension(sp[f], &spDim);
347:     if (!funcs[f]) {
348:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
349:       continue;
350:     }
351:     PetscDualSpaceGetDM(sp[f], &dm);
352:     PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
353:     PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
354:     DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
355:     for (q = 0; q < numPoints; ++q, ++tp) {
356:       if (isAffine) {
357:         CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
358:       } else {
359:         fegeom.v    = &fgeom->v[tp * dE];
360:         fegeom.J    = &fgeom->J[tp * dE * dE];
361:         fegeom.invJ = &fgeom->invJ[tp * dE * dE];
362:         fegeom.detJ = &fgeom->detJ[tp];
363:         fegeom.n    = &fgeom->n[tp * dE];

365:         cgeom.J    = &fgeom->suppJ[0][tp * dE * dE];
366:         cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
367:         cgeom.detJ = &fgeom->suppDetJ[0][tp];
368:       }
369:       /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
370:       PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
371:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
372:       (*funcs[f])(dE, Nf, 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]);
373:     }
374:     PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
375:     DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
376:     v += spDim;
377:   }
378:   DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
379:   if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
380:   return 0;
381: }

383: 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[])
384: {
385:   PetscFVCellGeom fvgeom;
386:   PetscInt        dim, dimEmbed;

389:   DMGetDimension(dm, &dim);
390:   DMGetCoordinateDim(dm, &dimEmbed);
391:   if (hasFV) DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);
392:   switch (type) {
393:   case DM_BC_ESSENTIAL:
394:   case DM_BC_NATURAL:
395:     DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode(**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values);
396:     break;
397:   case DM_BC_ESSENTIAL_FIELD:
398:   case DM_BC_NATURAL_FIELD:
399:     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);
400:     break;
401:   case DM_BC_ESSENTIAL_BD_FIELD:
402:     DMProjectPoint_BdField_Private(dm, ds, dmIn, 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);
403:     break;
404:   default:
405:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
406:   }
407:   return 0;
408: }

410: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
411: {
412:   PetscReal *points;
413:   PetscInt   f, numPoints;

415:   if (!dim) {
416:     PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
417:     return 0;
418:   }
419:   numPoints = 0;
420:   for (f = 0; f < Nf; ++f) {
421:     if (funcs[f]) {
422:       PetscQuadrature fAllPoints;
423:       PetscInt        fNumPoints;

425:       PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
426:       PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
427:       numPoints += fNumPoints;
428:     }
429:   }
430:   PetscMalloc1(dim * numPoints, &points);
431:   numPoints = 0;
432:   for (f = 0; f < Nf; ++f) {
433:     if (funcs[f]) {
434:       PetscQuadrature  fAllPoints;
435:       PetscInt         qdim, fNumPoints, q;
436:       const PetscReal *fPoints;

438:       PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
439:       PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
441:       for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
442:       numPoints += fNumPoints;
443:     }
444:   }
445:   PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
446:   PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL);
447:   return 0;
448: }

450: /*@C
451:   DMGetFirstLabeledPoint - Find first labeled point p_o in odm such that the corresponding point p in dm has the specified height. Return p and the corresponding ds.

453:   Input Parameters:
454:   dm - the `DM`
455:   odm - the enclosing `DM`
456:   label - label for `DM` domain, or NULL for whole domain
457:   numIds - the number of ids
458:   ids - An array of the label ids in sequence for the domain
459:   height - Height of target cells in `DMPLEX` topology

461:   Output Parameters:
462:   point - the first labeled point
463:   ds - the ds corresponding to the first labeled point

465:   Level: developer

467: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
468: @*/
469: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
470: {
471:   DM              plex;
472:   DMEnclosureType enc;
473:   PetscInt        ls = -1;

475:   if (point) *point = -1;
476:   if (!label) return 0;
477:   DMGetEnclosureRelation(dm, odm, &enc);
478:   DMConvert(dm, DMPLEX, &plex);
479:   for (PetscInt i = 0; i < numIds; ++i) {
480:     IS       labelIS;
481:     PetscInt num_points, pStart, pEnd;
482:     DMLabelGetStratumIS(label, ids[i], &labelIS);
483:     if (!labelIS) continue; /* No points with that id on this process */
484:     DMPlexGetHeightStratum(plex, height, &pStart, &pEnd);
485:     ISGetSize(labelIS, &num_points);
486:     if (num_points) {
487:       const PetscInt *points;
488:       ISGetIndices(labelIS, &points);
489:       for (PetscInt i = 0; i < num_points; i++) {
490:         PetscInt point;
491:         DMGetEnclosurePoint(dm, odm, enc, points[i], &point);
492:         if (pStart <= point && point < pEnd) {
493:           ls = point;
494:           if (ds) DMGetCellDS(dm, ls, ds);
495:         }
496:       }
497:       ISRestoreIndices(labelIS, &points);
498:     }
499:     ISDestroy(&labelIS);
500:     if (ls >= 0) break;
501:   }
502:   if (point) *point = ls;
503:   DMDestroy(&plex);
504:   return 0;
505: }

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

510:   There are several different scenarios:

512:   1) Volumetric mesh with volumetric auxiliary data

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

516:   2) Boundary mesh with boundary auxiliary data

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

520:   3) Volumetric mesh with boundary auxiliary data

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

524:   4) Volumetric input mesh with boundary output mesh

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

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

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

532:   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.
533:     - 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
534:       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
535:       dual spaces for the boundary from our input spaces.
536:     - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.

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

540:   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.
541: */
542: 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)
543: {
544:   DM               plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
545:   DMEnclosureType  encIn, encAux;
546:   PetscDS          ds = NULL, dsIn = NULL, dsAux = NULL;
547:   Vec              localA = NULL, tv;
548:   IS               fieldIS;
549:   PetscSection     section;
550:   PetscDualSpace  *sp, *cellsp, *spIn, *cellspIn;
551:   PetscTabulation *T = NULL, *TAux = NULL;
552:   PetscInt        *Nc;
553:   PetscInt         dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
554:   PetscBool       *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, transform;
555:   DMField          coordField;
556:   DMLabel          depthLabel;
557:   PetscQuadrature  allPoints = NULL;

559:   if (localU) VecGetDM(localU, &dmIn);
560:   else dmIn = dm;
561:   DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA);
562:   if (localA) VecGetDM(localA, &dmAux);
563:   else dmAux = NULL;
564:   DMConvert(dm, DMPLEX, &plex);
565:   DMConvert(dmIn, DMPLEX, &plexIn);
566:   DMGetEnclosureRelation(dmIn, dm, &encIn);
567:   DMGetEnclosureRelation(dmAux, dm, &encAux);
568:   DMGetDimension(dm, &dim);
569:   DMPlexGetVTKCellHeight(plex, &minHeight);
570:   DMGetBasisTransformDM_Internal(dm, &tdm);
571:   DMGetBasisTransformVec_Internal(dm, &tv);
572:   DMHasBasisTransform(dm, &transform);
573:   /* Auxiliary information can only be used with interpolation of field functions */
574:   if (dmAux) {
575:     DMConvert(dmAux, DMPLEX, &plexAux);
577:   }
578:   if (localU && localU != localX) DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);
579:   DMGetCoordinateField(dm, &coordField);
580:   /**** No collective calls below this point ****/
581:   /* Determine height for iteration of all meshes */
582:   {
583:     DMPolytopeType ct, ctIn, ctAux;
584:     PetscInt       minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
585:     PetscInt       dim = -1, dimIn = -1, dimAux = -1;

587:     DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
588:     if (pEnd > pStart) {
589:       DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL);
590:       p = lStart < 0 ? pStart : lStart;
591:       DMPlexGetCellType(plex, p, &ct);
592:       dim = DMPolytopeTypeGetDim(ct);
593:       DMPlexGetVTKCellHeight(plexIn, &minHeightIn);
594:       DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL);
595:       DMPlexGetCellType(plexIn, pStartIn, &ctIn);
596:       dimIn = DMPolytopeTypeGetDim(ctIn);
597:       if (dmAux) {
598:         DMPlexGetVTKCellHeight(plexAux, &minHeightAux);
599:         DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux);
600:         if (pStartAux < pEndAux) {
601:           DMPlexGetCellType(plexAux, pStartAux, &ctAux);
602:           dimAux = DMPolytopeTypeGetDim(ctAux);
603:         }
604:       } else dimAux = dim;
605:     } else {
606:       DMDestroy(&plex);
607:       DMDestroy(&plexIn);
608:       if (dmAux) DMDestroy(&plexAux);
609:       return 0;
610:     }
611:     if (dim < 0) {
612:       DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;

614:       /* Fall back to determination based on being a submesh */
615:       DMPlexGetSubpointMap(plex, &spmap);
616:       DMPlexGetSubpointMap(plexIn, &spmapIn);
617:       if (plexAux) DMPlexGetSubpointMap(plexAux, &spmapAux);
618:       dim    = spmap ? 1 : 0;
619:       dimIn  = spmapIn ? 1 : 0;
620:       dimAux = spmapAux ? 1 : 0;
621:     }
622:     {
623:       PetscInt dimProj   = PetscMin(PetscMin(dim, dimIn), (dimAux < 0 ? PETSC_MAX_INT : dimAux));
624:       PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;

627:       if (dimProj < dim) minHeight = 1;
628:       htInc    = dim - dimProj;
629:       htIncIn  = dimIn - dimProj;
630:       htIncAux = dimAuxEff - dimProj;
631:     }
632:   }
633:   DMPlexGetDepth(plex, &depth);
634:   DMPlexGetDepthLabel(plex, &depthLabel);
635:   DMPlexGetMaxProjectionHeight(plex, &maxHeight);
636:   maxHeight = PetscMax(maxHeight, minHeight);
638:   DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, 0, NULL, &ds);
639:   if (!ds) DMGetDS(dm, &ds);
640:   DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
641:   if (!dsIn) DMGetDS(dmIn, &dsIn);
642:   PetscDSGetNumFields(ds, &Nf);
643:   PetscDSGetNumFields(dsIn, &NfIn);
644:   DMGetNumFields(dm, &NfTot);
645:   DMFindRegionNum(dm, ds, &regionNum);
646:   DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
647:   PetscDSIsCohesive(ds, &isCohesive);
648:   DMGetCoordinateDim(dm, &dimEmbed);
649:   DMGetLocalSection(dm, &section);
650:   if (dmAux) {
651:     DMGetDS(dmAux, &dsAux);
652:     PetscDSGetNumFields(dsAux, &NfAux);
653:   }
654:   PetscDSGetComponents(ds, &Nc);
655:   PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
656:   if (maxHeight > 0) PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);
657:   else {
658:     cellsp   = sp;
659:     cellspIn = spIn;
660:   }
661:   /* Get cell dual spaces */
662:   for (f = 0; f < Nf; ++f) {
663:     PetscDiscType disctype;

665:     PetscDSGetDiscType_Internal(ds, f, &disctype);
666:     if (disctype == PETSC_DISC_FE) {
667:       PetscFE fe;

669:       isFE[f] = PETSC_TRUE;
670:       hasFE   = PETSC_TRUE;
671:       PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
672:       PetscFEGetDualSpace(fe, &cellsp[f]);
673:     } else if (disctype == PETSC_DISC_FV) {
674:       PetscFV fv;

676:       isFE[f] = PETSC_FALSE;
677:       hasFV   = PETSC_TRUE;
678:       PetscDSGetDiscretization(ds, f, (PetscObject *)&fv);
679:       PetscFVGetDualSpace(fv, &cellsp[f]);
680:     } else {
681:       isFE[f]   = PETSC_FALSE;
682:       cellsp[f] = NULL;
683:     }
684:   }
685:   for (f = 0; f < NfIn; ++f) {
686:     PetscDiscType disctype;

688:     PetscDSGetDiscType_Internal(dsIn, f, &disctype);
689:     if (disctype == PETSC_DISC_FE) {
690:       PetscFE fe;

692:       PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe);
693:       PetscFEGetDualSpace(fe, &cellspIn[f]);
694:     } else if (disctype == PETSC_DISC_FV) {
695:       PetscFV fv;

697:       PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv);
698:       PetscFVGetDualSpace(fv, &cellspIn[f]);
699:     } else {
700:       cellspIn[f] = NULL;
701:     }
702:   }
703:   for (f = 0; f < Nf; ++f) {
704:     if (!htInc) {
705:       sp[f] = cellsp[f];
706:     } else PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);
707:   }
708:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
709:     PetscFE          fem, subfem;
710:     PetscDiscType    disctype;
711:     const PetscReal *points;
712:     PetscInt         numPoints;

715:     PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints);
716:     PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
717:     PetscMalloc2(NfIn, &T, NfAux, &TAux);
718:     for (f = 0; f < NfIn; ++f) {
719:       if (!htIncIn) {
720:         spIn[f] = cellspIn[f];
721:       } else PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);

723:       PetscDSGetDiscType_Internal(dsIn, f, &disctype);
724:       if (disctype != PETSC_DISC_FE) continue;
725:       PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem);
726:       if (!htIncIn) {
727:         subfem = fem;
728:       } else PetscFEGetHeightSubspace(fem, htIncIn, &subfem);
729:       PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
730:     }
731:     for (f = 0; f < NfAux; ++f) {
732:       PetscDSGetDiscType_Internal(dsAux, f, &disctype);
733:       if (disctype != PETSC_DISC_FE) continue;
734:       PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem);
735:       if (!htIncAux) {
736:         subfem = fem;
737:       } else PetscFEGetHeightSubspace(fem, htIncAux, &subfem);
738:       PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
739:     }
740:   }
741:   /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
742:   for (h = minHeight; h <= maxHeight; h++) {
743:     PetscInt     hEff     = h - minHeight + htInc;
744:     PetscInt     hEffIn   = h - minHeight + htIncIn;
745:     PetscInt     hEffAux  = h - minHeight + htIncAux;
746:     PetscDS      dsEff    = ds;
747:     PetscDS      dsEffIn  = dsIn;
748:     PetscDS      dsEffAux = dsAux;
749:     PetscScalar *values;
750:     PetscBool   *fieldActive;
751:     PetscInt     maxDegree;
752:     PetscInt     pStart, pEnd, p, lStart, spDim, totDim, numValues;
753:     IS           heightIS;

755:     if (h > minHeight) {
756:       for (f = 0; f < Nf; ++f) PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);
757:     }
758:     DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
759:     DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL);
760:     DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
761:     if (pEnd <= pStart) {
762:       ISDestroy(&heightIS);
763:       continue;
764:     }
765:     /* Compute totDim, the number of dofs in the closure of a point at this height */
766:     totDim = 0;
767:     for (f = 0; f < Nf; ++f) {
768:       PetscBool cohesive;

770:       if (!sp[f]) continue;
771:       PetscDSGetCohesive(ds, f, &cohesive);
772:       PetscDualSpaceGetDimension(sp[f], &spDim);
773:       totDim += spDim;
774:       if (isCohesive && !cohesive) totDim += spDim;
775:     }
776:     p = lStart < 0 ? pStart : lStart;
777:     DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
779:     if (!totDim) {
780:       ISDestroy(&heightIS);
781:       continue;
782:     }
783:     if (htInc) PetscDSGetHeightSubspace(ds, hEff, &dsEff);
784:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
785:     if (localU) {
786:       PetscInt totDimIn, pIn, numValuesIn;

788:       totDimIn = 0;
789:       for (f = 0; f < NfIn; ++f) {
790:         PetscBool cohesive;

792:         if (!spIn[f]) continue;
793:         PetscDSGetCohesive(dsIn, f, &cohesive);
794:         PetscDualSpaceGetDimension(spIn[f], &spDim);
795:         totDimIn += spDim;
796:         if (isCohesive && !cohesive) totDimIn += spDim;
797:       }
798:       DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
799:       DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
801:       if (htIncIn) PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);
802:     }
803:     if (htIncAux) PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);
804:     /* Loop over points at this height */
805:     DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
806:     DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
807:     {
808:       const PetscInt *fields;

810:       ISGetIndices(fieldIS, &fields);
811:       for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
812:       for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
813:       ISRestoreIndices(fieldIS, &fields);
814:     }
815:     if (label) {
816:       PetscInt i;

818:       for (i = 0; i < numIds; ++i) {
819:         IS              pointIS, isectIS;
820:         const PetscInt *points;
821:         PetscInt        n;
822:         PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
823:         PetscQuadrature quad = NULL;

825:         DMLabelGetStratumIS(label, ids[i], &pointIS);
826:         if (!pointIS) continue; /* No points with that id on this process */
827:         ISIntersect(pointIS, heightIS, &isectIS);
828:         ISDestroy(&pointIS);
829:         if (!isectIS) continue;
830:         ISGetLocalSize(isectIS, &n);
831:         ISGetIndices(isectIS, &points);
832:         DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree);
833:         if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad);
834:         if (!quad) {
835:           if (!h && allPoints) {
836:             quad      = allPoints;
837:             allPoints = NULL;
838:           } else {
839:             PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad);
840:           }
841:         }
842:         DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
843:         for (p = 0; p < n; ++p) {
844:           const PetscInt point = points[p];

846:           PetscArrayzero(values, numValues);
847:           PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom);
848:           DMPlexSetActivePoint(dm, point);
849:           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);
850:           if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);
851:           DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
852:         }
853:         PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom);
854:         PetscFEGeomDestroy(&fegeom);
855:         PetscQuadratureDestroy(&quad);
856:         ISRestoreIndices(isectIS, &points);
857:         ISDestroy(&isectIS);
858:       }
859:     } else {
860:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
861:       PetscQuadrature quad = NULL;
862:       IS              pointIS;

864:       ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS);
865:       DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
866:       if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
867:       if (!quad) {
868:         if (!h && allPoints) {
869:           quad      = allPoints;
870:           allPoints = NULL;
871:         } else {
872:           PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad);
873:         }
874:       }
875:       DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
876:       for (p = pStart; p < pEnd; ++p) {
877:         PetscArrayzero(values, numValues);
878:         PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom);
879:         DMPlexSetActivePoint(dm, p);
880:         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);
881:         if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);
882:         DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
883:       }
884:       PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom);
885:       PetscFEGeomDestroy(&fegeom);
886:       PetscQuadratureDestroy(&quad);
887:       ISDestroy(&pointIS);
888:     }
889:     ISDestroy(&heightIS);
890:     DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
891:     DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
892:   }
893:   /* Cleanup */
894:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
895:     for (f = 0; f < NfIn; ++f) PetscTabulationDestroy(&T[f]);
896:     for (f = 0; f < NfAux; ++f) PetscTabulationDestroy(&TAux[f]);
897:     PetscFree2(T, TAux);
898:   }
899:   PetscQuadratureDestroy(&allPoints);
900:   PetscFree3(isFE, sp, spIn);
901:   if (maxHeight > 0) PetscFree2(cellsp, cellspIn);
902:   DMDestroy(&plex);
903:   DMDestroy(&plexIn);
904:   if (dmAux) DMDestroy(&plexAux);
905:   return 0;
906: }

908: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
909: {
910:   DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
911:   return 0;
912: }

914: 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)
915: {
916:   DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
917:   return 0;
918: }

920: 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)
921: {
922:   DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
923:   return 0;
924: }

926: 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)
927: {
928:   DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
929:   return 0;
930: }

932: 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)
933: {
934:   DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX);
935:   return 0;
936: }