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: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`
19: @*/
20: PetscErrorCode DMPlexGetActivePoint(DM dm, PetscInt *point)
21: {
22: PetscFunctionBeginHot;
23: *point = ((DM_Plex *)dm->data)->activePoint;
24: PetscFunctionReturn(PETSC_SUCCESS);
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: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetActivePoint()`
39: @*/
40: PetscErrorCode DMPlexSetActivePoint(DM dm, PetscInt point)
41: {
42: PetscFunctionBeginHot;
43: ((DM_Plex *)dm->data)->activePoint = point;
44: PetscFunctionReturn(PETSC_SUCCESS);
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:[](ch_unstructured), `DM`, `DMPLEX`, `PetscDS`, `PetscFEGeom`, `PetscFVCellGeom`, `PetscDualSpace`
69: */
70: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, PetscScalar values[])
71: {
72: PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp;
73: PetscBool isAffine, isCohesive, transform;
75: PetscFunctionBeginHot;
76: PetscCall(DMGetCoordinateDim(dmIn, &coordDim));
77: PetscCall(DMHasBasisTransform(dmIn, &transform));
78: PetscCall(PetscDSGetNumFields(ds, &Nf));
79: PetscCall(PetscDSGetComponents(ds, &Nc));
80: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
81: /* Get values for closure */
82: isAffine = fegeom->isAffine;
83: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
84: void *const ctx = ctxs ? ctxs[f] : NULL;
85: PetscBool cohesive;
87: if (!sp[f]) continue;
88: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
89: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
90: if (funcs[f]) {
91: if (isFE[f]) {
92: PetscQuadrature allPoints;
93: PetscInt q, dim, numPoints;
94: const PetscReal *points;
95: PetscScalar *pointEval;
96: PetscReal *x;
97: DM rdm;
99: PetscCall(PetscDualSpaceGetDM(sp[f], &rdm));
100: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
101: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
102: PetscCall(DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
103: PetscCall(DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x));
104: PetscCall(PetscArrayzero(pointEval, numPoints * Nc[f]));
105: for (q = 0; q < numPoints; q++, tp++) {
106: const PetscReal *v0;
108: if (isAffine) {
109: const PetscReal *refpoint = &points[q * dim];
110: PetscReal injpoint[3] = {0., 0., 0.};
112: if (dim != fegeom->dim) {
113: if (isCohesive) {
114: /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */
115: for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d];
116: refpoint = injpoint;
117: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %" PetscInt_FMT " != %" PetscInt_FMT " dual basis spatial dimension", fegeom->dim, dim);
118: }
119: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x);
120: v0 = x;
121: } else {
122: v0 = &fegeom->v[tp * coordDim];
123: }
124: if (transform) {
125: PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx));
126: v0 = x;
127: }
128: PetscCall((*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx));
129: }
130: /* Transform point evaluations pointEval[q,c] */
131: PetscCall(PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval));
132: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
133: PetscCall(DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x));
134: PetscCall(DMRestoreWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
135: v += spDim;
136: if (isCohesive && !cohesive) {
137: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
138: }
139: } else {
140: for (d = 0; d < spDim; ++d, ++v) PetscCall(PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]));
141: }
142: } else {
143: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
144: if (isCohesive && !cohesive) {
145: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
146: }
147: }
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
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: Level: developer
178: Note:
179: Not supported for FV
181: .seealso: `DMProjectPoint_Field_Private()`
182: */
183: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
184: {
185: PetscSection section, sectionAux = NULL;
186: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
187: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
188: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
189: const PetscScalar *constants;
190: PetscReal *x;
191: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
192: PetscFEGeom fegeom;
193: const PetscInt dE = cgeom->dimEmbed, *cone, *ornt;
194: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
195: PetscBool isAffine, isCohesive, isCohesiveIn, transform;
196: DMPolytopeType qct;
198: PetscFunctionBeginHot;
199: PetscCall(PetscDSGetNumFields(ds, &Nf));
200: PetscCall(PetscDSGetComponents(ds, &Nc));
201: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
202: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
203: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
204: PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
205: PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
206: PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
207: PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
208: PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
209: PetscCall(DMHasBasisTransform(dmIn, &transform));
210: PetscCall(DMGetLocalSection(dmIn, §ion));
211: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
212: // Get cohesive cell hanging off face
213: if (isCohesiveIn) {
214: PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
215: if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
216: DMPolytopeType ct;
217: const PetscInt *support;
218: PetscInt Ns, s;
220: PetscCall(DMPlexGetSupport(dmIn, inp, &support));
221: PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
222: for (s = 0; s < Ns; ++s) {
223: PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
224: if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
225: inp = support[s];
226: break;
227: }
228: }
229: PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
230: PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
231: PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
232: face[0] = 0;
233: face[1] = 0;
234: }
235: }
236: if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
237: if (dmAux) {
238: PetscInt subp;
240: PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
241: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
242: PetscCall(DMGetLocalSection(dmAux, §ionAux));
243: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
244: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
245: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
246: PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
247: }
248: /* Get values for closure */
249: isAffine = cgeom->isAffine;
250: fegeom.dim = cgeom->dim;
251: fegeom.dimEmbed = cgeom->dimEmbed;
252: if (isAffine) {
253: fegeom.v = x;
254: fegeom.xi = cgeom->xi;
255: fegeom.J = cgeom->J;
256: fegeom.invJ = cgeom->invJ;
257: fegeom.detJ = cgeom->detJ;
258: }
259: for (f = 0, v = 0; f < Nf; ++f) {
260: PetscQuadrature allPoints;
261: PetscInt q, dim, numPoints;
262: const PetscReal *points;
263: PetscScalar *pointEval;
264: PetscBool cohesive;
265: DM dm;
267: if (!sp[f]) continue;
268: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
269: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
270: if (!funcs[f]) {
271: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
272: if (isCohesive && !cohesive) {
273: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
274: }
275: continue;
276: }
277: const PetscInt ***perms;
278: PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
279: PetscCall(PetscDualSpaceGetSymmetries(sp[f], &perms, NULL));
280: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
281: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
282: PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
283: for (q = 0; q < numPoints; ++q, ++tp) {
284: PetscInt qpt[2];
286: if (isCohesiveIn) {
287: qpt[0] = perms ? perms[0][ornt[0]][q] : q;
288: qpt[1] = perms ? perms[0][DMPolytopeTypeComposeOrientationInv(qct, ornt[1], 0)][q] : q;
289: }
290: if (isAffine) {
291: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
292: } else {
293: fegeom.v = &cgeom->v[tp * dE];
294: fegeom.J = &cgeom->J[tp * dE * dE];
295: fegeom.invJ = &cgeom->invJ[tp * dE * dE];
296: fegeom.detJ = &cgeom->detJ[tp];
297: }
298: if (coefficients) {
299: if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
300: else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
301: }
302: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
303: if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
304: (*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]);
305: }
306: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
307: PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
308: v += spDim;
309: /* TODO: For now, set both sides equal, but this should use info from other support cell */
310: if (isCohesive && !cohesive) {
311: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
312: }
313: }
314: if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
315: if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
316: if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, PetscScalar values[])
321: {
322: PetscSection section, sectionAux = NULL;
323: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
324: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
325: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
326: const PetscScalar *constants;
327: PetscReal *x;
328: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
329: PetscFEGeom fegeom, cgeom;
330: const PetscInt dE = fgeom->dimEmbed, *cone, *ornt;
331: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
332: PetscBool isAffine, isCohesive, isCohesiveIn, transform;
333: DMPolytopeType qct;
335: PetscFunctionBeginHot;
336: PetscCall(PetscDSGetNumFields(ds, &Nf));
337: PetscCall(PetscDSGetComponents(ds, &Nc));
338: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
339: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
340: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
341: PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
342: PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
343: PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
344: PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
345: PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
346: PetscCall(DMHasBasisTransform(dmIn, &transform));
347: PetscCall(DMGetLocalSection(dmIn, §ion));
348: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
349: // Get cohesive cell hanging off face
350: if (isCohesiveIn) {
351: PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
352: if ((qct != DM_POLYTOPE_POINT_PRISM_TENSOR) && (qct != DM_POLYTOPE_SEG_PRISM_TENSOR) && (qct != DM_POLYTOPE_TRI_PRISM_TENSOR) && (qct != DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
353: DMPolytopeType ct;
354: const PetscInt *support;
355: PetscInt Ns, s;
357: PetscCall(DMPlexGetSupport(dmIn, inp, &support));
358: PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
359: for (s = 0; s < Ns; ++s) {
360: PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
361: if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
362: inp = support[s];
363: break;
364: }
365: }
366: PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
367: PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
368: PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
369: face[0] = 0;
370: face[1] = 0;
371: }
372: }
373: if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
374: if (dmAux) {
375: PetscInt subp;
377: PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
378: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
379: PetscCall(DMGetLocalSection(dmAux, §ionAux));
380: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
381: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
382: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
383: PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
384: }
385: /* Get values for closure */
386: isAffine = fgeom->isAffine;
387: fegeom.n = NULL;
388: fegeom.J = NULL;
389: fegeom.v = NULL;
390: fegeom.xi = NULL;
391: cgeom.dim = fgeom->dim;
392: cgeom.dimEmbed = fgeom->dimEmbed;
393: if (isAffine) {
394: fegeom.v = x;
395: fegeom.xi = fgeom->xi;
396: fegeom.J = fgeom->J;
397: fegeom.invJ = fgeom->invJ;
398: fegeom.detJ = fgeom->detJ;
399: fegeom.n = fgeom->n;
401: cgeom.J = fgeom->suppJ[0];
402: cgeom.invJ = fgeom->suppInvJ[0];
403: cgeom.detJ = fgeom->suppDetJ[0];
404: }
405: for (f = 0, v = 0; f < Nf; ++f) {
406: PetscQuadrature allPoints;
407: PetscInt q, dim, numPoints;
408: const PetscReal *points;
409: PetscScalar *pointEval;
410: PetscBool cohesive;
411: DM dm;
413: if (!sp[f]) continue;
414: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
415: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
416: if (!funcs[f]) {
417: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
418: if (isCohesive && !cohesive) {
419: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
420: }
421: continue;
422: }
423: PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
424: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
425: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
426: PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
427: for (q = 0; q < numPoints; ++q, ++tp) {
428: PetscInt qpt[2];
430: if (isCohesiveIn) {
431: // These points are not integration quadratures, but dual space quadratures
432: // If they had multiple points we should match them from both sides, similar to hybrid residual eval
433: qpt[0] = qpt[1] = q;
434: }
435: if (isAffine) {
436: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
437: } else {
438: fegeom.v = &fgeom->v[tp * dE];
439: fegeom.J = &fgeom->J[tp * dE * dE];
440: fegeom.invJ = &fgeom->invJ[tp * dE * dE];
441: fegeom.detJ = &fgeom->detJ[tp];
442: fegeom.n = &fgeom->n[tp * dE];
444: cgeom.J = &fgeom->suppJ[0][tp * dE * dE];
445: cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
446: cgeom.detJ = &fgeom->suppDetJ[0][tp];
447: }
448: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
449: if (coefficients) {
450: if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
451: else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
452: }
453: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
454: if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
455: (*funcs[f])(dE, NfIn, 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]);
456: }
457: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
458: PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
459: v += spDim;
460: /* TODO: For now, set both sides equal, but this should use info from other support cell */
461: if (isCohesive && !cohesive) {
462: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
463: }
464: }
465: if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
466: if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
467: if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscFEGeom *fegeom, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
472: {
473: PetscFVCellGeom fvgeom;
474: PetscInt dim, dimEmbed;
476: PetscFunctionBeginHot;
477: PetscCall(DMGetDimension(dm, &dim));
478: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
479: if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
480: switch (type) {
481: case DM_BC_ESSENTIAL:
482: case DM_BC_NATURAL:
483: PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
484: break;
485: case DM_BC_ESSENTIAL_FIELD:
486: case DM_BC_NATURAL_FIELD:
487: PetscCall(DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
488: break;
489: case DM_BC_ESSENTIAL_BD_FIELD:
490: PetscCall(DMProjectPoint_BdField_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))funcs, ctxs, values));
491: break;
492: default:
493: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
494: }
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
499: {
500: PetscReal *points;
501: PetscInt f, numPoints;
503: PetscFunctionBegin;
504: if (!dim) {
505: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
508: numPoints = 0;
509: for (f = 0; f < Nf; ++f) {
510: if (funcs[f]) {
511: PetscQuadrature fAllPoints;
512: PetscInt fNumPoints;
514: PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
515: PetscCall(PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL));
516: numPoints += fNumPoints;
517: }
518: }
519: PetscCall(PetscMalloc1(dim * numPoints, &points));
520: numPoints = 0;
521: for (f = 0; f < Nf; ++f) {
522: if (funcs[f]) {
523: PetscQuadrature fAllPoints;
524: PetscInt qdim, fNumPoints, q;
525: const PetscReal *fPoints;
527: PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
528: PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
529: PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %" PetscInt_FMT " for dual basis does not match input dimension %" PetscInt_FMT, qdim, dim);
530: for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
531: numPoints += fNumPoints;
532: }
533: }
534: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
535: PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@C
540: DMGetFirstLabeledPoint - Find first labeled `point` in `odm` such that the corresponding point in `dm` has the specified `height`. Return `point` and the corresponding `ds`.
542: Input Parameters:
543: + dm - the `DM`
544: . odm - the enclosing `DM`
545: . label - label for `DM` domain, or `NULL` for whole domain
546: . numIds - the number of `ids`
547: . ids - An array of the label ids in sequence for the domain
548: - height - Height of target cells in `DMPLEX` topology
550: Output Parameters:
551: + point - the first labeled point
552: - ds - the `PetscDS` corresponding to the first labeled point
554: Level: developer
556: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
557: @*/
558: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
559: {
560: DM plex;
561: DMEnclosureType enc;
562: PetscInt ls = -1;
564: PetscFunctionBegin;
565: if (point) *point = -1;
566: if (!label) PetscFunctionReturn(PETSC_SUCCESS);
567: PetscCall(DMGetEnclosureRelation(dm, odm, &enc));
568: PetscCall(DMConvert(dm, DMPLEX, &plex));
569: for (PetscInt i = 0; i < numIds; ++i) {
570: IS labelIS;
571: PetscInt num_points, pStart, pEnd;
572: PetscCall(DMLabelGetStratumIS(label, ids[i], &labelIS));
573: if (!labelIS) continue; /* No points with that id on this process */
574: PetscCall(DMPlexGetHeightStratum(plex, height, &pStart, &pEnd));
575: PetscCall(ISGetSize(labelIS, &num_points));
576: if (num_points) {
577: const PetscInt *points;
578: PetscCall(ISGetIndices(labelIS, &points));
579: for (PetscInt i = 0; i < num_points; i++) {
580: PetscInt point;
581: PetscCall(DMGetEnclosurePoint(dm, odm, enc, points[i], &point));
582: if (pStart <= point && point < pEnd) {
583: ls = point;
584: if (ds) {
585: // If this is a face of a cohesive cell, then prefer that DS
586: if (height == 1) {
587: const PetscInt *supp;
588: PetscInt suppSize;
589: DMPolytopeType ct;
591: PetscCall(DMPlexGetSupport(dm, ls, &supp));
592: PetscCall(DMPlexGetSupportSize(dm, ls, &suppSize));
593: for (PetscInt s = 0; s < suppSize; ++s) {
594: PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
595: if ((ct == DM_POLYTOPE_POINT_PRISM_TENSOR) || (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) || (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) || (ct == DM_POLYTOPE_QUAD_PRISM_TENSOR)) {
596: ls = supp[s];
597: break;
598: }
599: }
600: }
601: PetscCall(DMGetCellDS(dm, ls, ds, NULL));
602: }
603: if (ls >= 0) break;
604: }
605: }
606: PetscCall(ISRestoreIndices(labelIS, &points));
607: }
608: PetscCall(ISDestroy(&labelIS));
609: if (ls >= 0) break;
610: }
611: if (point) *point = ls;
612: PetscCall(DMDestroy(&plex));
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: /*
617: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
619: There are several different scenarios:
621: 1) Volumetric mesh with volumetric auxiliary data
623: Here minHeight=0 since we loop over cells.
625: 2) Boundary mesh with boundary auxiliary data
627: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
629: 3) Volumetric mesh with boundary auxiliary data
631: 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.
633: 4) Volumetric input mesh with boundary output mesh
635: Here we must get a subspace for the input DS
637: The maxHeight is used to support enforcement of constraints in DMForest.
639: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
641: 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.
642: - 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
643: 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
644: dual spaces for the boundary from our input spaces.
645: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
647: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
649: 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.
650: */
651: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, InsertMode mode, Vec localX)
652: {
653: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
654: DMEnclosureType encIn, encAux;
655: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
656: Vec localA = NULL, tv;
657: IS fieldIS;
658: PetscSection section;
659: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
660: PetscTabulation *T = NULL, *TAux = NULL;
661: PetscInt *Nc;
662: PetscInt dim, dimEmbed, depth, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, minHeightIn, minHeightAux, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
663: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, isCohesiveIn = PETSC_FALSE, transform;
664: DMField coordField;
665: DMLabel depthLabel;
666: PetscQuadrature allPoints = NULL;
668: PetscFunctionBegin;
669: if (localU) PetscCall(VecGetDM(localU, &dmIn));
670: else dmIn = dm;
671: PetscCall(DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA));
672: if (localA) PetscCall(VecGetDM(localA, &dmAux));
673: else dmAux = NULL;
674: PetscCall(DMConvert(dm, DMPLEX, &plex));
675: PetscCall(DMConvert(dmIn, DMPLEX, &plexIn));
676: PetscCall(DMGetEnclosureRelation(dmIn, dm, &encIn));
677: PetscCall(DMGetEnclosureRelation(dmAux, dm, &encAux));
678: PetscCall(DMGetDimension(dm, &dim));
679: PetscCall(DMPlexGetVTKCellHeight(plex, &minHeight));
680: PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
681: PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));
682: PetscCall(DMHasBasisTransform(dm, &transform));
683: /* Auxiliary information can only be used with interpolation of field functions */
684: if (dmAux) {
685: PetscCall(DMConvert(dmAux, DMPLEX, &plexAux));
686: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) PetscCheck(localA, PETSC_COMM_SELF, PETSC_ERR_USER, "Missing localA vector");
687: }
688: if (localU && localU != localX) PetscCall(DMPlexInsertBoundaryValues(plexIn, PETSC_TRUE, localU, time, NULL, NULL, NULL));
689: PetscCall(DMGetCoordinateField(dm, &coordField));
690: PetscCheck(coordField, PETSC_COMM_SELF, PETSC_ERR_USER, "DM must have a coordinate field");
691: /**** No collective calls below this point ****/
692: /* Determine height for iteration of all meshes */
693: {
694: DMPolytopeType ct, ctIn, ctAux;
695: PetscInt lStart, pStart, pEnd, p, pStartIn, pStartAux, pEndAux;
696: PetscInt dim = -1, dimIn = -1, dimAux = -1;
698: PetscCall(DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd));
699: if (pEnd > pStart) {
700: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL));
701: p = lStart < 0 ? pStart : lStart;
702: PetscCall(DMPlexGetCellType(plex, p, &ct));
703: dim = DMPolytopeTypeGetDim(ct);
704: PetscCall(DMPlexGetVTKCellHeight(plexIn, &minHeightIn));
705: PetscCall(DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL));
706: PetscCall(DMPlexGetCellType(plexIn, pStartIn, &ctIn));
707: dimIn = DMPolytopeTypeGetDim(ctIn);
708: if (dmAux) {
709: PetscCall(DMPlexGetVTKCellHeight(plexAux, &minHeightAux));
710: PetscCall(DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux));
711: if (pStartAux < pEndAux) {
712: PetscCall(DMPlexGetCellType(plexAux, pStartAux, &ctAux));
713: dimAux = DMPolytopeTypeGetDim(ctAux);
714: }
715: } else dimAux = dim;
716: } else {
717: PetscCall(DMDestroy(&plex));
718: PetscCall(DMDestroy(&plexIn));
719: if (dmAux) PetscCall(DMDestroy(&plexAux));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
722: if (dim < 0) {
723: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
725: /* Fall back to determination based on being a submesh */
726: PetscCall(DMPlexGetSubpointMap(plex, &spmap));
727: PetscCall(DMPlexGetSubpointMap(plexIn, &spmapIn));
728: if (plexAux) PetscCall(DMPlexGetSubpointMap(plexAux, &spmapAux));
729: dim = spmap ? 1 : 0;
730: dimIn = spmapIn ? 1 : 0;
731: dimAux = spmapAux ? 1 : 0;
732: }
733: {
734: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), dimAux < 0 ? PETSC_INT_MAX : dimAux);
735: PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;
737: PetscCheck(PetscAbsInt(dimProj - dim) <= 1 && PetscAbsInt(dimProj - dimIn) <= 1 && PetscAbsInt(dimProj - dimAuxEff) <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not currently support differences of more than 1 in dimension");
738: if (dimProj < dim) minHeight = 1;
739: htInc = dim - dimProj;
740: htIncIn = dimIn - dimProj;
741: htIncAux = dimAuxEff - dimProj;
742: }
743: }
744: PetscCall(DMPlexGetDepth(plex, &depth));
745: PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
746: PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
747: maxHeight = PetscMax(maxHeight, minHeight);
748: PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
749: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, NULL, &ds));
750: if (!ds) PetscCall(DMGetDS(dm, &ds));
751: PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, minHeight, NULL, &dsIn));
752: if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
753: PetscCall(PetscDSGetNumFields(ds, &Nf));
754: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
755: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
756: if (isCohesiveIn) --htIncIn; // Should be rearranged
757: PetscCall(DMGetNumFields(dm, &NfTot));
758: PetscCall(DMFindRegionNum(dm, ds, ®ionNum));
759: PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
760: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
761: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
762: PetscCall(DMGetLocalSection(dm, §ion));
763: if (dmAux) {
764: PetscCall(DMGetDS(dmAux, &dsAux));
765: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
766: }
767: PetscCall(PetscDSGetComponents(ds, &Nc));
768: PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
769: if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
770: else {
771: cellsp = sp;
772: cellspIn = spIn;
773: }
774: /* Get cell dual spaces */
775: for (f = 0; f < Nf; ++f) {
776: PetscDiscType disctype;
778: PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
779: if (disctype == PETSC_DISC_FE) {
780: PetscFE fe;
782: isFE[f] = PETSC_TRUE;
783: hasFE = PETSC_TRUE;
784: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
785: PetscCall(PetscFEGetDualSpace(fe, &cellsp[f]));
786: } else if (disctype == PETSC_DISC_FV) {
787: PetscFV fv;
789: isFE[f] = PETSC_FALSE;
790: hasFV = PETSC_TRUE;
791: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fv));
792: PetscCall(PetscFVGetDualSpace(fv, &cellsp[f]));
793: } else {
794: isFE[f] = PETSC_FALSE;
795: cellsp[f] = NULL;
796: }
797: }
798: for (f = 0; f < NfIn; ++f) {
799: PetscDiscType disctype;
801: PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
802: if (disctype == PETSC_DISC_FE) {
803: PetscFE fe;
805: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
806: PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
807: } else if (disctype == PETSC_DISC_FV) {
808: PetscFV fv;
810: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv));
811: PetscCall(PetscFVGetDualSpace(fv, &cellspIn[f]));
812: } else {
813: cellspIn[f] = NULL;
814: }
815: }
816: for (f = 0; f < Nf; ++f) {
817: if (!htInc) {
818: sp[f] = cellsp[f];
819: } else PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]));
820: }
821: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
822: PetscFE fem, subfem;
823: PetscDiscType disctype;
824: const PetscReal *points;
825: PetscInt numPoints;
827: PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
828: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
829: PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
830: PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
831: for (f = 0; f < NfIn; ++f) {
832: if (!htIncIn) {
833: spIn[f] = cellspIn[f];
834: } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));
836: PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
837: if (disctype != PETSC_DISC_FE) continue;
838: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem));
839: if (!htIncIn) {
840: subfem = fem;
841: } else PetscCall(PetscFEGetHeightSubspace(fem, htIncIn, &subfem));
842: PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]));
843: }
844: for (f = 0; f < NfAux; ++f) {
845: PetscCall(PetscDSGetDiscType_Internal(dsAux, f, &disctype));
846: if (disctype != PETSC_DISC_FE) continue;
847: PetscCall(PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem));
848: if (!htIncAux) {
849: subfem = fem;
850: } else PetscCall(PetscFEGetHeightSubspace(fem, htIncAux, &subfem));
851: PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]));
852: }
853: }
854: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
855: for (h = minHeight; h <= maxHeight; h++) {
856: PetscInt hEff = h - minHeight + htInc;
857: PetscInt hEffIn = h - minHeight + htIncIn;
858: PetscInt hEffAux = h - minHeight + htIncAux;
859: PetscDS dsEff = ds;
860: PetscDS dsEffIn = dsIn;
861: PetscDS dsEffAux = dsAux;
862: PetscScalar *values;
863: PetscBool *fieldActive;
864: PetscInt maxDegree;
865: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
866: IS heightIS;
868: if (h > minHeight) {
869: for (f = 0; f < Nf; ++f) PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]));
870: }
871: PetscCall(DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd));
872: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL));
873: PetscCall(DMLabelGetStratumIS(depthLabel, depth - h, &heightIS));
874: if (pEnd <= pStart) {
875: PetscCall(ISDestroy(&heightIS));
876: continue;
877: }
878: /* Compute totDim, the number of dofs in the closure of a point at this height */
879: totDim = 0;
880: for (f = 0; f < Nf; ++f) {
881: PetscBool cohesive;
883: if (!sp[f]) continue;
884: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
885: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
886: totDim += spDim;
887: if (isCohesive && !cohesive) totDim += spDim;
888: }
889: p = lStart < 0 ? pStart : lStart;
890: PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
891: PetscCheck(numValues == totDim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The output section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT " in [%" PetscInt_FMT ", %" PetscInt_FMT "]", p, numValues, totDim, h, minHeight, maxHeight);
892: if (!totDim) {
893: PetscCall(ISDestroy(&heightIS));
894: continue;
895: }
896: if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
897: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
898: if (localU) {
899: PetscInt totDimIn, pIn, numValuesIn;
901: totDimIn = 0;
902: for (f = 0; f < NfIn; ++f) {
903: PetscBool cohesive;
905: if (!spIn[f]) continue;
906: PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
907: PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
908: totDimIn += spDim;
909: if (isCohesiveIn && !cohesive) totDimIn += spDim;
910: }
911: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
912: PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
913: // TODO We could check that pIn is a cohesive cell for this check
914: PetscCheck(isCohesiveIn || (numValuesIn == totDimIn), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The input section point (%" PetscInt_FMT ") closure size %" PetscInt_FMT " != dual space dimension %" PetscInt_FMT " at height %" PetscInt_FMT, pIn, numValuesIn, totDimIn, htIncIn);
915: if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
916: }
917: if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
918: /* Loop over points at this height */
919: PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
920: PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
921: {
922: const PetscInt *fields;
924: PetscCall(ISGetIndices(fieldIS, &fields));
925: for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
926: for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
927: PetscCall(ISRestoreIndices(fieldIS, &fields));
928: }
929: if (label) {
930: PetscInt i;
932: for (i = 0; i < numIds; ++i) {
933: IS pointIS, isectIS;
934: const PetscInt *points;
935: PetscInt n;
936: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
937: PetscQuadrature quad = NULL;
939: PetscCall(DMLabelGetStratumIS(label, ids[i], &pointIS));
940: if (!pointIS) continue; /* No points with that id on this process */
941: PetscCall(ISIntersect(pointIS, heightIS, &isectIS));
942: PetscCall(ISDestroy(&pointIS));
943: if (!isectIS) continue;
944: PetscCall(ISGetLocalSize(isectIS, &n));
945: PetscCall(ISGetIndices(isectIS, &points));
946: PetscCall(DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree));
947: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad));
948: if (!quad) {
949: if (!h && allPoints) {
950: quad = allPoints;
951: allPoints = NULL;
952: } else {
953: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad));
954: }
955: }
956: PetscBool computeFaceGeom = htInc && h == minHeight ? PETSC_TRUE : PETSC_FALSE;
958: if (n) {
959: PetscInt depth, dep;
961: PetscCall(DMPlexGetDepth(dm, &depth));
962: PetscCall(DMPlexGetPointDepth(dm, points[0], &dep));
963: if (dep < depth && h == minHeight) computeFaceGeom = PETSC_TRUE;
964: }
965: PetscCall(DMFieldCreateFEGeom(coordField, isectIS, quad, computeFaceGeom, &fegeom));
966: for (p = 0; p < n; ++p) {
967: const PetscInt point = points[p];
969: PetscCall(PetscArrayzero(values, numValues));
970: PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
971: PetscCall(DMPlexSetActivePoint(dm, point));
972: PetscCall(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));
973: if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
974: PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
975: }
976: PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
977: PetscCall(PetscFEGeomDestroy(&fegeom));
978: PetscCall(PetscQuadratureDestroy(&quad));
979: PetscCall(ISRestoreIndices(isectIS, &points));
980: PetscCall(ISDestroy(&isectIS));
981: }
982: } else {
983: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
984: PetscQuadrature quad = NULL;
985: IS pointIS;
987: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS));
988: PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
989: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
990: if (!quad) {
991: if (!h && allPoints) {
992: quad = allPoints;
993: allPoints = NULL;
994: } else {
995: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad));
996: }
997: }
998: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_TRUE : PETSC_FALSE, &fegeom));
999: for (p = pStart; p < pEnd; ++p) {
1000: PetscCall(PetscArrayzero(values, numValues));
1001: PetscCall(PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom));
1002: PetscCall(DMPlexSetActivePoint(dm, p));
1003: PetscCall(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));
1004: if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values));
1005: PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode));
1006: }
1007: PetscCall(PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom));
1008: PetscCall(PetscFEGeomDestroy(&fegeom));
1009: PetscCall(PetscQuadratureDestroy(&quad));
1010: PetscCall(ISDestroy(&pointIS));
1011: }
1012: PetscCall(ISDestroy(&heightIS));
1013: PetscCall(DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values));
1014: PetscCall(DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive));
1015: }
1016: /* Cleanup */
1017: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
1018: for (f = 0; f < NfIn; ++f) PetscCall(PetscTabulationDestroy(&T[f]));
1019: for (f = 0; f < NfAux; ++f) PetscCall(PetscTabulationDestroy(&TAux[f]));
1020: PetscCall(PetscFree2(T, TAux));
1021: }
1022: PetscCall(PetscQuadratureDestroy(&allPoints));
1023: PetscCall(PetscFree3(isFE, sp, spIn));
1024: if (maxHeight > 0) PetscCall(PetscFree2(cellsp, cellspIn));
1025: PetscCall(DMDestroy(&plex));
1026: PetscCall(DMDestroy(&plexIn));
1027: if (dmAux) PetscCall(DMDestroy(&plexAux));
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1032: {
1033: PetscFunctionBegin;
1034: PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: 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)
1039: {
1040: PetscFunctionBegin;
1041: PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
1046: {
1047: PetscFunctionBegin;
1048: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1049: PetscFunctionReturn(PETSC_SUCCESS);
1050: }
1052: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
1053: {
1054: PetscFunctionBegin;
1055: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1056: PetscFunctionReturn(PETSC_SUCCESS);
1057: }
1059: PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
1060: {
1061: PetscFunctionBegin;
1062: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }