Actual source code: plexproject.c
petsc-3.10.5 2019-03-28
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;
14: DMGetCoordinateDim(dm,&coordDim);
15: PetscDSGetNumFields(prob, &Nf);
16: PetscDSGetComponents(prob, &Nc);
17: PetscDSGetTotalDimension(prob, &totDim);
18: /* Get values for closure */
19: isAffine = fegeom->isAffine;
20: for (f = 0, v = 0, tp = 0; f < Nf; ++f) {
21: void * const ctx = ctxs ? ctxs[f] : NULL;
23: if (!sp[f]) continue;
24: PetscDualSpaceGetDimension(sp[f], &spDim);
25: if (funcs[f]) {
26: if (isFE[f]) {
27: PetscQuadrature allPoints;
28: PetscInt q, dim, numPoints;
29: const PetscReal *points;
30: PetscScalar *pointEval;
31: PetscReal *x;
32: DM dm;
34: PetscDualSpaceGetDM(sp[f],&dm);
35: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
36: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
37: PetscDualSpaceGetDM(sp[f],&dm);
38: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
39: DMGetWorkArray(dm,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: (*funcs[f])(coordDim,time,v0, Nc[f], &pointEval[Nc[f]*q], ctx);
50: }
51: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
52: DMRestoreWorkArray(dm,coordDim,MPIU_REAL,&x);
53: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
54: v += spDim;
55: } else {
56: for (d = 0; d < spDim; ++d, ++v) {
57: PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);
58: }
59: }
60: } else {
61: for (d = 0; d < spDim; d++, v++) {values[v] = 0.0;}
62: }
63: }
64: return(0);
65: }
67: static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS prob, DM dmAux, PetscDS probAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fegeom, PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
68: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
69: void (**funcs)(PetscInt, PetscInt, PetscInt,
70: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
71: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
72: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs,
73: PetscScalar values[])
74: {
75: PetscSection section, sectionAux = NULL;
76: PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL, *bc;
77: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
78: PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL;
79: const PetscScalar *constants;
80: PetscReal *x;
81: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
82: const PetscInt dE = fegeom->dimEmbed;
83: PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0;
84: PetscBool isAffine;
85: PetscErrorCode ierr;
88: PetscDSGetNumFields(prob, &Nf);
89: PetscDSGetDimensions(prob, &Nb);
90: PetscDSGetComponents(prob, &Nc);
91: PetscDSGetComponentOffsets(prob, &uOff);
92: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
93: PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);
94: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
95: PetscDSGetConstants(prob, &numConstants, &constants);
96: DMGetSection(dm, §ion);
97: DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);
98: if (dmAux) {
99: DMLabel spmap;
100: PetscInt subp = p;
102: /* If dm is a submesh, do not get subpoint */
103: DMPlexGetSubpointMap(dm, &spmap);
104: if (!spmap) {DMPlexGetSubpoint(dmAux, p, &subp);}
105: PetscDSGetSpatialDimension(probAux, &dimAux);
106: PetscDSGetNumFields(probAux, &NfAux);
107: PetscDSGetDimensions(probAux, &NbAux);
108: PetscDSGetComponents(probAux, &NcAux);
109: DMGetSection(dmAux, §ionAux);
110: PetscDSGetComponentOffsets(probAux, &aOff);
111: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
112: PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);
113: PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
114: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
115: }
116: /* Get values for closure */
117: isAffine = fegeom->isAffine;
118: for (f = 0, v = 0; f < Nf; ++f) {
119: PetscQuadrature allPoints;
120: PetscInt q, dim, numPoints;
121: const PetscReal *points;
122: PetscScalar *pointEval;
123: DM dm;
125: if (!sp[f]) continue;
126: PetscDualSpaceGetDimension(sp[f], &spDim);
127: if (!funcs[f]) {
128: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
129: continue;
130: }
131: PetscDualSpaceGetDM(sp[f],&dm);
132: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
133: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
134: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
135: for (q = 0; q < numPoints; ++q, ++tp) {
136: const PetscReal *v0;
137: const PetscReal *invJ;
139: if (isAffine) {
140: CoordinatesRefToReal(dE, dim, fegeom->xi, fegeom->v, fegeom->J, &points[q*dim], x);
141: v0 = x;
142: invJ = fegeom->invJ;
143: } else {
144: v0 = &fegeom->v[tp*dE];
145: invJ = &fegeom->invJ[tp*dE*dE];
146: }
147: EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, invJ, coefficients, coefficients_t, u, u_x, u_t);
148: if (probAux) {EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);}
149: (*funcs[f])(dE, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, v0, numConstants, constants, &pointEval[Nc[f]*q]);
150: }
151: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
152: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
153: v += spDim;
154: }
155: DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);
156: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
157: return(0);
158: }
160: 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[],
161: PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
162: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
163: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
164: {
165: PetscFVCellGeom fvgeom;
166: PetscInt dim, dimEmbed;
167: PetscErrorCode ierr;
170: DMGetDimension(dm, &dim);
171: DMGetCoordinateDim(dm, &dimEmbed);
172: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
173: switch (type) {
174: case DM_BC_ESSENTIAL:
175: case DM_BC_NATURAL:
176: DMProjectPoint_Func_Private(dm, prob, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
177: case DM_BC_ESSENTIAL_FIELD:
178: case DM_BC_NATURAL_FIELD:
179: DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps,
180: basisTab, basisDerTab, basisTabAux, basisDerTabAux,
181: (void (**)(PetscInt, PetscInt, PetscInt,
182: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
183: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
184: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
185: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
186: }
187: return(0);
188: }
190: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
191: {
192: PetscReal *points;
193: PetscInt f, numPoints;
197: numPoints = 0;
198: for (f = 0; f < Nf; ++f) {
199: if (funcs[f]) {
200: PetscQuadrature fAllPoints;
201: PetscInt fNumPoints;
203: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
204: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
205: numPoints += fNumPoints;
206: }
207: }
208: PetscMalloc1(dim*numPoints,&points);
209: numPoints = 0;
210: for (f = 0; f < Nf; ++f) {
211: PetscInt spDim;
213: PetscDualSpaceGetDimension(sp[f], &spDim);
214: if (funcs[f]) {
215: PetscQuadrature fAllPoints;
216: PetscInt qdim, fNumPoints, q;
217: const PetscReal *fPoints;
219: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
220: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
221: 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);
222: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
223: numPoints += fNumPoints;
224: }
225: }
226: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
227: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
228: return(0);
229: }
231: /*
232: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
234: There are several different scenarios:
236: 1) Volumetric mesh with volumetric auxiliary data
238: Here minHeight=0 since we loop over cells.
240: 2) Boundary mesh with boundary auxiliary data
242: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
244: 3) Volumetric mesh with boundary auxiliary data
246: 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.
248: The maxHeight is used to support enforcement of constraints in DMForest.
250: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
252: 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.
253: - 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
254: 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
255: dual spaces for the boundary from our input spaces.
256: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
258: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
260: 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.
261: */
262: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
263: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
264: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
265: InsertMode mode, Vec localX)
266: {
267: DM dmAux = NULL;
268: PetscDS prob, probAux = NULL;
269: Vec localA = NULL;
270: PetscSection section;
271: PetscDualSpace *sp, *cellsp;
272: PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL;
273: PetscInt *Nc;
274: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfAux = 0, f;
275: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE;
276: DMField coordField;
277: DMLabel depthLabel;
278: PetscQuadrature allPoints = NULL;
279: PetscErrorCode ierr;
282: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
283: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
284: DMGetDimension(dm, &dim);
285: DMPlexGetVTKCellHeight(dm, &minHeight);
286: /* Auxiliary information can only be used with interpolation of field functions */
287: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
288: if (!minHeight && dmAux) {
289: DMLabel spmap;
291: /* If dmAux is a surface, then force the projection to take place over a surface */
292: DMPlexGetSubpointMap(dmAux, &spmap);
293: if (spmap) {minHeight = 1; auxBd = PETSC_TRUE;}
294: }
295: }
296: DMPlexGetMaxProjectionHeight(dm, &maxHeight);
297: maxHeight = PetscMax(maxHeight, minHeight);
298: if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);}
299: DMGetDS(dm, &prob);
300: DMGetCoordinateDim(dm, &dimEmbed);
301: DMGetSection(dm, §ion);
302: PetscSectionGetNumFields(section, &Nf);
303: if (dmAux) {
304: DMGetDS(dmAux, &probAux);
305: PetscDSGetNumFields(probAux, &NfAux);
306: }
307: PetscDSGetComponents(prob, &Nc);
308: PetscMalloc2(Nf, &isFE, Nf, &sp);
309: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
310: else {cellsp = sp;}
311: if (localU && localU != localX) {DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
312: /* Get cell dual spaces */
313: for (f = 0; f < Nf; ++f) {
314: PetscObject obj;
315: PetscClassId id;
317: DMGetField(dm, f, &obj);
318: PetscObjectGetClassId(obj, &id);
319: if (id == PETSCFE_CLASSID) {
320: PetscFE fe = (PetscFE) obj;
322: hasFE = PETSC_TRUE;
323: isFE[f] = PETSC_TRUE;
324: PetscFEGetDualSpace(fe, &cellsp[f]);
325: } else if (id == PETSCFV_CLASSID) {
326: PetscFV fv = (PetscFV) obj;
328: hasFV = PETSC_TRUE;
329: isFE[f] = PETSC_FALSE;
330: PetscFVGetDualSpace(fv, &cellsp[f]);
331: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
332: }
333: DMGetCoordinateField(dm,&coordField);
334: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
335: PetscInt effectiveHeight = auxBd ? minHeight : 0;
336: PetscFE fem, subfem;
337: const PetscReal *points;
338: PetscInt numPoints;
340: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
341: for (f = 0; f < Nf; ++f) {
342: if (!effectiveHeight) {sp[f] = cellsp[f];}
343: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
344: }
345: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
346: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
347: PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);
348: for (f = 0; f < Nf; ++f) {
349: if (!isFE[f]) continue;
350: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
351: if (!effectiveHeight) {subfem = fem;}
352: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
353: PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);
354: }
355: for (f = 0; f < NfAux; ++f) {
356: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
357: if (!effectiveHeight || auxBd) {subfem = fem;}
358: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
359: PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);
360: }
361: }
362: DMPlexGetDepth(dm,&depth);
363: DMPlexGetDepthLabel(dm,&depthLabel);
364: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
365: for (h = minHeight; h <= maxHeight; h++) {
366: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
367: PetscDS probEff = prob;
368: PetscScalar *values;
369: PetscBool *fieldActive;
370: PetscInt maxDegree;
371: PetscInt pStart, pEnd, p, spDim, totDim, numValues;
372: IS heightIS;
374: if (effectiveHeight) {PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);}
375: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
376: DMLabelGetStratumIS(depthLabel,depth - h,&heightIS);
377: if (!h) {
378: PetscInt cEndInterior;
380: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
381: pEnd = cEndInterior < 0 ? pEnd : cEndInterior;
382: }
383: if (pEnd <= pStart) {
384: ISDestroy(&heightIS);
385: continue;
386: }
387: /* Compute totDim, the number of dofs in the closure of a point at this height */
388: totDim = 0;
389: for (f = 0; f < Nf; ++f) {
390: if (!effectiveHeight) {
391: sp[f] = cellsp[f];
392: } else {
393: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
394: if (!sp[f]) continue;
395: }
396: PetscDualSpaceGetDimension(sp[f], &spDim);
397: totDim += spDim;
398: }
399: DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
400: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim);
401: if (!totDim) {
402: ISDestroy(&heightIS);
403: continue;
404: }
405: /* Loop over points at this height */
406: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
407: DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);
408: for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
409: if (label) {
410: PetscInt i;
412: for (i = 0; i < numIds; ++i) {
413: IS pointIS, isectIS;
414: const PetscInt *points;
415: PetscInt n;
416: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
417: PetscQuadrature quad = NULL;
419: DMLabelGetStratumIS(label, ids[i], &pointIS);
420: if (!pointIS) continue; /* No points with that id on this process */
421: ISIntersect(pointIS,heightIS,&isectIS);
422: ISDestroy(&pointIS);
423: if (!isectIS) continue;
424: ISGetLocalSize(isectIS, &n);
425: ISGetIndices(isectIS, &points);
426: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
427: if (maxDegree <= 1) {
428: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
429: }
430: if (!quad) {
431: if (!h && allPoints) {
432: quad = allPoints;
433: allPoints = NULL;
434: } else {
435: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
436: }
437: }
438: DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);
439: for (p = 0; p < n; ++p) {
440: const PetscInt point = points[p];
442: PetscMemzero(values, numValues * sizeof(PetscScalar));
443: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
444: 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);
445: if (ierr) {
446: PetscErrorCode ierr2;
447: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
448: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
449:
450: }
451: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);
452: }
453: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
454: PetscFEGeomDestroy(&fegeom);
455: PetscQuadratureDestroy(&quad);
456: ISRestoreIndices(isectIS, &points);
457: ISDestroy(&isectIS);
458: }
459: } else {
460: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
461: PetscQuadrature quad = NULL;
462: IS pointIS;
464: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
465: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
466: if (maxDegree <= 1) {
467: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
468: }
469: if (!quad) {
470: if (!h && allPoints) {
471: quad = allPoints;
472: allPoints = NULL;
473: } else {
474: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
475: }
476: }
477: DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);
478: for (p = pStart; p < pEnd; ++p) {
479: PetscMemzero(values, numValues * sizeof(PetscScalar));
480: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
481: 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);
482: if (ierr) {
483: PetscErrorCode ierr2;
484: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
485: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
486:
487: }
488: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);
489: }
490: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
491: PetscFEGeomDestroy(&fegeom);
492: PetscQuadratureDestroy(&quad);
493: ISDestroy(&pointIS);
494: }
495: ISDestroy(&heightIS);
496: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
497: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
498: }
499: /* Cleanup */
500: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
501: PetscInt effectiveHeight = auxBd ? minHeight : 0;
502: PetscFE fem, subfem;
504: for (f = 0; f < Nf; ++f) {
505: if (!isFE[f]) continue;
506: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
507: if (!effectiveHeight) {subfem = fem;}
508: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
509: PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);
510: }
511: for (f = 0; f < NfAux; ++f) {
512: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
513: if (!effectiveHeight || auxBd) {subfem = fem;}
514: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
515: PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);
516: }
517: PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);
518: }
519: PetscQuadratureDestroy(&allPoints);
520: PetscFree2(isFE, sp);
521: if (maxHeight > 0) {PetscFree(cellsp);}
522: return(0);
523: }
525: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
526: {
530: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
531: return(0);
532: }
534: 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)
535: {
539: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
540: return(0);
541: }
543: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
544: void (**funcs)(PetscInt, PetscInt, PetscInt,
545: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
546: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
547: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
548: InsertMode mode, Vec localX)
549: {
553: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
554: return(0);
555: }
557: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
558: void (**funcs)(PetscInt, PetscInt, PetscInt,
559: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
560: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
561: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
562: InsertMode mode, Vec localX)
563: {
567: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
568: return(0);
569: }