Actual source code: plexproject.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS prob, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[],
6: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs,
7: PetscScalar values[])
8: {
9: PetscInt coordDim, Nf, *Nc, f, totDim, spDim, d, v, tp;
10: PetscBool isAffine, transform;
14: DMGetCoordinateDim(dm,&coordDim);
15: DMHasBasisTransform(dm, &transform);
16: PetscDSGetNumFields(prob, &Nf);
17: PetscDSGetComponents(prob, &Nc);
18: PetscDSGetTotalDimension(prob, &totDim);
19: /* Get values for closure */
20: isAffine = fegeom->isAffine;
21: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
22: void * const ctx = ctxs ? ctxs[f] : NULL;
24: if (!sp[f]) continue;
25: PetscDualSpaceGetDimension(sp[f], &spDim);
26: if (funcs[f]) {
27: if (isFE[f]) {
28: PetscQuadrature allPoints;
29: PetscInt q, dim, numPoints;
30: const PetscReal *points;
31: PetscScalar *pointEval;
32: PetscReal *x;
33: DM rdm;
35: PetscDualSpaceGetDM(sp[f],&rdm);
36: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
37: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
38: DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
39: DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);
40: for (q = 0; q < numPoints; q++, tp++) {
41: const PetscReal *v0;
43: if (isAffine) {
44: CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, &points[q*dim], x);
45: v0 = x;
46: } else {
47: v0 = &fegeom->v[tp*coordDim];
48: }
49: if (transform) {DMPlexBasisTransformApplyReal_Internal(dm, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx); v0 = x;}
50: (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);
51: }
52: /* Transform point evaluations pointEval[q,c] */
53: PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);
54: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
55: DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);
56: DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
57: v += spDim;
58: } else {
59: for (d = 0; d < spDim; ++d, ++v) {
60: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
61: }
62: }
63: } else {
64: for (d = 0; d < spDim; d++, v++) {values[v] = 0.0;}
65: }
66: }
67: return(0);
68: }
70: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS prob, DM dmAux, PetscDS probAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
71: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
72: void (**funcs)(PetscInt, PetscInt, PetscInt,
73: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
74: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
75: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
76: PetscScalar values[])
77: {
78: PetscSection section, sectionAux = NULL;
79: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
80: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
81: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
82: const PetscScalar *constants;
83: PetscReal *x;
84: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
85: PetscFEGeom fegeom;
86: const PetscInt dE = cgeom->dimEmbed;
87: PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
88: PetscBool isAffine, transform;
89: PetscErrorCode ierr;
92: PetscDSGetNumFields(prob, &Nf);
93: PetscDSGetDimensions(prob, &Nb);
94: PetscDSGetComponents(prob, &Nc);
95: PetscDSGetComponentOffsets(prob, &uOff);
96: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
97: PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);
98: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
99: PetscDSGetConstants(prob, &numConstants, &constants);
100: DMHasBasisTransform(dm, &transform);
101: DMGetLocalSection(dm, §ion);
102: DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);
103: if (dmAux) {
104: PetscInt subp;
106: DMPlexGetAuxiliaryPoint(dm, dmAux, p, &subp);
107: PetscDSGetSpatialDimension(probAux, &dimAux);
108: PetscDSGetNumFields(probAux, &NfAux);
109: PetscDSGetDimensions(probAux, &NbAux);
110: PetscDSGetComponents(probAux, &NcAux);
111: DMGetLocalSection(dmAux, §ionAux);
112: PetscDSGetComponentOffsets(probAux, &aOff);
113: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
114: PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);
115: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
116: }
117: /* Get values for closure */
118: isAffine = cgeom->isAffine;
119: if (isAffine) {
120: fegeom.v = x;
121: fegeom.xi = cgeom->xi;
122: fegeom.J = cgeom->J;
123: fegeom.invJ = cgeom->invJ;
124: fegeom.detJ = cgeom->detJ;
125: }
126: for (f = 0, v = 0; f < Nf; ++f) {
127: PetscQuadrature allPoints;
128: PetscInt q, dim, numPoints;
129: const PetscReal *points;
130: PetscScalar *pointEval;
131: DM dm;
133: if (!sp[f]) continue;
134: PetscDualSpaceGetDimension(sp[f], &spDim);
135: if (!funcs[f]) {
136: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
137: continue;
138: }
139: PetscDualSpaceGetDM(sp[f],&dm);
140: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
141: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
142: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
143: for (q = 0; q < numPoints; ++q, ++tp) {
144: if (isAffine) {
145: CoordinatesRefToReal(dE, dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x);
146: } else {
147: fegeom.v = &cgeom->v[tp*dE];
148: fegeom.J = &cgeom->J[tp*dE*dE];
149: fegeom.invJ = &cgeom->invJ[tp*dE*dE];
150: fegeom.detJ = &cgeom->detJ[tp];
151: }
152: PetscFEEvaluateFieldJets_Internal(prob, dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, &fegeom, coefficients, coefficients_t, u, u_x, u_t);
153: if (probAux) {PetscFEEvaluateFieldJets_Internal(probAux, dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
154: if (transform) {DMPlexBasisTransformApplyReal_Internal(dm, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);}
155: (*funcs[f])(dE, Nf, 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]);
156: }
157: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
158: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
159: v += spDim;
160: }
161: DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);
162: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
163: return(0);
164: }
166: static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS prob, DM dmAux, PetscDS probAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
167: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
168: void (**funcs)(PetscInt, PetscInt, PetscInt,
169: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
170: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
171: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
172: PetscScalar values[])
173: {
174: PetscSection section, sectionAux = NULL;
175: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc;
176: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
177: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
178: const PetscScalar *constants;
179: PetscReal *x;
180: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
181: PetscFEGeom fegeom, cgeom;
182: const PetscInt dE = fgeom->dimEmbed;
183: PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
184: PetscBool isAffine;
185: PetscErrorCode ierr;
188: PetscDSGetNumFields(prob, &Nf);
189: PetscDSGetDimensions(prob, &Nb);
190: PetscDSGetComponents(prob, &Nc);
191: PetscDSGetComponentOffsets(prob, &uOff);
192: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
193: PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);
194: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
195: PetscDSGetConstants(prob, &numConstants, &constants);
196: DMGetLocalSection(dm, §ion);
197: DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);
198: if (dmAux) {
199: DMLabel spmap;
200: PetscInt subp = p;
202: /* If dm is a submesh, do not get subpoint */
203: DMPlexGetSubpointMap(dm, &spmap);
204: if (!spmap) {DMPlexGetSubpoint(dmAux, p, &subp);}
205: PetscDSGetSpatialDimension(probAux, &dimAux);
206: PetscDSGetNumFields(probAux, &NfAux);
207: PetscDSGetDimensions(probAux, &NbAux);
208: PetscDSGetComponents(probAux, &NcAux);
209: DMGetLocalSection(dmAux, §ionAux);
210: PetscDSGetComponentOffsets(probAux, &aOff);
211: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
212: PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);
213: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
214: }
215: /* Get values for closure */
216: isAffine = fgeom->isAffine;
217: fegeom.n = 0;
218: fegeom.J = 0;
219: fegeom.v = 0;
220: fegeom.xi = 0;
221: if (isAffine) {
222: fegeom.v = x;
223: fegeom.xi = fgeom->xi;
224: fegeom.J = fgeom->J;
225: fegeom.invJ = fgeom->invJ;
226: fegeom.detJ = fgeom->detJ;
227: fegeom.n = fgeom->n;
229: cgeom.J = fgeom->suppJ[0];
230: cgeom.invJ = fgeom->suppInvJ[0];
231: cgeom.detJ = fgeom->suppDetJ[0];
232: }
233: for (f = 0, v = 0; f < Nf; ++f) {
234: PetscQuadrature allPoints;
235: PetscInt q, dim, numPoints;
236: const PetscReal *points;
237: PetscScalar *pointEval;
238: DM dm;
240: if (!sp[f]) continue;
241: PetscDualSpaceGetDimension(sp[f], &spDim);
242: if (!funcs[f]) {
243: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
244: continue;
245: }
246: PetscDualSpaceGetDM(sp[f],&dm);
247: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
248: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
249: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
250: for (q = 0; q < numPoints; ++q, ++tp) {
251: if (isAffine) {
252: CoordinatesRefToReal(dE, dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x);
253: } else {
254: fegeom.v = &fgeom->v[tp*dE];
255: fegeom.J = &fgeom->J[tp*dE*dE];
256: fegeom.invJ = &fgeom->invJ[tp*dE*dE];
257: fegeom.detJ = &fgeom->detJ[tp];
258: fegeom.n = &fgeom->n[tp*dE];
260: cgeom.J = &fgeom->suppJ[0][tp*dE*dE];
261: cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE];
262: cgeom.detJ = &fgeom->suppDetJ[0][tp];
263: }
264: PetscFEEvaluateFieldJets_Internal(prob, dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, &cgeom, coefficients, coefficients_t, u, u_x, u_t);
265: if (probAux) {PetscFEEvaluateFieldJets_Internal(probAux, dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
266: (*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]);
267: }
268: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
269: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
270: v += spDim;
271: }
272: DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);
273: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
274: return(0);
275: }
277: static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS prob, PetscFEGeom *fegeom, DM dmAux, PetscDS probAux, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[],
278: PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
279: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
280: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
281: {
282: PetscFVCellGeom fvgeom;
283: PetscInt dim, dimEmbed;
284: PetscErrorCode ierr;
287: DMGetDimension(dm, &dim);
288: DMGetCoordinateDim(dm, &dimEmbed);
289: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
290: switch (type) {
291: case DM_BC_ESSENTIAL:
292: case DM_BC_NATURAL:
293: DMProjectPoint_Func_Private(dm, prob, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
294: case DM_BC_ESSENTIAL_FIELD:
295: case DM_BC_NATURAL_FIELD:
296: DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps,
297: basisTab, basisDerTab, basisTabAux, basisDerTabAux,
298: (void (**)(PetscInt, PetscInt, PetscInt,
299: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
300: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
301: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
302: case DM_BC_ESSENTIAL_BD_FIELD:
303: DMProjectPoint_BdField_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps,
304: basisTab, basisDerTab, basisTabAux, basisDerTabAux,
305: (void (**)(PetscInt, PetscInt, PetscInt,
306: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
307: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
308: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
309: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
310: }
311: return(0);
312: }
314: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
315: {
316: PetscReal *points;
317: PetscInt f, numPoints;
321: numPoints = 0;
322: for (f = 0; f < Nf; ++f) {
323: if (funcs[f]) {
324: PetscQuadrature fAllPoints;
325: PetscInt fNumPoints;
327: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
328: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
329: numPoints += fNumPoints;
330: }
331: }
332: PetscMalloc1(dim*numPoints,&points);
333: numPoints = 0;
334: for (f = 0; f < Nf; ++f) {
335: PetscInt spDim;
337: PetscDualSpaceGetDimension(sp[f], &spDim);
338: if (funcs[f]) {
339: PetscQuadrature fAllPoints;
340: PetscInt qdim, fNumPoints, q;
341: const PetscReal *fPoints;
343: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
344: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
345: 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);
346: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
347: numPoints += fNumPoints;
348: }
349: }
350: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
351: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
352: return(0);
353: }
355: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *prob)
356: {
357: DMLabel depthLabel;
358: PetscInt dim, cdepth, ls = -1, i;
362: if (lStart) *lStart = -1;
363: if (!label) return(0);
364: DMGetDimension(dm, &dim);
365: DMPlexGetDepthLabel(dm, &depthLabel);
366: cdepth = dim - height;
367: for (i = 0; i < numIds; ++i) {
368: IS pointIS;
369: const PetscInt *points;
370: PetscInt pdepth;
372: DMLabelGetStratumIS(label, ids[i], &pointIS);
373: if (!pointIS) continue; /* No points with that id on this process */
374: ISGetIndices(pointIS, &points);
375: DMLabelGetValue(depthLabel, points[0], &pdepth);
376: if (pdepth == cdepth) {
377: ls = points[0];
378: if (prob) {DMGetCellDS(dm, ls, prob);}
379: }
380: ISRestoreIndices(pointIS, &points);
381: ISDestroy(&pointIS);
382: if (ls >= 0) break;
383: }
384: if (lStart) *lStart = ls;
385: return(0);
386: }
388: /*
389: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
391: There are several different scenarios:
393: 1) Volumetric mesh with volumetric auxiliary data
395: Here minHeight=0 since we loop over cells.
397: 2) Boundary mesh with boundary auxiliary data
399: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
401: 3) Volumetric mesh with boundary auxiliary data
403: 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.
405: The maxHeight is used to support enforcement of constraints in DMForest.
407: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
409: 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.
410: - 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
411: 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
412: dual spaces for the boundary from our input spaces.
413: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
415: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
417: 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.
418: */
419: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
420: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
421: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
422: InsertMode mode, Vec localX)
423: {
424: DM dmAux = NULL, tdm;
425: PetscDS prob = NULL, probAux = NULL;
426: Vec localA = NULL, tv;
427: PetscSection section;
428: PetscDualSpace *sp, *cellsp;
429: PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL;
430: PetscInt *Nc;
431: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfAux = 0, f;
432: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, transform;
433: DMField coordField;
434: DMLabel depthLabel;
435: PetscQuadrature allPoints = NULL;
436: PetscErrorCode ierr;
439: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
440: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
441: DMGetDimension(dm, &dim);
442: DMPlexGetVTKCellHeight(dm, &minHeight);
443: DMGetBasisTransformDM_Internal(dm, &tdm);
444: DMGetBasisTransformVec_Internal(dm, &tv);
445: DMHasBasisTransform(dm, &transform);
446: /* Auxiliary information can only be used with interpolation of field functions */
447: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
448: if (dmAux && !localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
449: if (!minHeight && dmAux) {
450: DMLabel spmap;
452: /* If dmAux is a surface, then force the projection to take place over a surface */
453: DMPlexGetSubpointMap(dmAux, &spmap);
454: if (spmap) {
455: DMPlexGetVTKCellHeight(dmAux, &minHeight);
456: auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
457: }
458: }
459: }
460: DMPlexGetDepth(dm,&depth);
461: DMPlexGetDepthLabel(dm,&depthLabel);
462: DMPlexGetMaxProjectionHeight(dm, &maxHeight);
463: maxHeight = PetscMax(maxHeight, minHeight);
464: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
465: DMGetFirstLabelEntry_Private(dm, label, numIds, ids, 0, NULL, &prob);
466: if (!prob) {DMGetDS(dm, &prob);}
467: PetscDSGetNumFields(prob, &Nf);
468: DMGetCoordinateDim(dm, &dimEmbed);
469: DMGetLocalSection(dm, §ion);
470: if (dmAux) {
471: DMGetDS(dmAux, &probAux);
472: PetscDSGetNumFields(probAux, &NfAux);
473: }
474: PetscDSGetComponents(prob, &Nc);
475: PetscMalloc2(Nf, &isFE, Nf, &sp);
476: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
477: else {cellsp = sp;}
478: if (localU && localU != localX) {DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
479: /* Get cell dual spaces */
480: for (f = 0; f < Nf; ++f) {
481: PetscObject obj;
482: PetscClassId id;
484: PetscDSGetDiscretization(prob, f, &obj);
485: PetscObjectGetClassId(obj, &id);
486: if (id == PETSCFE_CLASSID) {
487: PetscFE fe = (PetscFE) obj;
489: hasFE = PETSC_TRUE;
490: isFE[f] = PETSC_TRUE;
491: PetscFEGetDualSpace(fe, &cellsp[f]);
492: } else if (id == PETSCFV_CLASSID) {
493: PetscFV fv = (PetscFV) obj;
495: hasFV = PETSC_TRUE;
496: isFE[f] = PETSC_FALSE;
497: PetscFVGetDualSpace(fv, &cellsp[f]);
498: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
499: }
500: DMGetCoordinateField(dm,&coordField);
501: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
502: PetscInt effectiveHeight = auxBd ? minHeight : 0;
503: PetscFE fem, subfem;
504: const PetscReal *points;
505: PetscInt numPoints;
507: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
508: for (f = 0; f < Nf; ++f) {
509: if (!effectiveHeight) {sp[f] = cellsp[f];}
510: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
511: }
512: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
513: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
514: PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);
515: for (f = 0; f < Nf; ++f) {
516: if (!isFE[f]) continue;
517: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
518: if (!effectiveHeight) {subfem = fem;}
519: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
520: PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);
521: }
522: for (f = 0; f < NfAux; ++f) {
523: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
524: if (!effectiveHeight || auxBd) {subfem = fem;}
525: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
526: PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);
527: }
528: }
529: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
530: for (h = minHeight; h <= maxHeight; h++) {
531: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
532: PetscDS probEff = prob;
533: PetscScalar *values;
534: PetscBool *fieldActive;
535: PetscInt maxDegree;
536: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
537: IS heightIS;
539: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
540: DMGetFirstLabelEntry_Private(dm, label, numIds, ids, h, &lStart, NULL);
541: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
542: if (!h) {
543: PetscInt cEndInterior;
545: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
546: pEnd = cEndInterior < 0 ? pEnd : cEndInterior;
547: }
548: if (pEnd <= pStart) {
549: ISDestroy(&heightIS);
550: continue;
551: }
552: /* Compute totDim, the number of dofs in the closure of a point at this height */
553: totDim = 0;
554: for (f = 0; f < Nf; ++f) {
555: if (!effectiveHeight) {
556: sp[f] = cellsp[f];
557: } else {
558: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
559: if (!sp[f]) continue;
560: }
561: PetscDualSpaceGetDimension(sp[f], &spDim);
562: totDim += spDim;
563: }
564: DMPlexVecGetClosure(dm, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
565: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim);
566: if (!totDim) {
567: ISDestroy(&heightIS);
568: continue;
569: }
570: if (effectiveHeight) {PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);}
571: /* Loop over points at this height */
572: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
573: DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);
574: for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
575: if (label) {
576: PetscInt i;
578: for (i = 0; i < numIds; ++i) {
579: IS pointIS, isectIS;
580: const PetscInt *points;
581: PetscInt n;
582: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
583: PetscQuadrature quad = NULL;
585: DMLabelGetStratumIS(label, ids[i], &pointIS);
586: if (!pointIS) continue; /* No points with that id on this process */
587: ISIntersect(pointIS,heightIS,&isectIS);
588: ISDestroy(&pointIS);
589: if (!isectIS) continue;
590: ISGetLocalSize(isectIS, &n);
591: ISGetIndices(isectIS, &points);
592: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
593: if (maxDegree <= 1) {
594: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
595: }
596: if (!quad) {
597: if (!h && allPoints) {
598: quad = allPoints;
599: allPoints = NULL;
600: } else {
601: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
602: }
603: }
604: DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);
605: for (p = 0; p < n; ++p) {
606: const PetscInt point = points[p];
608: PetscArrayzero(values, numValues);
609: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
610: DMProjectPoint_Private(dm, probEff, chunkgeom, dmAux, probAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values);
611: if (ierr) {
612: PetscErrorCode ierr2;
613: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
614: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
615:
616: }
617: if (transform) {DMPlexBasisTransformPoint_Internal(dm, tdm, tv, point, fieldActive, PETSC_FALSE, values);}
618: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);
619: }
620: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
621: PetscFEGeomDestroy(&fegeom);
622: PetscQuadratureDestroy(&quad);
623: ISRestoreIndices(isectIS, &points);
624: ISDestroy(&isectIS);
625: }
626: } else {
627: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
628: PetscQuadrature quad = NULL;
629: IS pointIS;
631: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
632: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
633: if (maxDegree <= 1) {
634: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
635: }
636: if (!quad) {
637: if (!h && allPoints) {
638: quad = allPoints;
639: allPoints = NULL;
640: } else {
641: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
642: }
643: }
644: DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);
645: for (p = pStart; p < pEnd; ++p) {
646: PetscArrayzero(values, numValues);
647: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
648: DMProjectPoint_Private(dm, probEff, chunkgeom, dmAux, probAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values);
649: if (ierr) {
650: PetscErrorCode ierr2;
651: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
652: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
653:
654: }
655: if (transform) {DMPlexBasisTransformPoint_Internal(dm, tdm, tv, p, fieldActive, PETSC_FALSE, values);}
656: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);
657: }
658: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
659: PetscFEGeomDestroy(&fegeom);
660: PetscQuadratureDestroy(&quad);
661: ISDestroy(&pointIS);
662: }
663: ISDestroy(&heightIS);
664: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
665: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
666: }
667: /* Cleanup */
668: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
669: PetscInt effectiveHeight = auxBd ? minHeight : 0;
670: PetscFE fem, subfem;
672: for (f = 0; f < Nf; ++f) {
673: if (!isFE[f]) continue;
674: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
675: if (!effectiveHeight) {subfem = fem;}
676: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
677: PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);
678: }
679: for (f = 0; f < NfAux; ++f) {
680: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
681: if (!effectiveHeight || auxBd) {subfem = fem;}
682: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
683: PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);
684: }
685: PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);
686: }
687: PetscQuadratureDestroy(&allPoints);
688: PetscFree2(isFE, sp);
689: if (maxHeight > 0) {PetscFree(cellsp);}
690: return(0);
691: }
693: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
694: {
698: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
699: return(0);
700: }
702: 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)
703: {
707: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
708: return(0);
709: }
711: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
712: void (**funcs)(PetscInt, PetscInt, PetscInt,
713: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
714: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
715: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
716: InsertMode mode, Vec localX)
717: {
721: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
722: return(0);
723: }
725: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
726: void (**funcs)(PetscInt, PetscInt, PetscInt,
727: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
728: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
729: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
730: InsertMode mode, Vec localX)
731: {
735: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
736: return(0);
737: }