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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`
19: @*/
20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
21: {
23: *point = ((DM_Plex *)dm->data)->activePoint;
24: return 0;
25: }
27: /*@
28: DMPlexSetActivePoint - Set the point on which projection is currently working
30: Not collective
32: Input Parameters:
33: + dm - the `DM`
34: - point - The mesh point involved in the current projection
36: Level: developer
38: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexGetActivePoint()`
39: @*/
40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
41: {
43: ((DM_Plex *)dm->data)->activePoint = point;
44: return 0;
45: }
47: /*
48: DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point
50: Input Parameters:
51: + dm - The output `DM`
52: . ds - The output `DS`
53: . dmIn - The input `DM`
54: . dsIn - The input `DS`
55: . time - The time for this evaluation
56: . fegeom - The FE geometry for this point
57: . fvgeom - The FV geometry for this point
58: . isFE - Flag indicating whether each output field has an FE discretization
59: . sp - The output `PetscDualSpace` for each field
60: . funcs - The evaluation function for each field
61: - ctxs - The user context for each field
63: Output Parameter:
64: . values - The value for each dual basis vector in the output dual space
66: Level: developer
68: .seealso:[](chapter_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;
76: DMGetCoordinateDim(dmIn, &coordDim);
77: DMHasBasisTransform(dmIn, &transform);
78: PetscDSGetNumFields(ds, &Nf);
79: PetscDSGetComponents(ds, &Nc);
80: PetscDSIsCohesive(ds, &isCohesive);
81: /* Get values for closure */
82: isAffine = fegeom->isAffine;
83: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
84: void *const ctx = ctxs ? ctxs[f] : NULL;
85: PetscBool cohesive;
87: if (!sp[f]) continue;
88: PetscDSGetCohesive(ds, f, &cohesive);
89: PetscDualSpaceGetDimension(sp[f], &spDim);
90: if (funcs[f]) {
91: if (isFE[f]) {
92: PetscQuadrature allPoints;
93: PetscInt q, dim, numPoints;
94: const PetscReal *points;
95: PetscScalar *pointEval;
96: PetscReal *x;
97: DM rdm;
99: PetscDualSpaceGetDM(sp[f], &rdm);
100: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
101: PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
102: DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
103: DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x);
104: PetscArrayzero(pointEval, numPoints * Nc[f]);
105: for (q = 0; q < numPoints; q++, tp++) {
106: const PetscReal *v0;
108: if (isAffine) {
109: const PetscReal *refpoint = &points[q * dim];
110: PetscReal injpoint[3] = {0., 0., 0.};
112: if (dim != fegeom->dim) {
113: if (isCohesive) {
114: /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
115: for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
116: refpoint = injpoint;
117: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %" PetscInt_FMT " != %" PetscInt_FMT " dual basis spatial dimension", fegeom->dim, dim);
118: }
119: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
120: v0 = x;
121: } else {
122: v0 = &fegeom->v[tp * coordDim];
123: }
124: if (transform) {
125: DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx);
126: v0 = x;
127: }
128: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx);
129: }
130: /* Transform point evaluations pointEval[q,c] */
131: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
132: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
133: DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x);
134: DMRestoreWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
135: v += spDim;
136: if (isCohesive && !cohesive) {
137: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
138: }
139: } else {
140: for (d = 0; d < spDim; ++d, ++v) PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
141: }
142: } else {
143: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
144: if (isCohesive && !cohesive) {
145: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
146: }
147: }
148: }
149: return 0;
150: }
152: /*
153: DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point
155: Input Parameters:
156: + dm - The output DM
157: . ds - The output DS
158: . dmIn - The input DM
159: . dsIn - The input DS
160: . dmAux - The auxiliary DM, which is always for the input space
161: . dsAux - The auxiliary DS, which is always for the input space
162: . time - The time for this evaluation
163: . localU - The local solution
164: . localA - The local auziliary fields
165: . cgeom - The FE geometry for this point
166: . sp - The output PetscDualSpace for each field
167: . p - The point in the output DM
168: . T - Input basis and derivatives for each field tabulated on the quadrature points
169: . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
170: . funcs - The evaluation function for each field
171: - ctxs - The user context for each field
173: Output Parameter:
174: . values - The value for each dual basis vector in the output dual space
176: 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;
198: PetscDSGetNumFields(ds, &Nf);
199: PetscDSGetComponents(ds, &Nc);
200: PetscDSIsCohesive(ds, &isCohesive);
201: PetscDSGetNumFields(dsIn, &NfIn);
202: PetscDSGetComponentOffsets(dsIn, &uOff);
203: PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
204: PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
205: PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
206: PetscDSGetConstants(dsIn, &numConstants, &constants);
207: DMHasBasisTransform(dmIn, &transform);
208: DMGetLocalSection(dmIn, §ion);
209: DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
210: DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
211: if (dmAux) {
212: PetscInt subp;
214: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
215: PetscDSGetNumFields(dsAux, &NfAux);
216: DMGetLocalSection(dmAux, §ionAux);
217: PetscDSGetComponentOffsets(dsAux, &aOff);
218: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
219: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
220: 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: PetscDSGetCohesive(ds, f, &cohesive);
243: 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: PetscDualSpaceGetDM(sp[f], &dm);
252: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
253: PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
254: 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: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
265: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
266: if (transform) 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: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
270: 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: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
278: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
279: return 0;
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;
298: PetscDSGetNumFields(ds, &Nf);
299: PetscDSGetComponents(ds, &Nc);
300: PetscDSGetComponentOffsets(ds, &uOff);
301: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
302: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
303: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
304: PetscDSGetConstants(ds, &numConstants, &constants);
305: DMGetLocalSection(dm, §ion);
306: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
307: if (dmAux) {
308: PetscInt subp;
310: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
311: PetscDSGetNumFields(dsAux, &NfAux);
312: DMGetLocalSection(dmAux, §ionAux);
313: PetscDSGetComponentOffsets(dsAux, &aOff);
314: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
315: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
316: 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: PetscDualSpaceGetDimension(sp[f], &spDim);
347: if (!funcs[f]) {
348: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
349: continue;
350: }
351: PetscDualSpaceGetDM(sp[f], &dm);
352: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
353: PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
354: 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: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
371: if (dsAux) 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: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
375: DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
376: v += spDim;
377: }
378: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
379: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
380: return 0;
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;
389: DMGetDimension(dm, &dim);
390: DMGetCoordinateDim(dm, &dimEmbed);
391: if (hasFV) DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);
392: switch (type) {
393: case DM_BC_ESSENTIAL:
394: case DM_BC_NATURAL:
395: 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: 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: 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: return 0;
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: if (!dim) {
416: PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
417: return 0;
418: }
419: numPoints = 0;
420: for (f = 0; f < Nf; ++f) {
421: if (funcs[f]) {
422: PetscQuadrature fAllPoints;
423: PetscInt fNumPoints;
425: PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
426: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
427: numPoints += fNumPoints;
428: }
429: }
430: PetscMalloc1(dim * numPoints, &points);
431: numPoints = 0;
432: for (f = 0; f < Nf; ++f) {
433: if (funcs[f]) {
434: PetscQuadrature fAllPoints;
435: PetscInt qdim, fNumPoints, q;
436: const PetscReal *fPoints;
438: PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
439: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
441: for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
442: numPoints += fNumPoints;
443: }
444: }
445: PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
446: PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL);
447: return 0;
448: }
450: /*@C
451: DMGetFirstLabeledPoint - Find first labeled point p_o in odm such that the corresponding point p in dm has the specified height. Return p and the corresponding ds.
453: Input Parameters:
454: dm - the `DM`
455: odm - the enclosing `DM`
456: label - label for `DM` domain, or NULL for whole domain
457: numIds - the number of ids
458: ids - An array of the label ids in sequence for the domain
459: height - Height of target cells in `DMPLEX` topology
461: Output Parameters:
462: point - the first labeled point
463: ds - the ds corresponding to the first labeled point
465: Level: developer
467: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
468: @*/
469: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
470: {
471: DM plex;
472: DMEnclosureType enc;
473: PetscInt ls = -1;
475: if (point) *point = -1;
476: if (!label) return 0;
477: DMGetEnclosureRelation(dm, odm, &enc);
478: DMConvert(dm, DMPLEX, &plex);
479: for (PetscInt i = 0; i < numIds; ++i) {
480: IS labelIS;
481: PetscInt num_points, pStart, pEnd;
482: DMLabelGetStratumIS(label, ids[i], &labelIS);
483: if (!labelIS) continue; /* No points with that id on this process */
484: DMPlexGetHeightStratum(plex, height, &pStart, &pEnd);
485: ISGetSize(labelIS, &num_points);
486: if (num_points) {
487: const PetscInt *points;
488: ISGetIndices(labelIS, &points);
489: for (PetscInt i = 0; i < num_points; i++) {
490: PetscInt point;
491: DMGetEnclosurePoint(dm, odm, enc, points[i], &point);
492: if (pStart <= point && point < pEnd) {
493: ls = point;
494: if (ds) DMGetCellDS(dm, ls, ds);
495: }
496: }
497: ISRestoreIndices(labelIS, &points);
498: }
499: ISDestroy(&labelIS);
500: if (ls >= 0) break;
501: }
502: if (point) *point = ls;
503: DMDestroy(&plex);
504: return 0;
505: }
507: /*
508: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
510: There are several different scenarios:
512: 1) Volumetric mesh with volumetric auxiliary data
514: Here minHeight=0 since we loop over cells.
516: 2) Boundary mesh with boundary auxiliary data
518: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
520: 3) Volumetric mesh with boundary auxiliary data
522: 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.
524: 4) Volumetric input mesh with boundary output mesh
526: Here we must get a subspace for the input DS
528: The maxHeight is used to support enforcement of constraints in DMForest.
530: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
532: 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.
533: - 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
534: 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
535: dual spaces for the boundary from our input spaces.
536: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
538: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
540: 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.
541: */
542: 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)
543: {
544: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
545: DMEnclosureType encIn, encAux;
546: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
547: Vec localA = NULL, tv;
548: IS fieldIS;
549: PetscSection section;
550: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
551: PetscTabulation *T = NULL, *TAux = NULL;
552: PetscInt *Nc;
553: PetscInt dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
554: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, transform;
555: DMField coordField;
556: DMLabel depthLabel;
557: PetscQuadrature allPoints = NULL;
559: if (localU) VecGetDM(localU, &dmIn);
560: else dmIn = dm;
561: DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA);
562: if (localA) VecGetDM(localA, &dmAux);
563: else dmAux = NULL;
564: DMConvert(dm, DMPLEX, &plex);
565: DMConvert(dmIn, DMPLEX, &plexIn);
566: DMGetEnclosureRelation(dmIn, dm, &encIn);
567: DMGetEnclosureRelation(dmAux, dm, &encAux);
568: DMGetDimension(dm, &dim);
569: DMPlexGetVTKCellHeight(plex, &minHeight);
570: DMGetBasisTransformDM_Internal(dm, &tdm);
571: DMGetBasisTransformVec_Internal(dm, &tv);
572: DMHasBasisTransform(dm, &transform);
573: /* Auxiliary information can only be used with interpolation of field functions */
574: if (dmAux) {
575: DMConvert(dmAux, DMPLEX, &plexAux);
577: }
578: if (localU && localU != localX) DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);
579: DMGetCoordinateField(dm, &coordField);
580: /**** No collective calls below this point ****/
581: /* Determine height for iteration of all meshes */
582: {
583: DMPolytopeType ct, ctIn, ctAux;
584: PetscInt minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
585: PetscInt dim = -1, dimIn = -1, dimAux = -1;
587: DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
588: if (pEnd > pStart) {
589: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL);
590: p = lStart < 0 ? pStart : lStart;
591: DMPlexGetCellType(plex, p, &ct);
592: dim = DMPolytopeTypeGetDim(ct);
593: DMPlexGetVTKCellHeight(plexIn, &minHeightIn);
594: DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL);
595: DMPlexGetCellType(plexIn, pStartIn, &ctIn);
596: dimIn = DMPolytopeTypeGetDim(ctIn);
597: if (dmAux) {
598: DMPlexGetVTKCellHeight(plexAux, &minHeightAux);
599: DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux);
600: if (pStartAux < pEndAux) {
601: DMPlexGetCellType(plexAux, pStartAux, &ctAux);
602: dimAux = DMPolytopeTypeGetDim(ctAux);
603: }
604: } else dimAux = dim;
605: } else {
606: DMDestroy(&plex);
607: DMDestroy(&plexIn);
608: if (dmAux) DMDestroy(&plexAux);
609: return 0;
610: }
611: if (dim < 0) {
612: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
614: /* Fall back to determination based on being a submesh */
615: DMPlexGetSubpointMap(plex, &spmap);
616: DMPlexGetSubpointMap(plexIn, &spmapIn);
617: if (plexAux) DMPlexGetSubpointMap(plexAux, &spmapAux);
618: dim = spmap ? 1 : 0;
619: dimIn = spmapIn ? 1 : 0;
620: dimAux = spmapAux ? 1 : 0;
621: }
622: {
623: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), (dimAux < 0 ? PETSC_MAX_INT : dimAux));
624: PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;
627: if (dimProj < dim) minHeight = 1;
628: htInc = dim - dimProj;
629: htIncIn = dimIn - dimProj;
630: htIncAux = dimAuxEff - dimProj;
631: }
632: }
633: DMPlexGetDepth(plex, &depth);
634: DMPlexGetDepthLabel(plex, &depthLabel);
635: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
636: maxHeight = PetscMax(maxHeight, minHeight);
638: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, 0, NULL, &ds);
639: if (!ds) DMGetDS(dm, &ds);
640: DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
641: if (!dsIn) DMGetDS(dmIn, &dsIn);
642: PetscDSGetNumFields(ds, &Nf);
643: PetscDSGetNumFields(dsIn, &NfIn);
644: DMGetNumFields(dm, &NfTot);
645: DMFindRegionNum(dm, ds, ®ionNum);
646: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
647: PetscDSIsCohesive(ds, &isCohesive);
648: DMGetCoordinateDim(dm, &dimEmbed);
649: DMGetLocalSection(dm, §ion);
650: if (dmAux) {
651: DMGetDS(dmAux, &dsAux);
652: PetscDSGetNumFields(dsAux, &NfAux);
653: }
654: PetscDSGetComponents(ds, &Nc);
655: PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
656: if (maxHeight > 0) PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);
657: else {
658: cellsp = sp;
659: cellspIn = spIn;
660: }
661: /* Get cell dual spaces */
662: for (f = 0; f < Nf; ++f) {
663: PetscDiscType disctype;
665: PetscDSGetDiscType_Internal(ds, f, &disctype);
666: if (disctype == PETSC_DISC_FE) {
667: PetscFE fe;
669: isFE[f] = PETSC_TRUE;
670: hasFE = PETSC_TRUE;
671: PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
672: PetscFEGetDualSpace(fe, &cellsp[f]);
673: } else if (disctype == PETSC_DISC_FV) {
674: PetscFV fv;
676: isFE[f] = PETSC_FALSE;
677: hasFV = PETSC_TRUE;
678: PetscDSGetDiscretization(ds, f, (PetscObject *)&fv);
679: PetscFVGetDualSpace(fv, &cellsp[f]);
680: } else {
681: isFE[f] = PETSC_FALSE;
682: cellsp[f] = NULL;
683: }
684: }
685: for (f = 0; f < NfIn; ++f) {
686: PetscDiscType disctype;
688: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
689: if (disctype == PETSC_DISC_FE) {
690: PetscFE fe;
692: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe);
693: PetscFEGetDualSpace(fe, &cellspIn[f]);
694: } else if (disctype == PETSC_DISC_FV) {
695: PetscFV fv;
697: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv);
698: PetscFVGetDualSpace(fv, &cellspIn[f]);
699: } else {
700: cellspIn[f] = NULL;
701: }
702: }
703: for (f = 0; f < Nf; ++f) {
704: if (!htInc) {
705: sp[f] = cellsp[f];
706: } else PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);
707: }
708: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
709: PetscFE fem, subfem;
710: PetscDiscType disctype;
711: const PetscReal *points;
712: PetscInt numPoints;
715: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints);
716: PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
717: PetscMalloc2(NfIn, &T, NfAux, &TAux);
718: for (f = 0; f < NfIn; ++f) {
719: if (!htIncIn) {
720: spIn[f] = cellspIn[f];
721: } else PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);
723: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
724: if (disctype != PETSC_DISC_FE) continue;
725: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem);
726: if (!htIncIn) {
727: subfem = fem;
728: } else PetscFEGetHeightSubspace(fem, htIncIn, &subfem);
729: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
730: }
731: for (f = 0; f < NfAux; ++f) {
732: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
733: if (disctype != PETSC_DISC_FE) continue;
734: PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem);
735: if (!htIncAux) {
736: subfem = fem;
737: } else PetscFEGetHeightSubspace(fem, htIncAux, &subfem);
738: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
739: }
740: }
741: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
742: for (h = minHeight; h <= maxHeight; h++) {
743: PetscInt hEff = h - minHeight + htInc;
744: PetscInt hEffIn = h - minHeight + htIncIn;
745: PetscInt hEffAux = h - minHeight + htIncAux;
746: PetscDS dsEff = ds;
747: PetscDS dsEffIn = dsIn;
748: PetscDS dsEffAux = dsAux;
749: PetscScalar *values;
750: PetscBool *fieldActive;
751: PetscInt maxDegree;
752: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
753: IS heightIS;
755: if (h > minHeight) {
756: for (f = 0; f < Nf; ++f) PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);
757: }
758: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
759: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL);
760: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
761: if (pEnd <= pStart) {
762: ISDestroy(&heightIS);
763: continue;
764: }
765: /* Compute totDim, the number of dofs in the closure of a point at this height */
766: totDim = 0;
767: for (f = 0; f < Nf; ++f) {
768: PetscBool cohesive;
770: if (!sp[f]) continue;
771: PetscDSGetCohesive(ds, f, &cohesive);
772: PetscDualSpaceGetDimension(sp[f], &spDim);
773: totDim += spDim;
774: if (isCohesive && !cohesive) totDim += spDim;
775: }
776: p = lStart < 0 ? pStart : lStart;
777: DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
779: if (!totDim) {
780: ISDestroy(&heightIS);
781: continue;
782: }
783: if (htInc) PetscDSGetHeightSubspace(ds, hEff, &dsEff);
784: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
785: if (localU) {
786: PetscInt totDimIn, pIn, numValuesIn;
788: totDimIn = 0;
789: for (f = 0; f < NfIn; ++f) {
790: PetscBool cohesive;
792: if (!spIn[f]) continue;
793: PetscDSGetCohesive(dsIn, f, &cohesive);
794: PetscDualSpaceGetDimension(spIn[f], &spDim);
795: totDimIn += spDim;
796: if (isCohesive && !cohesive) totDimIn += spDim;
797: }
798: DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
799: DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
801: if (htIncIn) PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);
802: }
803: if (htIncAux) PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);
804: /* Loop over points at this height */
805: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
806: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
807: {
808: const PetscInt *fields;
810: ISGetIndices(fieldIS, &fields);
811: for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
812: for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
813: ISRestoreIndices(fieldIS, &fields);
814: }
815: if (label) {
816: PetscInt i;
818: for (i = 0; i < numIds; ++i) {
819: IS pointIS, isectIS;
820: const PetscInt *points;
821: PetscInt n;
822: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
823: PetscQuadrature quad = NULL;
825: DMLabelGetStratumIS(label, ids[i], &pointIS);
826: if (!pointIS) continue; /* No points with that id on this process */
827: ISIntersect(pointIS, heightIS, &isectIS);
828: ISDestroy(&pointIS);
829: if (!isectIS) continue;
830: ISGetLocalSize(isectIS, &n);
831: ISGetIndices(isectIS, &points);
832: DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree);
833: if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad);
834: if (!quad) {
835: if (!h && allPoints) {
836: quad = allPoints;
837: allPoints = NULL;
838: } else {
839: PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad);
840: }
841: }
842: DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
843: for (p = 0; p < n; ++p) {
844: const PetscInt point = points[p];
846: PetscArrayzero(values, numValues);
847: PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom);
848: DMPlexSetActivePoint(dm, point);
849: 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);
850: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);
851: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
852: }
853: PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom);
854: PetscFEGeomDestroy(&fegeom);
855: PetscQuadratureDestroy(&quad);
856: ISRestoreIndices(isectIS, &points);
857: ISDestroy(&isectIS);
858: }
859: } else {
860: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
861: PetscQuadrature quad = NULL;
862: IS pointIS;
864: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS);
865: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
866: if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
867: if (!quad) {
868: if (!h && allPoints) {
869: quad = allPoints;
870: allPoints = NULL;
871: } else {
872: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad);
873: }
874: }
875: DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
876: for (p = pStart; p < pEnd; ++p) {
877: PetscArrayzero(values, numValues);
878: PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom);
879: DMPlexSetActivePoint(dm, p);
880: 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);
881: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);
882: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
883: }
884: PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom);
885: PetscFEGeomDestroy(&fegeom);
886: PetscQuadratureDestroy(&quad);
887: ISDestroy(&pointIS);
888: }
889: ISDestroy(&heightIS);
890: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
891: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
892: }
893: /* Cleanup */
894: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
895: for (f = 0; f < NfIn; ++f) PetscTabulationDestroy(&T[f]);
896: for (f = 0; f < NfAux; ++f) PetscTabulationDestroy(&TAux[f]);
897: PetscFree2(T, TAux);
898: }
899: PetscQuadratureDestroy(&allPoints);
900: PetscFree3(isFE, sp, spIn);
901: if (maxHeight > 0) PetscFree2(cellsp, cellspIn);
902: DMDestroy(&plex);
903: DMDestroy(&plexIn);
904: if (dmAux) DMDestroy(&plexAux);
905: return 0;
906: }
908: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
909: {
910: DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
911: return 0;
912: }
914: PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
915: {
916: DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
917: return 0;
918: }
920: 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)
921: {
922: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
923: return 0;
924: }
926: 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)
927: {
928: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
929: return 0;
930: }
932: 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)
933: {
934: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX);
935: return 0;
936: }