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[],
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, isCohesive, transform;
78: DMGetCoordinateDim(dmIn, &coordDim);
79: DMHasBasisTransform(dmIn, &transform);
80: PetscDSGetNumFields(ds, &Nf);
81: PetscDSGetComponents(ds, &Nc);
82: PetscDSIsCohesive(ds, &isCohesive);
83: /* Get values for closure */
84: isAffine = fegeom->isAffine;
85: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
86: void * const ctx = ctxs ? ctxs[f] : NULL;
87: PetscBool cohesive;
89: if (!sp[f]) continue;
90: PetscDSGetCohesive(ds, f, &cohesive);
91: PetscDualSpaceGetDimension(sp[f], &spDim);
92: if (funcs[f]) {
93: if (isFE[f]) {
94: PetscQuadrature allPoints;
95: PetscInt q, dim, numPoints;
96: const PetscReal *points;
97: PetscScalar *pointEval;
98: PetscReal *x;
99: DM rdm;
101: PetscDualSpaceGetDM(sp[f],&rdm);
102: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
103: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
104: DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
105: DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
106: for (q = 0; q < numPoints; q++, tp++) {
107: const PetscReal *v0;
109: if (isAffine) {
110: const PetscReal *refpoint = &points[q*dim];
111: PetscReal injpoint[3] = {0., 0., 0.};
113: if (dim != fegeom->dim) {
114: if (isCohesive) {
115: /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
116: for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
117: refpoint = injpoint;
118: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D dual basis spatial dimension", fegeom->dim, dim);
119: }
120: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
121: v0 = x;
122: } else {
123: v0 = &fegeom->v[tp*coordDim];
124: }
125: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
126: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
127: }
128: /* Transform point evaluations pointEval[q,c] */
129: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
130: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
131: DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
132: DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
133: v += spDim;
134: if (isCohesive && !cohesive) {
135: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
136: }
137: } else {
138: for (d = 0; d < spDim; ++d, ++v) {
139: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
140: }
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,
183: PetscTabulation *T, PetscTabulation *TAux,
184: void (**funcs)(PetscInt, PetscInt, PetscInt,
185: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
186: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
187: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
188: PetscScalar values[])
189: {
190: PetscSection section, sectionAux = NULL;
191: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
192: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
193: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
194: const PetscScalar *constants;
195: PetscReal *x;
196: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
197: PetscFEGeom fegeom;
198: const PetscInt dE = cgeom->dimEmbed;
199: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
200: PetscBool isAffine, isCohesive, transform;
203: PetscDSGetNumFields(ds, &Nf);
204: PetscDSGetComponents(ds, &Nc);
205: PetscDSIsCohesive(ds, &isCohesive);
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: PetscBool cohesive;
244: DM dm;
246: if (!sp[f]) continue;
247: PetscDSGetCohesive(ds, f, &cohesive);
248: PetscDualSpaceGetDimension(sp[f], &spDim);
249: if (!funcs[f]) {
250: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
251: if (isCohesive && !cohesive) {
252: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
253: }
254: continue;
255: }
256: PetscDualSpaceGetDM(sp[f],&dm);
257: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
258: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
259: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
260: for (q = 0; q < numPoints; ++q, ++tp) {
261: if (isAffine) {
262: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x);
263: } else {
264: fegeom.v = &cgeom->v[tp*dE];
265: fegeom.J = &cgeom->J[tp*dE*dE];
266: fegeom.invJ = &cgeom->invJ[tp*dE*dE];
267: fegeom.detJ = &cgeom->detJ[tp];
268: }
269: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
270: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
271: if (transform) DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);
272: (*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]);
273: }
274: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
275: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
276: v += spDim;
277: /* TODO: For now, set both sides equal, but this should use info from other support cell */
278: if (isCohesive && !cohesive) {
279: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
280: }
281: }
282: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
283: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
284: return 0;
285: }
287: 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,
288: PetscTabulation *T, PetscTabulation *TAux,
289: void (**funcs)(PetscInt, PetscInt, PetscInt,
290: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
291: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
292: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
293: PetscScalar values[])
294: {
295: PetscSection section, sectionAux = NULL;
296: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
297: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
298: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
299: const PetscScalar *constants;
300: PetscReal *x;
301: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
302: PetscFEGeom fegeom, cgeom;
303: const PetscInt dE = fgeom->dimEmbed;
304: PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
305: PetscBool isAffine;
309: PetscDSGetNumFields(ds, &Nf);
310: PetscDSGetComponents(ds, &Nc);
311: PetscDSGetComponentOffsets(ds, &uOff);
312: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
313: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
314: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
315: PetscDSGetConstants(ds, &numConstants, &constants);
316: DMGetLocalSection(dm, §ion);
317: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
318: if (dmAux) {
319: PetscInt subp;
321: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
322: PetscDSGetNumFields(dsAux, &NfAux);
323: DMGetLocalSection(dmAux, §ionAux);
324: PetscDSGetComponentOffsets(dsAux, &aOff);
325: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
326: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
327: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
328: }
329: /* Get values for closure */
330: isAffine = fgeom->isAffine;
331: fegeom.n = NULL;
332: fegeom.J = NULL;
333: fegeom.v = NULL;
334: fegeom.xi = NULL;
335: cgeom.dim = fgeom->dim;
336: cgeom.dimEmbed = fgeom->dimEmbed;
337: if (isAffine) {
338: fegeom.v = x;
339: fegeom.xi = fgeom->xi;
340: fegeom.J = fgeom->J;
341: fegeom.invJ = fgeom->invJ;
342: fegeom.detJ = fgeom->detJ;
343: fegeom.n = fgeom->n;
345: cgeom.J = fgeom->suppJ[0];
346: cgeom.invJ = fgeom->suppInvJ[0];
347: cgeom.detJ = fgeom->suppDetJ[0];
348: }
349: for (f = 0, v = 0; f < Nf; ++f) {
350: PetscQuadrature allPoints;
351: PetscInt q, dim, numPoints;
352: const PetscReal *points;
353: PetscScalar *pointEval;
354: DM dm;
356: if (!sp[f]) continue;
357: PetscDualSpaceGetDimension(sp[f], &spDim);
358: if (!funcs[f]) {
359: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
360: continue;
361: }
362: PetscDualSpaceGetDM(sp[f],&dm);
363: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
364: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
365: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
366: for (q = 0; q < numPoints; ++q, ++tp) {
367: if (isAffine) {
368: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
369: } else {
370: fegeom.v = &fgeom->v[tp*dE];
371: fegeom.J = &fgeom->J[tp*dE*dE];
372: fegeom.invJ = &fgeom->invJ[tp*dE*dE];
373: fegeom.detJ = &fgeom->detJ[tp];
374: fegeom.n = &fgeom->n[tp*dE];
376: cgeom.J = &fgeom->suppJ[0][tp*dE*dE];
377: cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE];
378: cgeom.detJ = &fgeom->suppDetJ[0][tp];
379: }
380: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
381: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
382: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
383: (*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]);
384: }
385: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
386: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
387: v += spDim;
388: }
389: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
390: if (dmAux) DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);
391: return 0;
392: }
394: 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[],
395: PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux,
396: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
397: {
398: PetscFVCellGeom fvgeom;
399: PetscInt dim, dimEmbed;
400: PetscErrorCode ierr;
403: DMGetDimension(dm, &dim);
404: DMGetCoordinateDim(dm, &dimEmbed);
405: if (hasFV) DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);
406: switch (type) {
407: case DM_BC_ESSENTIAL:
408: case DM_BC_NATURAL:
409: DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
410: case DM_BC_ESSENTIAL_FIELD:
411: case DM_BC_NATURAL_FIELD:
412: DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
413: (void (**)(PetscInt, PetscInt, PetscInt,
414: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
415: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
416: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
417: case DM_BC_ESSENTIAL_BD_FIELD:
418: DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
419: (void (**)(PetscInt, PetscInt, PetscInt,
420: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
421: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
422: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
423: default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
424: }
425: return 0;
426: }
428: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
429: {
430: PetscReal *points;
431: PetscInt f, numPoints;
433: numPoints = 0;
434: for (f = 0; f < Nf; ++f) {
435: if (funcs[f]) {
436: PetscQuadrature fAllPoints;
437: PetscInt fNumPoints;
439: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
440: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
441: numPoints += fNumPoints;
442: }
443: }
444: PetscMalloc1(dim*numPoints,&points);
445: numPoints = 0;
446: for (f = 0; f < Nf; ++f) {
447: if (funcs[f]) {
448: PetscQuadrature fAllPoints;
449: PetscInt qdim, fNumPoints, q;
450: const PetscReal *fPoints;
452: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
453: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
455: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
456: numPoints += fNumPoints;
457: }
458: }
459: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
460: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
461: return 0;
462: }
464: /*@C
465: 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.
467: Input Parameters:
468: dm - the DM
469: odm - the enclosing DM
470: label - label for DM domain, or NULL for whole domain
471: numIds - the number of ids
472: ids - An array of the label ids in sequence for the domain
473: height - Height of target cells in DMPlex topology
475: Output Parameters:
476: point - the first labeled point
477: ds - the ds corresponding to the first labeled point
479: Level: developer
480: @*/
481: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
482: {
483: DM plex;
484: DMEnclosureType enc;
485: PetscInt ls = -1;
487: if (point) *point = -1;
488: if (!label) return 0;
489: DMGetEnclosureRelation(dm, odm, &enc);
490: DMConvert(dm, DMPLEX, &plex);
491: for (PetscInt i = 0; i < numIds; ++i) {
492: IS labelIS;
493: PetscInt num_points, pStart, pEnd;
494: DMLabelGetStratumIS(label, ids[i], &labelIS);
495: if (!labelIS) continue; /* No points with that id on this process */
496: DMPlexGetHeightStratum(plex, height, &pStart, &pEnd);
497: ISGetSize(labelIS, &num_points);
498: if (num_points) {
499: const PetscInt *points;
500: ISGetIndices(labelIS, &points);
501: for (PetscInt i=0; i<num_points; i++) {
502: PetscInt point;
503: DMGetEnclosurePoint(dm, odm, enc, points[i], &point);
504: if (pStart <= point && point < pEnd) {
505: ls = point;
506: if (ds) DMGetCellDS(dm, ls, ds);
507: }
508: }
509: ISRestoreIndices(labelIS, &points);
510: }
511: ISDestroy(&labelIS);
512: if (ls >= 0) break;
513: }
514: if (point) *point = ls;
515: DMDestroy(&plex);
516: return 0;
517: }
519: /*
520: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
522: There are several different scenarios:
524: 1) Volumetric mesh with volumetric auxiliary data
526: Here minHeight=0 since we loop over cells.
528: 2) Boundary mesh with boundary auxiliary data
530: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
532: 3) Volumetric mesh with boundary auxiliary data
534: 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.
536: 4) Volumetric input mesh with boundary output mesh
538: Here we must get a subspace for the input DS
540: The maxHeight is used to support enforcement of constraints in DMForest.
542: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
544: 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.
545: - 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
546: 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
547: dual spaces for the boundary from our input spaces.
548: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
550: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
552: 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.
553: */
554: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
555: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
556: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
557: InsertMode mode, Vec localX)
558: {
559: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
560: DMEnclosureType encIn, encAux;
561: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
562: Vec localA = NULL, tv;
563: IS fieldIS;
564: PetscSection section;
565: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
566: PetscTabulation *T = NULL, *TAux = NULL;
567: PetscInt *Nc;
568: PetscInt dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
569: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, transform;
570: DMField coordField;
571: DMLabel depthLabel;
572: PetscQuadrature allPoints = NULL;
573: PetscErrorCode ierr;
575: if (localU) VecGetDM(localU, &dmIn);
576: else {dmIn = dm;}
577: DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA);
578: if (localA) VecGetDM(localA, &dmAux); else {dmAux = NULL;}
579: DMConvert(dm, DMPLEX, &plex);
580: DMConvert(dmIn, DMPLEX, &plexIn);
581: DMGetEnclosureRelation(dmIn, dm, &encIn);
582: DMGetEnclosureRelation(dmAux, dm, &encAux);
583: DMGetDimension(dm, &dim);
584: DMPlexGetVTKCellHeight(plex, &minHeight);
585: DMGetBasisTransformDM_Internal(dm, &tdm);
586: DMGetBasisTransformVec_Internal(dm, &tv);
587: DMHasBasisTransform(dm, &transform);
588: /* Auxiliary information can only be used with interpolation of field functions */
589: if (dmAux) {
590: DMConvert(dmAux, DMPLEX, &plexAux);
591: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
593: }
594: }
595: /* Determine height for iteration of all meshes */
596: {
597: DMPolytopeType ct, ctIn, ctAux;
598: PetscInt minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux;
599: PetscInt dim = -1, dimIn, dimAux;
601: DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
602: if (pEnd > pStart) {
603: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL);
604: p = lStart < 0 ? pStart : lStart;
605: DMPlexGetCellType(plex, p, &ct);
606: dim = DMPolytopeTypeGetDim(ct);
607: DMPlexGetVTKCellHeight(plexIn, &minHeightIn);
608: DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL);
609: DMPlexGetCellType(plexIn, pStartIn, &ctIn);
610: dimIn = DMPolytopeTypeGetDim(ctIn);
611: if (dmAux) {
612: DMPlexGetVTKCellHeight(plexAux, &minHeightAux);
613: DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, NULL);
614: DMPlexGetCellType(plexAux, pStartAux, &ctAux);
615: dimAux = DMPolytopeTypeGetDim(ctAux);
616: } else dimAux = dim;
617: }
618: if (dim < 0) {
619: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
621: /* Fall back to determination based on being a submesh */
622: DMPlexGetSubpointMap(plex, &spmap);
623: DMPlexGetSubpointMap(plexIn, &spmapIn);
624: if (plexAux) DMPlexGetSubpointMap(plexAux, &spmapAux);
625: dim = spmap ? 1 : 0;
626: dimIn = spmapIn ? 1 : 0;
627: dimAux = spmapAux ? 1 : 0;
628: }
629: {
630: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), dimAux);
633: if (dimProj < dim) minHeight = 1;
634: htInc = dim - dimProj;
635: htIncIn = dimIn - dimProj;
636: htIncAux = dimAux - dimProj;
637: }
638: }
639: DMPlexGetDepth(plex, &depth);
640: DMPlexGetDepthLabel(plex, &depthLabel);
641: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
642: maxHeight = PetscMax(maxHeight, minHeight);
644: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, 0, NULL, &ds);
645: if (!ds) DMGetDS(dm, &ds);
646: DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
647: if (!dsIn) DMGetDS(dmIn, &dsIn);
648: PetscDSGetNumFields(ds, &Nf);
649: PetscDSGetNumFields(dsIn, &NfIn);
650: DMGetNumFields(dm, &NfTot);
651: DMFindRegionNum(dm, ds, ®ionNum);
652: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
653: PetscDSIsCohesive(ds, &isCohesive);
654: DMGetCoordinateDim(dm, &dimEmbed);
655: DMGetLocalSection(dm, §ion);
656: if (dmAux) {
657: DMGetDS(dmAux, &dsAux);
658: PetscDSGetNumFields(dsAux, &NfAux);
659: }
660: PetscDSGetComponents(ds, &Nc);
661: PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
662: if (maxHeight > 0) PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);
663: else {cellsp = sp; cellspIn = spIn;}
664: if (localU && localU != localX) DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);
665: /* Get cell dual spaces */
666: for (f = 0; f < Nf; ++f) {
667: PetscDiscType disctype;
669: PetscDSGetDiscType_Internal(ds, f, &disctype);
670: if (disctype == PETSC_DISC_FE) {
671: PetscFE fe;
673: isFE[f] = PETSC_TRUE;
674: hasFE = PETSC_TRUE;
675: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
676: PetscFEGetDualSpace(fe, &cellsp[f]);
677: } else if (disctype == PETSC_DISC_FV) {
678: PetscFV fv;
680: isFE[f] = PETSC_FALSE;
681: hasFV = PETSC_TRUE;
682: PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
683: PetscFVGetDualSpace(fv, &cellsp[f]);
684: } else {
685: isFE[f] = PETSC_FALSE;
686: cellsp[f] = NULL;
687: }
688: }
689: for (f = 0; f < NfIn; ++f) {
690: PetscDiscType disctype;
692: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
693: if (disctype == PETSC_DISC_FE) {
694: PetscFE fe;
696: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fe);
697: PetscFEGetDualSpace(fe, &cellspIn[f]);
698: } else if (disctype == PETSC_DISC_FV) {
699: PetscFV fv;
701: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fv);
702: PetscFVGetDualSpace(fv, &cellspIn[f]);
703: } else {
704: cellspIn[f] = NULL;
705: }
706: }
707: DMGetCoordinateField(dm,&coordField);
708: for (f = 0; f < Nf; ++f) {
709: if (!htInc) {sp[f] = cellsp[f];}
710: else PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);
711: }
712: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
713: PetscFE fem, subfem;
714: PetscDiscType disctype;
715: const PetscReal *points;
716: PetscInt numPoints;
719: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &allPoints);
720: PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
721: PetscMalloc2(NfIn, &T, NfAux, &TAux);
722: for (f = 0; f < NfIn; ++f) {
723: if (!htIncIn) {spIn[f] = cellspIn[f];}
724: else PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);
726: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
727: if (disctype != PETSC_DISC_FE) continue;
728: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
729: if (!htIncIn) {subfem = fem;}
730: else PetscFEGetHeightSubspace(fem, htIncIn, &subfem);
731: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
732: }
733: for (f = 0; f < NfAux; ++f) {
734: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
735: if (disctype != PETSC_DISC_FE) continue;
736: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
737: if (!htIncAux) {subfem = fem;}
738: else PetscFEGetHeightSubspace(fem, htIncAux, &subfem);
739: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
740: }
741: }
742: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
743: for (h = minHeight; h <= maxHeight; h++) {
744: PetscInt hEff = h - minHeight + htInc;
745: PetscInt hEffIn = h - minHeight + htIncIn;
746: PetscInt hEffAux = h - minHeight + htIncAux;
747: PetscDS dsEff = ds;
748: PetscDS dsEffIn = dsIn;
749: PetscDS dsEffAux = dsAux;
750: PetscScalar *values;
751: PetscBool *fieldActive;
752: PetscInt maxDegree;
753: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
754: IS heightIS;
756: if (h > minHeight) {
757: for (f = 0; f < Nf; ++f) PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);
758: }
759: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
760: DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL);
761: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
762: if (pEnd <= pStart) {
763: ISDestroy(&heightIS);
764: continue;
765: }
766: /* Compute totDim, the number of dofs in the closure of a point at this height */
767: totDim = 0;
768: for (f = 0; f < Nf; ++f) {
769: PetscBool cohesive;
771: if (!sp[f]) continue;
772: PetscDSGetCohesive(ds, f, &cohesive);
773: PetscDualSpaceGetDimension(sp[f], &spDim);
774: totDim += spDim;
775: if (isCohesive && !cohesive) totDim += spDim;
776: }
777: p = lStart < 0 ? pStart : lStart;
778: DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
780: if (!totDim) {
781: ISDestroy(&heightIS);
782: continue;
783: }
784: if (htInc) PetscDSGetHeightSubspace(ds, hEff, &dsEff);
785: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
786: if (localU) {
787: PetscInt totDimIn, pIn, numValuesIn;
789: totDimIn = 0;
790: for (f = 0; f < NfIn; ++f) {
791: PetscBool cohesive;
793: if (!spIn[f]) continue;
794: PetscDSGetCohesive(dsIn, f, &cohesive);
795: PetscDualSpaceGetDimension(spIn[f], &spDim);
796: totDimIn += spDim;
797: if (isCohesive && !cohesive) totDimIn += spDim;
798: }
799: DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
800: DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
802: if (htIncIn) PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);
803: }
804: if (htIncAux) PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);
805: /* Loop over points at this height */
806: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
807: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
808: {
809: const PetscInt *fields;
811: ISGetIndices(fieldIS, &fields);
812: for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
813: for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
814: ISRestoreIndices(fieldIS, &fields);
815: }
816: if (label) {
817: PetscInt i;
819: for (i = 0; i < numIds; ++i) {
820: IS pointIS, isectIS;
821: const PetscInt *points;
822: PetscInt n;
823: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
824: PetscQuadrature quad = NULL;
826: DMLabelGetStratumIS(label, ids[i], &pointIS);
827: if (!pointIS) continue; /* No points with that id on this process */
828: ISIntersect(pointIS,heightIS,&isectIS);
829: ISDestroy(&pointIS);
830: if (!isectIS) continue;
831: ISGetLocalSize(isectIS, &n);
832: ISGetIndices(isectIS, &points);
833: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
834: if (maxDegree <= 1) {
835: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
836: }
837: if (!quad) {
838: if (!h && allPoints) {
839: quad = allPoints;
840: allPoints = NULL;
841: } else {
842: PetscDualSpaceGetAllPointsUnion(Nf,sp,isCohesive ? dim-htInc-1 : dim-htInc,funcs,&quad);
843: }
844: }
845: DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
846: for (p = 0; p < n; ++p) {
847: const PetscInt point = points[p];
849: PetscArrayzero(values, numValues);
850: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
851: DMPlexSetActivePoint(dm, point);
852: 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);
853: if (ierr) {
854: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
855: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
856: ierr;
857: }
858: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);
859: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
860: }
861: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
862: PetscFEGeomDestroy(&fegeom);
863: PetscQuadratureDestroy(&quad);
864: ISRestoreIndices(isectIS, &points);
865: ISDestroy(&isectIS);
866: }
867: } else {
868: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
869: PetscQuadrature quad = NULL;
870: IS pointIS;
872: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
873: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
874: if (maxDegree <= 1) {
875: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
876: }
877: if (!quad) {
878: if (!h && allPoints) {
879: quad = allPoints;
880: allPoints = NULL;
881: } else {
882: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &quad);
883: }
884: }
885: DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
886: for (p = pStart; p < pEnd; ++p) {
887: PetscArrayzero(values, numValues);
888: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
889: DMPlexSetActivePoint(dm, p);
890: 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);
891: if (ierr) {
892: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
893: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
894: ierr;
895: }
896: if (transform) DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);
897: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
898: }
899: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
900: PetscFEGeomDestroy(&fegeom);
901: PetscQuadratureDestroy(&quad);
902: ISDestroy(&pointIS);
903: }
904: ISDestroy(&heightIS);
905: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
906: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
907: }
908: /* Cleanup */
909: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
910: for (f = 0; f < NfIn; ++f) PetscTabulationDestroy(&T[f]);
911: for (f = 0; f < NfAux; ++f) PetscTabulationDestroy(&TAux[f]);
912: PetscFree2(T, TAux);
913: }
914: PetscQuadratureDestroy(&allPoints);
915: PetscFree3(isFE, sp, spIn);
916: if (maxHeight > 0) PetscFree2(cellsp, cellspIn);
917: DMDestroy(&plex);
918: DMDestroy(&plexIn);
919: if (dmAux) DMDestroy(&plexAux);
920: return 0;
921: }
923: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
924: {
925: DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
926: return 0;
927: }
929: 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)
930: {
931: DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
932: return 0;
933: }
935: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
936: void (**funcs)(PetscInt, PetscInt, PetscInt,
937: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
938: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
939: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
940: InsertMode mode, Vec localX)
941: {
942: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
943: return 0;
944: }
946: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
947: void (**funcs)(PetscInt, PetscInt, PetscInt,
948: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
949: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
950: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
951: InsertMode mode, Vec localX)
952: {
953: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
954: return 0;
955: }
957: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
958: void (**funcs)(PetscInt, PetscInt, PetscInt,
959: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
960: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
961: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
962: InsertMode mode, Vec localX)
963: {
964: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
965: return 0;
966: }