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, 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 derivatives 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: 4) Volumetric input mesh with boundary output mesh
523: Here we must get a subspace for the input DS
525: The maxHeight is used to support enforcement of constraints in DMForest.
527: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
529: If we are using an input field (DM_BC_ESSENTIAL_FIELD or DM_BC_NATURAL_FIELD), we need to evaluate it at all the quadrature points of the dual basis functionals.
530: - We use effectiveHeight to mean the height above our incoming DS. For example, if the DS is for a submesh then the effective height is zero, whereas if the DS
531: is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When the effective height is nonzero, we need to extract
532: dual spaces for the boundary from our input spaces.
533: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
535: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
537: If we have a label, we iterate over those points. This will probably break the maxHeight functionality since we do not check the height of those points.
538: */
539: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
540: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
541: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
542: InsertMode mode, Vec localX)
543: {
544: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
545: DMEnclosureType encIn, encAux;
546: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
547: Vec localA = NULL, tv;
548: IS fieldIS;
549: PetscSection section;
550: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
551: PetscTabulation *T = NULL, *TAux = NULL;
552: PetscInt *Nc;
553: PetscInt dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
554: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isHybrid = PETSC_FALSE, transform;
555: DMField coordField;
556: DMLabel depthLabel;
557: PetscQuadrature allPoints = NULL;
558: PetscErrorCode ierr;
561: if (localU) {VecGetDM(localU, &dmIn);}
562: else {dmIn = dm;}
563: DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, &localA);
564: if (localA) {VecGetDM(localA, &dmAux);} else {dmAux = NULL;}
565: DMConvert(dm, DMPLEX, &plex);
566: DMConvert(dmIn, DMPLEX, &plexIn);
567: DMGetEnclosureRelation(dmIn, dm, &encIn);
568: DMGetEnclosureRelation(dmAux, dm, &encAux);
569: DMGetDimension(dm, &dim);
570: DMPlexGetVTKCellHeight(plex, &minHeight);
571: DMGetBasisTransformDM_Internal(dm, &tdm);
572: DMGetBasisTransformVec_Internal(dm, &tv);
573: DMHasBasisTransform(dm, &transform);
574: /* Auxiliary information can only be used with interpolation of field functions */
575: if (dmAux) {
576: DMConvert(dmAux, DMPLEX, &plexAux);
577: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
578: if (!localA) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Missing localA vector");
579: }
580: }
581: /* Determine height for iteration of all meshes */
582: {
583: DMPolytopeType ct, ctIn, ctAux;
584: PetscInt minHeightIn, minHeightAux, lStart, pStart, pEnd, p, pStartIn, pStartAux;
585: PetscInt dim = -1, dimIn, dimAux;
587: DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd);
588: if (pEnd > pStart) {
589: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, minHeight, &lStart, NULL);
590: p = lStart < 0 ? pStart : lStart;
591: DMPlexGetCellType(plex, p, &ct);
592: dim = DMPolytopeTypeGetDim(ct);
593: DMPlexGetVTKCellHeight(plexIn, &minHeightIn);
594: DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL);
595: DMPlexGetCellType(plexIn, pStartIn, &ctIn);
596: dimIn = DMPolytopeTypeGetDim(ctIn);
597: if (dmAux) {
598: DMPlexGetVTKCellHeight(plexAux, &minHeightAux);
599: DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, NULL);
600: DMPlexGetCellType(plexAux, pStartAux, &ctAux);
601: dimAux = DMPolytopeTypeGetDim(ctAux);
602: } else dimAux = dim;
603: }
604: if (dim < 0) {
605: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
607: /* Fall back to determination based on being a submesh */
608: DMPlexGetSubpointMap(plex, &spmap);
609: DMPlexGetSubpointMap(plexIn, &spmapIn);
610: if (plexAux) {DMPlexGetSubpointMap(plexAux, &spmapAux);}
611: dim = spmap ? 1 : 0;
612: dimIn = spmapIn ? 1 : 0;
613: dimAux = spmapAux ? 1 : 0;
614: }
615: {
616: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), dimAux);
618: if (PetscAbsInt(dimProj - dim) > 1 || PetscAbsInt(dimProj - dimIn) > 1 || PetscAbsInt(dimProj - dimAux) > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not currently support differences of more than 1 in dimension");
619: if (dimProj < dim) minHeight = 1;
620: htInc = dim - dimProj;
621: htIncIn = dimIn - dimProj;
622: htIncAux = dimAux - dimProj;
623: }
624: }
625: DMPlexGetDepth(plex, &depth);
626: DMPlexGetDepthLabel(plex, &depthLabel);
627: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
628: maxHeight = PetscMax(maxHeight, minHeight);
629: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
630: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
631: if (!ds) {DMGetDS(dm, &ds);}
632: DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
633: if (!dsIn) {DMGetDS(dmIn, &dsIn);}
634: PetscDSGetNumFields(ds, &Nf);
635: PetscDSGetNumFields(dsIn, &NfIn);
636: DMGetNumFields(dm, &NfTot);
637: DMFindRegionNum(dm, ds, ®ionNum);
638: DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL);
639: PetscDSGetHybrid(ds, &isHybrid);
640: DMGetCoordinateDim(dm, &dimEmbed);
641: DMGetLocalSection(dm, §ion);
642: if (dmAux) {
643: DMGetDS(dmAux, &dsAux);
644: PetscDSGetNumFields(dsAux, &NfAux);
645: }
646: PetscDSGetComponents(ds, &Nc);
647: PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn);
648: if (maxHeight > 0) {PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn);}
649: else {cellsp = sp; cellspIn = spIn;}
650: if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
651: /* Get cell dual spaces */
652: for (f = 0; f < Nf; ++f) {
653: PetscDiscType disctype;
655: PetscDSGetDiscType_Internal(ds, f, &disctype);
656: if (disctype == PETSC_DISC_FE) {
657: PetscFE fe;
659: isFE[f] = PETSC_TRUE;
660: hasFE = PETSC_TRUE;
661: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
662: PetscFEGetDualSpace(fe, &cellsp[f]);
663: } else if (disctype == PETSC_DISC_FV) {
664: PetscFV fv;
666: isFE[f] = PETSC_FALSE;
667: hasFV = PETSC_TRUE;
668: PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);
669: PetscFVGetDualSpace(fv, &cellsp[f]);
670: } else {
671: isFE[f] = PETSC_FALSE;
672: cellsp[f] = NULL;
673: }
674: }
675: for (f = 0; f < NfIn; ++f) {
676: PetscDiscType disctype;
678: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
679: if (disctype == PETSC_DISC_FE) {
680: PetscFE fe;
682: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fe);
683: PetscFEGetDualSpace(fe, &cellspIn[f]);
684: } else if (disctype == PETSC_DISC_FV) {
685: PetscFV fv;
687: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fv);
688: PetscFVGetDualSpace(fv, &cellspIn[f]);
689: } else {
690: cellspIn[f] = NULL;
691: }
692: }
693: DMGetCoordinateField(dm,&coordField);
694: for (f = 0; f < Nf; ++f) {
695: if (!htInc) {sp[f] = cellsp[f];}
696: else {PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]);}
697: }
698: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
699: PetscFE fem, subfem;
700: PetscDiscType disctype;
701: const PetscReal *points;
702: PetscInt numPoints;
704: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
705: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &allPoints);
706: PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL);
707: PetscMalloc2(NfIn, &T, NfAux, &TAux);
708: for (f = 0; f < NfIn; ++f) {
709: if (!htIncIn) {spIn[f] = cellspIn[f];}
710: else {PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]);}
712: PetscDSGetDiscType_Internal(dsIn, f, &disctype);
713: if (disctype != PETSC_DISC_FE) continue;
714: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
715: if (!htIncIn) {subfem = fem;}
716: else {PetscFEGetHeightSubspace(fem, htIncIn, &subfem);}
717: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
718: }
719: for (f = 0; f < NfAux; ++f) {
720: PetscDSGetDiscType_Internal(dsAux, f, &disctype);
721: if (disctype != PETSC_DISC_FE) continue;
722: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
723: if (!htIncAux) {subfem = fem;}
724: else {PetscFEGetHeightSubspace(fem, htIncAux, &subfem);}
725: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
726: }
727: }
728: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
729: for (h = minHeight; h <= maxHeight; h++) {
730: PetscInt hEff = h - minHeight + htInc;
731: PetscInt hEffIn = h - minHeight + htIncIn;
732: PetscInt hEffAux = h - minHeight + htIncAux;
733: PetscDS dsEff = ds;
734: PetscDS dsEffIn = dsIn;
735: PetscDS dsEffAux = dsAux;
736: PetscScalar *values;
737: PetscBool *fieldActive;
738: PetscInt maxDegree;
739: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
740: IS heightIS;
742: if (h > minHeight) {
743: for (f = 0; f < Nf; ++f) {PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]);}
744: }
745: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
746: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
747: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
748: if (pEnd <= pStart) {
749: ISDestroy(&heightIS);
750: continue;
751: }
752: /* Compute totDim, the number of dofs in the closure of a point at this height */
753: totDim = 0;
754: for (f = 0; f < Nf; ++f) {
755: if (!sp[f]) continue;
756: PetscDualSpaceGetDimension(sp[f], &spDim);
757: totDim += spDim;
758: if (isHybrid && (f < Nf-1)) totDim += spDim;
759: }
760: p = lStart < 0 ? pStart : lStart;
761: DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL);
762: if (numValues != totDim) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The output section point (%D) closure size %D != dual space dimension %D at height %D in [%D, %D]", p, numValues, totDim, h, minHeight, maxHeight);
763: if (!totDim) {
764: ISDestroy(&heightIS);
765: continue;
766: }
767: if (htInc) {PetscDSGetHeightSubspace(ds, hEff, &dsEff);}
768: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
769: if (localU) {
770: PetscInt totDimIn, pIn, numValuesIn;
772: totDimIn = 0;
773: for (f = 0; f < NfIn; ++f) {
774: if (!spIn[f]) continue;
775: PetscDualSpaceGetDimension(spIn[f], &spDim);
776: totDimIn += spDim;
777: if (isHybrid && (f < Nf-1)) totDimIn += spDim;
778: }
779: DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn);
780: DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL);
781: if (numValuesIn != totDimIn) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The input section point (%D) closure size %D != dual space dimension %D at height %D", pIn, numValuesIn, totDimIn, htIncIn);
782: if (htIncIn) {PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn);}
783: }
784: if (htIncAux) {PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux);}
785: /* Loop over points at this height */
786: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
787: DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive);
788: {
789: const PetscInt *fields;
791: ISGetIndices(fieldIS, &fields);
792: for (f = 0; f < NfTot; ++f) {fieldActive[f] = PETSC_FALSE;}
793: for (f = 0; f < Nf; ++f) {fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;}
794: ISRestoreIndices(fieldIS, &fields);
795: }
796: if (label) {
797: PetscInt i;
799: for (i = 0; i < numIds; ++i) {
800: IS pointIS, isectIS;
801: const PetscInt *points;
802: PetscInt n;
803: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
804: PetscQuadrature quad = NULL;
806: DMLabelGetStratumIS(label, ids[i], &pointIS);
807: if (!pointIS) continue; /* No points with that id on this process */
808: ISIntersect(pointIS,heightIS,&isectIS);
809: ISDestroy(&pointIS);
810: if (!isectIS) continue;
811: ISGetLocalSize(isectIS, &n);
812: ISGetIndices(isectIS, &points);
813: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
814: if (maxDegree <= 1) {
815: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
816: }
817: if (!quad) {
818: if (!h && allPoints) {
819: quad = allPoints;
820: allPoints = NULL;
821: } else {
822: PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-htInc-1 : dim-htInc,funcs,&quad);
823: }
824: }
825: DMFieldCreateFEGeom(coordField, isectIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
826: for (p = 0; p < n; ++p) {
827: const PetscInt point = points[p];
829: PetscArrayzero(values, numValues);
830: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
831: DMPlexSetActivePoint(dm, point);
832: 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);
833: if (ierr) {
834: PetscErrorCode ierr2;
835: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
836: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
837:
838: }
839: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
840: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);
841: }
842: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
843: PetscFEGeomDestroy(&fegeom);
844: PetscQuadratureDestroy(&quad);
845: ISRestoreIndices(isectIS, &points);
846: ISDestroy(&isectIS);
847: }
848: } else {
849: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
850: PetscQuadrature quad = NULL;
851: IS pointIS;
853: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
854: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
855: if (maxDegree <= 1) {
856: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
857: }
858: if (!quad) {
859: if (!h && allPoints) {
860: quad = allPoints;
861: allPoints = NULL;
862: } else {
863: PetscDualSpaceGetAllPointsUnion(Nf, sp, dim-htInc, funcs, &quad);
864: }
865: }
866: DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom);
867: for (p = pStart; p < pEnd; ++p) {
868: PetscArrayzero(values, numValues);
869: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
870: DMPlexSetActivePoint(dm, p);
871: 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);
872: if (ierr) {
873: PetscErrorCode ierr2;
874: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
875: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
876:
877: }
878: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
879: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);
880: }
881: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
882: PetscFEGeomDestroy(&fegeom);
883: PetscQuadratureDestroy(&quad);
884: ISDestroy(&pointIS);
885: }
886: ISDestroy(&heightIS);
887: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
888: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
889: }
890: /* Cleanup */
891: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
892: for (f = 0; f < NfIn; ++f) {PetscTabulationDestroy(&T[f]);}
893: for (f = 0; f < NfAux; ++f) {PetscTabulationDestroy(&TAux[f]);}
894: PetscFree2(T, TAux);
895: }
896: PetscQuadratureDestroy(&allPoints);
897: PetscFree3(isFE, sp, spIn);
898: if (maxHeight > 0) {PetscFree2(cellsp, cellspIn);}
899: DMDestroy(&plex);
900: DMDestroy(&plexIn);
901: if (dmAux) {DMDestroy(&plexAux);}
902: return(0);
903: }
905: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
906: {
910: DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
911: return(0);
912: }
914: PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
915: {
919: DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
920: return(0);
921: }
923: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
924: void (**funcs)(PetscInt, PetscInt, PetscInt,
925: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
926: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
927: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
928: InsertMode mode, Vec localX)
929: {
933: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
934: return(0);
935: }
937: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
938: void (**funcs)(PetscInt, PetscInt, PetscInt,
939: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
940: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
941: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
942: InsertMode mode, Vec localX)
943: {
947: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
948: return(0);
949: }
951: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
952: void (**funcs)(PetscInt, PetscInt, PetscInt,
953: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
954: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
955: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
956: InsertMode mode, Vec localX)
957: {
961: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
962: return(0);
963: }