Actual source code: plexproject.c
1: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: /*@
6: DMPlexGetActivePoint - Get the point on which projection is currently working
8: Not collective
10: Input Parameter:
11: . dm - the DM
13: Output Parameter:
14: . point - The mesh point involved in the current projection
16: Level: developer
18: .seealso: `DMPlexSetActivePoint()`
19: @*/
20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
21: {
23: *point = ((DM_Plex *)dm->data)->activePoint;
24: return 0;
25: }
27: /*@
28: DMPlexSetActivePoint - Set the point on which projection is currently working
30: Not collective
32: Input Parameters:
33: + dm - the DM
34: - point - The mesh point involved in the current projection
36: Level: developer
38: .seealso: `DMPlexGetActivePoint()`
39: @*/
40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
41: {
43: ((DM_Plex *)dm->data)->activePoint = point;
44: return 0;
45: }
47: /*
48: DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point
50: Input Parameters:
51: + dm - The output DM
52: . ds - The output DS
53: . dmIn - The input DM
54: . dsIn - The input DS
55: . time - The time for this evaluation
56: . fegeom - The FE geometry for this point
57: . fvgeom - The FV geometry for this point
58: . isFE - Flag indicating whether each output field has an FE discretization
59: . sp - The output PetscDualSpace for each field
60: . funcs - The evaluation function for each field
61: - ctxs - The user context for each field
63: Output Parameter:
64: . values - The value for each dual basis vector in the output dual space
66: Level: developer
68: .seealso: `DMProjectPoint_Field_Private()`
69: */
70: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], 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: Note: Not supported for FV
178: Level: developer
180: .seealso: `DMProjectPoint_Field_Private()`
181: */
182: 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[])
183: {
184: PetscSection section, sectionAux = NULL;
185: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
186: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
187: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
188: const PetscScalar *constants;
189: PetscReal *x;
190: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
191: PetscFEGeom fegeom;
192: const PetscInt dE = cgeom->dimEmbed;
193: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
194: PetscBool isAffine, isCohesive, transform;
197: PetscDSGetNumFields(ds, &Nf);
198: PetscDSGetComponents(ds, &Nc);
199: PetscDSIsCohesive(ds, &isCohesive);
200: PetscDSGetNumFields(dsIn, &NfIn);
201: PetscDSGetComponentOffsets(dsIn, &uOff);
202: PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
203: PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
204: PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
205: PetscDSGetConstants(dsIn, &numConstants, &constants);
206: DMHasBasisTransform(dmIn, &transform);
207: DMGetLocalSection(dmIn, §ion);
208: DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
209: DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
210: if (dmAux) {
211: PetscInt subp;
213: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
214: PetscDSGetNumFields(dsAux, &NfAux);
215: DMGetLocalSection(dmAux, §ionAux);
216: PetscDSGetComponentOffsets(dsAux, &aOff);
217: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
218: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
219: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
220: }
221: /* Get values for closure */
222: isAffine = cgeom->isAffine;
223: fegeom.dim = cgeom->dim;
224: fegeom.dimEmbed = cgeom->dimEmbed;
225: if (isAffine) {
226: fegeom.v = x;
227: fegeom.xi = cgeom->xi;
228: fegeom.J = cgeom->J;
229: fegeom.invJ = cgeom->invJ;
230: fegeom.detJ = cgeom->detJ;
231: }
232: for (f = 0, v = 0; f < Nf; ++f) {
233: PetscQuadrature allPoints;
234: PetscInt q, dim, numPoints;
235: const PetscReal *points;
236: PetscScalar *pointEval;
237: PetscBool cohesive;
238: DM dm;
240: if (!sp[f]) continue;
241: PetscDSGetCohesive(ds, f, &cohesive);
242: PetscDualSpaceGetDimension(sp[f], &spDim);
243: if (!funcs[f]) {
244: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
245: if (isCohesive && !cohesive) {
246: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
247: }
248: continue;
249: }
250: PetscDualSpaceGetDM(sp[f], &dm);
251: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
252: PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
253: DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
254: for (q = 0; q < numPoints; ++q, ++tp) {
255: if (isAffine) {
256: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
257: } else {
258: fegeom.v = &cgeom->v[tp * dE];
259: fegeom.J = &cgeom->J[tp * dE * dE];
260: fegeom.invJ = &cgeom->invJ[tp * dE * dE];
261: fegeom.detJ = &cgeom->detJ[tp];
262: }
263: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
264: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
265: if (transform) DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);
266: (*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]);
267: }
268: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
269: DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
270: v += spDim;
271: /* TODO: For now, set both sides equal, but this should use info from other support cell */
272: if (isCohesive && !cohesive) {
273: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
274: }
275: }
276: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
277: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
278: return 0;
279: }
281: 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[])
282: {
283: PetscSection section, sectionAux = NULL;
284: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
285: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
286: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
287: const PetscScalar *constants;
288: PetscReal *x;
289: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
290: PetscFEGeom fegeom, cgeom;
291: const PetscInt dE = fgeom->dimEmbed;
292: PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
293: PetscBool isAffine;
297: PetscDSGetNumFields(ds, &Nf);
298: PetscDSGetComponents(ds, &Nc);
299: PetscDSGetComponentOffsets(ds, &uOff);
300: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
301: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
302: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
303: PetscDSGetConstants(ds, &numConstants, &constants);
304: DMGetLocalSection(dm, §ion);
305: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
306: if (dmAux) {
307: PetscInt subp;
309: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
310: PetscDSGetNumFields(dsAux, &NfAux);
311: DMGetLocalSection(dmAux, §ionAux);
312: PetscDSGetComponentOffsets(dsAux, &aOff);
313: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
314: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
315: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
316: }
317: /* Get values for closure */
318: isAffine = fgeom->isAffine;
319: fegeom.n = NULL;
320: fegeom.J = NULL;
321: fegeom.v = NULL;
322: fegeom.xi = NULL;
323: cgeom.dim = fgeom->dim;
324: cgeom.dimEmbed = fgeom->dimEmbed;
325: if (isAffine) {
326: fegeom.v = x;
327: fegeom.xi = fgeom->xi;
328: fegeom.J = fgeom->J;
329: fegeom.invJ = fgeom->invJ;
330: fegeom.detJ = fgeom->detJ;
331: fegeom.n = fgeom->n;
333: cgeom.J = fgeom->suppJ[0];
334: cgeom.invJ = fgeom->suppInvJ[0];
335: cgeom.detJ = fgeom->suppDetJ[0];
336: }
337: for (f = 0, v = 0; f < Nf; ++f) {
338: PetscQuadrature allPoints;
339: PetscInt q, dim, numPoints;
340: const PetscReal *points;
341: PetscScalar *pointEval;
342: DM dm;
344: if (!sp[f]) continue;
345: PetscDualSpaceGetDimension(sp[f], &spDim);
346: if (!funcs[f]) {
347: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
348: continue;
349: }
350: PetscDualSpaceGetDM(sp[f], &dm);
351: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
352: PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL);
353: DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
354: for (q = 0; q < numPoints; ++q, ++tp) {
355: if (isAffine) {
356: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
357: } else {
358: fegeom.v = &fgeom->v[tp * dE];
359: fegeom.J = &fgeom->J[tp * dE * dE];
360: fegeom.invJ = &fgeom->invJ[tp * dE * dE];
361: fegeom.detJ = &fgeom->detJ[tp];
362: fegeom.n = &fgeom->n[tp * dE];
364: cgeom.J = &fgeom->suppJ[0][tp * dE * dE];
365: cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
366: cgeom.detJ = &fgeom->suppDetJ[0][tp];
367: }
368: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
369: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
370: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
371: (*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]);
372: }
373: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
374: DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval);
375: v += spDim;
376: }
377: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
378: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
379: return 0;
380: }
382: 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[])
383: {
384: PetscFVCellGeom fvgeom;
385: PetscInt dim, dimEmbed;
388: DMGetDimension(dm, &dim);
389: DMGetCoordinateDim(dm, &dimEmbed);
390: if (hasFV) DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);
391: switch (type) {
392: case DM_BC_ESSENTIAL:
393: case DM_BC_NATURAL:
394: DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode(**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values);
395: break;
396: case DM_BC_ESSENTIAL_FIELD:
397: case DM_BC_NATURAL_FIELD:
398: 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);
399: break;
400: case DM_BC_ESSENTIAL_BD_FIELD:
401: 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);
402: break;
403: default:
404: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
405: }
406: return 0;
407: }
409: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
410: {
411: PetscReal *points;
412: PetscInt f, numPoints;
414: if (!dim) {
415: PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
416: return 0;
417: }
418: numPoints = 0;
419: for (f = 0; f < Nf; ++f) {
420: if (funcs[f]) {
421: PetscQuadrature fAllPoints;
422: PetscInt fNumPoints;
424: PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
425: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
426: numPoints += fNumPoints;
427: }
428: }
429: PetscMalloc1(dim * numPoints, &points);
430: numPoints = 0;
431: for (f = 0; f < Nf; ++f) {
432: if (funcs[f]) {
433: PetscQuadrature fAllPoints;
434: PetscInt qdim, fNumPoints, q;
435: const PetscReal *fPoints;
437: PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL);
438: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
440: for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
441: numPoints += fNumPoints;
442: }
443: }
444: PetscQuadratureCreate(PETSC_COMM_SELF, allPoints);
445: PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL);
446: return 0;
447: }
449: /*@C
450: 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.
452: Input Parameters:
453: dm - the DM
454: odm - the enclosing DM
455: label - label for DM domain, or NULL for whole domain
456: numIds - the number of ids
457: ids - An array of the label ids in sequence for the domain
458: height - Height of target cells in DMPlex topology
460: Output Parameters:
461: point - the first labeled point
462: ds - the ds corresponding to the first labeled point
464: Level: developer
465: @*/
466: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
467: {
468: DM plex;
469: DMEnclosureType enc;
470: PetscInt ls = -1;
472: if (point) *point = -1;
473: if (!label) return 0;
474: DMGetEnclosureRelation(dm, odm, &enc);
475: DMConvert(dm, DMPLEX, &plex);
476: for (PetscInt i = 0; i < numIds; ++i) {
477: IS labelIS;
478: PetscInt num_points, pStart, pEnd;
479: DMLabelGetStratumIS(label, ids[i], &labelIS);
480: if (!labelIS) continue; /* No points with that id on this process */
481: DMPlexGetHeightStratum(plex, height, &pStart, &pEnd);
482: ISGetSize(labelIS, &num_points);
483: if (num_points) {
484: const PetscInt *points;
485: ISGetIndices(labelIS, &points);
486: for (PetscInt i = 0; i < num_points; i++) {
487: PetscInt point;
488: DMGetEnclosurePoint(dm, odm, enc, points[i], &point);
489: if (pStart <= point && point < pEnd) {
490: ls = point;
491: if (ds) DMGetCellDS(dm, ls, ds);
492: }
493: }
494: ISRestoreIndices(labelIS, &points);
495: }
496: ISDestroy(&labelIS);
497: if (ls >= 0) break;
498: }
499: if (point) *point = ls;
500: DMDestroy(&plex);
501: return 0;
502: }
504: /*
505: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
507: There are several different scenarios:
509: 1) Volumetric mesh with volumetric auxiliary data
511: Here minHeight=0 since we loop over cells.
513: 2) Boundary mesh with boundary auxiliary data
515: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
517: 3) Volumetric mesh with boundary auxiliary data
519: Here minHeight=1 and auxbd=PETSC_TRUE since we loop over faces and use data only supported on those faces. This is common when imposing Dirichlet boundary conditions.
521: 4) Volumetric input mesh with boundary output mesh
523: Here we must get a subspace for the input DS
525: The maxHeight is used to support enforcement of constraints in DMForest.
527: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
529: If we are using an input field (DM_BC_ESSENTIAL_FIELD or DM_BC_NATURAL_FIELD), we need to evaluate it at all the quadrature points of the dual basis functionals.
530: - We use effectiveHeight to mean the height above our incoming DS. For example, if the DS is for a submesh then the effective height is zero, whereas if the DS
531: is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When the effective height is nonzero, we need to extract
532: dual spaces for the boundary from our input spaces.
533: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
535: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
537: If we have a label, we iterate over those points. This will probably break the maxHeight functionality since we do not check the height of those points.
538: */
539: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, InsertMode mode, Vec localX)
540: {
541: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
542: DMEnclosureType encIn, encAux;
543: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
544: Vec localA = NULL, tv;
545: IS fieldIS;
546: PetscSection section;
547: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
548: PetscTabulation *T = NULL, *TAux = NULL;
549: PetscInt *Nc;
550: PetscInt dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
551: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, transform;
552: DMField coordField;
553: DMLabel depthLabel;
554: PetscQuadrature allPoints = NULL;
556: if (localU) VecGetDM(localU, &dmIn);
557: else dmIn = dm;
558: DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA);
559: if (localA) VecGetDM(localA, &dmAux);
560: else dmAux = NULL;
561: DMConvert(dm, DMPLEX, &plex);
562: DMConvert(dmIn, DMPLEX, &plexIn);
563: DMGetEnclosureRelation(dmIn, dm, &encIn);
564: DMGetEnclosureRelation(dmAux, dm, &encAux);
565: DMGetDimension(dm, &dim);
566: DMPlexGetVTKCellHeight(plex, &minHeight);
567: DMGetBasisTransformDM_Internal(dm, &tdm);
568: DMGetBasisTransformVec_Internal(dm, &tv);
569: DMHasBasisTransform(dm, &transform);
570: /* Auxiliary information can only be used with interpolation of field functions */
571: if (dmAux) {
572: DMConvert(dmAux, DMPLEX, &plexAux);
574: }
575: if (localU && localU != localX) DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);
576: DMGetCoordinateField(dm, &coordField);
577: /**** No collective calls below this point ****/
578: /* Determine height for iteration of all meshes */
579: {
580: DMPolytopeType ct, ctIn, ctAux;
581: PetscInt minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
582: PetscInt dim = -1, dimIn = -1, dimAux = -1;
584: DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
585: if (pEnd > pStart) {
586: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL);
587: p = lStart < 0 ? pStart : lStart;
588: DMPlexGetCellType(plex, p, &ct);
589: dim = DMPolytopeTypeGetDim(ct);
590: DMPlexGetVTKCellHeight(plexIn, &minHeightIn);
591: DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL);
592: DMPlexGetCellType(plexIn, pStartIn, &ctIn);
593: dimIn = DMPolytopeTypeGetDim(ctIn);
594: if (dmAux) {
595: DMPlexGetVTKCellHeight(plexAux, &minHeightAux);
596: DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux);
597: if (pStartAux < pEndAux) {
598: DMPlexGetCellType(plexAux, pStartAux, &ctAux);
599: dimAux = DMPolytopeTypeGetDim(ctAux);
600: }
601: } else dimAux = dim;
602: } else {
603: DMDestroy(&plex);
604: DMDestroy(&plexIn);
605: if (dmAux) DMDestroy(&plexAux);
606: return 0;
607: }
608: if (dim < 0) {
609: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
611: /* Fall back to determination based on being a submesh */
612: DMPlexGetSubpointMap(plex, &spmap);
613: DMPlexGetSubpointMap(plexIn, &spmapIn);
614: if (plexAux) DMPlexGetSubpointMap(plexAux, &spmapAux);
615: dim = spmap ? 1 : 0;
616: dimIn = spmapIn ? 1 : 0;
617: dimAux = spmapAux ? 1 : 0;
618: }
619: {
620: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), (dimAux < 0 ? PETSC_MAX_INT : dimAux));
621: PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;
624: if (dimProj < dim) minHeight = 1;
625: htInc = dim - dimProj;
626: htIncIn = dimIn - dimProj;
627: htIncAux = dimAuxEff - dimProj;
628: }
629: }
630: DMPlexGetDepth(plex, &depth);
631: DMPlexGetDepthLabel(plex, &depthLabel);
632: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
633: maxHeight = PetscMax(maxHeight, minHeight);
635: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, 0, NULL, &ds);
636: if (!ds) DMGetDS(dm, &ds);
637: DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
638: if (!dsIn) DMGetDS(dmIn, &dsIn);
639: PetscDSGetNumFields(ds, &Nf);
640: PetscDSGetNumFields(dsIn, &NfIn);
641: DMGetNumFields(dm, &NfTot);
642: DMFindRegionNum(dm, ds, ®ionNum);
643: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
644: PetscDSIsCohesive(ds, &isCohesive);
645: DMGetCoordinateDim(dm, &dimEmbed);
646: DMGetLocalSection(dm, §ion);
647: if (dmAux) {
648: DMGetDS(dmAux, &dsAux);
649: PetscDSGetNumFields(dsAux, &NfAux);
650: }
651: PetscDSGetComponents(ds, &Nc);
652: PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
653: if (maxHeight > 0) PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);
654: else {
655: cellsp = sp;
656: cellspIn = spIn;
657: }
658: /* Get cell dual spaces */
659: for (f = 0; f < Nf; ++f) {
660: PetscDiscType disctype;
662: PetscDSGetDiscType_Internal(ds, f, &disctype);
663: if (disctype == PETSC_DISC_FE) {
664: PetscFE fe;
666: isFE[f] = PETSC_TRUE;
667: hasFE = PETSC_TRUE;
668: PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
669: PetscFEGetDualSpace(fe, &cellsp[f]);
670: } else if (disctype == PETSC_DISC_FV) {
671: PetscFV fv;
673: isFE[f] = PETSC_FALSE;
674: hasFV = PETSC_TRUE;
675: PetscDSGetDiscretization(ds, f, (PetscObject *)&fv);
676: PetscFVGetDualSpace(fv, &cellsp[f]);
677: } else {
678: isFE[f] = PETSC_FALSE;
679: cellsp[f] = NULL;
680: }
681: }
682: for (f = 0; f < NfIn; ++f) {
683: PetscDiscType disctype;
685: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
686: if (disctype == PETSC_DISC_FE) {
687: PetscFE fe;
689: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe);
690: PetscFEGetDualSpace(fe, &cellspIn[f]);
691: } else if (disctype == PETSC_DISC_FV) {
692: PetscFV fv;
694: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv);
695: PetscFVGetDualSpace(fv, &cellspIn[f]);
696: } else {
697: cellspIn[f] = NULL;
698: }
699: }
700: for (f = 0; f < Nf; ++f) {
701: if (!htInc) {
702: sp[f] = cellsp[f];
703: } else PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);
704: }
705: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
706: PetscFE fem, subfem;
707: PetscDiscType disctype;
708: const PetscReal *points;
709: PetscInt numPoints;
712: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints);
713: PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
714: PetscMalloc2(NfIn, &T, NfAux, &TAux);
715: for (f = 0; f < NfIn; ++f) {
716: if (!htIncIn) {
717: spIn[f] = cellspIn[f];
718: } else PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);
720: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
721: if (disctype != PETSC_DISC_FE) continue;
722: PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem);
723: if (!htIncIn) {
724: subfem = fem;
725: } else PetscFEGetHeightSubspace(fem, htIncIn, &subfem);
726: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
727: }
728: for (f = 0; f < NfAux; ++f) {
729: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
730: if (disctype != PETSC_DISC_FE) continue;
731: PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem);
732: if (!htIncAux) {
733: subfem = fem;
734: } else PetscFEGetHeightSubspace(fem, htIncAux, &subfem);
735: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
736: }
737: }
738: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
739: for (h = minHeight; h <= maxHeight; h++) {
740: PetscInt hEff = h - minHeight + htInc;
741: PetscInt hEffIn = h - minHeight + htIncIn;
742: PetscInt hEffAux = h - minHeight + htIncAux;
743: PetscDS dsEff = ds;
744: PetscDS dsEffIn = dsIn;
745: PetscDS dsEffAux = dsAux;
746: PetscScalar *values;
747: PetscBool *fieldActive;
748: PetscInt maxDegree;
749: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
750: IS heightIS;
752: if (h > minHeight) {
753: for (f = 0; f < Nf; ++f) PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);
754: }
755: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
756: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL);
757: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
758: if (pEnd <= pStart) {
759: ISDestroy(&heightIS);
760: continue;
761: }
762: /* Compute totDim, the number of dofs in the closure of a point at this height */
763: totDim = 0;
764: for (f = 0; f < Nf; ++f) {
765: PetscBool cohesive;
767: if (!sp[f]) continue;
768: PetscDSGetCohesive(ds, f, &cohesive);
769: PetscDualSpaceGetDimension(sp[f], &spDim);
770: totDim += spDim;
771: if (isCohesive && !cohesive) totDim += spDim;
772: }
773: p = lStart < 0 ? pStart : lStart;
774: DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
776: if (!totDim) {
777: ISDestroy(&heightIS);
778: continue;
779: }
780: if (htInc) PetscDSGetHeightSubspace(ds, hEff, &dsEff);
781: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
782: if (localU) {
783: PetscInt totDimIn, pIn, numValuesIn;
785: totDimIn = 0;
786: for (f = 0; f < NfIn; ++f) {
787: PetscBool cohesive;
789: if (!spIn[f]) continue;
790: PetscDSGetCohesive(dsIn, f, &cohesive);
791: PetscDualSpaceGetDimension(spIn[f], &spDim);
792: totDimIn += spDim;
793: if (isCohesive && !cohesive) totDimIn += spDim;
794: }
795: DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
796: DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
798: if (htIncIn) PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);
799: }
800: if (htIncAux) PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);
801: /* Loop over points at this height */
802: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
803: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
804: {
805: const PetscInt *fields;
807: ISGetIndices(fieldIS, &fields);
808: for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
809: for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
810: ISRestoreIndices(fieldIS, &fields);
811: }
812: if (label) {
813: PetscInt i;
815: for (i = 0; i < numIds; ++i) {
816: IS pointIS, isectIS;
817: const PetscInt *points;
818: PetscInt n;
819: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
820: PetscQuadrature quad = NULL;
822: DMLabelGetStratumIS(label, ids[i], &pointIS);
823: if (!pointIS) continue; /* No points with that id on this process */
824: ISIntersect(pointIS, heightIS, &isectIS);
825: ISDestroy(&pointIS);
826: if (!isectIS) continue;
827: ISGetLocalSize(isectIS, &n);
828: ISGetIndices(isectIS, &points);
829: DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree);
830: if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad);
831: if (!quad) {
832: if (!h && allPoints) {
833: quad = allPoints;
834: allPoints = NULL;
835: } else {
836: PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad);
837: }
838: }
839: DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
840: for (p = 0; p < n; ++p) {
841: const PetscInt point = points[p];
843: PetscArrayzero(values, numValues);
844: PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom);
845: DMPlexSetActivePoint(dm, point);
846: 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);
847: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);
848: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
849: }
850: PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom);
851: PetscFEGeomDestroy(&fegeom);
852: PetscQuadratureDestroy(&quad);
853: ISRestoreIndices(isectIS, &points);
854: ISDestroy(&isectIS);
855: }
856: } else {
857: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
858: PetscQuadrature quad = NULL;
859: IS pointIS;
861: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS);
862: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
863: if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
864: if (!quad) {
865: if (!h && allPoints) {
866: quad = allPoints;
867: allPoints = NULL;
868: } else {
869: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad);
870: }
871: }
872: DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
873: for (p = pStart; p < pEnd; ++p) {
874: PetscArrayzero(values, numValues);
875: PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom);
876: DMPlexSetActivePoint(dm, p);
877: 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);
878: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);
879: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
880: }
881: PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom);
882: PetscFEGeomDestroy(&fegeom);
883: PetscQuadratureDestroy(&quad);
884: ISDestroy(&pointIS);
885: }
886: ISDestroy(&heightIS);
887: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
888: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
889: }
890: /* Cleanup */
891: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
892: for (f = 0; f < NfIn; ++f) PetscTabulationDestroy(&T[f]);
893: for (f = 0; f < NfAux; ++f) PetscTabulationDestroy(&TAux[f]);
894: PetscFree2(T, TAux);
895: }
896: PetscQuadratureDestroy(&allPoints);
897: PetscFree3(isFE, sp, spIn);
898: if (maxHeight > 0) PetscFree2(cellsp, cellspIn);
899: DMDestroy(&plex);
900: DMDestroy(&plexIn);
901: if (dmAux) DMDestroy(&plexAux);
902: return 0;
903: }
905: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
906: {
907: DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
908: return 0;
909: }
911: 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)
912: {
913: DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX);
914: return 0;
915: }
917: 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)
918: {
919: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
920: return 0;
921: }
923: 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)
924: {
925: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX);
926: return 0;
927: }
929: 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)
930: {
931: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX);
932: return 0;
933: }