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, face[2];
192:   PetscFEGeom        fegeom;
193:   const PetscInt     dE = cgeom->dimEmbed, *cone, *ornt;
194:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
195:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
196:   DMPolytopeType     qct;

198:   PetscFunctionBeginHot;
199:   PetscCall(PetscDSGetNumFields(ds, &Nf));
200:   PetscCall(PetscDSGetComponents(ds, &Nc));
201:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
202:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
203:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
204:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
205:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
206:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
207:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
208:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
209:   PetscCall(DMHasBasisTransform(dmIn, &transform));
210:   PetscCall(DMGetLocalSection(dmIn, &section));
211:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
212:   // Get cohesive cell hanging off face
213:   if (isCohesiveIn) {
214:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
215:     if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
216:       DMPolytopeType  ct;
217:       const PetscInt *support;
218:       PetscInt        Ns, s;

220:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
221:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
222:       for (s = 0; s < Ns; ++s) {
223:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
224:         if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
225:           inp = support[s];
226:           break;
227:         }
228:       }
229:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
230:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
231:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
232:       face[0] = 0;
233:       face[1] = 0;
234:     }
235:   }
236:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
237:   if (dmAux) {
238:     PetscInt subp;

240:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
241:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
242:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
243:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
244:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
245:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
246:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
247:   }
248:   /* Get values for closure */
249:   isAffine        = cgeom->isAffine;
250:   fegeom.dim      = cgeom->dim;
251:   fegeom.dimEmbed = cgeom->dimEmbed;
252:   if (isAffine) {
253:     fegeom.v    = x;
254:     fegeom.xi   = cgeom->xi;
255:     fegeom.J    = cgeom->J;
256:     fegeom.invJ = cgeom->invJ;
257:     fegeom.detJ = cgeom->detJ;
258:   }
259:   for (f = 0, v = 0; f < Nf; ++f) {
260:     PetscQuadrature  allPoints;
261:     PetscInt         q, dim, numPoints;
262:     const PetscReal *points;
263:     PetscScalar     *pointEval;
264:     PetscBool        cohesive;
265:     DM               dm;

267:     if (!sp[f]) continue;
268:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
269:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
270:     if (!funcs[f]) {
271:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
272:       if (isCohesive && !cohesive) {
273:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
274:       }
275:       continue;
276:     }
277:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
278:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
279:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
280:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
281:     for (q = 0; q < numPoints; ++q, ++tp) {
282:       PetscInt qpt[2];

284:       if (isCohesiveIn) {
285:         PetscCall(PetscDSPermuteQuadPoint(dsIn, ornt[0], f, q, &qpt[0]));
286:         PetscCall(PetscDSPermuteQuadPoint(dsIn, DMPolytopeTypeComposeOrientationInv(qct, ornt[1], 0), f, q, &qpt[1]));
287:       }
288:       if (isAffine) {
289:         CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
290:       } else {
291:         fegeom.v    = &cgeom->v[tp * dE];
292:         fegeom.J    = &cgeom->J[tp * dE * dE];
293:         fegeom.invJ = &cgeom->invJ[tp * dE * dE];
294:         fegeom.detJ = &cgeom->detJ[tp];
295:       }
296:       if (coefficients) {
297:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
298:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
299:       }
300:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
301:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
302:       (*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]);
303:     }
304:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
305:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
306:     v += spDim;
307:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
308:     if (isCohesive && !cohesive) {
309:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
310:     }
311:   }
312:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
313:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
314:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
315:   PetscFunctionReturn(PETSC_SUCCESS);
316: }

