Actual source code: plexproject.c

  1: #include <petsc/private/dmpleximpl.h>

  3: #include <petsc/private/petscfeimpl.h>

  5: /*@
  6:   DMPlexGetActivePoint - Get the point on which projection is currently working

  8:   Not Collective

 10:   Input Parameter:
 11: . dm   - the `DM`

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

 16:   Level: developer

 18: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`
 19: @*/
 20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
 21: {
 22:   PetscFunctionBeginHot;
 23:   *point = ((DM_Plex *)dm->data)->activePoint;
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

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

 30:   Not Collective

 32:   Input Parameters:
 33: + dm   - the `DM`
 34: - point - The mesh point involved in the current projection

 36:   Level: developer

 38: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetActivePoint()`
 39: @*/
 40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
 41: {
 42:   PetscFunctionBeginHot;
 43:   ((DM_Plex *)dm->data)->activePoint = point;
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: /*
 48:   DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point

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

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

 66:   Level: developer

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

 75:   PetscFunctionBeginHot;
 76:   PetscCall(DMGetCoordinateDim(dmIn, &coordDim));
 77:   PetscCall(DMHasBasisTransform(dmIn, &transform));
 78:   PetscCall(PetscDSGetNumFields(ds, &Nf));
 79:   PetscCall(PetscDSGetComponents(ds, &Nc));
 80:   PetscCall(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:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
 89:     PetscCall(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:         PetscCall(PetscDualSpaceGetDM(sp[f], &rdm));
100:         PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
101:         PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
102:         PetscCall(DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
103:         PetscCall(DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x));
104:         PetscCall(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:             PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx));
126:             v0 = x;
127:           }
128:           PetscCall((*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx));
129:         }
130:         /* Transform point evaluations pointEval[q,c] */
131:         PetscCall(PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval));
132:         PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
133:         PetscCall(DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x));
134:         PetscCall(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) PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

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

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

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

176:   Level: developer

178:   Note:
179:   Not supported for FV

181: .seealso: `DMProjectPoint_Field_Private()`
182: */
183: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
184: {
185:   PetscSection       section, sectionAux = NULL;
186:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
187:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
188:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
189:   const PetscScalar *constants;
190:   PetscReal         *x;
191:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
192:   PetscFEGeom        fegeom;
193:   const PetscInt     dE = cgeom->dimEmbed;
194:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
195:   PetscBool          isAffine, isCohesive, transform;

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

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

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

282: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
283: {
284:   PetscSection       section, sectionAux = NULL;
285:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
286:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
287:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
288:   const PetscScalar *constants;
289:   PetscReal         *x;
290:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
291:   PetscFEGeom        fegeom, cgeom;
292:   const PetscInt     dE = fgeom->dimEmbed;
293:   PetscInt           numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
294:   PetscBool          isAffine;

296:   PetscFunctionBeginHot;
297:   PetscCheck(dm == dmIn, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM");
298:   PetscCall(PetscDSGetNumFields(ds, &Nf));
299:   PetscCall(PetscDSGetComponents(ds, &Nc));
300:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
301:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
302:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x));
303:   PetscCall(PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL));
304:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
305:   PetscCall(DMGetLocalSection(dm, &section));
306:   PetscCall(DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients));
307:   if (dmAux) {
308:     PetscInt subp;

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

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

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

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

383: static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscFEGeom *fegeom, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
384: {
385:   PetscFVCellGeom fvgeom;
386:   PetscInt        dim, dimEmbed;

388:   PetscFunctionBeginHot;
389:   PetscCall(DMGetDimension(dm, &dim));
390:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
391:   if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
392:   switch (type) {
393:   case DM_BC_ESSENTIAL:
394:   case DM_BC_NATURAL:
395:     PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode(**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
396:     break;
397:   case DM_BC_ESSENTIAL_FIELD:
398:   case DM_BC_NATURAL_FIELD:
399:     PetscCall(DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
400:     break;
401:   case DM_BC_ESSENTIAL_BD_FIELD:
402:     PetscCall(DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
403:     break;
404:   default:
405:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
406:   }
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

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

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

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

439:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
440:       PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
441:       PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %" PetscInt_FMT " for dual basis does not match input dimension %" PetscInt_FMT, qdim, dim);
442:       for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
443:       numPoints += fNumPoints;
444:     }
445:   }
446:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
447:   PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

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

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

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

466:   Level: developer

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

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

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

512:   There are several different scenarios:

514:   1) Volumetric mesh with volumetric auxiliary data

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

518:   2) Boundary mesh with boundary auxiliary data

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

522:   3) Volumetric mesh with boundary auxiliary data

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

526:   4) Volumetric input mesh with boundary output mesh

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

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

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

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

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

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

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

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

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

629:       PetscCheck(PetscAbsInt(dimProj - dim) <= 1 && PetscAbsInt(dimProj - dimIn) <= 1 && PetscAbsInt(dimProj - dimAuxEff) <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not currently support differences of more than 1 in dimension");
630:       if (dimProj < dim) minHeight = 1;
631:       htInc    = dim - dimProj;
632:       htIncIn  = dimIn - dimProj;
633:       htIncAux = dimAuxEff - dimProj;
634:     }
635:   }
636:   PetscCall(DMPlexGetDepth(plex, &depth));
637:   PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
638:   PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
639:   maxHeight = PetscMax(maxHeight, minHeight);
640:   PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
641:   PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, 0, NULL, &ds));
642:   if (!ds) PetscCall(DMGetDS(dm, &ds));
643:   PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn));
644:   if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
645:   PetscCall(PetscDSGetNumFields(ds, &Nf));
646:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
647:   PetscCall(DMGetNumFields(dm, &NfTot));
648:   PetscCall(DMFindRegionNum(dm, ds, &regionNum));
649:   PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
650:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
651:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
652:   PetscCall(DMGetLocalSection(dm, &section));
653:   if (dmAux) {
654:     PetscCall(DMGetDS(dmAux, &dsAux));
655:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
656:   }
657:   PetscCall(PetscDSGetComponents(ds, &Nc));
658:   PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
659:   if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
660:   else {
661:     cellsp   = sp;
662:     cellspIn = spIn;
663:   }
664:   /* Get cell dual spaces */
665:   for (f = 0; f < Nf; ++f) {
666:     PetscDiscType disctype;

668:     PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
669:     if (disctype == PETSC_DISC_FE) {
670:       PetscFE fe;

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

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

691:     PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
692:     if (disctype == PETSC_DISC_FE) {
693:       PetscFE fe;

695:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
696:       PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
697:     } else if (disctype == PETSC_DISC_FV) {
698:       PetscFV fv;

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

717:     PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
718:     PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
719:     PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
720:     PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
721:     for (f = 0; f < NfIn; ++f) {
722:       if (!htIncIn) {
723:         spIn[f] = cellspIn[f];
724:       } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));

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

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

773:       if (!sp[f]) continue;
774:       PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
775:       PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
776:       totDim += spDim;
777:       if (isCohesive && !cohesive) totDim += spDim;
778:     }
779:     p = lStart < 0 ? pStart : lStart;
780:     PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
781:     PetscCheck(numValues == totDim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The output section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT " in [%" PetscInt_FMT ", %" PetscInt_FMT "]", p, numValues, totDim, h, minHeight, maxHeight);
782:     if (!totDim) {
783:       PetscCall(ISDestroy(&heightIS));
784:       continue;
785:     }
786:     if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
787:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
788:     if (localU) {
789:       PetscInt totDimIn, pIn, numValuesIn;

791:       totDimIn = 0;
792:       for (f = 0; f < NfIn; ++f) {
793:         PetscBool cohesive;

795:         if (!spIn[f]) continue;
796:         PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
797:         PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
798:         totDimIn += spDim;
799:         if (isCohesive && !cohesive) totDimIn += spDim;
800:       }
801:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
802:       PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
803:       PetscCheck(numValuesIn == totDimIn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The input section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT, pIn, numValuesIn, totDimIn, htIncIn);
804:       if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
805:     }
806:     if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
807:     /* Loop over points at this height */
808:     PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
809:     PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
810:     {
811:       const PetscInt *fields;

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

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

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

849:           PetscCall(PetscArrayzero(values, numValues));
850:           PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
851:           PetscCall(DMPlexSetActivePoint(dm, point));
852:           PetscCall(DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsEffIn, plexAux, encAux, dsEffAux, chunkgeom, htInc, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values));
853:           if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
854:           PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
855:         }
856:         PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
857:         PetscCall(PetscFEGeomDestroy(&fegeom));
858:         PetscCall(PetscQuadratureDestroy(&quad));
859:         PetscCall(ISRestoreIndices(isectIS, &points));
860:         PetscCall(ISDestroy(&isectIS));
861:       }
862:     } else {
863:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
864:       PetscQuadrature quad = NULL;
865:       IS              pointIS;

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

911: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
912: {
913:   PetscFunctionBegin;
914:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
915:   PetscFunctionReturn(PETSC_SUCCESS);
916: }

918: 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)
919: {
920:   PetscFunctionBegin;
921:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
922:   PetscFunctionReturn(PETSC_SUCCESS);
923: }

925: 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)
926: {
927:   PetscFunctionBegin;
928:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
929:   PetscFunctionReturn(PETSC_SUCCESS);
930: }

932: 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)
933: {
934:   PetscFunctionBegin;
935:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
936:   PetscFunctionReturn(PETSC_SUCCESS);
937: }

939: 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)
940: {
941:   PetscFunctionBegin;
942:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX));
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }