Actual source code: plexproject.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: /*
6: DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point
8: Input Parameters:
9: + dm - The output DM
10: . ds - The output DS
11: . dmIn - The input DM
12: . dsIn - The input DS
13: . time - The time for this evaluation
14: . fegeom - The FE geometry for this point
15: . fvgeom - The FV geometry for this point
16: . isFE - Flag indicating whether each output field has an FE discretization
17: . sp - The output PetscDualSpace for each field
18: . funcs - The evaluation function for each field
19: - ctxs - The user context for each field
21: Output Parameter:
22: . values - The value for each dual basis vector in the output dual space
24: Level: developer
26: .seealso: DMProjectPoint_Field_Private()
27: */
28: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[],
29: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
30: PetscScalar values[])
31: {
32: PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp;
33: PetscBool isAffine, transform;
37: DMGetCoordinateDim(dmIn, &coordDim);
38: DMHasBasisTransform(dmIn, &transform);
39: PetscDSGetNumFields(ds, &Nf);
40: PetscDSGetComponents(ds, &Nc);
41: /* Get values for closure */
42: isAffine = fegeom->isAffine;
43: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
44: void * const ctx = ctxs ? ctxs[f] : NULL;
46: if (!sp[f]) continue;
47: PetscDualSpaceGetDimension(sp[f], &spDim);
48: if (funcs[f]) {
49: if (isFE[f]) {
50: PetscQuadrature allPoints;
51: PetscInt q, dim, numPoints;
52: const PetscReal *points;
53: PetscScalar *pointEval;
54: PetscReal *x;
55: DM rdm;
57: PetscDualSpaceGetDM(sp[f],&rdm);
58: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
59: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
60: DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
61: DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
62: for (q = 0; q < numPoints; q++, tp++) {
63: const PetscReal *v0;
65: if (isAffine) {
66: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, &points[q*dim], x);
67: v0 = x;
68: } else {
69: v0 = &fegeom->v[tp*coordDim];
70: }
71: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
72: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
73: }
74: /* Transform point evaluations pointEval[q,c] */
75: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
76: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
77: DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
78: DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
79: v += spDim;
80: } else {
81: for (d = 0; d < spDim; ++d, ++v) {
82: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
83: }
84: }
85: } else {
86: for (d = 0; d < spDim; d++, v++) {values[v] = 0.0;}
87: }
88: }
89: return(0);
90: }
92: /*
93: DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point
95: Input Parameters:
96: + dm - The output DM
97: . ds - The output DS
98: . dmIn - The input DM
99: . dsIn - The input DS
100: . dmAux - The auxiliary DM, which is always for the input space
101: . dsAux - The auxiliary DS, which is always for the input space
102: . time - The time for this evaluation
103: . localU - The local solution
104: . localA - The local auziliary fields
105: . cgeom - The FE geometry for this point
106: . sp - The output PetscDualSpace for each field
107: . p - The point in the output DM
108: . T - Input basis and derviatives for each field tabulated on the quadrature points
109: . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points
110: . funcs - The evaluation function for each field
111: - ctxs - The user context for each field
113: Output Parameter:
114: . values - The value for each dual basis vector in the output dual space
116: Note: Not supported for FV
118: Level: developer
120: .seealso: DMProjectPoint_Field_Private()
121: */
122: 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,
123: PetscTabulation *T, PetscTabulation *TAux,
124: void (**funcs)(PetscInt, PetscInt, PetscInt,
125: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
126: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
127: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
128: PetscScalar values[])
129: {
130: PetscSection section, sectionAux = NULL;
131: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
132: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
133: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
134: const PetscScalar *constants;
135: PetscReal *x;
136: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
137: PetscFEGeom fegeom;
138: const PetscInt dE = cgeom->dimEmbed;
139: PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0;
140: PetscBool isAffine, transform;
141: PetscErrorCode ierr;
144: PetscDSGetNumFields(ds, &Nf);
145: PetscDSGetComponents(ds, &Nc);
146: PetscDSGetNumFields(dsIn, &NfIn);
147: PetscDSGetComponentOffsets(dsIn, &uOff);
148: PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);
149: PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);
150: PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);
151: PetscDSGetConstants(dsIn, &numConstants, &constants);
152: DMHasBasisTransform(dmIn, &transform);
153: DMGetLocalSection(dmIn, §ion);
154: DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);
155: DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);
156: if (dmAux) {
157: PetscInt subp;
159: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
160: PetscDSGetNumFields(dsAux, &NfAux);
161: DMGetLocalSection(dmAux, §ionAux);
162: PetscDSGetComponentOffsets(dsAux, &aOff);
163: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
164: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
165: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
166: }
167: /* Get values for closure */
168: isAffine = cgeom->isAffine;
169: if (isAffine) {
170: fegeom.v = x;
171: fegeom.xi = cgeom->xi;
172: fegeom.J = cgeom->J;
173: fegeom.invJ = cgeom->invJ;
174: fegeom.detJ = cgeom->detJ;
175: }
176: for (f = 0, v = 0; f < Nf; ++f) {
177: PetscQuadrature allPoints;
178: PetscInt q, dim, numPoints;
179: const PetscReal *points;
180: PetscScalar *pointEval;
181: DM dm;
183: if (!sp[f]) continue;
184: PetscDualSpaceGetDimension(sp[f], &spDim);
185: if (!funcs[f]) {
186: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
187: continue;
188: }
189: PetscDualSpaceGetDM(sp[f],&dm);
190: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
191: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
192: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
193: for (q = 0; q < numPoints; ++q, ++tp) {
194: if (isAffine) {
195: CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x);
196: } else {
197: fegeom.v = &cgeom->v[tp*dE];
198: fegeom.J = &cgeom->J[tp*dE*dE];
199: fegeom.invJ = &cgeom->invJ[tp*dE*dE];
200: fegeom.detJ = &cgeom->detJ[tp];
201: }
202: PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
203: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
204: if (transform) {DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);}
205: (*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]);
206: }
207: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
208: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
209: v += spDim;
210: }
211: DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);
212: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
213: return(0);
214: }
216: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p,
217: PetscTabulation *T, PetscTabulation *TAux,
218: void (**funcs)(PetscInt, PetscInt, PetscInt,
219: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
220: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
221: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
222: PetscScalar values[])
223: {
224: PetscSection section, sectionAux = NULL;
225: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
226: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
227: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
228: const PetscScalar *constants;
229: PetscReal *x;
230: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc;
231: PetscFEGeom fegeom, cgeom;
232: const PetscInt dE = fgeom->dimEmbed;
233: PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
234: PetscBool isAffine;
235: PetscErrorCode ierr;
238: if (dm != dmIn) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM");
239: PetscDSGetNumFields(ds, &Nf);
240: PetscDSGetComponents(ds, &Nc);
241: PetscDSGetComponentOffsets(ds, &uOff);
242: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
243: PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);
244: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
245: PetscDSGetConstants(ds, &numConstants, &constants);
246: DMGetLocalSection(dm, §ion);
247: DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);
248: if (dmAux) {
249: PetscInt subp;
251: DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);
252: PetscDSGetNumFields(dsAux, &NfAux);
253: DMGetLocalSection(dmAux, §ionAux);
254: PetscDSGetComponentOffsets(dsAux, &aOff);
255: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
256: PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);
257: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
258: }
259: /* Get values for closure */
260: isAffine = fgeom->isAffine;
261: fegeom.n = 0;
262: fegeom.J = 0;
263: fegeom.v = 0;
264: fegeom.xi = 0;
265: cgeom.dim = fgeom->dim;
266: cgeom.dimEmbed = fgeom->dimEmbed;
267: if (isAffine) {
268: fegeom.v = x;
269: fegeom.xi = fgeom->xi;
270: fegeom.J = fgeom->J;
271: fegeom.invJ = fgeom->invJ;
272: fegeom.detJ = fgeom->detJ;
273: fegeom.n = fgeom->n;
275: cgeom.J = fgeom->suppJ[0];
276: cgeom.invJ = fgeom->suppInvJ[0];
277: cgeom.detJ = fgeom->suppDetJ[0];
278: }
279: for (f = 0, v = 0; f < Nf; ++f) {
280: PetscQuadrature allPoints;
281: PetscInt q, dim, numPoints;
282: const PetscReal *points;
283: PetscScalar *pointEval;
284: DM dm;
286: if (!sp[f]) continue;
287: PetscDualSpaceGetDimension(sp[f], &spDim);
288: if (!funcs[f]) {
289: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
290: continue;
291: }
292: PetscDualSpaceGetDM(sp[f],&dm);
293: PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);
294: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
295: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
296: for (q = 0; q < numPoints; ++q, ++tp) {
297: if (isAffine) {
298: CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
299: } else {
300: fegeom.v = &fgeom->v[tp*dE];
301: fegeom.J = &fgeom->J[tp*dE*dE];
302: fegeom.invJ = &fgeom->invJ[tp*dE*dE];
303: fegeom.detJ = &fgeom->detJ[tp];
304: fegeom.n = &fgeom->n[tp*dE];
306: cgeom.J = &fgeom->suppJ[0][tp*dE*dE];
307: cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE];
308: cgeom.detJ = &fgeom->suppDetJ[0][tp];
309: }
310: /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */
311: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
312: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
313: (*funcs[f])(dE, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, fegeom.v, fegeom.n, numConstants, constants, &pointEval[Nc[f]*q]);
314: }
315: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
316: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
317: v += spDim;
318: }
319: DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);
320: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
321: return(0);
322: }
324: 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[],
325: PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux,
326: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
327: {
328: PetscFVCellGeom fvgeom;
329: PetscInt dim, dimEmbed;
330: PetscErrorCode ierr;
333: DMGetDimension(dm, &dim);
334: DMGetCoordinateDim(dm, &dimEmbed);
335: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
336: switch (type) {
337: case DM_BC_ESSENTIAL:
338: case DM_BC_NATURAL:
339: DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
340: case DM_BC_ESSENTIAL_FIELD:
341: case DM_BC_NATURAL_FIELD:
342: DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
343: (void (**)(PetscInt, PetscInt, PetscInt,
344: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
345: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
346: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
347: case DM_BC_ESSENTIAL_BD_FIELD:
348: DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux,
349: (void (**)(PetscInt, PetscInt, PetscInt,
350: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
351: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
352: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
353: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
354: }
355: return(0);
356: }
358: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
359: {
360: PetscReal *points;
361: PetscInt f, numPoints;
365: numPoints = 0;
366: for (f = 0; f < Nf; ++f) {
367: if (funcs[f]) {
368: PetscQuadrature fAllPoints;
369: PetscInt fNumPoints;
371: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
372: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
373: numPoints += fNumPoints;
374: }
375: }
376: PetscMalloc1(dim*numPoints,&points);
377: numPoints = 0;
378: for (f = 0; f < Nf; ++f) {
379: PetscInt spDim;
381: PetscDualSpaceGetDimension(sp[f], &spDim);
382: if (funcs[f]) {
383: PetscQuadrature fAllPoints;
384: PetscInt qdim, fNumPoints, q;
385: const PetscReal *fPoints;
387: PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);
388: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
389: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Spatial dimension %D for dual basis does not match input dimension %D", qdim, dim);
390: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
391: numPoints += fNumPoints;
392: }
393: }
394: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
395: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
396: return(0);
397: }
399: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *ds)
400: {
401: DM plex;
402: DMEnclosureType enc;
403: DMLabel depthLabel;
404: PetscInt dim, cdepth, ls = -1, i;
405: PetscErrorCode ierr;
408: if (lStart) *lStart = -1;
409: if (!label) return(0);
410: DMGetEnclosureRelation(dm, odm, &enc);
411: DMGetDimension(dm, &dim);
412: DMConvert(dm, DMPLEX, &plex);
413: DMPlexGetDepthLabel(plex, &depthLabel);
414: cdepth = dim - height;
415: for (i = 0; i < numIds; ++i) {
416: IS pointIS;
417: const PetscInt *points;
418: PetscInt pdepth, point;
420: DMLabelGetStratumIS(label, ids[i], &pointIS);
421: if (!pointIS) continue; /* No points with that id on this process */
422: ISGetIndices(pointIS, &points);
423: DMGetEnclosurePoint(dm, odm, enc, points[0], &point);
424: DMLabelGetValue(depthLabel, point, &pdepth);
425: if (pdepth == cdepth) {
426: ls = point;
427: if (ds) {DMGetCellDS(dm, ls, ds);}
428: }
429: ISRestoreIndices(pointIS, &points);
430: ISDestroy(&pointIS);
431: if (ls >= 0) break;
432: }
433: if (lStart) *lStart = ls;
434: DMDestroy(&plex);
435: return(0);
436: }
438: /*
439: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
441: There are several different scenarios:
443: 1) Volumetric mesh with volumetric auxiliary data
445: Here minHeight=0 since we loop over cells.
447: 2) Boundary mesh with boundary auxiliary data
449: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
451: 3) Volumetric mesh with boundary auxiliary data
453: 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.
455: The maxHeight is used to support enforcement of constraints in DMForest.
457: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
459: 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.
460: - 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
461: is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When th effective height is nonzero, we need to extract
462: dual spaces for the boundary from our input spaces.
463: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
465: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
467: 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.
468: */
469: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
470: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
471: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
472: InsertMode mode, Vec localX)
473: {
474: DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm;
475: DMEnclosureType encIn, encAux;
476: PetscDS ds = NULL, dsIn = NULL, dsAux = NULL;
477: Vec localA = NULL, tv;
478: PetscSection section;
479: PetscDualSpace *sp, *cellsp;
480: PetscTabulation *T = NULL, *TAux = NULL;
481: PetscInt *Nc;
482: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfIn, NfAux = 0, f;
483: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, transform;
484: DMField coordField;
485: DMLabel depthLabel;
486: PetscQuadrature allPoints = NULL;
487: PetscErrorCode ierr;
490: if (localU) {VecGetDM(localU, &dmIn);}
491: else {dmIn = dm;}
492: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
493: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
494: DMConvert(dm, DMPLEX, &plex);
495: DMConvert(dmIn, DMPLEX, &plexIn);
496: DMGetEnclosureRelation(dmIn, dm, &encIn);
497: DMGetEnclosureRelation(dmAux, dm, &encAux);
498: DMGetDimension(dm, &dim);
499: DMPlexGetVTKCellHeight(plex, &minHeight);
500: DMGetBasisTransformDM_Internal(dm, &tdm);
501: DMGetBasisTransformVec_Internal(dm, &tv);
502: DMHasBasisTransform(dm, &transform);
503: /* Auxiliary information can only be used with interpolation of field functions */
504: if (dmAux) {
505: DMConvert(dmAux, DMPLEX, &plexAux);
506: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) {
507: if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
508: if (!minHeight) {
509: DMLabel spmap;
511: /* If dmAux is a surface, then force the projection to take place over a surface */
512: DMPlexGetSubpointMap(plexAux, &spmap);
513: if (spmap) {
514: DMPlexGetVTKCellHeight(plexAux, &minHeight);
515: auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
516: }
517: }
518: }
519: }
520: DMPlexGetDepth(plex, &depth);
521: DMPlexGetDepthLabel(plex, &depthLabel);
522: DMPlexGetMaxProjectionHeight(plex, &maxHeight);
523: maxHeight = PetscMax(maxHeight, minHeight);
524: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
525: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);
526: if (!ds) {DMGetDS(dm, &ds);}
527: DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);
528: if (!dsIn) {DMGetDS(dmIn, &dsIn);}
529: PetscDSGetNumFields(ds, &Nf);
530: PetscDSGetNumFields(dsIn, &NfIn);
531: DMGetCoordinateDim(dm, &dimEmbed);
532: DMGetLocalSection(dm, §ion);
533: if (dmAux) {
534: DMGetDS(dmAux, &dsAux);
535: PetscDSGetNumFields(dsAux, &NfAux);
536: }
537: PetscDSGetComponents(ds, &Nc);
538: PetscMalloc2(Nf, &isFE, Nf, &sp);
539: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
540: else {cellsp = sp;}
541: if (localU && localU != localX) {DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
542: /* Get cell dual spaces */
543: for (f = 0; f < Nf; ++f) {
544: PetscObject obj;
545: PetscClassId id;
547: PetscDSGetDiscretization(ds, f, &obj);
548: PetscObjectGetClassId(obj, &id);
549: if (id == PETSCFE_CLASSID) {
550: PetscFE fe = (PetscFE) obj;
552: hasFE = PETSC_TRUE;
553: isFE[f] = PETSC_TRUE;
554: PetscFEGetDualSpace(fe, &cellsp[f]);
555: } else if (id == PETSCFV_CLASSID) {
556: PetscFV fv = (PetscFV) obj;
558: hasFV = PETSC_TRUE;
559: isFE[f] = PETSC_FALSE;
560: PetscFVGetDualSpace(fv, &cellsp[f]);
561: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
562: }
563: DMGetCoordinateField(dm,&coordField);
564: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
565: PetscInt effectiveHeight = auxBd ? minHeight : 0;
566: PetscFE fem, subfem;
567: PetscBool isfe;
568: const PetscReal *points;
569: PetscInt numPoints;
571: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
572: for (f = 0; f < Nf; ++f) {
573: if (!effectiveHeight) {sp[f] = cellsp[f];}
574: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
575: }
576: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
577: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
578: PetscMalloc2(NfIn, &T, NfAux, &TAux);
579: for (f = 0; f < NfIn; ++f) {
580: PetscDSIsFE_Internal(dsIn, f, &isfe);
581: if (!isfe) continue;
582: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
583: if (!effectiveHeight) {subfem = fem;}
584: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
585: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);
586: }
587: for (f = 0; f < NfAux; ++f) {
588: PetscDSIsFE_Internal(dsAux, f, &isfe);
589: if (!isfe) continue;
590: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
591: if (!effectiveHeight || auxBd) {subfem = fem;}
592: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
593: PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);
594: }
595: }
596: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
597: for (h = minHeight; h <= maxHeight; h++) {
598: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
599: PetscDS dsEff = ds;
600: PetscScalar *values;
601: PetscBool *fieldActive;
602: PetscInt maxDegree;
603: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
604: IS heightIS;
606: /* Note we assume that dm and dmIn share the same topology */
607: DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);
608: DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);
609: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
610: if (pEnd <= pStart) {
611: ISDestroy(&heightIS);
612: continue;
613: }
614: /* Compute totDim, the number of dofs in the closure of a point at this height */
615: totDim = 0;
616: for (f = 0; f < Nf; ++f) {
617: if (!effectiveHeight) {
618: sp[f] = cellsp[f];
619: } else {
620: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
621: if (!sp[f]) continue;
622: }
623: PetscDualSpaceGetDimension(sp[f], &spDim);
624: totDim += spDim;
625: }
626: DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
627: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim);
628: if (!totDim) {
629: ISDestroy(&heightIS);
630: continue;
631: }
632: if (effectiveHeight) {PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);}
633: /* Loop over points at this height */
634: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
635: DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);
636: for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
637: if (label) {
638: PetscInt i;
640: for (i = 0; i < numIds; ++i) {
641: IS pointIS, isectIS;
642: const PetscInt *points;
643: PetscInt n;
644: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
645: PetscQuadrature quad = NULL;
647: DMLabelGetStratumIS(label, ids[i], &pointIS);
648: if (!pointIS) continue; /* No points with that id on this process */
649: ISIntersect(pointIS,heightIS,&isectIS);
650: ISDestroy(&pointIS);
651: if (!isectIS) continue;
652: ISGetLocalSize(isectIS, &n);
653: ISGetIndices(isectIS, &points);
654: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
655: if (maxDegree <= 1) {
656: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
657: }
658: if (!quad) {
659: if (!h && allPoints) {
660: quad = allPoints;
661: allPoints = NULL;
662: } else {
663: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
664: }
665: }
666: DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
667: for (p = 0; p < n; ++p) {
668: const PetscInt point = points[p];
670: PetscArrayzero(values, numValues);
671: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
672: DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values);
673: if (ierr) {
674: PetscErrorCode ierr2;
675: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
676: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
677:
678: }
679: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
680: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, values, mode);
681: }
682: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
683: PetscFEGeomDestroy(&fegeom);
684: PetscQuadratureDestroy(&quad);
685: ISRestoreIndices(isectIS, &points);
686: ISDestroy(&isectIS);
687: }
688: } else {
689: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
690: PetscQuadrature quad = NULL;
691: IS pointIS;
693: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
694: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
695: if (maxDegree <= 1) {
696: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
697: }
698: if (!quad) {
699: if (!h && allPoints) {
700: quad = allPoints;
701: allPoints = NULL;
702: } else {
703: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
704: }
705: }
706: DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);
707: for (p = pStart; p < pEnd; ++p) {
708: PetscArrayzero(values, numValues);
709: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
710: DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, T, TAux, type, funcs, ctxs, fieldActive, values);
711: if (ierr) {
712: PetscErrorCode ierr2;
713: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
714: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
715:
716: }
717: if (transform) {DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
718: DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, values, mode);
719: }
720: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
721: PetscFEGeomDestroy(&fegeom);
722: PetscQuadratureDestroy(&quad);
723: ISDestroy(&pointIS);
724: }
725: ISDestroy(&heightIS);
726: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
727: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
728: }
729: /* Cleanup */
730: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
731: PetscInt effectiveHeight = auxBd ? minHeight : 0;
732: PetscFE fem, subfem;
733: PetscBool isfe;
735: for (f = 0; f < NfIn; ++f) {
736: PetscDSIsFE_Internal(dsIn, f, &isfe);
737: PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);
738: if (!effectiveHeight) {subfem = fem;}
739: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
740: PetscTabulationDestroy(&T[f]);
741: }
742: for (f = 0; f < NfAux; ++f) {
743: PetscDSIsFE_Internal(dsAux, f, &isfe);
744: PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);
745: if (!effectiveHeight || auxBd) {subfem = fem;}
746: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
747: PetscTabulationDestroy(&TAux[f]);
748: }
749: PetscFree2(T, TAux);
750: }
751: PetscQuadratureDestroy(&allPoints);
752: PetscFree2(isFE, sp);
753: if (maxHeight > 0) {PetscFree(cellsp);}
754: DMDestroy(&plex);
755: DMDestroy(&plexIn);
756: if (dmAux) {DMDestroy(&plexAux);}
757: return(0);
758: }
760: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
761: {
765: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
766: return(0);
767: }
769: 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)
770: {
774: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
775: return(0);
776: }
778: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
779: void (**funcs)(PetscInt, PetscInt, PetscInt,
780: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
781: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
782: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
783: InsertMode mode, Vec localX)
784: {
788: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
789: return(0);
790: }
792: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
793: void (**funcs)(PetscInt, PetscInt, PetscInt,
794: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
795: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
796: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
797: InsertMode mode, Vec localX)
798: {
802: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
803: return(0);
804: }