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[],
 71:                                                   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
 72:                                                   PetscScalar values[])
 73: {
 74:   PetscInt       coordDim, Nf, *Nc, f, spDim, d, v, tp;
 75:   PetscBool      isAffine, isHybrid, transform;

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

 89:     if (!sp[f]) continue;
 90:     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:         PetscDualSpaceGetDM(sp[f],&rdm);
101:         PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
102:         PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
103:         DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
104:         DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
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 (isHybrid) {
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 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D 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) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
125:           (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
126:         }
127:         /* Transform point evaluations pointEval[q,c] */
128:         PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
129:         PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
130:         DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
131:         DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
132:         v += spDim;
133:         if (isHybrid && (f < Nf-1)) {
134:           for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
135:         }
136:       } else {
137:         for (d = 0; d < spDim; ++d, ++v) {
138:           PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
139:         }
140:       }
141:     } else {
142:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
143:       if (isHybrid && (f < Nf-1)) {
144:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
145:       }
146:     }
147:   }
148:   return(0);
149: }

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

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

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

175:   Note: Not supported for FV

177:   Level: developer

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

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

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

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

285: 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,
286:                                                      PetscTabulation *T, PetscTabulation *TAux,
287:                                                      void (**funcs)(PetscInt, PetscInt, PetscInt,
288:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
289:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
290:                                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
291:                                                      PetscScalar values[])
292: {
293:   PetscSection       section, sectionAux = NULL;
294:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
295:   PetscScalar       *coefficients   = NULL, *coefficientsAux   = NULL;
296:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
297:   const PetscScalar *constants;
298:   PetscReal         *x;
299:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
300:   PetscFEGeom        fegeom, cgeom;
301:   const PetscInt     dE = fgeom->dimEmbed;
302:   PetscInt           numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
303:   PetscBool          isAffine;
304:   PetscErrorCode     ierr;

307:   if (dm != dmIn) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM");
308:   PetscDSGetNumFields(ds, &Nf);
309:   PetscDSGetComponents(ds, &Nc);
310:   PetscDSGetComponentOffsets(ds, &uOff);
311:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
312:   PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
313:   PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
314:   PetscDSGetConstants(ds, &numConstants, &constants);
315:   DMGetLocalSection(dm, &section);
316:   DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
317:   if (dmAux) {
318:     PetscInt subp;

320:     DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
321:     PetscDSGetNumFields(dsAux, &NfAux);
322:     DMGetLocalSection(dmAux, &sectionAux);
323:     PetscDSGetComponentOffsets(dsAux, &aOff);
324:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
325:     PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
326:     DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
327:   }
328:   /* Get values for closure */
329:   isAffine = fgeom->isAffine;
330:   fegeom.n  = NULL;
331:   fegeom.J  = NULL;
332:   fegeom.v  = NULL;
333:   fegeom.xi = NULL;
334:   cgeom.dim      = fgeom->dim;
335:   cgeom.dimEmbed = fgeom->dimEmbed;
336:   if (isAffine) {
337:     fegeom.v    = x;
338:     fegeom.xi   = fgeom->xi;
339:     fegeom.J    = fgeom->J;
340:     fegeom.invJ = fgeom->invJ;
341:     fegeom.detJ = fgeom->detJ;
342:     fegeom.n    = fgeom->n;

344:     cgeom.J     = fgeom->suppJ[0];
345:     cgeom.invJ  = fgeom->suppInvJ[0];
346:     cgeom.detJ  = fgeom->suppDetJ[0];
347:   }
348:   for (f = 0, v = 0; f < Nf; ++f) {
349:     PetscQuadrature   allPoints;
350:     PetscInt          q, dim, numPoints;
351:     const PetscReal   *points;
352:     PetscScalar       *pointEval;
353:     DM                dm;

355:     if (!sp[f]) continue;
356:     PetscDualSpaceGetDimension(sp[f], &spDim);
357:     if (!funcs[f]) {
358:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
359:       continue;
360:     }
361:     PetscDualSpaceGetDM(sp[f],&dm);
362:     PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
363:     PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
364:     DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
365:     for (q = 0; q < numPoints; ++q, ++tp) {
366:       if (isAffine) {
367:         CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
368:       } else {
369:         fegeom.v    = &fgeom->v[tp*dE];
370:         fegeom.J    = &fgeom->J[tp*dE*dE];
371:         fegeom.invJ = &fgeom->invJ[tp*dE*dE];
372:         fegeom.detJ = &fgeom->detJ[tp];
373:         fegeom.n    = &fgeom->n[tp*dE];

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

393: 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[],
394:                                              PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux,
395:                                              DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
396: {
397:   PetscFVCellGeom fvgeom;
398:   PetscInt        dim, dimEmbed;
399:   PetscErrorCode  ierr;

402:   DMGetDimension(dm, &dim);
403:   DMGetCoordinateDim(dm, &dimEmbed);
404:   if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
405:   switch (type) {
406:   case DM_BC_ESSENTIAL:
407:   case DM_BC_NATURAL:
408:     DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
409:   case DM_BC_ESSENTIAL_FIELD:
410:   case DM_BC_NATURAL_FIELD:
411:     DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
412:                                         (void (**)(PetscInt, PetscInt, PetscInt,
413:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
414:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
415:                                                    PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
416:   case DM_BC_ESSENTIAL_BD_FIELD:
417:     DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
418:                                           (void (**)(PetscInt, PetscInt, PetscInt,
419:                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
420:                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
421:                                                      PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
422:   default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
423:   }
424:   return(0);
425: }

427: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
428: {
429:   PetscReal      *points;
430:   PetscInt       f, numPoints;

434:   numPoints = 0;
435:   for (f = 0; f < Nf; ++f) {
436:     if (funcs[f]) {
437:       PetscQuadrature fAllPoints;
438:       PetscInt        fNumPoints;

440:       PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
441:       PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
442:       numPoints += fNumPoints;
443:     }
444:   }
445:   PetscMalloc1(dim*numPoints,&points);
446:   numPoints = 0;
447:   for (f = 0; f < Nf; ++f) {
448:     if (funcs[f]) {
449:       PetscQuadrature fAllPoints;
450:       PetscInt        qdim, fNumPoints, q;
451:       const PetscReal *fPoints;

453:       PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
454:       PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
455:       if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %D for dual basis does not match input dimension %D", qdim, dim);
456:       for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
457:       numPoints += fNumPoints;
458:     }
459:   }
460:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
461:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
462:   return(0);
463: }

465: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *ds)
466: {
467:   DM              plex;
468:   DMEnclosureType enc;
469:   DMLabel         depthLabel;
470:   PetscInt        dim, cdepth, ls = -1, i;
471:   PetscErrorCode  ierr;

474:   if (lStart) *lStart = -1;
475:   if (!label) return(0);
476:   DMGetEnclosureRelation(dm, odm, &enc);
477:   DMGetDimension(dm, &dim);
478:   DMConvert(dm, DMPLEX, &plex);
479:   DMPlexGetDepthLabel(plex, &depthLabel);
480:   cdepth = dim - height;
481:   for (i = 0; i < numIds; ++i) {
482:     IS              pointIS;
483:     const PetscInt *points;
484:     PetscInt        pdepth, point;

486:     DMLabelGetStratumIS(label, ids[i], &pointIS);
487:     if (!pointIS) continue; /* No points with that id on this process */
488:     ISGetIndices(pointIS, &points);
489:     DMGetEnclosurePoint(dm, odm, enc, points[0], &point);
490:     DMLabelGetValue(depthLabel, point, &pdepth);
491:     if (pdepth == cdepth) {
492:       ls = point;
493:       if (ds) {DMGetCellDS(dm, ls, ds);}
494:     }
495:     ISRestoreIndices(pointIS, &points);
496:     ISDestroy(&pointIS);
497:     if (ls >= 0) break;
498:   }
499:   if (lStart) *lStart = 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,
540:                                                   PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
541:                                                   DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
542:                                                   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, isHybrid = PETSC_FALSE, transform;
555:   DMField            coordField;
556:   DMLabel            depthLabel;
557:   PetscQuadrature    allPoints = NULL;
558:   PetscErrorCode     ierr;

561:   if (localU) {VecGetDM(localU, &dmIn);}
562:   else        {dmIn = dm;}
563:   DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, &localA);
564:   if (localA) {VecGetDM(localA, &dmAux);} else {dmAux = NULL;}
565:   DMConvert(dm, DMPLEX, &plex);
566:   DMConvert(dmIn, DMPLEX, &plexIn);
567:   DMGetEnclosureRelation(dmIn, dm, &encIn);
568:   DMGetEnclosureRelation(dmAux, dm, &encAux);
569:   DMGetDimension(dm, &dim);
570:   DMPlexGetVTKCellHeight(plex, &minHeight);
571:   DMGetBasisTransformDM_Internal(dm, &tdm);
572:   DMGetBasisTransformVec_Internal(dm, &tv);
573:   DMHasBasisTransform(dm, &transform);
574:   /* Auxiliary information can only be used with interpolation of field functions */
575:   if (dmAux) {
576:     DMConvert(dmAux, DMPLEX, &plexAux);
577:     if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
578:       if (!localA) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Missing localA vector");
579:     }
580:   }
581:   /* Determine height for iteration of all meshes */
582:   {
583:     DMPolytopeType ct, ctIn, ctAux;
584:     PetscInt       minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux;
585:     PetscInt       dim = -1, dimIn, dimAux;

587:     DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
588:     if (pEnd > pStart) {
589:       DMGetFirstLabelEntry_Private(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, NULL);
600:         DMPlexGetCellType(plexAux, pStartAux, &ctAux);
601:         dimAux = DMPolytopeTypeGetDim(ctAux);
602:       } else dimAux = dim;
603:     }
604:     if (dim < 0) {
605:       DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;

607:       /* Fall back to determination based on being a submesh */
608:       DMPlexGetSubpointMap(plex,   &spmap);
609:       DMPlexGetSubpointMap(plexIn, &spmapIn);
610:       if (plexAux) {DMPlexGetSubpointMap(plexAux, &spmapAux);}
611:       dim    = spmap    ? 1 : 0;
612:       dimIn  = spmapIn  ? 1 : 0;
613:       dimAux = spmapAux ? 1 : 0;
614:     }
615:     {
616:       PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), dimAux);

618:       if (PetscAbsInt(dimProj - dim) > 1 || PetscAbsInt(dimProj - dimIn) > 1 || PetscAbsInt(dimProj - dimAux) > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not currently support differences of more than 1 in dimension");
619:       if (dimProj < dim) minHeight = 1;
620:       htInc    =  dim    - dimProj;
621:       htIncIn  =  dimIn  - dimProj;
622:       htIncAux =  dimAux - dimProj;
623:     }
624:   }
625:   DMPlexGetDepth(plex, &depth);
626:   DMPlexGetDepthLabel(plex, &depthLabel);
627:   DMPlexGetMaxProjectionHeight(plex, &maxHeight);
628:   maxHeight = PetscMax(maxHeight, minHeight);
629:   if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
630:   DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
631:   if (!ds) {DMGetDS(dm, &ds);}
632:   DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
633:   if (!dsIn) {DMGetDS(dmIn, &dsIn);}
634:   PetscDSGetNumFields(ds, &Nf);
635:   PetscDSGetNumFields(dsIn, &NfIn);
636:   DMGetNumFields(dm, &NfTot);
637:   DMFindRegionNum(dm, ds, &regionNum);
638:   DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
639:   PetscDSGetHybrid(ds, &isHybrid);
640:   DMGetCoordinateDim(dm, &dimEmbed);
641:   DMGetLocalSection(dm, &section);
642:   if (dmAux) {
643:     DMGetDS(dmAux, &dsAux);
644:     PetscDSGetNumFields(dsAux, &NfAux);
645:   }
646:   PetscDSGetComponents(ds, &Nc);
647:   PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
648:   if (maxHeight > 0) {PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);}
649:   else               {cellsp = sp; cellspIn = spIn;}
650:   if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
651:   /* Get cell dual spaces */
652:   for (f = 0; f < Nf; ++f) {
653:     PetscDiscType disctype;

655:     PetscDSGetDiscType_Internal(ds, f, &disctype);
656:     if (disctype == PETSC_DISC_FE) {
657:       PetscFE fe;

659:       isFE[f] = PETSC_TRUE;
660:       hasFE   = PETSC_TRUE;
661:       PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
662:       PetscFEGetDualSpace(fe, &cellsp[f]);
663:     } else if (disctype == PETSC_DISC_FV) {
664:       PetscFV fv;

666:       isFE[f] = PETSC_FALSE;
667:       hasFV   = PETSC_TRUE;
668:       PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
669:       PetscFVGetDualSpace(fv, &cellsp[f]);
670:     } else {
671:       isFE[f]   = PETSC_FALSE;
672:       cellsp[f] = NULL;
673:     }
674:   }
675:   for (f = 0; f < NfIn; ++f) {
676:     PetscDiscType disctype;

678:     PetscDSGetDiscType_Internal(dsIn, f, &disctype);
679:     if (disctype == PETSC_DISC_FE) {
680:       PetscFE fe;

682:       PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fe);
683:       PetscFEGetDualSpace(fe, &cellspIn[f]);
684:     } else if (disctype == PETSC_DISC_FV) {
685:       PetscFV fv;

687:       PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fv);
688:       PetscFVGetDualSpace(fv, &cellspIn[f]);
689:     } else {
690:       cellspIn[f] = NULL;
691:     }
692:   }
693:   DMGetCoordinateField(dm,&coordField);
694:   for (f = 0; f < Nf; ++f) {
695:     if (!htInc) {sp[f] = cellsp[f];}
696:     else        {PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);}
697:   }
698:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
699:     PetscFE          fem, subfem;
700:     PetscDiscType    disctype;
701:     const PetscReal *points;
702:     PetscInt         numPoints;

