Actual source code: plexproject.c

petsc-master 2020-11-28
Report Typos and Errors
  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 Argument:
 11: . dm   - the DM

 13:   Output Argument:
 14: . point - The mesh point involved in the current projection

 16:   Level: developer

 18: .seealso: DMPlexSetActivePoint()
 19: @*/
 20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point) {
 22:   *point = ((DM_Plex *) dm->data)->activePoint;
 23:   return(0);
 24: }

 26: /*@
 27:   DMPlexSetActivePoint - Set the point on which projection is currently working

 29:   Not collective

 31:   Input Arguments:
 32: + dm   - the DM
 33: - point - The mesh point involved in the current projection

 35:   Level: developer

 37: .seealso: DMPlexGetActivePoint()
 38: @*/
 39: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point) {
 41:   ((DM_Plex *) dm->data)->activePoint = point;
 42:   return(0);
 43: }

 45: /*
 46:   DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point

 48:   Input Parameters:
 49: + dm     - The output DM
 50: . ds     - The output DS
 51: . dmIn   - The input DM
 52: . dsIn   - The input DS
 53: . time   - The time for this evaluation
 54: . fegeom - The FE geometry for this point
 55: . fvgeom - The FV geometry for this point
 56: . isFE   - Flag indicating whether each output field has an FE discretization
 57: . sp     - The output PetscDualSpace for each field
 58: . funcs  - The evaluation function for each field
 59: - ctxs   - The user context for each field

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

 64:   Level: developer

 66: .seealso: DMProjectPoint_Field_Private()
 67: */
 68: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[],
 69:                                                   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
 70:                                                   PetscScalar values[])
 71: {
 72:   PetscInt       coordDim, Nf, *Nc, f, spDim, d, v, tp;
 73:   PetscBool      isAffine, isHybrid, transform;

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

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

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

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

110:             if (dim != fegeom->dim) {
111:               if (isHybrid) {
112:                 /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
113:                 for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
114:                 refpoint = injpoint;
115:               } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D dual basis spatial dimension", fegeom->dim, dim);
116:             }
117:             CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
118:             v0 = x;
119:           } else {
120:             v0 = &fegeom->v[tp*coordDim];
121:           }
122:           if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
123:           (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
124:         }
125:         /* Transform point evaluations pointEval[q,c] */
126:         PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
127:         PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
128:         DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
129:         DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
130:         v += spDim;
131:         if (isHybrid && (f < Nf-1)) {
132:           for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
133:         }
134:       } else {
135:         for (d = 0; d < spDim; ++d, ++v) {
136:           PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
137:         }
138:       }
139:     } else {
140:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
141:       if (isHybrid && (f < Nf-1)) {
142:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
143:       }
144:     }
145:   }
146:   return(0);
147: }

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

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

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

173:   Note: Not supported for FV

175:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

451:       PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
452:       PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
453:       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);
454:       for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
455:       numPoints += fNumPoints;
456:     }
457:   }
458:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
459:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
460:   return(0);
461: }

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

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

484:     DMLabelGetStratumIS(label, ids[i], &pointIS);
485:     if (!pointIS) continue; /* No points with that id on this process */
486:     ISGetIndices(pointIS, &points);
487:     DMGetEnclosurePoint(dm, odm, enc, points[0], &point);
488:     DMLabelGetValue(depthLabel, point, &pdepth);
489:     if (pdepth == cdepth) {
490:       ls = point;
491:       if (ds) {DMGetCellDS(dm, ls, ds);}
492:     }
493:     ISRestoreIndices(pointIS, &points);
494:     ISDestroy(&pointIS);
495:     if (ls >= 0) break;
496:   }
497:   if (lStart) *lStart = ls;
498:   DMDestroy(&plex);
499:   return(0);
500: }

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

505:   There are several different scenarios:

507:   1) Volumetric mesh with volumetric auxiliary data

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

511:   2) Boundary mesh with boundary auxiliary data

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

515:   3) Volumetric mesh with boundary auxiliary data

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

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

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