318: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, 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[])
319: {
320:   PetscSection       section, sectionAux = NULL;
321:   PetscScalar       *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
322:   PetscScalar       *coefficients = NULL, *coefficientsAux = NULL;
323:   PetscScalar       *coefficients_t = NULL, *coefficientsAux_t = NULL;
324:   const PetscScalar *constants;
325:   PetscReal         *x;
326:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
327:   PetscFEGeom        fegeom, cgeom;
328:   const PetscInt     dE = fgeom->dimEmbed, *cone, *ornt;
329:   PetscInt           numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
330:   PetscBool          isAffine, isCohesive, isCohesiveIn, transform;
331:   DMPolytopeType     qct;

333:   PetscFunctionBeginHot;
334:   PetscCall(PetscDSGetNumFields(ds, &Nf));
335:   PetscCall(PetscDSGetComponents(ds, &Nc));
336:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
337:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
338:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
339:   PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
340:   PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
341:   PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
342:   PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
343:   PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
344:   PetscCall(DMHasBasisTransform(dmIn, &transform));
345:   PetscCall(DMGetLocalSection(dmIn, &section));
346:   PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
347:   // Get cohesive cell hanging off face
348:   if (isCohesiveIn) {
349:     PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
350:     if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
351:       DMPolytopeType  ct;
352:       const PetscInt *support;
353:       PetscInt        Ns, s;

355:       PetscCall(DMPlexGetSupport(dmIn, inp, &support));
356:       PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
357:       for (s = 0; s < Ns; ++s) {
358:         PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
359:         if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
360:           inp = support[s];
361:           break;
362:         }
363:       }
364:       PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
365:       PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
366:       PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
367:       face[0] = 0;
368:       face[1] = 0;
369:     }
370:   }
371:   if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
372:   if (dmAux) {
373:     PetscInt subp;

375:     PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
376:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
377:     PetscCall(DMGetLocalSection(dmAux, &sectionAux));
378:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
379:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
380:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
381:     PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
382:   }
383:   /* Get values for closure */
384:   isAffine       = fgeom->isAffine;
385:   fegeom.n       = NULL;
386:   fegeom.J       = NULL;
387:   fegeom.v       = NULL;
388:   fegeom.xi      = NULL;
389:   cgeom.dim      = fgeom->dim;
390:   cgeom.dimEmbed = fgeom->dimEmbed;
391:   if (isAffine) {
392:     fegeom.v    = x;
393:     fegeom.xi   = fgeom->xi;
394:     fegeom.J    = fgeom->J;
395:     fegeom.invJ = fgeom->invJ;
396:     fegeom.detJ = fgeom->detJ;
397:     fegeom.n    = fgeom->n;

399:     cgeom.J    = fgeom->suppJ[0];
400:     cgeom.invJ = fgeom->suppInvJ[0];
401:     cgeom.detJ = fgeom->suppDetJ[0];
402:   }
403:   for (f = 0, v = 0; f < Nf; ++f) {
404:     PetscQuadrature  allPoints;
405:     PetscInt         q, dim, numPoints;
406:     const PetscReal *points;
407:     PetscScalar     *pointEval;
408:     PetscBool        cohesive;
409:     DM               dm;

411:     if (!sp[f]) continue;
412:     PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
413:     PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
414:     if (!funcs[f]) {
415:       for (d = 0; d < spDim; d++, v++) values[v] = 0.;
416:       if (isCohesive && !cohesive) {
417:         for (d = 0; d < spDim; d++, v++) values[v] = 0.;
418:       }
419:       continue;
420:     }
421:     PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
422:     PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
423:     PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
424:     PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
425:     for (q = 0; q < numPoints; ++q, ++tp) {
426:       PetscInt qpt[2];

428:       if (isCohesiveIn) {
429:         // These points are not integration quadratures, but dual space quadratures
430:         // If they had multiple points we should match them from both sides, similar to hybrid residual eval
431:         qpt[0] = qpt[1] = q;
432:       }
433:       if (isAffine) {
434:         CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
435:       } else {
436:         fegeom.v    = &fgeom->v[tp * dE];
437:         fegeom.J    = &fgeom->J[tp * dE * dE];
438:         fegeom.invJ = &fgeom->invJ[tp * dE * dE];
439:         fegeom.detJ = &fgeom->detJ[tp];
440:         fegeom.n    = &fgeom->n[tp * dE];

442:         cgeom.J    = &fgeom->suppJ[0][tp * dE * dE];
443:         cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
444:         cgeom.detJ = &fgeom->suppDetJ[0][tp];
445:       }
446:       /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
447:       if (coefficients) {
448:         if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
449:         else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
450:       }
451:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
452:       if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
453:       (*funcs[f])(dE, NfIn, 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]);
454:     }
455:     PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
456:     PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
457:     v += spDim;
458:     /* TODO: For now, set both sides equal, but this should use info from other support cell */
459:     if (isCohesive && !cohesive) {
460:       for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
461:     }
462:   }
463:   if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
464:   if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
465:   if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

