Actual source code: plexproject.c
petsc-3.11.4 2019-09-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: PetscInt subp;
101: DMPlexGetAuxiliaryPoint(dm, dmAux, p, &subp);
102: PetscDSGetSpatialDimension(probAux, &dimAux);
103: PetscDSGetNumFields(probAux, &NfAux);
104: PetscDSGetDimensions(probAux, &NbAux);
105: PetscDSGetComponents(probAux, &NcAux);
106: DMGetSection(dmAux, §ionAux);
107: PetscDSGetComponentOffsets(probAux, &aOff);
108: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
109: PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);
110: PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
111: DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);
112: }
113: /* Get values for closure */
114: isAffine = fegeom->isAffine;
115: for (f = 0, v = 0; f < Nf; ++f) {
116: PetscQuadrature allPoints;
117: PetscInt q, dim, numPoints;
118: const PetscReal *points;
119: PetscScalar *pointEval;
120: DM dm;
122: if (!sp[f]) continue;
123: PetscDualSpaceGetDimension(sp[f], &spDim);
124: if (!funcs[f]) {
125: for (d = 0; d < spDim; d++, v++) values[v] = 0.;
126: continue;
127: }
128: PetscDualSpaceGetDM(sp[f],&dm);
129: PetscDualSpaceGetAllPoints(sp[f], &allPoints);
130: PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);
131: DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
132: for (q = 0; q < numPoints; ++q, ++tp) {
133: const PetscReal *v0;
134: const PetscReal *invJ;
136: if (isAffine) {
137: CoordinatesRefToReal(dE, dim, fegeom->xi, fegeom->v, fegeom->J, &points[q*dim], x);
138: v0 = x;
139: invJ = fegeom->invJ;
140: } else {
141: v0 = &fegeom->v[tp*dE];
142: invJ = &fegeom->invJ[tp*dE*dE];
143: }
144: EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, invJ, coefficients, coefficients_t, u, u_x, u_t);
145: if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);
146: (*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]);
147: }
148: PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);
149: DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);
150: v += spDim;
151: }
152: DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);
153: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);}
154: return(0);
155: }
157: 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[],
158: PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[],
159: PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux,
160: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[])
161: {
162: PetscFVCellGeom fvgeom;
163: PetscInt dim, dimEmbed;
164: PetscErrorCode ierr;
167: DMGetDimension(dm, &dim);
168: DMGetCoordinateDim(dm, &dimEmbed);
169: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
170: switch (type) {
171: case DM_BC_ESSENTIAL:
172: case DM_BC_NATURAL:
173: DMProjectPoint_Func_Private(dm, prob, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);break;
174: case DM_BC_ESSENTIAL_FIELD:
175: case DM_BC_NATURAL_FIELD:
176: DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps,
177: basisTab, basisDerTab, basisTabAux, basisDerTabAux,
178: (void (**)(PetscInt, PetscInt, PetscInt,
179: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
180: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
181: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);break;
182: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type);
183: }
184: return(0);
185: }
187: static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints)
188: {
189: PetscReal *points;
190: PetscInt f, numPoints;
194: numPoints = 0;
195: for (f = 0; f < Nf; ++f) {
196: if (funcs[f]) {
197: PetscQuadrature fAllPoints;
198: PetscInt fNumPoints;
200: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
201: PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);
202: numPoints += fNumPoints;
203: }
204: }
205: PetscMalloc1(dim*numPoints,&points);
206: numPoints = 0;
207: for (f = 0; f < Nf; ++f) {
208: PetscInt spDim;
210: PetscDualSpaceGetDimension(sp[f], &spDim);
211: if (funcs[f]) {
212: PetscQuadrature fAllPoints;
213: PetscInt qdim, fNumPoints, q;
214: const PetscReal *fPoints;
216: PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);
217: PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);
218: 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);
219: for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q];
220: numPoints += fNumPoints;
221: }
222: }
223: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
224: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
225: return(0);
226: }
228: static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *prob)
229: {
230: DMLabel depthLabel;
231: PetscInt dim, cdepth, ls = -1, i;
235: if (lStart) *lStart = -1;
236: if (!label) return(0);
237: DMGetDimension(dm, &dim);
238: DMPlexGetDepthLabel(dm, &depthLabel);
239: cdepth = dim - height;
240: for (i = 0; i < numIds; ++i) {
241: IS pointIS;
242: const PetscInt *points;
243: PetscInt pdepth;
245: DMLabelGetStratumIS(label, ids[i], &pointIS);
246: if (!pointIS) continue; /* No points with that id on this process */
247: ISGetIndices(pointIS, &points);
248: DMLabelGetValue(depthLabel, points[0], &pdepth);
249: if (pdepth == cdepth) {
250: ls = points[0];
251: if (prob) {DMGetCellDS(dm, ls, prob);}
252: }
253: ISRestoreIndices(pointIS, &points);
254: ISDestroy(&pointIS);
255: if (ls >= 0) break;
256: }
257: if (lStart) *lStart = ls;
258: return(0);
259: }
261: /*
262: This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM
264: There are several different scenarios:
266: 1) Volumetric mesh with volumetric auxiliary data
268: Here minHeight=0 since we loop over cells.
270: 2) Boundary mesh with boundary auxiliary data
272: Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation.
274: 3) Volumetric mesh with boundary auxiliary data
276: 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.
278: The maxHeight is used to support enforcement of constraints in DMForest.
280: If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it.
282: 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.
283: - 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
284: 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
285: dual spaces for the boundary from our input spaces.
286: - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them.
288: We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration
290: 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.
291: */
292: static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU,
293: PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[],
294: DMBoundaryConditionType type, void (**funcs)(void), void **ctxs,
295: InsertMode mode, Vec localX)
296: {
297: DM dmAux = NULL;
298: PetscDS prob = NULL, probAux = NULL;
299: Vec localA = NULL;
300: PetscSection section;
301: PetscDualSpace *sp, *cellsp;
302: PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL;
303: PetscInt *Nc;
304: PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfAux = 0, f;
305: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE;
306: DMField coordField;
307: DMLabel depthLabel;
308: PetscQuadrature allPoints = NULL;
309: PetscErrorCode ierr;
312: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
313: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);
314: DMGetDimension(dm, &dim);
315: DMPlexGetVTKCellHeight(dm, &minHeight);
316: /* Auxiliary information can only be used with interpolation of field functions */
317: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
318: if (dmAux && !localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector");
319: if (!minHeight && dmAux) {
320: DMLabel spmap;
322: /* If dmAux is a surface, then force the projection to take place over a surface */
323: DMPlexGetSubpointMap(dmAux, &spmap);
324: if (spmap) {
325: DMPlexGetVTKCellHeight(dmAux, &minHeight);
326: auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE;
327: }
328: }
329: }
330: DMPlexGetDepth(dm,&depth);
331: DMPlexGetDepthLabel(dm,&depthLabel);
332: DMPlexGetMaxProjectionHeight(dm, &maxHeight);
333: maxHeight = PetscMax(maxHeight, minHeight);
334: if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim);
335: DMGetFirstLabelEntry_Private(dm, label, numIds, ids, 0, NULL, &prob);
336: if (!prob) {DMGetDS(dm, &prob);}
337: PetscDSGetNumFields(prob, &Nf);
338: DMGetCoordinateDim(dm, &dimEmbed);
339: DMGetSection(dm, §ion);
340: if (dmAux) {
341: DMGetDS(dmAux, &probAux);
342: PetscDSGetNumFields(probAux, &NfAux);
343: }
344: PetscDSGetComponents(prob, &Nc);
345: PetscMalloc2(Nf, &isFE, Nf, &sp);
346: if (maxHeight > 0) {PetscMalloc1(Nf, &cellsp);}
347: else {cellsp = sp;}
348: if (localU && localU != localX) {DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);}
349: /* Get cell dual spaces */
350: for (f = 0; f < Nf; ++f) {
351: PetscObject obj;
352: PetscClassId id;
354: PetscDSGetDiscretization(prob, f, &obj);
355: PetscObjectGetClassId(obj, &id);
356: if (id == PETSCFE_CLASSID) {
357: PetscFE fe = (PetscFE) obj;
359: hasFE = PETSC_TRUE;
360: isFE[f] = PETSC_TRUE;
361: PetscFEGetDualSpace(fe, &cellsp[f]);
362: } else if (id == PETSCFV_CLASSID) {
363: PetscFV fv = (PetscFV) obj;
365: hasFV = PETSC_TRUE;
366: isFE[f] = PETSC_FALSE;
367: PetscFVGetDualSpace(fv, &cellsp[f]);
368: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
369: }
370: DMGetCoordinateField(dm,&coordField);
371: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
372: PetscInt effectiveHeight = auxBd ? minHeight : 0;
373: PetscFE fem, subfem;
374: const PetscReal *points;
375: PetscInt numPoints;
377: if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation");
378: for (f = 0; f < Nf; ++f) {
379: if (!effectiveHeight) {sp[f] = cellsp[f];}
380: else {PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);}
381: }
382: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);
383: PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);
384: PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);
385: for (f = 0; f < Nf; ++f) {
386: if (!isFE[f]) continue;
387: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
388: if (!effectiveHeight) {subfem = fem;}
389: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
390: PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);
391: }
392: for (f = 0; f < NfAux; ++f) {
393: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
394: if (!effectiveHeight || auxBd) {subfem = fem;}
395: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
396: PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);
397: }
398: }
399: /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */
400: for (h = minHeight; h <= maxHeight; h++) {
401: PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight);
402: PetscDS probEff = prob;
403: PetscScalar *values;
404: PetscBool *fieldActive;
405: PetscInt maxDegree;
406: PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues;
407: IS heightIS;
409: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
410: DMGetFirstLabelEntry_Private(dm, label, numIds, ids, h, &lStart, NULL);
411: DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);
412: if (!h) {
413: PetscInt cEndInterior;
415: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
416: pEnd = cEndInterior < 0 ? pEnd : cEndInterior;
417: }
418: if (pEnd <= pStart) {
419: ISDestroy(&heightIS);
420: continue;
421: }
422: /* Compute totDim, the number of dofs in the closure of a point at this height */
423: totDim = 0;
424: for (f = 0; f < Nf; ++f) {
425: if (!effectiveHeight) {
426: sp[f] = cellsp[f];
427: } else {
428: PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);
429: if (!sp[f]) continue;
430: }
431: PetscDualSpaceGetDimension(sp[f], &spDim);
432: totDim += spDim;
433: }
434: DMPlexVecGetClosure(dm, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);
435: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim);
436: if (!totDim) {
437: ISDestroy(&heightIS);
438: continue;
439: }
440: if (effectiveHeight) {PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);}
441: /* Loop over points at this height */
442: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
443: DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);
444: for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
445: if (label) {
446: PetscInt i;
448: for (i = 0; i < numIds; ++i) {
449: IS pointIS, isectIS;
450: const PetscInt *points;
451: PetscInt n;
452: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
453: PetscQuadrature quad = NULL;
455: DMLabelGetStratumIS(label, ids[i], &pointIS);
456: if (!pointIS) continue; /* No points with that id on this process */
457: ISIntersect(pointIS,heightIS,&isectIS);
458: ISDestroy(&pointIS);
459: if (!isectIS) continue;
460: ISGetLocalSize(isectIS, &n);
461: ISGetIndices(isectIS, &points);
462: DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);
463: if (maxDegree <= 1) {
464: DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);
465: }
466: if (!quad) {
467: if (!h && allPoints) {
468: quad = allPoints;
469: allPoints = NULL;
470: } else {
471: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
472: }
473: }
474: DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);
475: for (p = 0; p < n; ++p) {
476: const PetscInt point = points[p];
478: PetscMemzero(values, numValues * sizeof(PetscScalar));
479: PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);
480: 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);
481: if (ierr) {
482: PetscErrorCode ierr2;
483: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
484: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
485:
486: }
487: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);
488: }
489: PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);
490: PetscFEGeomDestroy(&fegeom);
491: PetscQuadratureDestroy(&quad);
492: ISRestoreIndices(isectIS, &points);
493: ISDestroy(&isectIS);
494: }
495: } else {
496: PetscFEGeom *fegeom = NULL, *chunkgeom = NULL;
497: PetscQuadrature quad = NULL;
498: IS pointIS;
500: ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);
501: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
502: if (maxDegree <= 1) {
503: DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);
504: }
505: if (!quad) {
506: if (!h && allPoints) {
507: quad = allPoints;
508: allPoints = NULL;
509: } else {
510: PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);
511: }
512: }
513: DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);
514: for (p = pStart; p < pEnd; ++p) {
515: PetscMemzero(values, numValues * sizeof(PetscScalar));
516: PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);
517: 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);
518: if (ierr) {
519: PetscErrorCode ierr2;
520: ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2);
521: ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2);
522:
523: }
524: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);
525: }
526: PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);
527: PetscFEGeomDestroy(&fegeom);
528: PetscQuadratureDestroy(&quad);
529: ISDestroy(&pointIS);
530: }
531: ISDestroy(&heightIS);
532: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
533: DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);
534: }
535: /* Cleanup */
536: if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) {
537: PetscInt effectiveHeight = auxBd ? minHeight : 0;
538: PetscFE fem, subfem;
540: for (f = 0; f < Nf; ++f) {
541: if (!isFE[f]) continue;
542: PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);
543: if (!effectiveHeight) {subfem = fem;}
544: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
545: PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);
546: }
547: for (f = 0; f < NfAux; ++f) {
548: PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);
549: if (!effectiveHeight || auxBd) {subfem = fem;}
550: else {PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);}
551: PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);
552: }
553: PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);
554: }
555: PetscQuadratureDestroy(&allPoints);
556: PetscFree2(isFE, sp);
557: if (maxHeight > 0) {PetscFree(cellsp);}
558: return(0);
559: }
561: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
562: {
566: DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
567: return(0);
568: }
570: 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)
571: {
575: DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);
576: return(0);
577: }
579: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU,
580: void (**funcs)(PetscInt, PetscInt, PetscInt,
581: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
582: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
583: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
584: InsertMode mode, Vec localX)
585: {
589: DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
590: return(0);
591: }
593: PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU,
594: void (**funcs)(PetscInt, PetscInt, PetscInt,
595: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
596: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
597: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
598: InsertMode mode, Vec localX)
599: {
603: DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);
604: return(0);
605: }