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 debug = ((DM_Plex *)dm->data)->printProject;
73: PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp;
74: PetscBool isAffine, isCohesive, transform;
76: PetscFunctionBeginHot;
77: PetscCall(DMGetCoordinateDim(dmIn, &coordDim));
78: PetscCall(DMHasBasisTransform(dmIn, &transform));
79: PetscCall(PetscDSGetNumFields(ds, &Nf));
80: PetscCall(PetscDSGetComponents(ds, &Nc));
81: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
82: /* Get values for closure */
83: isAffine = fegeom->isAffine;
84: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
85: void *const ctx = ctxs ? ctxs[f] : NULL;
86: PetscBool cohesive;
88: if (!sp[f]) continue;
89: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
90: PetscCall(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: PetscCall(PetscDualSpaceGetDM(sp[f], &rdm));
101: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
102: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
103: PetscCall(DMGetWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
104: PetscCall(DMGetWorkArray(rdm, coordDim, MPIU_REAL, &x));
105: PetscCall(PetscArrayzero(pointEval, numPoints * Nc[f]));
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: PetscCheck(isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %" PetscInt_FMT " != %" PetscInt_FMT " dual basis spatial dimension", fegeom->dim, dim);
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: }
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: if (debug > 3) {
129: PetscInt ap;
130: PetscCall(DMPlexGetActivePoint(dm, &ap));
131: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Project point %" PetscInt_FMT ", analytic: ref (", ap));
132: for (PetscInt d = 0; d < dim; ++d) {
133: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
134: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)points[q * dim + d]));
135: }
136: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") real ("));
137: for (PetscInt d = 0; d < dim; ++d) {
138: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
139: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)v0[d]));
140: }
141: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
142: }
143: PetscCall((*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f] * q], ctx));
144: }
145: /* Transform point evaluations pointEval[q,c] */
146: PetscCall(PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval));
147: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
148: PetscCall(DMRestoreWorkArray(rdm, coordDim, MPIU_REAL, &x));
149: PetscCall(DMRestoreWorkArray(rdm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
150: v += spDim;
151: if (isCohesive && !cohesive) {
152: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
153: }
154: } else {
155: for (d = 0; d < spDim; ++d, ++v) PetscCall(PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]));
156: }
157: } else {
158: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
159: if (isCohesive && !cohesive) {
160: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
161: }
162: }
163: }
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*
168: DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point
170: Input Parameters:
171: + dm - The output DM
172: . ds - The output DS
173: . dmIn - The input DM
174: . dsIn - The input DS
175: . dmAux - The auxiliary DM, which is always for the input space
176: . dsAux - The auxiliary DS, which is always for the input space
177: . time - The time for this evaluation
178: . localU - The local solution
179: . localA - The local auziliary fields
180: . cgeom - The FE geometry for this point
181: . sp - The output PetscDualSpace for each field
182: . p - The point in the output DM
183: . T - Input basis and derivatives for each field tabulated on the quadrature points
184: . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
185: . funcs - The evaluation function for each field
186: - ctxs - The user context for each field
188: Output Parameter:
189: . values - The value for each dual basis vector in the output dual space
191: Level: developer
193: Note:
194: Not supported for FV
196: .seealso: `DMProjectPoint_Field_Private()`
197: */
198: 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[])
199: {
200: PetscSection section, sectionAux = NULL;
201: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
202: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
203: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
204: const PetscScalar *constants;
205: PetscReal *x;
206: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
207: PetscFEGeom fegeom, fgeomN[2];
208: const PetscInt dE = cgeom->dimEmbed, *cone, *ornt;
209: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
210: PetscBool isAffine, isCohesive, isCohesiveIn, transform;
211: DMPolytopeType qct;
213: PetscFunctionBeginHot;
214: PetscCall(PetscDSGetNumFields(ds, &Nf));
215: PetscCall(PetscDSGetComponents(ds, &Nc));
216: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
217: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
218: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
219: PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
220: PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
221: PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
222: PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
223: PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
224: PetscCall(DMHasBasisTransform(dmIn, &transform));
225: PetscCall(DMGetLocalSection(dmIn, §ion));
226: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
227: // Get cohesive cell hanging off face
228: if (isCohesiveIn) {
229: PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
230: 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)) {
231: DMPolytopeType ct;
232: const PetscInt *support;
233: PetscInt Ns, s;
235: PetscCall(DMPlexGetSupport(dmIn, inp, &support));
236: PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
237: for (s = 0; s < Ns; ++s) {
238: PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
239: 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)) {
240: inp = support[s];
241: break;
242: }
243: }
244: PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
245: PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
246: PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
247: face[0] = 0;
248: face[1] = 0;
249: }
250: }
251: if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
252: if (dmAux) {
253: PetscInt subp;
255: PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
256: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
257: PetscCall(DMGetLocalSection(dmAux, §ionAux));
258: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
259: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
260: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
261: PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
262: }
263: /* Get values for closure */
264: isAffine = cgeom->isAffine;
265: fegeom.dim = cgeom->dim;
266: fegeom.dimEmbed = cgeom->dimEmbed;
267: if (isCohesiveIn) {
268: fgeomN[0].dim = cgeom->dim;
269: fgeomN[0].dimEmbed = cgeom->dimEmbed;
270: fgeomN[1].dim = cgeom->dim;
271: fgeomN[1].dimEmbed = cgeom->dimEmbed;
272: }
273: if (isAffine) {
274: fegeom.v = x;
275: fegeom.xi = cgeom->xi;
276: fegeom.J = cgeom->J;
277: fegeom.invJ = cgeom->invJ;
278: fegeom.detJ = cgeom->detJ;
279: if (isCohesiveIn) {
280: fgeomN[0].J = cgeom->suppJ[0];
281: fgeomN[0].invJ = cgeom->suppInvJ[0];
282: fgeomN[0].detJ = cgeom->suppDetJ[0];
283: fgeomN[1].J = cgeom->suppJ[1];
284: fgeomN[1].invJ = cgeom->suppInvJ[1];
285: fgeomN[1].detJ = cgeom->suppDetJ[1];
286: }
287: }
288: for (f = 0, v = 0; f < Nf; ++f) {
289: PetscQuadrature allPoints;
290: PetscInt q, dim, numPoints;
291: const PetscReal *points;
292: PetscScalar *pointEval;
293: PetscBool cohesive;
294: DM dm;
296: if (!sp[f]) continue;
297: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
298: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
299: if (!funcs[f]) {
300: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
301: if (isCohesive && !cohesive) {
302: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
303: }
304: continue;
305: }
306: const PetscInt ***perms;
307: PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
308: PetscCall(PetscDualSpaceGetSymmetries(sp[f], &perms, NULL));
309: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
310: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
311: PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
312: for (q = 0; q < numPoints; ++q, ++tp) {
313: PetscInt qpt[2];
315: if (isCohesiveIn) {
316: qpt[0] = perms ? perms[0][ornt[0]][q] : q;
317: qpt[1] = perms ? perms[0][DMPolytopeTypeComposeOrientationInv(qct, ornt[1], 0)][q] : q;
318: }
319: if (isAffine) {
320: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q * dim], x);
321: } else {
322: fegeom.v = &cgeom->v[tp * dE];
323: fegeom.J = &cgeom->J[tp * dE * dE];
324: fegeom.invJ = &cgeom->invJ[tp * dE * dE];
325: fegeom.detJ = &cgeom->detJ[tp];
326: if (isCohesiveIn) {
327: fgeomN[0].J = &cgeom->suppJ[0][tp * dE * dE];
328: fgeomN[0].invJ = &cgeom->suppInvJ[0][tp * dE * dE];
329: fgeomN[0].detJ = &cgeom->suppDetJ[0][tp];
330: fgeomN[1].J = &cgeom->suppJ[1][tp * dE * dE];
331: fgeomN[1].invJ = &cgeom->suppInvJ[1][tp * dE * dE];
332: fgeomN[1].detJ = &cgeom->suppDetJ[1][tp];
333: }
334: }
335: if (coefficients) {
336: if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, fgeomN, coefficients, coefficients_t, u, u_x, u_t));
337: else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t));
338: }
339: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
340: if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
341: (*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]);
342: }
343: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
344: PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
345: v += spDim;
346: /* TODO: For now, set both sides equal, but this should use info from other support cell */
347: if (isCohesive && !cohesive) {
348: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
349: }
350: }
351: if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
352: if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
353: if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: 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[])
358: {
359: PetscSection section, sectionAux = NULL;
360: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
361: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
362: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
363: const PetscScalar *constants;
364: PetscReal *x;
365: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc, face[2];
366: PetscFEGeom fegeom, cgeom, fgeomN[2];
367: const PetscInt dE = fgeom->dimEmbed, *cone, *ornt;
368: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
369: PetscBool isAffine, isCohesive, isCohesiveIn, transform;
370: DMPolytopeType qct;
372: PetscFunctionBeginHot;
373: PetscCall(PetscDSGetNumFields(ds, &Nf));
374: PetscCall(PetscDSGetComponents(ds, &Nc));
375: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
376: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
377: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
378: PetscCall(PetscDSGetComponentOffsets(dsIn, &uOff));
379: PetscCall(PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x));
380: PetscCall(PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x));
381: PetscCall(PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL));
382: PetscCall(PetscDSGetConstants(dsIn, &numConstants, &constants));
383: PetscCall(DMHasBasisTransform(dmIn, &transform));
384: PetscCall(DMGetLocalSection(dmIn, §ion));
385: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp));
386: // Get cohesive cell hanging off face
387: if (isCohesiveIn) {
388: PetscCall(DMPlexGetCellType(dmIn, inp, &qct));
389: 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)) {
390: DMPolytopeType ct;
391: const PetscInt *support;
392: PetscInt Ns, s;
394: PetscCall(DMPlexGetSupport(dmIn, inp, &support));
395: PetscCall(DMPlexGetSupportSize(dmIn, inp, &Ns));
396: for (s = 0; s < Ns; ++s) {
397: PetscCall(DMPlexGetCellType(dmIn, support[s], &ct));
398: 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)) {
399: inp = support[s];
400: break;
401: }
402: }
403: PetscCheck(s < Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cell not found from face %" PetscInt_FMT, inp);
404: PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 2, &uOff));
405: PetscCall(DMPlexGetOrientedCone(dmIn, inp, &cone, &ornt));
406: face[0] = 0;
407: face[1] = 0;
408: }
409: }
410: if (localU) PetscCall(DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients));
411: if (dmAux) {
412: PetscInt subp;
414: PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp));
415: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
416: PetscCall(DMGetLocalSection(dmAux, §ionAux));
417: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
418: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
419: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x));
420: PetscCall(DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux));
421: }
422: /* Get values for closure */
423: isAffine = fgeom->isAffine;
424: fegeom.dim = fgeom->dim;
425: fegeom.dimEmbed = fgeom->dimEmbed;
426: fegeom.n = NULL;
427: fegeom.J = NULL;
428: fegeom.v = NULL;
429: fegeom.xi = NULL;
430: cgeom.dim = fgeom->dim;
431: cgeom.dimEmbed = fgeom->dimEmbed;
432: if (isCohesiveIn) {
433: fgeomN[0].dim = fgeom->dim;
434: fgeomN[0].dimEmbed = fgeom->dimEmbed;
435: fgeomN[1].dim = fgeom->dim;
436: fgeomN[1].dimEmbed = fgeom->dimEmbed;
437: }
438: if (isAffine) {
439: fegeom.v = x;
440: fegeom.xi = fgeom->xi;
441: fegeom.J = fgeom->J;
442: fegeom.invJ = fgeom->invJ;
443: fegeom.detJ = fgeom->detJ;
444: fegeom.n = fgeom->n;
446: cgeom.J = fgeom->suppJ[0];
447: cgeom.invJ = fgeom->suppInvJ[0];
448: cgeom.detJ = fgeom->suppDetJ[0];
450: if (isCohesiveIn) {
451: fgeomN[0].J = fgeom->suppJ[0];
452: fgeomN[0].invJ = fgeom->suppInvJ[0];
453: fgeomN[0].detJ = fgeom->suppDetJ[0];
454: fgeomN[1].J = fgeom->suppJ[1];
455: fgeomN[1].invJ = fgeom->suppInvJ[1];
456: fgeomN[1].detJ = fgeom->suppDetJ[1];
457: }
458: }
459: for (f = 0, v = 0; f < Nf; ++f) {
460: PetscQuadrature allPoints;
461: PetscInt q, dim, numPoints;
462: const PetscReal *points;
463: PetscScalar *pointEval;
464: PetscBool cohesive;
465: DM dm;
467: if (!sp[f]) continue;
468: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
469: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
470: if (!funcs[f]) {
471: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
472: if (isCohesive && !cohesive) {
473: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
474: }
475: continue;
476: }
477: PetscCall(PetscDualSpaceGetDM(sp[f], &dm));
478: PetscCall(PetscDualSpaceGetAllData(sp[f], &allPoints, NULL));
479: PetscCall(PetscQuadratureGetData(allPoints, &dim, NULL, &numPoints, &points, NULL));
480: PetscCall(DMGetWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
481: for (q = 0; q < numPoints; ++q, ++tp) {
482: PetscInt qpt[2];
484: if (isCohesiveIn) {
485: // These points are not integration quadratures, but dual space quadratures
486: // If they had multiple points we should match them from both sides, similar to hybrid residual eval
487: qpt[0] = qpt[1] = q;
488: }
489: if (isAffine) {
490: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q * dim], x);
491: } else {
492: fegeom.v = &fgeom->v[tp * dE];
493: fegeom.J = &fgeom->J[tp * dE * dE];
494: fegeom.invJ = &fgeom->invJ[tp * dE * dE];
495: fegeom.detJ = &fgeom->detJ[tp];
496: fegeom.n = &fgeom->n[tp * dE];
498: cgeom.J = &fgeom->suppJ[0][tp * dE * dE];
499: cgeom.invJ = &fgeom->suppInvJ[0][tp * dE * dE];
500: cgeom.detJ = &fgeom->suppDetJ[0][tp];
501: if (isCohesiveIn) {
502: fgeomN[0].J = &fgeom->suppJ[0][tp * dE * dE];
503: fgeomN[0].invJ = &fgeom->suppInvJ[0][tp * dE * dE];
504: fgeomN[0].detJ = &fgeom->suppDetJ[0][tp];
505: fgeomN[1].J = &fgeom->suppJ[1][tp * dE * dE];
506: fgeomN[1].invJ = &fgeom->suppInvJ[1][tp * dE * dE];
507: fgeomN[1].detJ = &fgeom->suppDetJ[1][tp];
508: }
509: }
510: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
511: if (coefficients) {
512: if (isCohesiveIn) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, NfIn, 0, tp, T, face, qpt, T, &fegeom, fgeomN, coefficients, coefficients_t, u, u_x, u_t));
513: else PetscCall(PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t));
514: }
515: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t));
516: if (transform) PetscCall(DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx));
517: (*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]);
518: }
519: PetscCall(PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]));
520: PetscCall(DMRestoreWorkArray(dm, numPoints * Nc[f], MPIU_SCALAR, &pointEval));
521: v += spDim;
522: /* TODO: For now, set both sides equal, but this should use info from other support cell */
523: if (isCohesive && !cohesive) {
524: for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim];
525: }
526: }
527: if (localU) PetscCall(DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients));
528: if (dmAux) PetscCall(DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux));
529: if (isCohesiveIn) PetscCall(DMPlexRestoreOrientedCone(dmIn, inp, &cone, &ornt));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: 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[])
534: {
535: PetscFVCellGeom fvgeom;
536: PetscInt dim, dimEmbed;
538: PetscFunctionBeginHot;
539: PetscCall(DMGetDimension(dm, &dim));
540: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
541: if (hasFV) PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL));
542: switch (type) {
543: case DM_BC_ESSENTIAL:
544: case DM_BC_NATURAL:
545: PetscCall(DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))funcs, ctxs, values));
546: break;
547: case DM_BC_ESSENTIAL_FIELD:
548: case DM_BC_NATURAL_FIELD:
549: 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));
550: break;
551: case DM_BC_ESSENTIAL_BD_FIELD:
552: 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));
553: break;
554: default:
555: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int)type);
556: }
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
561: {
562: PetscReal *points;
563: PetscInt f, numPoints;
565: PetscFunctionBegin;
566: if (!dim) {
567: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
570: numPoints = 0;
571: for (f = 0; f < Nf; ++f) {
572: if (funcs[f]) {
573: PetscQuadrature fAllPoints;
574: PetscInt fNumPoints;
576: PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
577: PetscCall(PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL));
578: numPoints += fNumPoints;
579: }
580: }
581: PetscCall(PetscMalloc1(dim * numPoints, &points));
582: numPoints = 0;
583: for (f = 0; f < Nf; ++f) {
584: if (funcs[f]) {
585: PetscQuadrature fAllPoints;
586: PetscInt qdim, fNumPoints, q;
587: const PetscReal *fPoints;
589: PetscCall(PetscDualSpaceGetAllData(sp[f], &fAllPoints, NULL));
590: PetscCall(PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL));
591: 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);
592: for (q = 0; q < fNumPoints * dim; ++q) points[numPoints * dim + q] = fPoints[q];
593: numPoints += fNumPoints;
594: }
595: }
596: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allPoints));
597: PetscCall(PetscQuadratureSetData(*allPoints, dim, 0, numPoints, points, NULL));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*@C
602: DMGetFirstLabeledPoint - Find first labeled `point` in `odm` such that the corresponding point in `dm` has the specified `height`. Return `point` and the corresponding `ds`.
604: Input Parameters:
605: + dm - the `DM`
606: . odm - the enclosing `DM`
607: . label - label for `DM` domain, or `NULL` for whole domain
608: . numIds - the number of `ids`
609: . ids - An array of the label ids in sequence for the domain
610: - height - Height of target cells in `DMPLEX` topology
612: Output Parameters:
613: + point - the first labeled point
614: - ds - the `PetscDS` corresponding to the first labeled point
616: Level: developer
618: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetActivePoint()`, `DMLabel`, `PetscDS`
619: @*/
620: PetscErrorCode DMGetFirstLabeledPoint(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *point, PetscDS *ds)
621: {
622: DM plex;
623: DMEnclosureType enc;
624: PetscInt ls = -1;
626: PetscFunctionBegin;
627: if (point) *point = -1;
628: if (!label) PetscFunctionReturn(PETSC_SUCCESS);
629: PetscCall(DMGetEnclosureRelation(dm, odm, &enc));
630: PetscCall(DMConvert(dm, DMPLEX, &plex));
631: for (PetscInt i = 0; i < numIds; ++i) {
632: IS labelIS;
633: PetscInt num_points, pStart, pEnd;
634: PetscCall(DMLabelGetStratumIS(label, ids[i], &labelIS));
635: if (!labelIS) continue; /* No points with that id on this process */
636: PetscCall(DMPlexGetHeightStratum(plex, height, &pStart, &pEnd));
637: PetscCall(ISGetSize(labelIS, &num_points));
638: if (num_points) {
639: const PetscInt *points;
640: PetscCall(ISGetIndices(labelIS, &points));
641: for (PetscInt i = 0; i < num_points; i++) {
642: PetscInt point;
643: PetscCall(DMGetEnclosurePoint(dm, odm, enc, points[i], &point));
644: if (pStart <= point && point < pEnd) {
645: ls = point;
646: if (ds) {
647: // If this is a face of a cohesive cell, then prefer that DS
648: if (height == 1) {
649: const PetscInt *supp;
650: PetscInt suppSize;
651: DMPolytopeType ct;
653: PetscCall(DMPlexGetSupport(dm, ls, &supp));
654: PetscCall(DMPlexGetSupportSize(dm, ls, &suppSize));
655: for (PetscInt s = 0; s < suppSize; ++s) {
656: PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
657: 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)) {
658: ls = supp[s];
659: break;
660: }
661: }
662: }
663: PetscCall(DMGetCellDS(dm, ls, ds, NULL));
664: }
665: if (ls >= 0) break;
666: }
667: }
668: PetscCall(ISRestoreIndices(labelIS, &points));
669: }
670: PetscCall(ISDestroy(&labelIS));
671: if (ls >= 0) break;
672: }
673: if (point) *point = ls;
674: PetscCall(DMDestroy(&plex));
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: /*
679: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
681: There are several different scenarios:
683: 1) Volumetric mesh with volumetric auxiliary data
685: Here minHeight=0 since we loop over cells.
687: 2) Boundary mesh with boundary auxiliary data
689: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
691: 3) Volumetric mesh with boundary auxiliary data
693: 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.
695: 4) Volumetric input mesh with boundary output mesh
697: Here we must get a subspace for the input DS
699: The maxHeight is used to support enforcement of constraints in DMForest.
701: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
703: 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.
704: - 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
705: 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
706: dual spaces for the boundary from our input spaces.
707: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
709: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
711: 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.
712: */
713: 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)
714: {
715: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
716: DMEnclosureType encIn, encAux;
717: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
718: Vec localA = NULL, tv;
719: IS fieldIS;
720: PetscSection section;
721: PetscDualSpace *sp, *cellsp, *spIn, *cellspIn;
722: PetscTabulation *T = NULL, *TAux = NULL;
723: PetscInt *Nc;
724: PetscInt dim, dimEmbed, depth, pStart, pEnd, lStart = PETSC_DETERMINE, htInc = 0, htIncIn = 0, htIncAux = 0, minHeight, maxHeight, minHeightIn, minHeightAux, h, regionNum, Nf, NfIn, NfAux = 0, NfTot, f;
725: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, isCohesive = PETSC_FALSE, isCohesiveIn = PETSC_FALSE, transform;
726: DMField coordField;
727: DMLabel depthLabel;
728: PetscQuadrature allPoints = NULL;
730: PetscFunctionBegin;
731: if (localU) PetscCall(VecGetDM(localU, &dmIn));
732: else dmIn = dm;
733: PetscCall(DMGetAuxiliaryVec(dm, label, numIds ? ids[0] : 0, 0, &localA));
734: if (localA) PetscCall(VecGetDM(localA, &dmAux));
735: else dmAux = NULL;
736: PetscCall(DMConvert(dm, DMPLEX, &plex));
737: PetscCall(DMConvert(dmIn, DMPLEX, &plexIn));
738: PetscCall(DMGetEnclosureRelation(dmIn, dm, &encIn));
739: PetscCall(DMGetEnclosureRelation(dmAux, dm, &encAux));
740: PetscCall(DMGetDimension(dm, &dim));
741: PetscCall(DMPlexGetVTKCellHeight(plex, &minHeight));
742: PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
743: PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));
744: PetscCall(DMHasBasisTransform(dm, &transform));
745: /* Auxiliary information can only be used with interpolation of field functions */
746: if (dmAux) {
747: PetscCall(DMConvert(dmAux, DMPLEX, &plexAux));
748: 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");
749: }
750: if (localU && localU != localX) PetscCall(DMPlexInsertBoundaryValues(plexIn, PETSC_TRUE, localU, time, NULL, NULL, NULL));
751: PetscCall(DMGetCoordinateField(dm, &coordField));
752: PetscCheck(coordField, PETSC_COMM_SELF, PETSC_ERR_USER, "DM must have a coordinate field");
753: /**** No collective calls below this point ****/
754: /* Determine height for iteration of all meshes */
755: {
756: DMPolytopeType ct, ctIn, ctAux;
757: PetscInt p, pStartIn, pStartAux, pEndAux;
758: PetscInt dim = -1, dimIn = -1, dimAux = -1;
760: PetscCall(DMPlexGetSimplexOrBoxCells(plex, minHeight, &pStart, &pEnd));
761: if (pEnd > pStart) {
762: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, &lStart, NULL));
763: p = lStart < 0 ? pStart : lStart;
764: PetscCall(DMPlexGetCellType(plex, p, &ct));
765: dim = DMPolytopeTypeGetDim(ct);
766: PetscCall(DMPlexGetVTKCellHeight(plexIn, &minHeightIn));
767: PetscCall(DMPlexGetSimplexOrBoxCells(plexIn, minHeightIn, &pStartIn, NULL));
768: PetscCall(DMPlexGetCellType(plexIn, pStartIn, &ctIn));
769: dimIn = DMPolytopeTypeGetDim(ctIn);
770: if (dmAux) {
771: PetscCall(DMPlexGetVTKCellHeight(plexAux, &minHeightAux));
772: PetscCall(DMPlexGetSimplexOrBoxCells(plexAux, minHeightAux, &pStartAux, &pEndAux));
773: if (pStartAux < pEndAux) {
774: PetscCall(DMPlexGetCellType(plexAux, pStartAux, &ctAux));
775: dimAux = DMPolytopeTypeGetDim(ctAux);
776: }
777: } else dimAux = dim;
778: } else {
779: PetscCall(DMDestroy(&plex));
780: PetscCall(DMDestroy(&plexIn));
781: if (dmAux) PetscCall(DMDestroy(&plexAux));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
784: if (dim < 0) {
785: DMLabel spmap = NULL, spmapIn = NULL, spmapAux = NULL;
787: /* Fall back to determination based on being a submesh */
788: PetscCall(DMPlexGetSubpointMap(plex, &spmap));
789: PetscCall(DMPlexGetSubpointMap(plexIn, &spmapIn));
790: if (plexAux) PetscCall(DMPlexGetSubpointMap(plexAux, &spmapAux));
791: dim = spmap ? 1 : 0;
792: dimIn = spmapIn ? 1 : 0;
793: dimAux = spmapAux ? 1 : 0;
794: }
795: {
796: PetscInt dimProj = PetscMin(PetscMin(dim, dimIn), dimAux < 0 ? PETSC_INT_MAX : dimAux);
797: PetscInt dimAuxEff = dimAux < 0 ? dimProj : dimAux;
799: 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");
800: if (dimProj < dim) minHeight = 1;
801: htInc = dim - dimProj;
802: htIncIn = dimIn - dimProj;
803: htIncAux = dimAuxEff - dimProj;
804: }
805: }
806: PetscCall(DMPlexGetDepth(plex, &depth));
807: PetscCall(DMPlexGetDepthLabel(plex, &depthLabel));
808: PetscCall(DMPlexGetMaxProjectionHeight(plex, &maxHeight));
809: maxHeight = PetscMax(maxHeight, minHeight);
810: PetscCheck(maxHeight >= 0 && maxHeight <= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", maxHeight, dim);
811: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, minHeight, NULL, &ds));
812: if (!ds) PetscCall(DMGetDS(dm, &ds));
813: PetscCall(DMGetFirstLabeledPoint(dmIn, dm, label, numIds, ids, minHeight, NULL, &dsIn));
814: if (!dsIn) PetscCall(DMGetDS(dmIn, &dsIn));
815: if (!dsIn) {
816: if (encIn == DM_ENC_SUPERMESH) {
817: PetscInt p = pStart, pIn;
819: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? p : lStart, &pIn));
820: // If the input mesh is higher dimensional than the output mesh, get a cell from the output mesh
821: if (htIncIn) PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &p, NULL));
822: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? p : lStart, &pIn));
823: PetscCall(DMGetCellDS(dmIn, pIn, &dsIn, NULL));
824: } else PetscCall(DMGetDS(dmIn, &dsIn));
825: }
826: PetscCall(PetscDSGetNumFields(ds, &Nf));
827: PetscCall(PetscDSGetNumFields(dsIn, &NfIn));
828: PetscCall(PetscDSIsCohesive(dsIn, &isCohesiveIn));
829: if (isCohesiveIn) --htIncIn; // Should be rearranged
830: PetscCall(DMGetNumFields(dm, &NfTot));
831: PetscCall(DMFindRegionNum(dm, ds, ®ionNum));
832: PetscCall(DMGetRegionNumDS(dm, regionNum, NULL, &fieldIS, NULL, NULL));
833: PetscCall(PetscDSIsCohesive(ds, &isCohesive));
834: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
835: PetscCall(DMGetLocalSection(dm, §ion));
836: if (dmAux) {
837: if (encAux == DM_ENC_SUPERMESH) {
838: PetscInt p = pStart, pAux;
840: // If the auxiliary mesh is higher dimensional than the output mesh, get a cell from the output mesh
841: if (htIncAux) PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &p, NULL));
842: PetscCall(DMGetEnclosurePoint(dmAux, dm, encAux, lStart < 0 ? p : lStart, &pAux));
843: PetscCall(DMGetCellDS(dmAux, pAux, &dsAux, NULL));
844: } else PetscCall(DMGetDS(dmAux, &dsAux));
845: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
846: }
847: PetscCall(PetscDSGetComponents(ds, &Nc));
848: PetscCall(PetscMalloc3(Nf, &isFE, Nf, &sp, NfIn, &spIn));
849: if (maxHeight > 0) PetscCall(PetscMalloc2(Nf, &cellsp, NfIn, &cellspIn));
850: else {
851: cellsp = sp;
852: cellspIn = spIn;
853: }
854: /* Get cell dual spaces */
855: for (f = 0; f < Nf; ++f) {
856: PetscDiscType disctype;
858: PetscCall(PetscDSGetDiscType_Internal(ds, f, &disctype));
859: if (disctype == PETSC_DISC_FE) {
860: PetscFE fe;
862: isFE[f] = PETSC_TRUE;
863: hasFE = PETSC_TRUE;
864: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
865: PetscCall(PetscFEGetDualSpace(fe, &cellsp[f]));
866: } else if (disctype == PETSC_DISC_FV) {
867: PetscFV fv;
869: isFE[f] = PETSC_FALSE;
870: hasFV = PETSC_TRUE;
871: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fv));
872: PetscCall(PetscFVGetDualSpace(fv, &cellsp[f]));
873: } else {
874: isFE[f] = PETSC_FALSE;
875: cellsp[f] = NULL;
876: }
877: }
878: for (f = 0; f < NfIn; ++f) {
879: PetscDiscType disctype;
881: PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
882: if (disctype == PETSC_DISC_FE) {
883: PetscFE fe;
885: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fe));
886: PetscCall(PetscFEGetDualSpace(fe, &cellspIn[f]));
887: } else if (disctype == PETSC_DISC_FV) {
888: PetscFV fv;
890: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fv));
891: PetscCall(PetscFVGetDualSpace(fv, &cellspIn[f]));
892: } else {
893: cellspIn[f] = NULL;
894: }
895: }
896: for (f = 0; f < Nf; ++f) {
897: if (!htInc) {
898: sp[f] = cellsp[f];
899: } else PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], htInc, &sp[f]));
900: }
901: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
902: PetscFE fem, subfem;
903: PetscDiscType disctype;
904: const PetscReal *points;
905: PetscInt numPoints, k;
907: PetscCheck(maxHeight <= minHeight, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
908: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &allPoints));
909: PetscCall(PetscQuadratureGetData(allPoints, NULL, NULL, &numPoints, &points, NULL));
910: PetscCall(PetscMalloc2(NfIn, &T, NfAux, &TAux));
911: for (f = 0; f < NfIn; ++f) {
912: if (!htIncIn) {
913: spIn[f] = cellspIn[f];
914: } else PetscCall(PetscDualSpaceGetHeightSubspace(cellspIn[f], htIncIn, &spIn[f]));
916: PetscCall(PetscDSGetDiscType_Internal(dsIn, f, &disctype));
917: if (disctype != PETSC_DISC_FE) continue;
918: PetscCall(PetscDSGetDiscretization(dsIn, f, (PetscObject *)&fem));
919: if (!htIncIn) {
920: subfem = fem;
921: } else PetscCall(PetscFEGetHeightSubspace(fem, htIncIn, &subfem));
922: PetscCall(PetscDSGetJetDegree(dsIn, f, &k));
923: PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, k, &T[f]));
924: }
925: for (f = 0; f < NfAux; ++f) {
926: PetscCall(PetscDSGetDiscType_Internal(dsAux, f, &disctype));
927: if (disctype != PETSC_DISC_FE) continue;
928: PetscCall(PetscDSGetDiscretization(dsAux, f, (PetscObject *)&fem));
929: if (!htIncAux) {
930: subfem = fem;
931: } else PetscCall(PetscFEGetHeightSubspace(fem, htIncAux, &subfem));
932: PetscCall(PetscDSGetJetDegree(dsAux, f, &k));
933: PetscCall(PetscFECreateTabulation(subfem, 1, numPoints, points, k, &TAux[f]));
934: }
935: }
936: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
937: for (h = minHeight; h <= maxHeight; h++) {
938: PetscInt hEff = h - minHeight + htInc;
939: PetscInt hEffIn = h - minHeight + htIncIn;
940: PetscInt hEffAux = h - minHeight + htIncAux;
941: PetscDS dsEff = ds;
942: PetscDS dsEffIn = dsIn;
943: PetscDS dsEffAux = dsAux;
944: PetscScalar *values;
945: PetscBool *fieldActive;
946: PetscInt maxDegree;
947: PetscInt p, spDim, totDim, numValues;
948: IS heightIS;
950: if (h > minHeight) {
951: for (f = 0; f < Nf; ++f) PetscCall(PetscDualSpaceGetHeightSubspace(cellsp[f], hEff, &sp[f]));
952: }
953: PetscCall(DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd));
954: PetscCall(DMGetFirstLabeledPoint(dm, dm, label, numIds, ids, h, &lStart, NULL));
955: PetscCall(DMLabelGetStratumIS(depthLabel, depth - h, &heightIS));
956: if (pEnd <= pStart) {
957: PetscCall(ISDestroy(&heightIS));
958: continue;
959: }
960: /* Compute totDim, the number of dofs in the closure of a point at this height */
961: totDim = 0;
962: for (f = 0; f < Nf; ++f) {
963: PetscBool cohesive;
965: if (!sp[f]) continue;
966: PetscCall(PetscDSGetCohesive(ds, f, &cohesive));
967: PetscCall(PetscDualSpaceGetDimension(sp[f], &spDim));
968: totDim += spDim;
969: if (isCohesive && !cohesive) totDim += spDim;
970: }
971: p = lStart < 0 ? pStart : lStart;
972: PetscCall(DMPlexVecGetClosure(plex, section, localX, p, &numValues, NULL));
973: 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);
974: if (!totDim) {
975: PetscCall(ISDestroy(&heightIS));
976: continue;
977: }
978: if (htInc) PetscCall(PetscDSGetHeightSubspace(ds, hEff, &dsEff));
979: /* Compute totDimIn, the number of dofs in the closure of a point at this height */
980: if (localU) {
981: PetscInt totDimIn, pIn, numValuesIn;
983: totDimIn = 0;
984: for (f = 0; f < NfIn; ++f) {
985: PetscBool cohesive;
987: if (!spIn[f]) continue;
988: PetscCall(PetscDSGetCohesive(dsIn, f, &cohesive));
989: PetscCall(PetscDualSpaceGetDimension(spIn[f], &spDim));
990: totDimIn += spDim;
991: if (isCohesiveIn && !cohesive) totDimIn += spDim;
992: }
993: PetscCall(DMGetEnclosurePoint(dmIn, dm, encIn, lStart < 0 ? pStart : lStart, &pIn));
994: PetscCall(DMPlexVecGetClosure(plexIn, NULL, localU, pIn, &numValuesIn, NULL));
995: // TODO We could check that pIn is a cohesive cell for this check
996: 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);
997: if (htIncIn) PetscCall(PetscDSGetHeightSubspace(dsIn, hEffIn, &dsEffIn));
998: }
999: if (htIncAux) PetscCall(PetscDSGetHeightSubspace(dsAux, hEffAux, &dsEffAux));
1000: /* Loop over points at this height */
1001: PetscCall(DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values));
1002: PetscCall(DMGetWorkArray(dm, NfTot, MPI_INT, &fieldActive));
1003: {
1004: const PetscInt *fields;
1006: PetscCall(ISGetIndices(fieldIS, &fields));
1007: for (f = 0; f < NfTot; ++f) fieldActive[f] = PETSC_FALSE;
1008: for (f = 0; f < Nf; ++f) fieldActive[fields[f]] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
1009: PetscCall(ISRestoreIndices(fieldIS, &fields));
1010: }
1011: if (label) {
1012: PetscInt i;
1014: for (i = 0; i < numIds; ++i) {
1015: IS pointIS, isectIS;
1016: const PetscInt *points;
1017: PetscInt n;
1018: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
1019: PetscQuadrature quad = NULL;
1021: PetscCall(DMLabelGetStratumIS(label, ids[i], &pointIS));
1022: if (!pointIS) continue; /* No points with that id on this process */
1023: PetscCall(ISIntersect(pointIS, heightIS, &isectIS));
1024: PetscCall(ISDestroy(&pointIS));
1025: if (!isectIS) continue;
1026: PetscCall(ISGetLocalSize(isectIS, &n));
1027: PetscCall(ISGetIndices(isectIS, &points));
1028: PetscCall(DMFieldGetDegree(coordField, isectIS, NULL, &maxDegree));
1029: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, isectIS, &quad));
1030: if (!quad) {
1031: if (!h && allPoints) {
1032: quad = allPoints;
1033: allPoints = NULL;
1034: } else {
1035: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, isCohesive ? dim - htInc - 1 : dim - htInc, funcs, &quad));
1036: }
1037: }
1038: PetscFEGeomMode geommode = htInc && h == minHeight ? PETSC_FEGEOM_BOUNDARY : PETSC_FEGEOM_BASIC;
1040: if (n) {
1041: PetscInt depth, dep;
1043: PetscCall(DMPlexGetDepth(dm, &depth));
1044: PetscCall(DMPlexGetPointDepth(dm, points[0], &dep));
1045: if (dep < depth && h == minHeight) geommode = PETSC_FEGEOM_BOUNDARY;
1046: }
1047: PetscCall(DMFieldCreateFEGeom(coordField, isectIS, quad, geommode, &fegeom));
1048: for (p = 0; p < n; ++p) {
1049: const PetscInt point = points[p];
1051: PetscCall(PetscArrayzero(values, numValues));
1052: PetscCall(PetscFEGeomGetChunk(fegeom, p, p + 1, &chunkgeom));
1053: PetscCall(DMPlexSetActivePoint(dm, point));
1054: 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));
1055: if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values));
1056: PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode));
1057: }
1058: PetscCall(PetscFEGeomRestoreChunk(fegeom, p, p + 1, &chunkgeom));
1059: PetscCall(PetscFEGeomDestroy(&fegeom));
1060: PetscCall(PetscQuadratureDestroy(&quad));
1061: PetscCall(ISRestoreIndices(isectIS, &points));
1062: PetscCall(ISDestroy(&isectIS));
1063: }
1064: } else {
1065: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
1066: PetscQuadrature quad = NULL;
1067: IS pointIS;
1069: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pointIS));
1070: PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
1071: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
1072: if (!quad) {
1073: if (!h && allPoints) {
1074: quad = allPoints;
1075: allPoints = NULL;
1076: } else {
1077: PetscCall(PetscDualSpaceGetAllPointsUnion(Nf, sp, dim - htInc, funcs, &quad));
1078: }
1079: }
1080: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, (htInc && h == minHeight) ? PETSC_FEGEOM_BOUNDARY : PETSC_FEGEOM_BASIC, &fegeom));
1081: for (p = pStart; p < pEnd; ++p) {
1082: PetscCall(PetscArrayzero(values, numValues));
1083: PetscCall(PetscFEGeomGetChunk(fegeom, p - pStart, p - pStart + 1, &chunkgeom));
1084: PetscCall(DMPlexSetActivePoint(dm, p));
1085: 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));
1086: if (transform) PetscCall(DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values));
1087: PetscCall(DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode));
1088: }
1089: PetscCall(PetscFEGeomRestoreChunk(fegeom, p - pStart, pStart - p + 1, &chunkgeom));
1090: PetscCall(PetscFEGeomDestroy(&fegeom));
1091: PetscCall(PetscQuadratureDestroy(&quad));
1092: PetscCall(ISDestroy(&pointIS));
1093: }
1094: PetscCall(ISDestroy(&heightIS));
1095: PetscCall(DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values));
1096: PetscCall(DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive));
1097: }
1098: /* Cleanup */
1099: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
1100: for (f = 0; f < NfIn; ++f) PetscCall(PetscTabulationDestroy(&T[f]));
1101: for (f = 0; f < NfAux; ++f) PetscCall(PetscTabulationDestroy(&TAux[f]));
1102: PetscCall(PetscFree2(T, TAux));
1103: }
1104: PetscCall(PetscQuadratureDestroy(&allPoints));
1105: PetscCall(PetscFree3(isFE, sp, spIn));
1106: if (maxHeight > 0) PetscCall(PetscFree2(cellsp, cellspIn));
1107: PetscCall(DMDestroy(&plex));
1108: PetscCall(DMDestroy(&plexIn));
1109: if (dmAux) PetscCall(DMDestroy(&plexAux));
1110: PetscFunctionReturn(PETSC_SUCCESS);
1111: }
1113: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1114: {
1115: PetscFunctionBegin;
1116: PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: 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)
1121: {
1122: PetscFunctionBegin;
1123: PetscCall(DMProjectLocal_Generic_Plex(dm, time, NULL, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void))funcs, ctxs, mode, localX));
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: 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)
1128: {
1129: PetscFunctionBegin;
1130: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: 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)
1135: {
1136: PetscFunctionBegin;
1137: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: 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)
1142: {
1143: PetscFunctionBegin;
1144: PetscCall(DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void))funcs, NULL, mode, localX));
1145: PetscFunctionReturn(PETSC_SUCCESS);
1146: }