469: 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[])
470: {
471:   PetscFVCellGeom fvgeom;
472:   PetscInt        dim, dimEmbed;

474:   PetscFunctionBeginHot;
475:   PetscCall(DMGetDimension(dm, &dim));
476:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
477:   if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
478:   switch (type) {
479:   case DM_BC_ESSENTIAL:
480:   case DM_BC_NATURAL:
481:     PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode(**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
482:     break;
483:   case DM_BC_ESSENTIAL_FIELD:
484:   case DM_BC_NATURAL_FIELD:
485:     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));
486:     break;
487:   case DM_BC_ESSENTIAL_BD_FIELD:
488:     PetscCall(DMProjectPoint_BdField_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[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
489:     break;
490:   default:
491:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
492:   }
493:   PetscFunctionReturn(PETSC_SUCCESS);
494: }

496: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
497: {
498:   PetscReal *points;
499:   PetscInt   f, numPoints;

501:   PetscFunctionBegin;
502:   if (!dim) {
503:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
504:     PetscFunctionReturn(PETSC_SUCCESS);
505:   }
506:   numPoints = 0;
507:   for (f = 0; f < Nf; ++f) {
508:     if (funcs[f]) {
509:       PetscQuadrature fAllPoints;
510:       PetscInt        fNumPoints;

512:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
513:       PetscCall(PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL));
514:       numPoints += fNumPoints;
515:     }
516:   }
517:   PetscCall(PetscMalloc1(dim * numPoints, &points));
518:   numPoints = 0;
519:   for (f = 0; f < Nf; ++f) {
520:     if (funcs[f]) {
521:       PetscQuadrature  fAllPoints;
522:       PetscInt         qdim, fNumPoints, q;
523:       const PetscReal *fPoints;

525:       PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
526:       PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
527:       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);
528:       for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
529:       numPoints += fNumPoints;
530:     }
531:   }
532:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
533:   PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

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

540:   Input Parameters:
541: + dm     - the `DM`
542: . odm    - the enclosing `DM`
543: . label  - label for `DM` domain, or `NULL` for whole domain
544: . numIds - the number of `ids`
545: . ids    - An array of the label ids in sequence for the domain
546: - height - Height of target cells in `DMPLEX` topology

548:   Output Parameters:
549: + point - the first labeled point
550: - ds    - the `PetscDS` corresponding to the first labeled point

552:   Level: developer

554: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
555: @*/
556: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
557: {
558:   DM              plex;
559:   DMEnclosureType enc;
560:   PetscInt        ls = -1;

562:   PetscFunctionBegin;
563:   if (point) *point = -1;
564:   if (!label) PetscFunctionReturn(PETSC_SUCCESS);
565:   PetscCall(DMGetEnclosureRelation(dm, odm, &enc));
566:   PetscCall(DMConvert(dm, DMPLEX, &plex));
567:   for (PetscInt i = 0; i < numIds; ++i) {
568:     IS       labelIS;
569:     PetscInt num_points, pStart, pEnd;
570:     PetscCall(DMLabelGetStratumIS(label, ids[i], &labelIS));
571:     if (!labelIS) continue; /* No points with that id on this process */
572:     PetscCall(DMPlexGetHeightStratum(plex, height, &pStart, &pEnd));
573:     PetscCall(ISGetSize(labelIS, &num_points));
574:     if (num_points) {
575:       const PetscInt *points;
576:       PetscCall(ISGetIndices(labelIS, &points));
577:       for (PetscInt i = 0; i < num_points; i++) {
578:         PetscInt point;
579:         PetscCall(DMGetEnclosurePoint(dm, odm, enc, points[i], &point));
580:         if (pStart <= point && point < pEnd) {
581:           ls = point;
582:           if (ds) {
583:             // If this is a face of a cohesive cell, then prefer that DS
584:             if (height == 1) {
585:               const PetscInt *supp;
586:               PetscInt        suppSize;
587:               DMPolytopeType  ct;

589:               PetscCall(DMPlexGetSupport(dm, ls, &supp));
590:               PetscCall(DMPlexGetSupportSize(dm, ls, &suppSize));
591:               for (PetscInt s = 0; s < suppSize; ++s) {
592:                 PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
593:                 if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
594:                   ls = supp[s];
595:                   break;
596:                 }
597:               }
598:             }
599:             PetscCall(DMGetCellDS(dm, ls, ds, NULL));
600:           }
601:           if (ls >= 0) break;
602:         }
603:       }
604:       PetscCall(ISRestoreIndices(labelIS, &points));
605:     }
606:     PetscCall(ISDestroy(&labelIS));
607:     if (ls >= 0) break;
608:   }
609:   if (point) *point = ls;
610:   PetscCall(DMDestroy(&plex));
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

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