523:   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.
524:     - 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
525:       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
526:       dual spaces for the boundary from our input spaces.
527:     - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.

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

531:   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.
532: */
533: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
534:                                                   PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
535:                                                   DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
536:                                                   InsertMode mode, Vec localX)
537: {
538:   DM                 plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
539:   DMEnclosureType    encIn, encAux;
540:   PetscDS            ds = NULL, dsIn = NULL, dsAux = NULL;
541:   Vec                localA = NULL, tv;
542:   IS                 fieldIS;
543:   PetscSection       section;
544:   PetscDualSpace    *sp, *cellsp;
545:   PetscTabulation *T = NULL, *TAux = NULL;
546:   PetscInt          *Nc;
547:   PetscInt           dim, dimEmbed, depth, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
548:   PetscBool         *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, isHybrid = PETSC_FALSE, transform;
549:   DMField            coordField;
550:   DMLabel            depthLabel;
551:   PetscQuadrature    allPoints = NULL;
552:   PetscErrorCode     ierr;

555:   if (localU) {VecGetDM(localU, &dmIn);}
556:   else        {dmIn = dm;}
557:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
558:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
559:   DMConvert(dm, DMPLEX, &plex);
560:   DMConvert(dmIn, DMPLEX, &plexIn);
561:   DMGetEnclosureRelation(dmIn, dm, &encIn);
562:   DMGetEnclosureRelation(dmAux, dm, &encAux);
563:   DMGetDimension(dm, &dim);
564:   DMPlexGetVTKCellHeight(plex, &minHeight);
565:   DMGetBasisTransformDM_Internal(dm, &tdm);
566:   DMGetBasisTransformVec_Internal(dm, &tv);
567:   DMHasBasisTransform(dm, &transform);
568:   /* Auxiliary information can only be used with interpolation of field functions */
569:   if (dmAux) {
570:     DMConvert(dmAux, DMPLEX, &plexAux);
571:     if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
572:       if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
573:       if (!minHeight) {
574:         DMLabel spmap;

576:         /* If dmAux is a surface, then force the projection to take place over a surface */
577:         DMPlexGetSubpointMap(plexAux, &spmap);
578:         if (spmap) {
579:           DMPlexGetVTKCellHeight(plexAux, &minHeight);
580:           auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
581:         }
582:       }
583:     }
584:   }
585:   DMPlexGetDepth(plex, &depth);
586:   DMPlexGetDepthLabel(plex, &depthLabel);
587:   DMPlexGetMaxProjectionHeight(plex, &maxHeight);
588:   maxHeight = PetscMax(maxHeight, minHeight);
589:   if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
590:   DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
591:   if (!ds) {DMGetDS(dm, &ds);}
592:   DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
593:   if (!dsIn) {DMGetDS(dmIn, &dsIn);}
594:   PetscDSGetNumFields(ds, &Nf);
595:   PetscDSGetNumFields(dsIn, &NfIn);
596:   DMGetNumFields(dm, &NfTot);
597:   DMFindRegionNum(dm, ds, &regionNum);
598:   DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
599:   PetscDSGetHybrid(ds, &isHybrid);
600:   DMGetCoordinateDim(dm, &dimEmbed);
601:   DMGetLocalSection(dm, &section);
602:   if (dmAux) {
603:     DMGetDS(dmAux, &dsAux);
604:     PetscDSGetNumFields(dsAux, &NfAux);
605:   }
606:   PetscDSGetComponents(ds, &Nc);
607:   PetscMalloc2(Nf, &isFE, Nf, &sp);
608:   if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
609:   else               {cellsp = sp;}
610:   if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
611:   /* Get cell dual spaces */
612:   for (f = 0; f < Nf; ++f) {
613:     PetscDiscType disctype;

615:     PetscDSGetDiscType_Internal(ds, f, &disctype);
616:     if (disctype == PETSC_DISC_FE) {
617:       PetscFE fe;

619:       isFE[f] = PETSC_TRUE;
620:       hasFE   = PETSC_TRUE;
621:       PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
622:       PetscFEGetDualSpace(fe, &cellsp[f]);
623:     } else if (disctype == PETSC_DISC_FV) {
624:       PetscFV fv;

626:       isFE[f] = PETSC_FALSE;
627:       hasFV   = PETSC_TRUE;
628:       PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
629:       PetscFVGetDualSpace(fv, &cellsp[f]);
630:     } else {
631:       isFE[f]   = PETSC_FALSE;
632:       cellsp[f] = NULL;
633:     }
634:   }
635:   DMGetCoordinateField(dm,&coordField);
636:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
637:     PetscInt         effectiveHeight = auxBd ? minHeight : 0;
638:     PetscFE          fem, subfem;
639:     PetscDiscType    disctype;
640:     const PetscReal *points;
641:     PetscInt         numPoints;

