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: `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: `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: `DMProjectPoint_Field_Private()`
 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:   Note: Not supported for FV

178:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

449: /*@C
450:   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.

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

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

464:   Level: developer
465: @*/
466: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
467: {
468:   DM              plex;
469:   DMEnclosureType enc;
470:   PetscInt        ls = -1;

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

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

507:   There are several different scenarios:

509:   1) Volumetric mesh with volumetric auxiliary data

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

513:   2) Boundary mesh with boundary auxiliary data

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

517:   3) Volumetric mesh with boundary auxiliary data

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

521:   4) Volumetric input mesh with boundary output mesh

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

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

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

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

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

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

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

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

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

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

662:     PetscDSGetDiscType_Internal(ds, f, &disctype);
663:     if (disctype == PETSC_DISC_FE) {
664:       PetscFE fe;

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

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

685:     PetscDSGetDiscType_Internal(dsIn, f, &disctype);
686:     if (disctype == PETSC_DISC_FE) {
687:       PetscFE fe;

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

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

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

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

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

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

785:       totDimIn = 0;
786:       for (f = 0; f < NfIn; ++f) {
787:         PetscBool cohesive;

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

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

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

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

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

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

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

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

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

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

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