617:   There are several different scenarios:

619:   1) Volumetric mesh with volumetric auxiliary data

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

623:   2) Boundary mesh with boundary auxiliary data

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

627:   3) Volumetric mesh with boundary auxiliary data

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

631:   4) Volumetric input mesh with boundary output mesh

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

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

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

639:   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.
640:     - 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
641:       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
642:       dual spaces for the boundary from our input spaces.
643:     - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.

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

647:   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.
648: */
649: 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)
650: {
651:   DM               plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
652:   DMEnclosureType  encIn, encAux;
653:   PetscDS          ds = NULL, dsIn = NULL, dsAux = NULL;
654:   Vec              localA = NULL, tv;
655:   IS               fieldIS;
656:   PetscSection     section;
657:   PetscDualSpace  *sp, *cellsp, *spIn, *cellspIn;
658:   PetscTabulation *T = NULL, *TAux = NULL;
659:   PetscInt        *Nc;
660:   PetscInt         dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, minHeightIn, minHeightAux, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
661:   PetscBool       *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, isCohesiveIn = PETSC_FALSE, transform;
662:   DMField          coordField;
663:   DMLabel          depthLabel;
664:   PetscQuadrature  allPoints = NULL;

666:   PetscFunctionBegin;
667:   if (localU) PetscCall(VecGetDM(localU, &dmIn));
668:   else dmIn = dm;
669:   PetscCall(DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA));
670:   if (localA) PetscCall(VecGetDM(localA, &dmAux));
671:   else dmAux = NULL;
672:   PetscCall(DMConvert(dm, DMPLEX, &plex));
673:   PetscCall(DMConvert(dmIn, DMPLEX, &plexIn));
674:   PetscCall(DMGetEnclosureRelation(dmIn, dm, &encIn));
675:   PetscCall(DMGetEnclosureRelation(dmAux, dm, &encAux));
676:   PetscCall(DMGetDimension(dm, &dim));
677:   PetscCall(DMPlexGetVTKCellHeight(plex, &minHeight));
678:   PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
679:   PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));
680:   PetscCall(DMHasBasisTransform(dm, &transform));
681:   /* Auxiliary information can only be used with interpolation of field functions */
682:   if (dmAux) {
683:     PetscCall(DMConvert(dmAux, DMPLEX, &plexAux));
684:     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");
685:   }
686:   if (localU && localU != localX) PetscCall(DMPlexInsertBoundaryValues(plexIn, PETSC_TRUE, localU, time, NULL, NULL, NULL));
687:   PetscCall(DMGetCoordinateField(dm, &coordField));
688:   PetscCheck(coordField, PETSC_COMM_SELF, PETSC_ERR_USER, "DM must have a coordinate field");
689:   /**** No collective calls below this point ****/
690:   /* Determine height for iteration of all meshes */
691:   {
692:     DMPolytopeType ct, ctIn, ctAux;
693:     PetscInt       lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
694:     PetscInt       dim = -1, dimIn = -1, dimAux = -1;

696:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd));
697:     if (pEnd > pStart) {
698:       PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL));
699:       p = lStart < 0 ? pStart : lStart;
700:       PetscCall(DMPlexGetCellType(plex, p, &ct));
701:       dim = DMPolytopeTypeGetDim(ct);
702:       PetscCall(DMPlexGetVTKCellHeight(plexIn, &minHeightIn));
703:       PetscCall(DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL));
704:       PetscCall(DMPlexGetCellType(plexIn, pStartIn, &ctIn));
705:       dimIn = DMPolytopeTypeGetDim(ctIn);
706:       if (dmAux) {
707:         PetscCall(DMPlexGetVTKCellHeight(plexAux, &minHeightAux));
708:         PetscCall(DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux));
709:         if (pStartAux < pEndAux) {
710:           PetscCall(DMPlexGetCellType(plexAux, pStartAux, &ctAux));
711:           dimAux = DMPolytopeTypeGetDim(ctAux);
712:         }
713:       } else dimAux = dim;
714:     } else {
715:       PetscCall(DMDestroy(&plex));
716:       PetscCall(DMDestroy(&plexIn));
717:       if (dmAux) PetscCall(DMDestroy(&plexAux));
718:       PetscFunctionReturn(PETSC_SUCCESS);
719:     }
720:     if (dim < 0) {
721:       DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;

723:       /* Fall back to determination based on being a submesh */
724:       PetscCall(DMPlexGetSubpointMap(plex, &spmap));
725:       PetscCall(DMPlexGetSubpointMap(plexIn, &spmapIn));
726:       if (plexAux) PetscCall(DMPlexGetSubpointMap(plexAux, &spmapAux));
727:       dim    = spmap ? 1 : 0;
728:       dimIn  = spmapIn ? 1 : 0;
729:       dimAux = spmapAux ? 1 : 0;
730:     }
731:     {
732:       PetscInt dimProj   = PetscMin(PetscMin(dim, dimIn), (dimAux < 0 ? PETSC_MAX_INT : dimAux));
733:       PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;

735:       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");
736:       if (dimProj < dim) minHeight = 1;
737:       htInc    = dim - dimProj;
738:       htIncIn  = dimIn - dimProj;
739:       htIncAux = dimAuxEff - dimProj;
740:     }
741:   }
742:   PetscCall(DMPlexGetDepth(plex, &depth));
743:   PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
744:   PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
745:   maxHeight = PetscMax(maxHeight, minHeight);
746:   PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
747:   PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, NULL, &ds));
748:   if (!ds) PetscCall(DMGetDS(dm, &ds));
749:   PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, minHeight, NULL, &dsIn));
750:   if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
751:   PetscCall(PetscDSGetNumFields(ds, &Nf));
752:   PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
753:   PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
754:   if (isCohesiveIn) --htIncIn; // Should be rearranged
755:   PetscCall(DMGetNumFields(dm, &NfTot));
756:   PetscCall(DMFindRegionNum(dm, ds, &regionNum));
757:   PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
758:   PetscCall(PetscDSIsCohesive(ds, &isCohesive));
759:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
760:   PetscCall(DMGetLocalSection(dm, &section));
761:   if (dmAux) {
762:     PetscCall(DMGetDS(dmAux, &dsAux));
763:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
764:   }
765:   PetscCall(PetscDSGetComponents(ds, &Nc));
766:   PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
767:   if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
768:   else {
769:     cellsp   = sp;
770:     cellspIn = spIn;
771:   }
772:   /* Get cell dual spaces */
773:   for (f = 0; f < Nf; ++f) {
774:     PetscDiscType disctype;

776:     PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
777:     if (disctype == PETSC_DISC_FE) {
778:       PetscFE fe;

780:       isFE[f] = PETSC_TRUE;
781:       hasFE   = PETSC_TRUE;
782:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
783:       PetscCall(PetscFEGetDualSpace(fe, &cellsp[f]));
784:     } else if (disctype == PETSC_DISC_FV) {
785:       PetscFV fv;

787:       isFE[f] = PETSC_FALSE;
788:       hasFV   = PETSC_TRUE;
789:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fv));
790:       PetscCall(PetscFVGetDualSpace(fv, &cellsp[f]));
791:     } else {
792:       isFE[f]   = PETSC_FALSE;
793:       cellsp[f] = NULL;
794:     }
795:   }
796:   for (f = 0; f < NfIn; ++f) {
797:     PetscDiscType disctype;

799:     PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
800:     if (disctype == PETSC_DISC_FE) {
801:       PetscFE fe;

803:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
804:       PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
805:     } else if (disctype == PETSC_DISC_FV) {
806:       PetscFV fv;

808:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv));
809:       PetscCall(PetscFVGetDualSpace(fv, &cellspIn[f]));
810:     } else {
811:       cellspIn[f] = NULL;
812:     }
813:   }
814:   for (f = 0; f < Nf; ++f) {
815:     if (!htInc) {
816:       sp[f] = cellsp[f];
817:     } else PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]));
818:   }
819:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
820:     PetscFE          fem, subfem;
821:     PetscDiscType    disctype;
822:     const PetscReal *points;
823:     PetscInt         numPoints;

