Actual source code: plexproject.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: /*@
6: DMPlexGetActivePoint - Get the point on which projection is currently working
8: Not collective
10: Input Argument:
11: . dm - the DM
13: Output Argument:
14: . point - The mesh point involved in the current projection
16: Level: developer
18: .seealso: DMPlexSetActivePoint()
19: @*/
20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point) {
22: *point = ((DM_Plex *) dm->data)->activePoint;
23: return(0);
24: }
26: /*@
27: DMPlexSetActivePoint - Set the point on which projection is currently working
29: Not collective
31: Input Arguments:
32: + dm - the DM
33: - point - The mesh point involved in the current projection
35: Level: developer
37: .seealso: DMPlexGetActivePoint()
38: @*/
39: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point) {
41: ((DM_Plex *) dm->data)->activePoint = point;
42: return(0);
43: }
45: /*
46: DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point
48: Input Parameters:
49: + dm - The output DM
50: . ds - The output DS
51: . dmIn - The input DM
52: . dsIn - The input DS
53: . time - The time for this evaluation
54: . fegeom - The FE geometry for this point
55: . fvgeom - The FV geometry for this point
56: . isFE - Flag indicating whether each output field has an FE discretization
57: . sp - The output PetscDualSpace for each field
58: . funcs - The evaluation function for each field
59: - ctxs - The user context for each field
61: Output Parameter:
62: . values - The value for each dual basis vector in the output dual space
64: Level: developer
66: .seealso: DMProjectPoint_Field_Private()
67: */
68: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[],
69: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
70: PetscScalar values[])
71: {
72: PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp;
73: PetscBool isAffine, isHybrid, transform;
77: DMGetCoordinateDim(dmIn, &coordDim);
78: DMHasBasisTransform(dmIn, &transform);
79: PetscDSGetNumFields(ds, &Nf);
80: PetscDSGetComponents(ds, &Nc);
81: PetscDSGetHybrid(ds, &isHybrid);
82: /* Get values for closure */
83: isAffine = fegeom->isAffine;
84: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
85: void * const ctx = ctxs ? ctxs[f] : NULL;
87: if (!sp[f]) continue;
88: PetscDualSpaceGetDimension(sp[f], &spDim);
89: if (funcs[f]) {
90: if (isFE[f]) {
91: PetscQuadrature allPoints;
92: PetscInt q, dim, numPoints;
93: const PetscReal *points;
94: PetscScalar *pointEval;
95: PetscReal *x;
96: DM rdm;
98: PetscDualSpaceGetDM(sp[f],&rdm);
99: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
100: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
101: DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
102: DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
103: for (q = 0; q < numPoints; q++, tp++) {
104: const PetscReal *v0;
106: if (isAffine) {
107: const PetscReal *refpoint = &points[q*dim];
108: PetscReal injpoint[3] = {0., 0., 0.};
110: if (dim != fegeom->dim) {
111: if (isHybrid) {
112: /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
113: for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
114: refpoint = injpoint;
115: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D dual basis spatial dimension", fegeom->dim, dim);
116: }
117: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
118: v0 = x;
119: } else {
120: v0 = &fegeom->v[tp*coordDim];
121: }
122: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
123: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
124: }
125: /* Transform point evaluations pointEval[q,c] */
126: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
127: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
128: DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
129: DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
130: v += spDim;
131: if (isHybrid && (f < Nf-1)) {
132: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
133: }
134: } else {
135: for (d = 0; d < spDim; ++d, ++v) {
136: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
137: }
138: }
139: } else {
140: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
141: if (isHybrid && (f < Nf-1)) {
142: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
143: }
144: }
145: }
146: return(0);
147: }
149: /*
150: DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point
152: Input Parameters:
153: + dm - The output DM
154: . ds - The output DS
155: . dmIn - The input DM
156: . dsIn - The input DS
157: . dmAux - The auxiliary DM, which is always for the input space
158: . dsAux - The auxiliary DS, which is always for the input space
159: . time - The time for this evaluation
160: . localU - The local solution
161: . localA - The local auziliary fields
162: . cgeom - The FE geometry for this point
163: . sp - The output PetscDualSpace for each field
164: . p - The point in the output DM
165: . T - Input basis and derviatives for each field tabulated on the quadrature points
166: . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
167: . funcs - The evaluation function for each field
168: - ctxs - The user context for each field
170: Output Parameter:
171: . values - The value for each dual basis vector in the output dual space
173: Note: Not supported for FV
175: Level: developer
177: .seealso: DMProjectPoint_Field_Private()
178: */
179: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p,
180: PetscTabulation *T, PetscTabulation *TAux,
181: void (**funcs)(PetscInt, PetscInt, PetscInt,
182: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
183: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
184: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
185: PetscScalar values[])
186: {
187: PetscSection section, sectionAux = NULL;
188: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
189: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
190: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
191: const PetscScalar *constants;
192: PetscReal *x;
193: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
194: PetscFEGeom fegeom;
195: const PetscInt dE = cgeom->dimEmbed;
196: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
197: PetscBool isAffine, isHybrid, transform;
198: PetscErrorCode ierr;
201: PetscDSGetNumFields(ds, &Nf);
202: PetscDSGetComponents(ds, &Nc);
203: PetscDSGetHybrid(ds, &isHybrid);
204: PetscDSGetNumFields(dsIn, &NfIn);
205: PetscDSGetComponentOffsets(dsIn, &uOff);
206: PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
207: PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
208: PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
209: PetscDSGetConstants(dsIn, &numConstants, &constants);
210: DMHasBasisTransform(dmIn, &transform);
211: DMGetLocalSection(dmIn, §ion);
212: DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
213: DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
214: if (dmAux) {
215: PetscInt subp;
217: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
218: PetscDSGetNumFields(dsAux, &NfAux);
219: DMGetLocalSection(dmAux, §ionAux);
220: PetscDSGetComponentOffsets(dsAux, &aOff);
221: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
222: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
223: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
224: }
225: /* Get values for closure */
226: isAffine = cgeom->isAffine;
227: fegeom.dim = cgeom->dim;
228: fegeom.dimEmbed = cgeom->dimEmbed;
229: if (isAffine) {
230: fegeom.v = x;
231: fegeom.xi = cgeom->xi;
232: fegeom.J = cgeom->J;
233: fegeom.invJ = cgeom->invJ;
234: fegeom.detJ = cgeom->detJ;
235: }
236: for (f = 0, v = 0; f < Nf; ++f) {
237: PetscQuadrature allPoints;
238: PetscInt q, dim, numPoints;
239: const PetscReal *points;
240: PetscScalar *pointEval;
241: DM dm;
243: if (!sp[f]) continue;
244: PetscDualSpaceGetDimension(sp[f], &spDim);
245: if (!funcs[f]) {
246: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
247: if (isHybrid && (f < Nf-1)) {
248: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
249: }
250: continue;
251: }
252: PetscDualSpaceGetDM(sp[f],&dm);
253: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
254: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
255: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
256: for (q = 0; q < numPoints; ++q, ++tp) {
257: if (isAffine) {
258: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x);
259: } else {
260: fegeom.v = &cgeom->v[tp*dE];
261: fegeom.J = &cgeom->J[tp*dE*dE];
262: fegeom.invJ = &cgeom->invJ[tp*dE*dE];
263: fegeom.detJ = &cgeom->detJ[tp];
264: }
265: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
266: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
267: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);}
268: (*funcs[f])(dE, NfIn, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, fegeom.v, numConstants, constants, &pointEval[Nc[f]*q]);
269: }
270: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
271: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
272: v += spDim;
273: /* TODO: For now, set both sides equal, but this should use info from other support cell */
274: if (isHybrid && (f < Nf-1)) {
275: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
276: }
277: }
278: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
279: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
280: return(0);
281: }
283: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p,
284: PetscTabulation *T, PetscTabulation *TAux,
285: void (**funcs)(PetscInt, PetscInt, PetscInt,
286: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
287: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
288: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
289: PetscScalar values[])
290: {
291: PetscSection section, sectionAux = NULL;
292: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
293: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
294: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
295: const PetscScalar *constants;
296: PetscReal *x;
297: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
298: PetscFEGeom fegeom, cgeom;
299: const PetscInt dE = fgeom->dimEmbed;
300: PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
301: PetscBool isAffine;
302: PetscErrorCode ierr;
305: if (dm != dmIn) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM");
306: PetscDSGetNumFields(ds, &Nf);
307: PetscDSGetComponents(ds, &Nc);
308: PetscDSGetComponentOffsets(ds, &uOff);
309: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
310: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
311: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
312: PetscDSGetConstants(ds, &numConstants, &constants);
313: DMGetLocalSection(dm, §ion);
314: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
315: if (dmAux) {
316: PetscInt subp;
318: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
319: PetscDSGetNumFields(dsAux, &NfAux);
320: DMGetLocalSection(dmAux, §ionAux);
321: PetscDSGetComponentOffsets(dsAux, &aOff);
322: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
323: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
324: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
325: }
326: /* Get values for closure */
327: isAffine = fgeom->isAffine;
328: fegeom.n = NULL;
329: fegeom.J = NULL;
330: fegeom.v = NULL;
331: fegeom.xi = NULL;
332: cgeom.dim = fgeom->dim;
333: cgeom.dimEmbed = fgeom->dimEmbed;
334: if (isAffine) {
335: fegeom.v = x;
336: fegeom.xi = fgeom->xi;
337: fegeom.J = fgeom->J;
338: fegeom.invJ = fgeom->invJ;
339: fegeom.detJ = fgeom->detJ;
340: fegeom.n = fgeom->n;
342: cgeom.J = fgeom->suppJ[0];
343: cgeom.invJ = fgeom->suppInvJ[0];
344: cgeom.detJ = fgeom->suppDetJ[0];
345: }
346: for (f = 0, v = 0; f < Nf; ++f) {
347: PetscQuadrature allPoints;
348: PetscInt q, dim, numPoints;
349: const PetscReal *points;
350: PetscScalar *pointEval;
351: DM dm;
353: if (!sp[f]) continue;
354: PetscDualSpaceGetDimension(sp[f], &spDim);
355: if (!funcs[f]) {
356: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
357: continue;
358: }
359: PetscDualSpaceGetDM(sp[f],&dm);
360: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
361: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
362: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
363: for (q = 0; q < numPoints; ++q, ++tp) {
364: if (isAffine) {
365: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
366: } else {
367: fegeom.v = &fgeom->v[tp*dE];
368: fegeom.J = &fgeom->J[tp*dE*dE];
369: fegeom.invJ = &fgeom->invJ[tp*dE*dE];
370: fegeom.detJ = &fgeom->detJ[tp];
371: fegeom.n = &fgeom->n[tp*dE];
373: cgeom.J = &fgeom->suppJ[0][tp*dE*dE];
374: cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE];
375: cgeom.detJ = &fgeom->suppDetJ[0][tp];
376: }
377: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
378: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
379: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
380: (*funcs[f])(dE, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, fegeom.v, fegeom.n, numConstants, constants, &pointEval[Nc[f]*q]);
381: }
382: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
383: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
384: v += spDim;
385: }
386: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
387: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
388: return(0);
389: }
391: static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscFEGeom *fegeom, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[],
392: PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux,
393: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
394: {
395: PetscFVCellGeom fvgeom;
396: PetscInt dim, dimEmbed;
397: PetscErrorCode ierr;
400: DMGetDimension(dm, &dim);
401: DMGetCoordinateDim(dm, &dimEmbed);
402: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
403: switch (type) {
404: case DM_BC_ESSENTIAL:
405: case DM_BC_NATURAL:
406: DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
407: case DM_BC_ESSENTIAL_FIELD:
408: case DM_BC_NATURAL_FIELD:
409: DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
410: (void (**)(PetscInt, PetscInt, PetscInt,
411: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
412: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
413: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
414: case DM_BC_ESSENTIAL_BD_FIELD:
415: DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
416: (void (**)(PetscInt, PetscInt, PetscInt,
417: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
418: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
419: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
420: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
421: }
422: return(0);
423: }
425: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
426: {
427: PetscReal *points;
428: PetscInt f, numPoints;
432: numPoints = 0;
433: for (f = 0; f < Nf; ++f) {
434: if (funcs[f]) {
435: PetscQuadrature fAllPoints;
436: PetscInt fNumPoints;
438: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
439: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
440: numPoints += fNumPoints;
441: }
442: }
443: PetscMalloc1(dim*numPoints,&points);
444: numPoints = 0;
445: for (f = 0; f < Nf; ++f) {
446: if (funcs[f]) {
447: PetscQuadrature fAllPoints;
448: PetscInt qdim, fNumPoints, q;
449: const PetscReal *fPoints;
451: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
452: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
453: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %D for dual basis does not match input dimension %D", qdim, dim);
454: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
455: numPoints += fNumPoints;
456: }
457: }
458: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
459: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
460: return(0);
461: }
463: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *ds)
464: {
465: DM plex;
466: DMEnclosureType enc;
467: DMLabel depthLabel;
468: PetscInt dim, cdepth, ls = -1, i;
469: PetscErrorCode ierr;
472: if (lStart) *lStart = -1;
473: if (!label) return(0);
474: DMGetEnclosureRelation(dm, odm, &enc);
475: DMGetDimension(dm, &dim);
476: DMConvert(dm, DMPLEX, &plex);
477: DMPlexGetDepthLabel(plex, &depthLabel);
478: cdepth = dim - height;
479: for (i = 0; i < numIds; ++i) {
480: IS pointIS;
481: const PetscInt *points;
482: PetscInt pdepth, point;
484: DMLabelGetStratumIS(label, ids[i], &pointIS);
485: if (!pointIS) continue; /* No points with that id on this process */
486: ISGetIndices(pointIS, &points);
487: DMGetEnclosurePoint(dm, odm, enc, points[0], &point);
488: DMLabelGetValue(depthLabel, point, &pdepth);
489: if (pdepth == cdepth) {
490: ls = point;
491: if (ds) {DMGetCellDS(dm, ls, ds);}
492: }
493: ISRestoreIndices(pointIS, &points);
494: ISDestroy(&pointIS);
495: if (ls >= 0) break;
496: }
497: if (lStart) *lStart = ls;
498: DMDestroy(&plex);
499: return(0);
500: }
502: /*
503: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
505: There are several different scenarios:
507: 1) Volumetric mesh with volumetric auxiliary data
509: Here minHeight=0 since we loop over cells.
511: 2) Boundary mesh with boundary auxiliary data
513: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
515: 3) Volumetric mesh with boundary auxiliary data
517: Here minHeight=1 and auxbd=PETSC_TRUE since we loop over faces and use data only supported on those faces. This is common when imposing Dirichlet boundary conditions.
519: The maxHeight is used to support enforcement of constraints in DMForest.
521: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
523: If we are using an input field (DM_BC_ESSENTIAL_FIELD or DM_BC_NATURAL_FIELD), we need to evaluate it at all the quadrature points of the dual basis functionals.
524: - We use effectiveHeight to mean the height above our incoming DS. For example, if the DS is for a submesh then the effective height is zero, whereas if the DS
525: is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When the effective height is nonzero, we need to extract
526: dual spaces for the boundary from our input spaces.
527: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
529: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
531: If we have a label, we iterate over those points. This will probably break the maxHeight functionality since we do not check the height of those points.
532: */
533: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
534: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
535: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
536: InsertMode mode, Vec localX)
537: {
538: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
539: DMEnclosureType encIn, encAux;
540: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
541: Vec localA = NULL, tv;
542: IS fieldIS;
543: PetscSection section;
544: PetscDualSpace *sp, *cellsp;
545: PetscTabulation *T = NULL, *TAux = NULL;
546: PetscInt *Nc;
547: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
548: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, isHybrid = PETSC_FALSE, transform;
549: DMField coordField;
550: DMLabel depthLabel;
551: PetscQuadrature allPoints = NULL;
552: PetscErrorCode ierr;
555: if (localU) {VecGetDM(localU, &dmIn);}
556: else {dmIn = dm;}
557: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
558: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
559: DMConvert(dm, DMPLEX, &plex);
560: DMConvert(dmIn, DMPLEX, &plexIn);
561: DMGetEnclosureRelation(dmIn, dm, &encIn);
562: DMGetEnclosureRelation(dmAux, dm, &encAux);
563: DMGetDimension(dm, &dim);
564: DMPlexGetVTKCellHeight(plex, &minHeight);
565: DMGetBasisTransformDM_Internal(dm, &tdm);
566: DMGetBasisTransformVec_Internal(dm, &tv);
567: DMHasBasisTransform(dm, &transform);
568: /* Auxiliary information can only be used with interpolation of field functions */
569: if (dmAux) {
570: DMConvert(dmAux, DMPLEX, &plexAux);
571: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
572: if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
573: if (!minHeight) {
574: DMLabel spmap;
576: /* If dmAux is a surface, then force the projection to take place over a surface */
577: DMPlexGetSubpointMap(plexAux, &spmap);
578: if (spmap) {
579: DMPlexGetVTKCellHeight(plexAux, &minHeight);
580: auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
581: }
582: }
583: }
584: }
585: DMPlexGetDepth(plex, &depth);
586: DMPlexGetDepthLabel(plex, &depthLabel);
587: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
588: maxHeight = PetscMax(maxHeight, minHeight);
589: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
590: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
591: if (!ds) {DMGetDS(dm, &ds);}
592: DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
593: if (!dsIn) {DMGetDS(dmIn, &dsIn);}
594: PetscDSGetNumFields(ds, &Nf);
595: PetscDSGetNumFields(dsIn, &NfIn);
596: DMGetNumFields(dm, &NfTot);
597: DMFindRegionNum(dm, ds, ®ionNum);
598: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
599: PetscDSGetHybrid(ds, &isHybrid);
600: DMGetCoordinateDim(dm, &dimEmbed);
601: DMGetLocalSection(dm, §ion);
602: if (dmAux) {
603: DMGetDS(dmAux, &dsAux);
604: PetscDSGetNumFields(dsAux, &NfAux);
605: }
606: PetscDSGetComponents(ds, &Nc);
607: PetscMalloc2(Nf, &isFE, Nf, &sp);
608: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
609: else {cellsp = sp;}
610: if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
611: /* Get cell dual spaces */
612: for (f = 0; f < Nf; ++f) {
613: PetscDiscType disctype;
615: PetscDSGetDiscType_Internal(ds, f, &disctype);
616: if (disctype == PETSC_DISC_FE) {
617: PetscFE fe;
619: isFE[f] = PETSC_TRUE;
620: hasFE = PETSC_TRUE;
621: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
622: PetscFEGetDualSpace(fe, &cellsp[f]);
623: } else if (disctype == PETSC_DISC_FV) {
624: PetscFV fv;
626: isFE[f] = PETSC_FALSE;
627: hasFV = PETSC_TRUE;
628: PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
629: PetscFVGetDualSpace(fv, &cellsp[f]);
630: } else {
631: isFE[f] = PETSC_FALSE;
632: cellsp[f] = NULL;
633: }
634: }
635: DMGetCoordinateField(dm,&coordField);
636: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
637: PetscInt effectiveHeight = auxBd ? minHeight : 0;
638: PetscFE fem, subfem;
639: PetscDiscType disctype;
640: const PetscReal *points;
641: PetscInt numPoints;
643: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
644: for (f = 0; f < Nf; ++f) {
645: if (!effectiveHeight) {sp[f] = cellsp[f];}
646: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
647: }
648: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
649: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
650: PetscMalloc2(NfIn, &T, NfAux, &TAux);
651: for (f = 0; f < NfIn; ++f) {
652: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
653: if (disctype != PETSC_DISC_FE) continue;
654: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
655: if (!effectiveHeight) {subfem = fem;}
656: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
657: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
658: }
659: for (f = 0; f < NfAux; ++f) {
660: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
661: if (disctype != PETSC_DISC_FE) continue;
662: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
663: if (!effectiveHeight || auxBd) {subfem = fem;}
664: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
665: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
666: }
667: }
668: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
669: for (h = minHeight; h <= maxHeight; h++) {
670: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
671: PetscDS dsEff = ds;
672: PetscScalar *values;
673: PetscBool *fieldActive;
674: PetscInt maxDegree;
675: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
676: IS heightIS;
678: /* Note we assume that dm and dmIn share the same topology */
679: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
680: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
681: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
682: if (pEnd <= pStart) {
683: ISDestroy(&heightIS);
684: continue;
685: }
686: /* Compute totDim, the number of dofs in the closure of a point at this height */
687: totDim = 0;
688: for (f = 0; f < Nf; ++f) {
689: if (!effectiveHeight) {
690: sp[f] = cellsp[f];
691: } else {
692: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
693: }
694: if (!sp[f]) continue;
695: PetscDualSpaceGetDimension(sp[f], &spDim);
696: totDim += spDim;
697: if (isHybrid && (f < Nf-1)) totDim += spDim;
698: }
699: DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
700: if (numValues != totDim) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point (%D) closure size %D != dual space dimension %D", lStart < 0 ? pStart : lStart, numValues, totDim);
701: if (!totDim) {
702: ISDestroy(&heightIS);
703: continue;
704: }
705: if (effectiveHeight) {PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);}
706: /* Loop over points at this height */
707: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
708: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
709: {
710: const PetscInt *fields;
712: ISGetIndices(fieldIS, &fields);
713: for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
714: for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
715: ISRestoreIndices(fieldIS, &fields);
716: }
717: if (label) {
718: PetscInt i;
720: for (i = 0; i < numIds; ++i) {
721: IS pointIS, isectIS;
722: const PetscInt *points;
723: PetscInt n;
724: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
725: PetscQuadrature quad = NULL;
727: DMLabelGetStratumIS(label, ids[i], &pointIS);
728: if (!pointIS) continue; /* No points with that id on this process */
729: ISIntersect(pointIS,heightIS,&isectIS);
730: ISDestroy(&pointIS);
731: if (!isectIS) continue;
732: ISGetLocalSize(isectIS, &n);
733: ISGetIndices(isectIS, &points);
734: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
735: if (maxDegree <= 1) {
736: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
737: }
738: if (!quad) {
739: if (!h && allPoints) {
740: quad = allPoints;
741: allPoints = NULL;
742: } else {
743: PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-effectiveHeight-1 : dim-effectiveHeight,funcs,&quad);
744: }
745: }
746: DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
747: for (p = 0; p < n; ++p) {
748: const PetscInt point = points[p];
750: PetscArrayzero(values, numValues);
751: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
752: DMPlexSetActivePoint(dm, point);
753: DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values);
754: if (ierr) {
755: PetscErrorCode ierr2;
756: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
757: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
758:
759: }
760: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
761: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
762: }
763: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
764: PetscFEGeomDestroy(&fegeom);
765: PetscQuadratureDestroy(&quad);
766: ISRestoreIndices(isectIS, &points);
767: ISDestroy(&isectIS);
768: }
769: } else {
770: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
771: PetscQuadrature quad = NULL;
772: IS pointIS;
774: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
775: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
776: if (maxDegree <= 1) {
777: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
778: }
779: if (!quad) {
780: if (!h && allPoints) {
781: quad = allPoints;
782: allPoints = NULL;
783: } else {
784: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&quad);
785: }
786: }
787: DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
788: for (p = pStart; p < pEnd; ++p) {
789: PetscArrayzero(values, numValues);
790: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
791: DMPlexSetActivePoint(dm, p);
792: DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, T, TAux, type, funcs, ctxs, fieldActive, values);
793: if (ierr) {
794: PetscErrorCode ierr2;
795: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
796: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
797:
798: }
799: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
800: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
801: }
802: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
803: PetscFEGeomDestroy(&fegeom);
804: PetscQuadratureDestroy(&quad);
805: ISDestroy(&pointIS);
806: }
807: ISDestroy(&heightIS);
808: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
809: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
810: }
811: /* Cleanup */
812: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
813: PetscInt effectiveHeight = auxBd ? minHeight : 0;
814: PetscFE fem, subfem;
816: for (f = 0; f < NfIn; ++f) {
817: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
818: if (!effectiveHeight) {subfem = fem;}
819: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
820: PetscTabulationDestroy(&T[f]);
821: }
822: for (f = 0; f < NfAux; ++f) {
823: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
824: if (!effectiveHeight || auxBd) {subfem = fem;}
825: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
826: PetscTabulationDestroy(&TAux[f]);
827: }
828: PetscFree2(T, TAux);
829: }
830: PetscQuadratureDestroy(&allPoints);
831: PetscFree2(isFE, sp);
832: if (maxHeight > 0) {PetscFree(cellsp);}
833: DMDestroy(&plex);
834: DMDestroy(&plexIn);
835: if (dmAux) {DMDestroy(&plexAux);}
836: return(0);
837: }
839: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
840: {
844: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
845: return(0);
846: }
848: PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
849: {
853: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
854: return(0);
855: }
857: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
858: void (**funcs)(PetscInt, PetscInt, PetscInt,
859: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
860: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
861: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
862: InsertMode mode, Vec localX)
863: {
867: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
868: return(0);
869: }
871: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
872: void (**funcs)(PetscInt, PetscInt, PetscInt,
873: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
874: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
875: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
876: InsertMode mode, Vec localX)
877: {
881: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
882: return(0);
883: }
885: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
886: void (**funcs)(PetscInt, PetscInt, PetscInt,
887: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
888: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
889: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
890: InsertMode mode, Vec localX)
891: {
895: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
896: return(0);
897: }