704:     if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
705:     PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &allPoints);
706:     PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
707:     PetscMalloc2(NfIn, &T, NfAux, &TAux);
708:     for (f = 0; f < NfIn; ++f) {
709:       if (!htIncIn) {spIn[f] = cellspIn[f];}
710:       else          {PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);}

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

742:     if (h > minHeight) {
743:       for (f = 0; f < Nf; ++f) {PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);}
744:     }
745:     DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
746:     DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
747:     DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
748:     if (pEnd <= pStart) {
749:       ISDestroy(&heightIS);
750:       continue;
751:     }
752:     /* Compute totDim, the number of dofs in the closure of a point at this height */
753:     totDim = 0;
754:     for (f = 0; f < Nf; ++f) {
755:       if (!sp[f]) continue;
756:       PetscDualSpaceGetDimension(sp[f], &spDim);
757:       totDim += spDim;
758:       if (isHybrid && (f < Nf-1)) totDim += spDim;
759:     }
760:     p    = lStart < 0 ? pStart : lStart;
761:     DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
762:     if (numValues != totDim) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The output section point (%D) closure size %D != dual space dimension %D at height %D in [%D, %D]", p, numValues, totDim, h, minHeight, maxHeight);
763:     if (!totDim) {
764:       ISDestroy(&heightIS);
765:       continue;
766:     }
767:     if (htInc) {PetscDSGetHeightSubspace(ds, hEff, &dsEff);}
768:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
769:     if (localU) {
770:       PetscInt totDimIn, pIn, numValuesIn;

772:       totDimIn = 0;
773:       for (f = 0; f < NfIn; ++f) {
774:         if (!spIn[f]) continue;
775:         PetscDualSpaceGetDimension(spIn[f], &spDim);
776:         totDimIn += spDim;
777:         if (isHybrid && (f < Nf-1)) totDimIn += spDim;
778:       }
779:       DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
780:       DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
781:       if (numValuesIn != totDimIn) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The input section point (%D) closure size %D != dual space dimension %D at height %D", pIn, numValuesIn, totDimIn, htIncIn);
782:       if (htIncIn) {PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);}
783:     }
784:     if (htIncAux) {PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);}
785:     /* Loop over points at this height */
786:     DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
787:     DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
788:     {
789:       const PetscInt *fields;

791:       ISGetIndices(fieldIS, &fields);
792:       for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
793:       for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
794:       ISRestoreIndices(fieldIS, &fields);
795:     }
796:     if (label) {
797:       PetscInt i;

799:       for (i = 0; i < numIds; ++i) {
800:         IS              pointIS, isectIS;
801:         const PetscInt *points;
802:         PetscInt        n;
803:         PetscFEGeom  *fegeom = NULL, *chunkgeom = NULL;
804:         PetscQuadrature quad = NULL;

806:         DMLabelGetStratumIS(label, ids[i], &pointIS);
807:         if (!pointIS) continue; /* No points with that id on this process */
808:         ISIntersect(pointIS,heightIS,&isectIS);
809:         ISDestroy(&pointIS);
810:         if (!isectIS) continue;
811:         ISGetLocalSize(isectIS, &n);
812:         ISGetIndices(isectIS, &points);
813:         DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
814:         if (maxDegree <= 1) {
815:           DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
816:         }
817:         if (!quad) {
818:           if (!h && allPoints) {
819:             quad = allPoints;
820:             allPoints = NULL;
821:           } else {
822:             PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-htInc-1 : dim-htInc,funcs,&quad);
823:           }
824:         }
825:         DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
826:         for (p = 0; p < n; ++p) {
827:           const PetscInt  point = points[p];

829:           PetscArrayzero(values, numValues);
830:           PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
831:           DMPlexSetActivePoint(dm, point);
832:           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);
833:           if (ierr) {
834:             PetscErrorCode ierr2;
835:             ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
836:             ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
837:             
838:           }
839:           if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
840:           DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
841:         }
842:         PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
843:         PetscFEGeomDestroy(&fegeom);
844:         PetscQuadratureDestroy(&quad);
845:         ISRestoreIndices(isectIS, &points);
846:         ISDestroy(&isectIS);
847:       }
848:     } else {
849:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
850:       PetscQuadrature quad = NULL;
851:       IS              pointIS;

853:       ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
854:       DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
855:       if (maxDegree <= 1) {
856:         DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
857:       }
858:       if (!quad) {
859:         if (!h && allPoints) {
860:           quad = allPoints;
861:           allPoints = NULL;
862:         } else {
863:           PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &quad);
864:         }
865:       }
866:       DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
867:       for (p = pStart; p < pEnd; ++p) {
868:         PetscArrayzero(values, numValues);
869:         PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
870:         DMPlexSetActivePoint(dm, p);
871:         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);
872:         if (ierr) {
873:           PetscErrorCode ierr2;
874:           ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
875:           ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
876:           
877:         }
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: {

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: {

919:   DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
920:   return(0);
921: }

923: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
924:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
925:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
926:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
927:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
928:                                         InsertMode mode, Vec localX)
929: {

933:   DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
934:   return(0);
935: }

937: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
938:                                              void (**funcs)(PetscInt, PetscInt, PetscInt,
939:                                                             const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
940:                                                             const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
941:                                                             PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
942:                                              InsertMode mode, Vec localX)
943: {

947:   DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
948:   return(0);
949: }

951: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
952:                                                void (**funcs)(PetscInt, PetscInt, PetscInt,
953:                                                               const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
954:                                                               const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
955:                                                               PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
956:                                                InsertMode mode, Vec localX)
957: {

961:   DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
962:   return(0);
963: }