825:     PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
826:     PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
827:     PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
828:     PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
829:     for (f = 0; f < NfIn; ++f) {
830:       if (!htIncIn) {
831:         spIn[f] = cellspIn[f];
832:       } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));

834:       PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
835:       if (disctype != PETSC_DISC_FE) continue;
836:       PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem));
837:       if (!htIncIn) {
838:         subfem = fem;
839:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncIn, &subfem));
840:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]));
841:     }
842:     for (f = 0; f < NfAux; ++f) {
843:       PetscCall(PetscDSGetDiscType_Internal(dsAux, f, &disctype));
844:       if (disctype != PETSC_DISC_FE) continue;
845:       PetscCall(PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem));
846:       if (!htIncAux) {
847:         subfem = fem;
848:       } else PetscCall(PetscFEGetHeightSubspace(fem, htIncAux, &subfem));
849:       PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]));
850:     }
851:   }
852:   /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
853:   for (h = minHeight; h <= maxHeight; h++) {
854:     PetscInt     hEff     = h - minHeight + htInc;
855:     PetscInt     hEffIn   = h - minHeight + htIncIn;
856:     PetscInt     hEffAux  = h - minHeight + htIncAux;
857:     PetscDS      dsEff    = ds;
858:     PetscDS      dsEffIn  = dsIn;
859:     PetscDS      dsEffAux = dsAux;
860:     PetscScalar *values;
861:     PetscBool   *fieldActive;
862:     PetscInt     maxDegree;
863:     PetscInt     pStart, pEnd, p, lStart, spDim, totDim, numValues;
864:     IS           heightIS;

866:     if (h > minHeight) {
867:       for (f = 0; f < Nf; ++f) PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]));
868:     }
869:     PetscCall(DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd));
870:     PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL));
871:     PetscCall(DMLabelGetStratumIS(depthLabel, depth - h, &heightIS));
872:     if (pEnd <= pStart) {
873:       PetscCall(ISDestroy(&heightIS));
874:       continue;
875:     }
876:     /* Compute totDim, the number of dofs in the closure of a point at this height */
877:     totDim = 0;
878:     for (f = 0; f < Nf; ++f) {
879:       PetscBool cohesive;

881:       if (!sp[f]) continue;
882:       PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
883:       PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
884:       totDim += spDim;
885:       if (isCohesive && !cohesive) totDim += spDim;
886:     }
887:     p = lStart < 0 ? pStart : lStart;
888:     PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
889:     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);
890:     if (!totDim) {
891:       PetscCall(ISDestroy(&heightIS));
892:       continue;
893:     }
894:     if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
895:     /* Compute totDimIn, the number of dofs in the closure of a point at this height */
896:     if (localU) {
897:       PetscInt totDimIn, pIn, numValuesIn;

899:       totDimIn = 0;
900:       for (f = 0; f < NfIn; ++f) {
901:         PetscBool cohesive;

903:         if (!spIn[f]) continue;
904:         PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
905:         PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
906:         totDimIn += spDim;
907:         if (isCohesiveIn && !cohesive) totDimIn += spDim;
908:       }
909:       PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
910:       PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
911:       // TODO We could check that pIn is a cohesive cell for this check
912:       PetscCheck(isCohesiveIn || (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);
913:       if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
914:     }
915:     if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
916:     /* Loop over points at this height */
917:     PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
918:     PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
919:     {
920:       const PetscInt *fields;

922:       PetscCall(ISGetIndices(fieldIS, &fields));
923:       for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
924:       for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
925:       PetscCall(ISRestoreIndices(fieldIS, &fields));
926:     }
927:     if (label) {
928:       PetscInt i;

930:       for (i = 0; i < numIds; ++i) {
931:         IS              pointIS, isectIS;
932:         const PetscInt *points;
933:         PetscInt        n;
934:         PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
935:         PetscQuadrature quad = NULL;

937:         PetscCall(DMLabelGetStratumIS(label, ids[i], &pointIS));
938:         if (!pointIS) continue; /* No points with that id on this process */
939:         PetscCall(ISIntersect(pointIS, heightIS, &isectIS));
940:         PetscCall(ISDestroy(&pointIS));
941:         if (!isectIS) continue;
942:         PetscCall(ISGetLocalSize(isectIS, &n));
943:         PetscCall(ISGetIndices(isectIS, &points));
944:         PetscCall(DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree));
945:         if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad));
946:         if (!quad) {
947:           if (!h && allPoints) {
948:             quad      = allPoints;
949:             allPoints = NULL;
950:           } else {
951:             PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad));
952:           }
953:         }
954:         PetscBool computeFaceGeom = htInc && h == minHeight ? PETSC_TRUE : PETSC_FALSE;

