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 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)
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 Arguments:
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[],
71: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
72: PetscScalar values[])
73: {
74: PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp;
75: PetscBool isAffine, isHybrid, transform;
79: DMGetCoordinateDim(dmIn, &coordDim);
80: DMHasBasisTransform(dmIn, &transform);
81: PetscDSGetNumFields(ds, &Nf);
82: PetscDSGetComponents(ds, &Nc);
83: PetscDSGetHybrid(ds, &isHybrid);
84: /* Get values for closure */
85: isAffine = fegeom->isAffine;
86: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
87: void * const ctx = ctxs ? ctxs[f] : NULL;
89: if (!sp[f]) continue;
90: PetscDualSpaceGetDimension(sp[f], &spDim);
91: if (funcs[f]) {
92: if (isFE[f]) {
93: PetscQuadrature allPoints;
94: PetscInt q, dim, numPoints;
95: const PetscReal *points;
96: PetscScalar *pointEval;
97: PetscReal *x;
98: DM rdm;
100: PetscDualSpaceGetDM(sp[f],&rdm);
101: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
102: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
103: DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
104: DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
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 (isHybrid) {
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 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D 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) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
125: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
126: }
127: /* Transform point evaluations pointEval[q,c] */
128: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
129: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
130: DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
131: DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
132: v += spDim;
133: if (isHybrid && (f < Nf-1)) {
134: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
135: }
136: } else {
137: for (d = 0; d < spDim; ++d, ++v) {
138: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
139: }
140: }
141: } else {
142: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
143: if (isHybrid && (f < Nf-1)) {
144: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
145: }
146: }
147: }
148: return(0);
149: }
151: /*
152: DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point
154: Input Parameters:
155: + dm - The output DM
156: . ds - The output DS
157: . dmIn - The input DM
158: . dsIn - The input DS
159: . dmAux - The auxiliary DM, which is always for the input space
160: . dsAux - The auxiliary DS, which is always for the input space
161: . time - The time for this evaluation
162: . localU - The local solution
163: . localA - The local auziliary fields
164: . cgeom - The FE geometry for this point
165: . sp - The output PetscDualSpace for each field
166: . p - The point in the output DM
167: . T - Input basis and derviatives for each field tabulated on the quadrature points
168: . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
169: . funcs - The evaluation function for each field
170: - ctxs - The user context for each field
172: Output Parameter:
173: . values - The value for each dual basis vector in the output dual space
175: Note: Not supported for FV
177: Level: developer
179: .seealso: DMProjectPoint_Field_Private()
180: */
181: 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,
182: PetscTabulation *T, PetscTabulation *TAux,
183: void (**funcs)(PetscInt, PetscInt, PetscInt,
184: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
185: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
186: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
187: PetscScalar values[])
188: {
189: PetscSection section, sectionAux = NULL;
190: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
191: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
192: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
193: const PetscScalar *constants;
194: PetscReal *x;
195: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
196: PetscFEGeom fegeom;
197: const PetscInt dE = cgeom->dimEmbed;
198: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
199: PetscBool isAffine, isHybrid, transform;
200: PetscErrorCode ierr;
203: PetscDSGetNumFields(ds, &Nf);
204: PetscDSGetComponents(ds, &Nc);
205: PetscDSGetHybrid(ds, &isHybrid);
206: PetscDSGetNumFields(dsIn, &NfIn);
207: PetscDSGetComponentOffsets(dsIn, &uOff);
208: PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
209: PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
210: PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
211: PetscDSGetConstants(dsIn, &numConstants, &constants);
212: DMHasBasisTransform(dmIn, &transform);
213: DMGetLocalSection(dmIn, §ion);
214: DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
215: DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
216: if (dmAux) {
217: PetscInt subp;
219: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
220: PetscDSGetNumFields(dsAux, &NfAux);
221: DMGetLocalSection(dmAux, §ionAux);
222: PetscDSGetComponentOffsets(dsAux, &aOff);
223: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
224: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
225: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
226: }
227: /* Get values for closure */
228: isAffine = cgeom->isAffine;
229: fegeom.dim = cgeom->dim;
230: fegeom.dimEmbed = cgeom->dimEmbed;
231: if (isAffine) {
232: fegeom.v = x;
233: fegeom.xi = cgeom->xi;
234: fegeom.J = cgeom->J;
235: fegeom.invJ = cgeom->invJ;
236: fegeom.detJ = cgeom->detJ;
237: }
238: for (f = 0, v = 0; f < Nf; ++f) {
239: PetscQuadrature allPoints;
240: PetscInt q, dim, numPoints;
241: const PetscReal *points;
242: PetscScalar *pointEval;
243: DM dm;
245: if (!sp[f]) continue;
246: PetscDualSpaceGetDimension(sp[f], &spDim);
247: if (!funcs[f]) {
248: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
249: if (isHybrid && (f < Nf-1)) {
250: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
251: }
252: continue;
253: }
254: PetscDualSpaceGetDM(sp[f],&dm);
255: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
256: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
257: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
258: for (q = 0; q < numPoints; ++q, ++tp) {
259: if (isAffine) {
260: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x);
261: } else {
262: fegeom.v = &cgeom->v[tp*dE];
263: fegeom.J = &cgeom->J[tp*dE*dE];
264: fegeom.invJ = &cgeom->invJ[tp*dE*dE];
265: fegeom.detJ = &cgeom->detJ[tp];
266: }
267: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
268: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
269: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);}
270: (*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]);
271: }
272: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
273: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
274: v += spDim;
275: /* TODO: For now, set both sides equal, but this should use info from other support cell */
276: if (isHybrid && (f < Nf-1)) {
277: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
278: }
279: }
280: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
281: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
282: return(0);
283: }
285: 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,
286: PetscTabulation *T, PetscTabulation *TAux,
287: void (**funcs)(PetscInt, PetscInt, PetscInt,
288: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
289: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
290: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
291: PetscScalar values[])
292: {
293: PetscSection section, sectionAux = NULL;
294: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
295: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
296: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
297: const PetscScalar *constants;
298: PetscReal *x;
299: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
300: PetscFEGeom fegeom, cgeom;
301: const PetscInt dE = fgeom->dimEmbed;
302: PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
303: PetscBool isAffine;
304: PetscErrorCode ierr;
307: if (dm != dmIn) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM");
308: PetscDSGetNumFields(ds, &Nf);
309: PetscDSGetComponents(ds, &Nc);
310: PetscDSGetComponentOffsets(ds, &uOff);
311: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
312: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
313: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
314: PetscDSGetConstants(ds, &numConstants, &constants);
315: DMGetLocalSection(dm, §ion);
316: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
317: if (dmAux) {
318: PetscInt subp;
320: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
321: PetscDSGetNumFields(dsAux, &NfAux);
322: DMGetLocalSection(dmAux, §ionAux);
323: PetscDSGetComponentOffsets(dsAux, &aOff);
324: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
325: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
326: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
327: }
328: /* Get values for closure */
329: isAffine = fgeom->isAffine;
330: fegeom.n = NULL;
331: fegeom.J = NULL;
332: fegeom.v = NULL;
333: fegeom.xi = NULL;
334: cgeom.dim = fgeom->dim;
335: cgeom.dimEmbed = fgeom->dimEmbed;
336: if (isAffine) {
337: fegeom.v = x;
338: fegeom.xi = fgeom->xi;
339: fegeom.J = fgeom->J;
340: fegeom.invJ = fgeom->invJ;
341: fegeom.detJ = fgeom->detJ;
342: fegeom.n = fgeom->n;
344: cgeom.J = fgeom->suppJ[0];
345: cgeom.invJ = fgeom->suppInvJ[0];
346: cgeom.detJ = fgeom->suppDetJ[0];
347: }
348: for (f = 0, v = 0; f < Nf; ++f) {
349: PetscQuadrature allPoints;
350: PetscInt q, dim, numPoints;
351: const PetscReal *points;
352: PetscScalar *pointEval;
353: DM dm;
355: if (!sp[f]) continue;
356: PetscDualSpaceGetDimension(sp[f], &spDim);
357: if (!funcs[f]) {
358: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
359: continue;
360: }
361: PetscDualSpaceGetDM(sp[f],&dm);
362: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
363: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
364: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
365: for (q = 0; q < numPoints; ++q, ++tp) {
366: if (isAffine) {
367: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
368: } else {
369: fegeom.v = &fgeom->v[tp*dE];
370: fegeom.J = &fgeom->J[tp*dE*dE];
371: fegeom.invJ = &fgeom->invJ[tp*dE*dE];
372: fegeom.detJ = &fgeom->detJ[tp];
373: fegeom.n = &fgeom->n[tp*dE];
375: cgeom.J = &fgeom->suppJ[0][tp*dE*dE];
376: cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE];
377: cgeom.detJ = &fgeom->suppDetJ[0][tp];
378: }
379: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
380: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
381: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
382: (*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]);
383: }
384: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
385: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
386: v += spDim;
387: }
388: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
389: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
390: return(0);
391: }
393: 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[],
394: PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux,
395: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
396: {
397: PetscFVCellGeom fvgeom;
398: PetscInt dim, dimEmbed;
399: PetscErrorCode ierr;
402: DMGetDimension(dm, &dim);
403: DMGetCoordinateDim(dm, &dimEmbed);
404: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
405: switch (type) {
406: case DM_BC_ESSENTIAL:
407: case DM_BC_NATURAL:
408: DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
409: case DM_BC_ESSENTIAL_FIELD:
410: case DM_BC_NATURAL_FIELD:
411: DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
412: (void (**)(PetscInt, PetscInt, PetscInt,
413: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
414: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
415: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
416: case DM_BC_ESSENTIAL_BD_FIELD:
417: DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
418: (void (**)(PetscInt, PetscInt, PetscInt,
419: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
420: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
421: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
422: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
423: }
424: return(0);
425: }
427: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
428: {
429: PetscReal *points;
430: PetscInt f, numPoints;
434: numPoints = 0;
435: for (f = 0; f < Nf; ++f) {
436: if (funcs[f]) {
437: PetscQuadrature fAllPoints;
438: PetscInt fNumPoints;
440: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
441: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
442: numPoints += fNumPoints;
443: }
444: }
445: PetscMalloc1(dim*numPoints,&points);
446: numPoints = 0;
447: for (f = 0; f < Nf; ++f) {
448: if (funcs[f]) {
449: PetscQuadrature fAllPoints;
450: PetscInt qdim, fNumPoints, q;
451: const PetscReal *fPoints;
453: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
454: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
455: 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);
456: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
457: numPoints += fNumPoints;
458: }
459: }
460: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
461: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
462: return(0);
463: }
465: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *ds)
466: {
467: DM plex;
468: DMEnclosureType enc;
469: DMLabel depthLabel;
470: PetscInt dim, cdepth, ls = -1, i;
471: PetscErrorCode ierr;
474: if (lStart) *lStart = -1;
475: if (!label) return(0);
476: DMGetEnclosureRelation(dm, odm, &enc);
477: DMGetDimension(dm, &dim);
478: DMConvert(dm, DMPLEX, &plex);
479: DMPlexGetDepthLabel(plex, &depthLabel);
480: cdepth = dim - height;
481: for (i = 0; i < numIds; ++i) {
482: IS pointIS;
483: const PetscInt *points;
484: PetscInt pdepth, point;
486: DMLabelGetStratumIS(label, ids[i], &pointIS);
487: if (!pointIS) continue; /* No points with that id on this process */
488: ISGetIndices(pointIS, &points);
489: DMGetEnclosurePoint(dm, odm, enc, points[0], &point);
490: DMLabelGetValue(depthLabel, point, &pdepth);
491: if (pdepth == cdepth) {
492: ls = point;
493: if (ds) {DMGetCellDS(dm, ls, ds);}
494: }
495: ISRestoreIndices(pointIS, &points);
496: ISDestroy(&pointIS);
497: if (ls >= 0) break;
498: }
499: if (lStart) *lStart = 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: The maxHeight is used to support enforcement of constraints in DMForest.
523: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
525: 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.
526: - 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
527: 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
528: dual spaces for the boundary from our input spaces.
529: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
531: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
533: 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.
534: */
535: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
536: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
537: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
538: InsertMode mode, Vec localX)
539: {
540: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
541: DMEnclosureType encIn, encAux;
542: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
543: Vec localA = NULL, tv;
544: IS fieldIS;
545: PetscSection section;
546: PetscDualSpace *sp, *cellsp;
547: PetscTabulation *T = NULL, *TAux = NULL;
548: PetscInt *Nc;
549: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
550: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, isHybrid = PETSC_FALSE, transform;
551: DMField coordField;
552: DMLabel depthLabel;
553: PetscQuadrature allPoints = NULL;
554: PetscErrorCode ierr;
557: if (localU) {VecGetDM(localU, &dmIn);}
558: else {dmIn = dm;}
559: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
560: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
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);
573: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
574: if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
575: if (!minHeight) {
576: DMLabel spmap;
578: /* If dmAux is a surface, then force the projection to take place over a surface */
579: DMPlexGetSubpointMap(plexAux, &spmap);
580: if (spmap) {
581: DMPlexGetVTKCellHeight(plexAux, &minHeight);
582: auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
583: }
584: }
585: }
586: }
587: DMPlexGetDepth(plex, &depth);
588: DMPlexGetDepthLabel(plex, &depthLabel);
589: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
590: maxHeight = PetscMax(maxHeight, minHeight);
591: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
592: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
593: if (!ds) {DMGetDS(dm, &ds);}
594: DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
595: if (!dsIn) {DMGetDS(dmIn, &dsIn);}
596: PetscDSGetNumFields(ds, &Nf);
597: PetscDSGetNumFields(dsIn, &NfIn);
598: DMGetNumFields(dm, &NfTot);
599: DMFindRegionNum(dm, ds, ®ionNum);
600: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
601: PetscDSGetHybrid(ds, &isHybrid);
602: DMGetCoordinateDim(dm, &dimEmbed);
603: DMGetLocalSection(dm, §ion);
604: if (dmAux) {
605: DMGetDS(dmAux, &dsAux);
606: PetscDSGetNumFields(dsAux, &NfAux);
607: }
608: PetscDSGetComponents(ds, &Nc);
609: PetscMalloc2(Nf, &isFE, Nf, &sp);
610: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
611: else {cellsp = sp;}
612: if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
613: /* Get cell dual spaces */
614: for (f = 0; f < Nf; ++f) {
615: PetscDiscType disctype;
617: PetscDSGetDiscType_Internal(ds, f, &disctype);
618: if (disctype == PETSC_DISC_FE) {
619: PetscFE fe;
621: isFE[f] = PETSC_TRUE;
622: hasFE = PETSC_TRUE;
623: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
624: PetscFEGetDualSpace(fe, &cellsp[f]);
625: } else if (disctype == PETSC_DISC_FV) {
626: PetscFV fv;
628: isFE[f] = PETSC_FALSE;
629: hasFV = PETSC_TRUE;
630: PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
631: PetscFVGetDualSpace(fv, &cellsp[f]);
632: } else {
633: isFE[f] = PETSC_FALSE;
634: cellsp[f] = NULL;
635: }
636: }
637: DMGetCoordinateField(dm,&coordField);
638: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
639: PetscInt effectiveHeight = auxBd ? minHeight : 0;
640: PetscFE fem, subfem;
641: PetscDiscType disctype;
642: const PetscReal *points;
643: PetscInt numPoints;
645: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
646: for (f = 0; f < Nf; ++f) {
647: if (!effectiveHeight) {sp[f] = cellsp[f];}
648: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
649: }
650: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
651: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
652: PetscMalloc2(NfIn, &T, NfAux, &TAux);
653: for (f = 0; f < NfIn; ++f) {
654: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
655: if (disctype != PETSC_DISC_FE) continue;
656: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
657: if (!effectiveHeight) {subfem = fem;}
658: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
659: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
660: }
661: for (f = 0; f < NfAux; ++f) {
662: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
663: if (disctype != PETSC_DISC_FE) continue;
664: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
665: if (!effectiveHeight || auxBd) {subfem = fem;}
666: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
667: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
668: }
669: }
670: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
671: for (h = minHeight; h <= maxHeight; h++) {
672: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
673: PetscDS dsEff = ds;
674: PetscScalar *values;
675: PetscBool *fieldActive;
676: PetscInt maxDegree;
677: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
678: IS heightIS;
680: /* Note we assume that dm and dmIn share the same topology */
681: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
682: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
683: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
684: if (pEnd <= pStart) {
685: ISDestroy(&heightIS);
686: continue;
687: }
688: /* Compute totDim, the number of dofs in the closure of a point at this height */
689: totDim = 0;
690: for (f = 0; f < Nf; ++f) {
691: if (!effectiveHeight) {
692: sp[f] = cellsp[f];
693: } else {
694: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
695: }
696: if (!sp[f]) continue;
697: PetscDualSpaceGetDimension(sp[f], &spDim);
698: totDim += spDim;
699: if (isHybrid && (f < Nf-1)) totDim += spDim;
700: }
701: DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
702: 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);
703: if (!totDim) {
704: ISDestroy(&heightIS);
705: continue;
706: }
707: if (effectiveHeight) {PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);}
708: /* Loop over points at this height */
709: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
710: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
711: {
712: const PetscInt *fields;
714: ISGetIndices(fieldIS, &fields);
715: for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
716: for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
717: ISRestoreIndices(fieldIS, &fields);
718: }
719: if (label) {
720: PetscInt i;
722: for (i = 0; i < numIds; ++i) {
723: IS pointIS, isectIS;
724: const PetscInt *points;
725: PetscInt n;
726: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
727: PetscQuadrature quad = NULL;
729: DMLabelGetStratumIS(label, ids[i], &pointIS);
730: if (!pointIS) continue; /* No points with that id on this process */
731: ISIntersect(pointIS,heightIS,&isectIS);
732: ISDestroy(&pointIS);
733: if (!isectIS) continue;
734: ISGetLocalSize(isectIS, &n);
735: ISGetIndices(isectIS, &points);
736: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
737: if (maxDegree <= 1) {
738: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
739: }
740: if (!quad) {
741: if (!h && allPoints) {
742: quad = allPoints;
743: allPoints = NULL;
744: } else {
745: PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-effectiveHeight-1 : dim-effectiveHeight,funcs,&quad);
746: }
747: }
748: DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
749: for (p = 0; p < n; ++p) {
750: const PetscInt point = points[p];
752: PetscArrayzero(values, numValues);
753: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
754: DMPlexSetActivePoint(dm, point);
755: 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);
756: if (ierr) {
757: PetscErrorCode ierr2;
758: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
759: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
760:
761: }
762: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
763: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
764: }
765: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
766: PetscFEGeomDestroy(&fegeom);
767: PetscQuadratureDestroy(&quad);
768: ISRestoreIndices(isectIS, &points);
769: ISDestroy(&isectIS);
770: }
771: } else {
772: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
773: PetscQuadrature quad = NULL;
774: IS pointIS;
776: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
777: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
778: if (maxDegree <= 1) {
779: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
780: }
781: if (!quad) {
782: if (!h && allPoints) {
783: quad = allPoints;
784: allPoints = NULL;
785: } else {
786: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&quad);
787: }
788: }
789: DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
790: for (p = pStart; p < pEnd; ++p) {
791: PetscArrayzero(values, numValues);
792: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
793: DMPlexSetActivePoint(dm, p);
794: 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);
795: if (ierr) {
796: PetscErrorCode ierr2;
797: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
798: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
799:
800: }
801: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
802: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
803: }
804: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
805: PetscFEGeomDestroy(&fegeom);
806: PetscQuadratureDestroy(&quad);
807: ISDestroy(&pointIS);
808: }
809: ISDestroy(&heightIS);
810: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
811: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
812: }
813: /* Cleanup */
814: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
815: PetscInt effectiveHeight = auxBd ? minHeight : 0;
816: PetscFE fem, subfem;
818: for (f = 0; f < NfIn; ++f) {
819: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
820: if (!effectiveHeight) {subfem = fem;}
821: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
822: PetscTabulationDestroy(&T[f]);
823: }
824: for (f = 0; f < NfAux; ++f) {
825: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
826: if (!effectiveHeight || auxBd) {subfem = fem;}
827: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
828: PetscTabulationDestroy(&TAux[f]);
829: }
830: PetscFree2(T, TAux);
831: }
832: PetscQuadratureDestroy(&allPoints);
833: PetscFree2(isFE, sp);
834: if (maxHeight > 0) {PetscFree(cellsp);}
835: DMDestroy(&plex);
836: DMDestroy(&plexIn);
837: if (dmAux) {DMDestroy(&plexAux);}
838: return(0);
839: }
841: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
842: {
846: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
847: return(0);
848: }
850: 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)
851: {
855: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
856: return(0);
857: }
859: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
860: void (**funcs)(PetscInt, PetscInt, PetscInt,
861: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
862: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
863: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
864: InsertMode mode, Vec localX)
865: {
869: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
870: return(0);
871: }
873: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
874: void (**funcs)(PetscInt, PetscInt, PetscInt,
875: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
876: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
877: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
878: InsertMode mode, Vec localX)
879: {
883: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
884: return(0);
885: }
887: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
888: void (**funcs)(PetscInt, PetscInt, PetscInt,
889: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
890: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
891: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
892: InsertMode mode, Vec localX)
893: {
897: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
898: return(0);
899: }