643:     if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
644:     for (f = 0; f < Nf; ++f) {
645:       if (!effectiveHeight) {sp[f] = cellsp[f];}
646:       else                  {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
647:     }
648:     PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
649:     PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
650:     PetscMalloc2(NfIn, &T, NfAux, &TAux);
651:     for (f = 0; f < NfIn; ++f) {
652:       PetscDSGetDiscType_Internal(dsIn, f, &disctype);
653:       if (disctype != PETSC_DISC_FE) continue;
654:       PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
655:       if (!effectiveHeight) {subfem = fem;}
656:       else                  {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
657:       PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
658:     }
659:     for (f = 0; f < NfAux; ++f) {
660:       PetscDSGetDiscType_Internal(dsAux, f, &disctype);
661:       if (disctype != PETSC_DISC_FE) continue;
662:       PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
663:       if (!effectiveHeight || auxBd) {subfem = fem;}
664:       else                           {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
665:       PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
666:     }
667:   }
668:   /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
669:   for (h = minHeight; h <= maxHeight; h++) {
670:     PetscInt     effectiveHeight = h - (auxBd ? 0 : minHeight);
671:     PetscDS      dsEff         = ds;
672:     PetscScalar *values;
673:     PetscBool   *fieldActive;
674:     PetscInt     maxDegree;
675:     PetscInt     pStart, pEnd, p, lStart, spDim, totDim, numValues;
676:     IS           heightIS;

678:     /* Note we assume that dm and dmIn share the same topology */
679:     DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
680:     DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
681:     DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
682:     if (pEnd <= pStart) {
683:       ISDestroy(&heightIS);
684:       continue;
685:     }
686:     /* Compute totDim, the number of dofs in the closure of a point at this height */
687:     totDim = 0;
688:     for (f = 0; f < Nf; ++f) {
689:       if (!effectiveHeight) {
690:         sp[f] = cellsp[f];
691:       } else {
692:         PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
693:       }
694:       if (!sp[f]) continue;
695:       PetscDualSpaceGetDimension(sp[f], &spDim);
696:       totDim += spDim;
697:       if (isHybrid && (f < Nf-1)) totDim += spDim;
698:     }
699:     DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
700:     if (numValues != totDim) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point (%D) closure size %D != dual space dimension %D", lStart < 0 ? pStart : lStart, numValues, totDim);
701:     if (!totDim) {
702:       ISDestroy(&heightIS);
703:       continue;
704:     }
705:     if (effectiveHeight) {PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);}
706:     /* Loop over points at this height */
707:     DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
708:     DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
709:     {
710:       const PetscInt *fields;

712:       ISGetIndices(fieldIS, &fields);
713:       for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
714:       for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
715:       ISRestoreIndices(fieldIS, &fields);
716:     }
717:     if (label) {
718:       PetscInt i;

720:       for (i = 0; i < numIds; ++i) {
721:         IS              pointIS, isectIS;
722:         const PetscInt *points;
723:         PetscInt        n;
724:         PetscFEGeom  *fegeom = NULL, *chunkgeom = NULL;
725:         PetscQuadrature quad = NULL;

727:         DMLabelGetStratumIS(label, ids[i], &pointIS);
728:         if (!pointIS) continue; /* No points with that id on this process */
729:         ISIntersect(pointIS,heightIS,&isectIS);
730:         ISDestroy(&pointIS);
731:         if (!isectIS) continue;
732:         ISGetLocalSize(isectIS, &n);
733:         ISGetIndices(isectIS, &points);
734:         DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
735:         if (maxDegree <= 1) {
736:           DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
737:         }
738:         if (!quad) {
739:           if (!h && allPoints) {
740:             quad = allPoints;
741:             allPoints = NULL;
742:           } else {
743:             PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-effectiveHeight-1 : dim-effectiveHeight,funcs,&quad);
744:           }
745:         }
746:         DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
747:         for (p = 0; p < n; ++p) {
748:           const PetscInt  point = points[p];

750:           PetscArrayzero(values, numValues);
751:           PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
752:           DMPlexSetActivePoint(dm, point);
753:           DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values);
754:           if (ierr) {
755:             PetscErrorCode ierr2;
756:             ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
757:             ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
758:             
759:           }
760:           if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
761:           DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
762:         }
763:         PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
764:         PetscFEGeomDestroy(&fegeom);
765:         PetscQuadratureDestroy(&quad);
766:         ISRestoreIndices(isectIS, &points);
767:         ISDestroy(&isectIS);
768:       }
769:     } else {
770:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
771:       PetscQuadrature quad = NULL;
772:       IS              pointIS;

774:       ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
775:       DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
776:       if (maxDegree <= 1) {
777:         DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
778:       }
779:       if (!quad) {
780:         if (!h && allPoints) {
781:           quad = allPoints;
782:           allPoints = NULL;
783:         } else {
784:           PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&quad);
785:         }
786:       }
787:       DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
788:       for (p = pStart; p < pEnd; ++p) {
789:         PetscArrayzero(values, numValues);
790:         PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
791:         DMPlexSetActivePoint(dm, p);
792:         DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, T, TAux, type, funcs, ctxs, fieldActive, values);
793:         if (ierr) {
794:           PetscErrorCode ierr2;
795:           ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
796:           ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
797:           
798:         }
799:         if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
800:         DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
801:       }
802:       PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
803:       PetscFEGeomDestroy(&fegeom);
804:       PetscQuadratureDestroy(&quad);
805:       ISDestroy(&pointIS);
806:     }
807:     ISDestroy(&heightIS);
808:     DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
809:     DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
810:   }
811:   /* Cleanup */
812:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
813:     PetscInt effectiveHeight = auxBd ? minHeight : 0;
814:     PetscFE  fem, subfem;