956:         if (n) {
957:           PetscInt depth, dep;

959:           PetscCall(DMPlexGetDepth(dm, &depth));
960:           PetscCall(DMPlexGetPointDepth(dm, points[0], &dep));
961:           if (dep < depth && h == minHeight) computeFaceGeom = PETSC_TRUE;
962:         }
963:         PetscCall(DMFieldCreateFEGeom(coordField, isectIS, quad, computeFaceGeom, &fegeom));
964:         for (p = 0; p < n; ++p) {
965:           const PetscInt point = points[p];

967:           PetscCall(PetscArrayzero(values, numValues));
968:           PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
969:           PetscCall(DMPlexSetActivePoint(dm, point));
970:           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));
971:           if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
972:           PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
973:         }
974:         PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
975:         PetscCall(PetscFEGeomDestroy(&fegeom));
976:         PetscCall(PetscQuadratureDestroy(&quad));
977:         PetscCall(ISRestoreIndices(isectIS, &points));
978:         PetscCall(ISDestroy(&isectIS));
979:       }
980:     } else {
981:       PetscFEGeom    *fegeom = NULL, *chunkgeom = NULL;
982:       PetscQuadrature quad = NULL;
983:       IS              pointIS;

985:       PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS));
986:       PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
987:       if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
988:       if (!quad) {
989:         if (!h && allPoints) {
990:           quad      = allPoints;
991:           allPoints = NULL;
992:         } else {
993:           PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad));
994:         }
995:       }
996:       PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom));
997:       for (p = pStart; p < pEnd; ++p) {
998:         PetscCall(PetscArrayzero(values, numValues));
999:         PetscCall(PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom));
1000:         PetscCall(DMPlexSetActivePoint(dm, p));
1001:         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));
1002:         if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values));
1003:         PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode));
1004:       }
1005:       PetscCall(PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom));
1006:       PetscCall(PetscFEGeomDestroy(&fegeom));
1007:       PetscCall(PetscQuadratureDestroy(&quad));
1008:       PetscCall(ISDestroy(&pointIS));
1009:     }
1010:     PetscCall(ISDestroy(&heightIS));
1011:     PetscCall(DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values));
1012:     PetscCall(DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive));
1013:   }
1014:   /* Cleanup */
1015:   if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
1016:     for (f = 0; f < NfIn; ++f) PetscCall(PetscTabulationDestroy(&T[f]));
1017:     for (f = 0; f < NfAux; ++f) PetscCall(PetscTabulationDestroy(&TAux[f]));
1018:     PetscCall(PetscFree2(T, TAux));
1019:   }
1020:   PetscCall(PetscQuadratureDestroy(&allPoints));
1021:   PetscCall(PetscFree3(isFE, sp, spIn));
1022:   if (maxHeight > 0) PetscCall(PetscFree2(cellsp, cellspIn));
1023:   PetscCall(DMDestroy(&plex));
1024:   PetscCall(DMDestroy(&plexIn));
1025:   if (dmAux) PetscCall(DMDestroy(&plexAux));
1026:   PetscFunctionReturn(PETSC_SUCCESS);
1027: }

1029: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1030: {
1031:   PetscFunctionBegin;
1032:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1033:   PetscFunctionReturn(PETSC_SUCCESS);
1034: }

1036: 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)
1037: {
1038:   PetscFunctionBegin;
1039:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1040:   PetscFunctionReturn(PETSC_SUCCESS);
1041: }

1043: 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)
1044: {
1045:   PetscFunctionBegin;
1046:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: 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)
1051: {
1052:   PetscFunctionBegin;
1053:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1054:   PetscFunctionReturn(PETSC_SUCCESS);
1055: }

1057: 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)
1058: {
1059:   PetscFunctionBegin;
1060:   PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1061:   PetscFunctionReturn(PETSC_SUCCESS);
1062: }