816:     for (f = 0; f < NfIn; ++f) {
817:       PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
818:       if (!effectiveHeight) {subfem = fem;}
819:       else                  {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
820:       PetscTabulationDestroy(&T[f]);
821:     }
822:     for (f = 0; f < NfAux; ++f) {
823:       PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
824:       if (!effectiveHeight || auxBd) {subfem = fem;}
825:       else                           {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
826:       PetscTabulationDestroy(&TAux[f]);
827:     }
828:     PetscFree2(T, TAux);
829:   }
830:   PetscQuadratureDestroy(&allPoints);
831:   PetscFree2(isFE, sp);
832:   if (maxHeight > 0) {PetscFree(cellsp);}
833:   DMDestroy(&plex);
834:   DMDestroy(&plexIn);
835:   if (dmAux) {DMDestroy(&plexAux);}
836:   return(0);
837: }

839: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
840: {

844:   DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
845:   return(0);
846: }

848: 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)
849: {

853:   DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
854:   return(0);
855: }

857: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
858:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
859:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
860:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
861:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
862:                                         InsertMode mode, Vec localX)
863: {

867:   DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
868:   return(0);
869: }

871: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
872:                                              void (**funcs)(PetscInt, PetscInt, PetscInt,
873:                                                             const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
874:                                                             const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
875:                                                             PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
876:                                              InsertMode mode, Vec localX)
877: {

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

885: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
886:                                                void (**funcs)(PetscInt, PetscInt, PetscInt,
887:                                                               const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
888:                                                               const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
889:                                                               PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
890:                                                InsertMode mode, Vec localX)
891: {

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