Actual source code: dmfieldds.c
1: #include <petsc/private/dmfieldimpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscfe.h>
4: #include <petscdmplex.h>
5: #include <petscds.h>
7: typedef struct _n_DMField_DS {
8: PetscBool multifieldVec;
9: PetscInt height; /* Point height at which we want values and number of discretizations */
10: PetscInt fieldNum; /* Number in DS of field which we evaluate */
11: PetscObject *disc; /* Discretizations of this field at each height */
12: Vec vec; /* Field values */
13: DM dmDG; /* DM for the DG values */
14: PetscObject *discDG; /* DG Discretizations of this field at each height */
15: Vec vecDG; /* DG Field values */
16: } DMField_DS;
18: static PetscErrorCode DMFieldDestroy_DS(DMField field)
19: {
20: DMField_DS *dsfield = (DMField_DS *)field->data;
21: PetscInt i;
23: PetscFunctionBegin;
24: PetscCall(VecDestroy(&dsfield->vec));
25: for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->disc[i]));
26: PetscCall(PetscFree(dsfield->disc));
27: PetscCall(VecDestroy(&dsfield->vecDG));
28: if (dsfield->discDG)
29: for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->discDG[i]));
30: PetscCall(PetscFree(dsfield->discDG));
31: PetscCall(PetscFree(dsfield));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode DMFieldView_DS(DMField field, PetscViewer viewer)
36: {
37: DMField_DS *dsfield = (DMField_DS *)field->data;
38: PetscObject disc;
39: PetscBool iascii;
41: PetscFunctionBegin;
42: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
43: disc = dsfield->disc[0];
44: if (iascii) {
45: PetscCall(PetscViewerASCIIPrintf(viewer, "PetscDS field %" PetscInt_FMT "\n", dsfield->fieldNum));
46: PetscCall(PetscViewerASCIIPushTab(viewer));
47: PetscCall(PetscObjectView(disc, viewer));
48: PetscCall(PetscViewerASCIIPopTab(viewer));
49: }
50: PetscCall(PetscViewerASCIIPushTab(viewer));
51: PetscCheck(!dsfield->multifieldVec, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "View of subfield not implemented yet");
52: PetscCall(VecView(dsfield->vec, viewer));
53: PetscCall(PetscViewerASCIIPopTab(viewer));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject discList[], PetscObject *disc)
58: {
59: PetscFunctionBegin;
60: PetscCheck(height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " must be non-negative", height);
61: if (!discList[height]) {
62: PetscClassId id;
64: PetscCall(PetscObjectGetClassId(discList[0], &id));
65: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height]));
66: }
67: *disc = discList[height];
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: /* y[m,c] = A[m,n,c] . b[n] */
72: #define DMFieldDSdot(y, A, b, m, n, c, cast) \
73: do { \
74: PetscInt _i, _j, _k; \
75: for (_i = 0; _i < (m); _i++) { \
76: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] = 0.; \
77: for (_j = 0; _j < (n); _j++) { \
78: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \
79: } \
80: } \
81: } while (0)
83: /*
84: Since this is used for coordinates, we need to allow for the possibility that values come from multiple sections/Vecs, so that we can have DG version of the coordinates for periodicity. This reproduces DMPlexGetCellCoordinates_Internal().
85: */
86: static PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
87: {
88: DMField_DS *dsfield = (DMField_DS *)field->data;
89: DM fdm = dsfield->dmDG;
90: PetscSection s = NULL;
91: const PetscScalar *cvalues;
92: PetscInt pStart, pEnd;
94: PetscFunctionBeginHot;
95: *isDG = PETSC_FALSE;
96: *Nc = 0;
97: *array = NULL;
98: *values = NULL;
99: /* Check for cellwise section */
100: if (fdm) PetscCall(DMGetLocalSection(fdm, &s));
101: if (!s) goto cg;
102: /* Check that the cell exists in the cellwise section */
103: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
104: if (cell < pStart || cell >= pEnd) goto cg;
105: /* Check for cellwise coordinates for this cell */
106: PetscCall(PetscSectionGetDof(s, cell, Nc));
107: if (!*Nc) goto cg;
108: /* Check for cellwise coordinates */
109: if (!dsfield->vecDG) goto cg;
110: /* Get cellwise coordinates */
111: PetscCall(VecGetArrayRead(dsfield->vecDG, array));
112: PetscCall(DMPlexPointLocalRead(fdm, cell, *array, &cvalues));
113: PetscCall(DMGetWorkArray(fdm, *Nc, MPIU_SCALAR, values));
114: PetscCall(PetscArraycpy(*values, cvalues, *Nc));
115: PetscCall(VecRestoreArrayRead(dsfield->vecDG, array));
116: *isDG = PETSC_TRUE;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: cg:
119: /* Use continuous values */
120: PetscCall(DMFieldGetDM(field, &fdm));
121: PetscCall(DMGetLocalSection(fdm, &s));
122: PetscCall(PetscSectionGetField(s, dsfield->fieldNum, &s));
123: PetscCall(DMPlexVecGetClosure(fdm, s, dsfield->vec, cell, Nc, values));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode DMFieldRestoreClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
128: {
129: DMField_DS *dsfield = (DMField_DS *)field->data;
130: DM fdm;
131: PetscSection s;
133: PetscFunctionBeginHot;
134: if (*isDG) {
135: PetscCall(DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values));
136: } else {
137: PetscCall(DMFieldGetDM(field, &fdm));
138: PetscCall(DMGetLocalSection(fdm, &s));
139: PetscCall(DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, values));
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
145: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
146: {
147: DMField_DS *dsfield = (DMField_DS *)field->data;
148: DM dm;
149: PetscObject disc;
150: PetscClassId classid;
151: PetscInt nq, nc, dim, meshDim, numCells;
152: PetscSection section;
153: const PetscReal *qpoints;
154: PetscBool isStride;
155: const PetscInt *points = NULL;
156: PetscInt sfirst = -1, stride = -1;
158: PetscFunctionBeginHot;
159: dm = field->dm;
160: nc = field->numComponents;
161: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL));
162: PetscCall(DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc));
163: PetscCall(DMGetDimension(dm, &meshDim));
164: PetscCall(DMGetLocalSection(dm, §ion));
165: PetscCall(PetscSectionGetField(section, dsfield->fieldNum, §ion));
166: PetscCall(PetscObjectGetClassId(disc, &classid));
167: /* TODO: batch */
168: PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
169: PetscCall(ISGetLocalSize(pointIS, &numCells));
170: if (isStride) PetscCall(ISStrideGetInfo(pointIS, &sfirst, &stride));
171: else PetscCall(ISGetIndices(pointIS, &points));
172: if (classid == PETSCFE_CLASSID) {
173: PetscFE fe = (PetscFE)disc;
174: PetscInt feDim, i;
175: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
176: PetscTabulation T;
178: PetscCall(PetscFEGetDimension(fe, &feDim));
179: PetscCall(PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T));
180: for (i = 0; i < numCells; i++) {
181: PetscInt c = isStride ? (sfirst + i * stride) : points[i];
182: PetscInt closureSize;
183: const PetscScalar *array;
184: PetscScalar *elem = NULL;
185: PetscBool isDG;
187: PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
188: if (B) {
189: /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
190: if (type == PETSC_SCALAR) {
191: PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i];
193: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
194: } else {
195: PetscReal *cB = &((PetscReal *)B)[nc * nq * i];
197: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
198: }
199: }
200: if (D) {
201: if (type == PETSC_SCALAR) {
202: PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i];
204: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
205: } else {
206: PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i];
208: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
209: }
210: }
211: if (H) {
212: if (type == PETSC_SCALAR) {
213: PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i];
215: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
216: } else {
217: PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i];
219: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
220: }
221: }
222: PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
223: }
224: PetscCall(PetscTabulationDestroy(&T));
225: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented");
226: if (!isStride) PetscCall(ISRestoreIndices(pointIS, &points));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
231: {
232: DMField_DS *dsfield = (DMField_DS *)field->data;
233: PetscSF cellSF = NULL;
234: const PetscSFNode *cells;
235: PetscInt c, nFound, numCells, feDim, nc;
236: const PetscInt *cellDegrees;
237: const PetscScalar *pointsArray;
238: PetscScalar *cellPoints;
239: PetscInt gatherSize, gatherMax;
240: PetscInt dim, dimR, offset;
241: MPI_Datatype pointType;
242: PetscObject cellDisc;
243: PetscFE cellFE;
244: PetscClassId discID;
245: PetscReal *coordsReal, *coordsRef;
246: PetscSection section;
247: PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
248: PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
249: PetscReal *v, *J, *invJ, *detJ;
250: PetscMPIInt idim;
252: PetscFunctionBegin;
253: nc = field->numComponents;
254: PetscCall(DMGetLocalSection(field->dm, §ion));
255: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
256: PetscCall(PetscObjectGetClassId(cellDisc, &discID));
257: PetscCheck(discID == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization type not supported");
258: cellFE = (PetscFE)cellDisc;
259: PetscCall(PetscFEGetDimension(cellFE, &feDim));
260: PetscCall(DMGetCoordinateDim(field->dm, &dim));
261: PetscCall(DMGetDimension(field->dm, &dimR));
262: PetscCall(DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF));
263: PetscCall(PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells));
264: for (c = 0; c < nFound; c++) PetscCheck(cells[c].index >= 0, PetscObjectComm((PetscObject)points), PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " could not be located", c);
265: PetscCall(PetscSFComputeDegreeBegin(cellSF, &cellDegrees));
266: PetscCall(PetscSFComputeDegreeEnd(cellSF, &cellDegrees));
267: for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
268: gatherMax = PetscMax(gatherMax, cellDegrees[c]);
269: gatherSize += cellDegrees[c];
270: }
271: PetscCall(PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef));
272: PetscCall(PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ));
273: if (datatype == PETSC_SCALAR) PetscCall(PetscMalloc3((B ? (size_t)nc * gatherSize : 0), &cellBs, (D ? (size_t)nc * dim * gatherSize : 0), &cellDs, (H ? (size_t)nc * dim * dim * gatherSize : 0), &cellHs));
274: else PetscCall(PetscMalloc3((B ? (size_t)nc * gatherSize : 0), &cellBr, (D ? (size_t)nc * dim * gatherSize : 0), &cellDr, (H ? (size_t)nc * dim * dim * gatherSize : 0), &cellHr));
276: PetscCall(PetscMPIIntCast(dim, &idim));
277: PetscCallMPI(MPI_Type_contiguous(idim, MPIU_SCALAR, &pointType));
278: PetscCallMPI(MPI_Type_commit(&pointType));
279: PetscCall(VecGetArrayRead(points, &pointsArray));
280: PetscCall(PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints));
281: PetscCall(PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints));
282: PetscCall(VecRestoreArrayRead(points, &pointsArray));
283: for (c = 0, offset = 0; c < numCells; c++) {
284: PetscInt nq = cellDegrees[c], p;
286: if (nq) {
287: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
288: PetscTabulation T;
289: PetscQuadrature quad;
290: const PetscScalar *array;
291: PetscScalar *elem = NULL;
292: PetscReal *quadPoints;
293: PetscBool isDG;
294: PetscInt closureSize, d, e, f, g;
296: for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
297: PetscCall(DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef));
298: PetscCall(PetscFECreateTabulation(cellFE, 1, nq, coordsRef, K, &T));
299: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
300: PetscCall(PetscMalloc1(dimR * nq, &quadPoints));
301: for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
302: PetscCall(PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL));
303: PetscCall(DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ));
304: PetscCall(PetscQuadratureDestroy(&quad));
305: PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
306: if (B) {
307: if (datatype == PETSC_SCALAR) {
308: PetscScalar *cB = &cellBs[nc * offset];
310: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
311: } else {
312: PetscReal *cB = &cellBr[nc * offset];
314: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
315: }
316: }
317: if (D) {
318: if (datatype == PETSC_SCALAR) {
319: PetscScalar *cD = &cellDs[nc * dim * offset];
321: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
322: for (p = 0; p < nq; p++) {
323: for (g = 0; g < nc; g++) {
324: PetscScalar vs[3];
326: for (d = 0; d < dimR; d++) {
327: vs[d] = 0.;
328: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
329: }
330: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = vs[d];
331: }
332: }
333: } else {
334: PetscReal *cD = &cellDr[nc * dim * offset];
336: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
337: for (p = 0; p < nq; p++) {
338: for (g = 0; g < nc; g++) {
339: for (d = 0; d < dimR; d++) {
340: v[d] = 0.;
341: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
342: }
343: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = v[d];
344: }
345: }
346: }
347: }
348: if (H) {
349: if (datatype == PETSC_SCALAR) {
350: PetscScalar *cH = &cellHs[nc * dim * dim * offset];
352: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
353: for (p = 0; p < nq; p++) {
354: for (g = 0; g < nc * dimR; g++) {
355: PetscScalar vs[3];
357: for (d = 0; d < dimR; d++) {
358: vs[d] = 0.;
359: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
360: }
361: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = vs[d];
362: }
363: for (g = 0; g < nc; g++) {
364: for (f = 0; f < dimR; f++) {
365: PetscScalar vs[3];
367: for (d = 0; d < dimR; d++) {
368: vs[d] = 0.;
369: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
370: }
371: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
372: }
373: }
374: }
375: } else {
376: PetscReal *cH = &cellHr[nc * dim * dim * offset];
378: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
379: for (p = 0; p < nq; p++) {
380: for (g = 0; g < nc * dimR; g++) {
381: for (d = 0; d < dimR; d++) {
382: v[d] = 0.;
383: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
384: }
385: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = v[d];
386: }
387: for (g = 0; g < nc; g++) {
388: for (f = 0; f < dimR; f++) {
389: for (d = 0; d < dimR; d++) {
390: v[d] = 0.;
391: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
392: }
393: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
394: }
395: }
396: }
397: }
398: }
399: PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
400: PetscCall(PetscTabulationDestroy(&T));
401: }
402: offset += nq;
403: }
404: {
405: MPI_Datatype origtype;
407: if (datatype == PETSC_SCALAR) {
408: origtype = MPIU_SCALAR;
409: } else {
410: origtype = MPIU_REAL;
411: }
412: if (B) {
413: MPI_Datatype Btype;
415: PetscCall(PetscMPIIntCast(nc, &idim));
416: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Btype));
417: PetscCallMPI(MPI_Type_commit(&Btype));
418: PetscCall(PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
419: PetscCall(PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
420: PetscCallMPI(MPI_Type_free(&Btype));
421: }
422: if (D) {
423: MPI_Datatype Dtype;
425: PetscCall(PetscMPIIntCast(nc * dim, &idim));
426: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Dtype));
427: PetscCallMPI(MPI_Type_commit(&Dtype));
428: PetscCall(PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
429: PetscCall(PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
430: PetscCallMPI(MPI_Type_free(&Dtype));
431: }
432: if (H) {
433: MPI_Datatype Htype;
435: PetscCall(PetscMPIIntCast(nc * dim * dim, &idim));
436: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Htype));
437: PetscCallMPI(MPI_Type_commit(&Htype));
438: PetscCall(PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
439: PetscCall(PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
440: PetscCallMPI(MPI_Type_free(&Htype));
441: }
442: }
443: PetscCall(PetscFree4(v, J, invJ, detJ));
444: PetscCall(PetscFree3(cellBr, cellDr, cellHr));
445: PetscCall(PetscFree3(cellBs, cellDs, cellHs));
446: PetscCall(PetscFree3(cellPoints, coordsReal, coordsRef));
447: PetscCallMPI(MPI_Type_free(&pointType));
448: PetscCall(PetscSFDestroy(&cellSF));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
453: {
454: DMField_DS *dsfield = (DMField_DS *)field->data;
455: PetscInt h, imin;
456: PetscInt dim;
457: PetscClassId id;
458: PetscQuadrature quad = NULL;
459: PetscInt maxDegree;
460: PetscFEGeom *geom;
461: PetscInt Nq, Nc, dimC, qNc, N;
462: PetscInt numPoints;
463: void *qB = NULL, *qD = NULL, *qH = NULL;
464: const PetscReal *weights;
465: MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
466: PetscObject disc;
467: DMField coordField;
469: PetscFunctionBegin;
470: Nc = field->numComponents;
471: PetscCall(DMGetCoordinateDim(field->dm, &dimC));
472: PetscCall(DMGetDimension(field->dm, &dim));
473: PetscCall(ISGetLocalSize(pointIS, &numPoints));
474: PetscCall(ISGetMinMax(pointIS, &imin, NULL));
475: for (h = 0; h < dsfield->height; h++) {
476: PetscInt hEnd;
478: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
479: if (imin < hEnd) break;
480: }
481: dim -= h;
482: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
483: PetscCall(PetscObjectGetClassId(disc, &id));
484: PetscCheck(id == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported");
485: PetscCall(DMGetCoordinateField(field->dm, &coordField));
486: PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
487: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
488: if (!quad) PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, &quad));
489: PetscCheck(quad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages");
490: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom));
491: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights));
492: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components");
493: N = numPoints * Nq * Nc;
494: if (B) PetscCall(DMGetWorkArray(field->dm, N, mpitype, &qB));
495: if (D) PetscCall(DMGetWorkArray(field->dm, N * dimC, mpitype, &qD));
496: if (H) PetscCall(DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
497: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH));
498: if (B) {
499: PetscInt i, j, k;
501: if (type == PETSC_SCALAR) {
502: PetscScalar *sB = (PetscScalar *)B;
503: PetscScalar *sqB = (PetscScalar *)qB;
505: for (i = 0; i < numPoints; i++) {
506: PetscReal vol = 0.;
508: for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.;
509: for (k = 0; k < Nq; k++) {
510: vol += geom->detJ[i * Nq + k] * weights[k];
511: for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j];
512: }
513: for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
514: }
515: } else {
516: PetscReal *rB = (PetscReal *)B;
517: PetscReal *rqB = (PetscReal *)qB;
519: for (i = 0; i < numPoints; i++) {
520: PetscReal vol = 0.;
522: for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.;
523: for (k = 0; k < Nq; k++) {
524: vol += geom->detJ[i * Nq + k] * weights[k];
525: for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j];
526: }
527: for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
528: }
529: }
530: }
531: if (D) {
532: PetscInt i, j, k, l, m;
534: if (type == PETSC_SCALAR) {
535: PetscScalar *sD = (PetscScalar *)D;
536: PetscScalar *sqD = (PetscScalar *)qD;
538: for (i = 0; i < numPoints; i++) {
539: PetscReal vol = 0.;
541: for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
542: for (k = 0; k < Nq; k++) {
543: vol += geom->detJ[i * Nq + k] * weights[k];
544: for (j = 0; j < Nc; j++) {
545: PetscScalar pD[3] = {0., 0., 0.};
547: for (l = 0; l < dimC; l++) {
548: for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m];
549: }
550: for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
551: }
552: }
553: for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
554: }
555: } else {
556: PetscReal *rD = (PetscReal *)D;
557: PetscReal *rqD = (PetscReal *)qD;
559: for (i = 0; i < numPoints; i++) {
560: PetscReal vol = 0.;
562: for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
563: for (k = 0; k < Nq; k++) {
564: vol += geom->detJ[i * Nq + k] * weights[k];
565: for (j = 0; j < Nc; j++) {
566: PetscReal pD[3] = {0., 0., 0.};
568: for (l = 0; l < dimC; l++) {
569: for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m];
570: }
571: for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
572: }
573: }
574: for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
575: }
576: }
577: }
578: if (H) {
579: PetscInt i, j, k, l, m, q, r;
581: if (type == PETSC_SCALAR) {
582: PetscScalar *sH = (PetscScalar *)H;
583: PetscScalar *sqH = (PetscScalar *)qH;
585: for (i = 0; i < numPoints; i++) {
586: PetscReal vol = 0.;
588: for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.;
589: for (k = 0; k < Nq; k++) {
590: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
592: vol += geom->detJ[i * Nq + k] * weights[k];
593: for (j = 0; j < Nc; j++) {
594: PetscScalar pH[3][3] = {
595: {0., 0., 0.},
596: {0., 0., 0.},
597: {0., 0., 0.}
598: };
599: const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];
601: for (l = 0; l < dimC; l++) {
602: for (m = 0; m < dimC; m++) {
603: for (q = 0; q < dim; q++) {
604: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
605: }
606: }
607: }
608: for (l = 0; l < dimC; l++) {
609: for (m = 0; m < dimC; m++) sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
610: }
611: }
612: }
613: for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
614: }
615: } else {
616: PetscReal *rH = (PetscReal *)H;
617: PetscReal *rqH = (PetscReal *)qH;
619: for (i = 0; i < numPoints; i++) {
620: PetscReal vol = 0.;
622: for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.;
623: for (k = 0; k < Nq; k++) {
624: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
626: vol += geom->detJ[i * Nq + k] * weights[k];
627: for (j = 0; j < Nc; j++) {
628: PetscReal pH[3][3] = {
629: {0., 0., 0.},
630: {0., 0., 0.},
631: {0., 0., 0.}
632: };
633: const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];
635: for (l = 0; l < dimC; l++) {
636: for (m = 0; m < dimC; m++) {
637: for (q = 0; q < dim; q++) {
638: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
639: }
640: }
641: }
642: for (l = 0; l < dimC; l++) {
643: for (m = 0; m < dimC; m++) rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
644: }
645: }
646: }
647: for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
648: }
649: }
650: }
651: if (B) PetscCall(DMRestoreWorkArray(field->dm, N, mpitype, &qB));
652: if (D) PetscCall(DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD));
653: if (H) PetscCall(DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
654: PetscCall(PetscFEGeomDestroy(&geom));
655: PetscCall(PetscQuadratureDestroy(&quad));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
660: {
661: DMField_DS *dsfield;
662: PetscObject disc;
663: PetscInt h, imin, imax;
664: PetscClassId id;
666: PetscFunctionBegin;
667: dsfield = (DMField_DS *)field->data;
668: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
669: if (imin >= imax) {
670: h = 0;
671: } else {
672: for (h = 0; h < dsfield->height; h++) {
673: PetscInt hEnd;
675: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
676: if (imin < hEnd) break;
677: }
678: }
679: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
680: PetscCall(PetscObjectGetClassId(disc, &id));
681: if (id == PETSCFE_CLASSID) {
682: PetscFE fe = (PetscFE)disc;
683: PetscSpace sp;
685: PetscCall(PetscFEGetBasisSpace(fe, &sp));
686: PetscCall(PetscSpaceGetDegree(sp, minDegree, maxDegree));
687: }
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
692: {
693: DM dm = field->dm;
694: const PetscInt *points;
695: DMPolytopeType ct;
696: PetscInt dim, n;
697: PetscBool isplex;
699: PetscFunctionBegin;
700: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
701: PetscCall(ISGetLocalSize(pointIS, &n));
702: if (isplex && n) {
703: PetscCall(DMGetDimension(dm, &dim));
704: PetscCall(ISGetIndices(pointIS, &points));
705: PetscCall(DMPlexGetCellType(dm, points[0], &ct));
706: switch (ct) {
707: case DM_POLYTOPE_TRIANGLE:
708: case DM_POLYTOPE_TETRAHEDRON:
709: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad));
710: break;
711: default:
712: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
713: }
714: PetscCall(ISRestoreIndices(pointIS, &points));
715: } else PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, quad));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
720: {
721: PetscInt h, dim, imax, imin, cellHeight;
722: DM dm;
723: DMField_DS *dsfield;
724: PetscObject disc;
725: PetscFE fe;
726: PetscClassId id;
728: PetscFunctionBegin;
729: dm = field->dm;
730: dsfield = (DMField_DS *)field->data;
731: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
732: PetscCall(DMGetDimension(dm, &dim));
733: for (h = 0; h <= dim; h++) {
734: PetscInt hStart, hEnd;
736: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
737: if (imax >= hStart && imin < hEnd) break;
738: }
739: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
740: h -= cellHeight;
741: *quad = NULL;
742: if (h < dsfield->height) {
743: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
744: PetscCall(PetscObjectGetClassId(disc, &id));
745: if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
746: fe = (PetscFE)disc;
747: PetscCall(PetscFEGetQuadrature(fe, quad));
748: PetscCall(PetscObjectReference((PetscObject)*quad));
749: }
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
754: {
755: const PetscInt *points;
756: PetscInt p, dim, dE, numFaces, Nq;
757: PetscInt maxDegree;
758: DMLabel depthLabel;
759: IS cellIS;
760: DM dm = field->dm;
762: PetscFunctionBegin;
763: dim = geom->dim;
764: dE = geom->dimEmbed;
765: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
766: PetscCall(DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS));
767: PetscCall(DMFieldGetDegree(field, cellIS, NULL, &maxDegree));
768: PetscCall(ISGetIndices(pointIS, &points));
769: numFaces = geom->numCells;
770: Nq = geom->numPoints;
771: /* First, set local faces and flip normals so that they are outward for the first supporting cell */
772: for (p = 0; p < numFaces; p++) {
773: PetscInt point = points[p];
774: PetscInt suppSize, s, coneSize, c, numChildren;
775: const PetscInt *supp;
777: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
778: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
779: PetscCall(DMPlexGetSupportSize(dm, point, &suppSize));
780: PetscCheck(suppSize <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, suppSize);
781: if (!suppSize) continue;
782: PetscCall(DMPlexGetSupport(dm, point, &supp));
783: for (s = 0; s < suppSize; ++s) {
784: const PetscInt *cone, *ornt;
786: PetscCall(DMPlexGetConeSize(dm, supp[s], &coneSize));
787: PetscCall(DMPlexGetOrientedCone(dm, supp[s], &cone, &ornt));
788: for (c = 0; c < coneSize; ++c)
789: if (cone[c] == point) break;
790: geom->face[p][s * 2 + 0] = c;
791: geom->face[p][s * 2 + 1] = ornt[c];
792: PetscCall(DMPlexRestoreOrientedCone(dm, supp[s], &cone, &ornt));
793: PetscCheck(c != coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid connectivity: point %" PetscInt_FMT " not found in cone of support point %" PetscInt_FMT, point, supp[s]);
794: }
795: if (geom->face[p][1] < 0) {
796: PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;
798: for (q = 0; q < Np; ++q)
799: for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
800: }
801: }
802: if (maxDegree <= 1) {
803: PetscQuadrature cellQuad = NULL;
804: PetscInt numCells, offset, *cells;
805: PetscFEGeom *cellGeom;
806: IS suppIS;
808: if (quad) {
809: DM dm;
810: PetscReal *points, *weights;
811: PetscInt tdim, Nc, Np;
813: PetscCall(DMFieldGetDM(field, &dm));
814: PetscCall(DMGetDimension(dm, &tdim));
815: if (tdim > dim) {
816: // Make a compatible cell quadrature (points don't matter since its affine)
817: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
818: PetscCall(PetscQuadratureGetData(quad, NULL, &Nc, &Np, NULL, NULL));
819: PetscCall(PetscCalloc1((dim + 1) * Np, &points));
820: PetscCall(PetscCalloc1(Nc * Np, &weights));
821: PetscCall(PetscQuadratureSetData(cellQuad, dim + 1, Nc, Np, points, weights));
822: } else {
823: // TODO J will be wrong here, but other things need to be fixed
824: // This path comes from calling DMProjectBdFieldLabelLocal() in Plex ex5
825: PetscCall(PetscObjectReference((PetscObject)quad));
826: cellQuad = quad;
827: }
828: }
829: for (p = 0, numCells = 0; p < numFaces; p++) {
830: PetscInt point = points[p];
831: PetscInt numSupp, numChildren;
833: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
834: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
835: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
836: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
837: numCells += numSupp;
838: }
839: PetscCall(PetscMalloc1(numCells, &cells));
840: for (p = 0, offset = 0; p < numFaces; p++) {
841: PetscInt point = points[p];
842: PetscInt numSupp, s;
843: const PetscInt *supp;
845: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
846: PetscCall(DMPlexGetSupport(dm, point, &supp));
847: for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
848: }
849: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
850: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom));
851: for (p = 0, offset = 0; p < numFaces; p++) {
852: PetscInt point = points[p];
853: PetscInt numSupp, s, q;
854: const PetscInt *supp;
856: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
857: PetscCall(DMPlexGetSupport(dm, point, &supp));
858: for (s = 0; s < numSupp; s++, offset++) {
859: for (q = 0; q < Nq * dE * dE; q++) {
860: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
861: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
862: }
863: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
864: }
865: }
866: PetscCall(PetscFEGeomDestroy(&cellGeom));
867: PetscCall(PetscQuadratureDestroy(&cellQuad));
868: PetscCall(ISDestroy(&suppIS));
869: PetscCall(PetscFree(cells));
870: } else {
871: DMField_DS *dsfield = (DMField_DS *)field->data;
872: PetscObject faceDisc, cellDisc;
873: PetscClassId faceId, cellId;
874: PetscDualSpace dsp;
875: DM K;
876: DMPolytopeType ct;
877: PetscInt(*co)[2][3];
878: PetscInt coneSize;
879: PetscInt **counts;
880: PetscInt f, i, o, q, s;
881: PetscBool found = PETSC_FALSE;
882: const PetscInt *coneK;
883: PetscInt eStart, minOrient, maxOrient, numOrient;
884: PetscInt *orients;
885: PetscReal **orientPoints;
886: PetscReal *cellPoints;
887: PetscReal *dummyWeights;
888: PetscQuadrature cellQuad = NULL;
890: PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc));
891: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
892: PetscCall(PetscObjectGetClassId(faceDisc, &faceId));
893: PetscCall(PetscObjectGetClassId(cellDisc, &cellId));
894: PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported");
895: PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp));
896: PetscCall(PetscDualSpaceGetDM(dsp, &K));
897: PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL));
898: PetscCall(DMPlexGetCellType(K, eStart, &ct));
899: PetscCall(DMPlexGetConeSize(K, 0, &coneSize));
900: PetscCall(DMPlexGetCone(K, 0, &coneK));
901: PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts));
902: PetscCall(PetscMalloc1(dE * Nq, &cellPoints));
903: PetscCall(PetscMalloc1(Nq, &dummyWeights));
904: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
905: PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights));
906: minOrient = PETSC_INT_MAX;
907: maxOrient = PETSC_INT_MIN;
908: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
909: PetscInt point = points[p];
910: PetscInt numSupp, numChildren;
911: const PetscInt *supp;
913: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
914: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
915: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
916: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
917: PetscCall(DMPlexGetSupport(dm, point, &supp));
918: for (s = 0; s < numSupp; s++) {
919: PetscInt cell = supp[s];
920: PetscInt numCone;
921: const PetscInt *cone, *orient;
923: PetscCall(DMPlexGetConeSize(dm, cell, &numCone));
924: // When we extract submeshes, we hang cells from the side that are not fully realized. We ignore these
925: if (numCone == 1) {
926: co[p][s][0] = -1;
927: co[p][s][1] = -1;
928: co[p][s][2] = -1;
929: continue;
930: }
931: PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
932: PetscCall(DMPlexGetCone(dm, cell, &cone));
933: PetscCall(DMPlexGetConeOrientation(dm, cell, &orient));
934: for (f = 0; f < coneSize; f++) {
935: if (cone[f] == point) break;
936: }
937: co[p][s][0] = f;
938: co[p][s][1] = orient[f];
939: co[p][s][2] = cell;
940: minOrient = PetscMin(minOrient, orient[f]);
941: maxOrient = PetscMax(maxOrient, orient[f]);
942: found = PETSC_TRUE;
943: }
944: for (; s < 2; s++) {
945: co[p][s][0] = -1;
946: co[p][s][1] = -1;
947: co[p][s][2] = -1;
948: }
949: }
950: numOrient = found ? maxOrient + 1 - minOrient : 0;
951: PetscCall(DMPlexGetCone(K, 0, &coneK));
952: /* count all (face,orientation) doubles that appear */
953: PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints));
954: for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f]));
955: for (p = 0; p < numFaces; p++) {
956: for (s = 0; s < 2; s++) {
957: if (co[p][s][0] >= 0) {
958: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
959: orients[co[p][s][1] - minOrient]++;
960: }
961: }
962: }
963: for (o = 0; o < numOrient; o++) {
964: if (orients[o]) {
965: PetscInt orient = o + minOrient;
966: PetscInt q;
968: PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o]));
969: /* rotate the quadrature points appropriately */
970: switch (ct) {
971: case DM_POLYTOPE_POINT:
972: break;
973: case DM_POLYTOPE_SEGMENT:
974: if (orient == -2 || orient == 1) {
975: for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
976: } else {
977: for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
978: }
979: break;
980: case DM_POLYTOPE_TRIANGLE:
981: for (q = 0; q < Nq; q++) {
982: PetscReal lambda[3];
983: PetscReal lambdao[3];
985: /* convert to barycentric */
986: lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
987: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
988: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
989: if (orient >= 0) {
990: for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
991: } else {
992: for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
993: }
994: /* convert to coordinates */
995: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
996: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
997: }
998: break;
999: case DM_POLYTOPE_QUADRILATERAL:
1000: for (q = 0; q < Nq; q++) {
1001: PetscReal xi[2], xio[2];
1002: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
1004: xi[0] = geom->xi[2 * q];
1005: xi[1] = geom->xi[2 * q + 1];
1006: switch (oabs) {
1007: case 1:
1008: xio[0] = xi[1];
1009: xio[1] = -xi[0];
1010: break;
1011: case 2:
1012: xio[0] = -xi[0];
1013: xio[1] = -xi[1];
1014: break;
1015: case 3:
1016: xio[0] = -xi[1];
1017: xio[1] = xi[0];
1018: break;
1019: case 0:
1020: default:
1021: xio[0] = xi[0];
1022: xio[1] = xi[1];
1023: break;
1024: }
1025: if (orient < 0) xio[0] = -xio[0];
1026: orientPoints[o][2 * q + 0] = xio[0];
1027: orientPoints[o][2 * q + 1] = xio[1];
1028: }
1029: break;
1030: default:
1031: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
1032: }
1033: }
1034: }
1035: for (f = 0; f < coneSize; f++) {
1036: PetscInt face = coneK[f];
1037: PetscReal v0[3];
1038: PetscReal J[9], detJ;
1039: PetscInt numCells, offset;
1040: PetscInt *cells;
1041: IS suppIS;
1043: PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ));
1044: for (o = 0; o <= numOrient; o++) {
1045: PetscFEGeom *cellGeom;
1047: if (!counts[f][o]) continue;
1048: /* If this (face,orientation) double appears,
1049: * convert the face quadrature points into volume quadrature points */
1050: for (q = 0; q < Nq; q++) {
1051: PetscReal xi0[3] = {-1., -1., -1.};
1053: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1054: }
1055: for (p = 0, numCells = 0; p < numFaces; p++) {
1056: for (s = 0; s < 2; s++) {
1057: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1058: }
1059: }
1060: PetscCall(PetscMalloc1(numCells, &cells));
1061: for (p = 0, offset = 0; p < numFaces; p++) {
1062: for (s = 0; s < 2; s++) {
1063: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1064: }
1065: }
1066: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
1067: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom));
1068: for (p = 0, offset = 0; p < numFaces; p++) {
1069: for (s = 0; s < 2; s++) {
1070: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1071: for (q = 0; q < Nq * dE * dE; q++) {
1072: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1073: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1074: }
1075: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1076: offset++;
1077: }
1078: }
1079: }
1080: PetscCall(PetscFEGeomDestroy(&cellGeom));
1081: PetscCall(ISDestroy(&suppIS));
1082: PetscCall(PetscFree(cells));
1083: }
1084: }
1085: for (o = 0; o < numOrient; o++) {
1086: if (orients[o]) PetscCall(PetscFree(orientPoints[o]));
1087: }
1088: PetscCall(PetscFree2(orients, orientPoints));
1089: PetscCall(PetscQuadratureDestroy(&cellQuad));
1090: for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f]));
1091: PetscCall(PetscFree2(co, counts));
1092: }
1093: PetscCall(ISRestoreIndices(pointIS, &points));
1094: PetscCall(ISDestroy(&cellIS));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1099: {
1100: PetscFunctionBegin;
1101: field->ops->destroy = DMFieldDestroy_DS;
1102: field->ops->evaluate = DMFieldEvaluate_DS;
1103: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1104: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1105: field->ops->getDegree = DMFieldGetDegree_DS;
1106: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1107: field->ops->view = DMFieldView_DS;
1108: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1113: {
1114: DMField_DS *dsfield;
1116: PetscFunctionBegin;
1117: PetscCall(PetscNew(&dsfield));
1118: field->data = dsfield;
1119: PetscCall(DMFieldInitialize_DS(field));
1120: PetscFunctionReturn(PETSC_SUCCESS);
1121: }
1123: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1124: {
1125: DMField b;
1126: DMField_DS *dsfield;
1127: PetscObject disc = NULL, discDG = NULL;
1128: PetscSection section;
1129: PetscBool isContainer = PETSC_FALSE;
1130: PetscClassId id = -1;
1131: PetscInt numComponents = -1, dsNumFields;
1133: PetscFunctionBegin;
1134: PetscCall(DMGetLocalSection(dm, §ion));
1135: PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents));
1136: PetscCall(DMGetNumFields(dm, &dsNumFields));
1137: if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc));
1138: if (dsNumFields && dmDG) {
1139: PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG));
1140: PetscCall(PetscObjectReference(discDG));
1141: }
1142: if (disc) {
1143: PetscCall(PetscObjectGetClassId(disc, &id));
1144: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1145: }
1146: if (!disc || isContainer) {
1147: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1148: PetscFE fe;
1149: DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1150: PetscInt dim, cStart, cEnd, cellHeight;
1152: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1153: PetscCall(DMGetDimension(dm, &dim));
1154: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1155: if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct));
1156: PetscCallMPI(MPIU_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm));
1157: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe));
1158: PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view"));
1159: disc = (PetscObject)fe;
1160: } else PetscCall(PetscObjectReference(disc));
1161: PetscCall(PetscObjectGetClassId(disc, &id));
1162: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents));
1163: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1164: PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b));
1165: PetscCall(DMFieldSetType(b, DMFIELDDS));
1166: dsfield = (DMField_DS *)b->data;
1167: dsfield->fieldNum = fieldNum;
1168: PetscCall(DMGetDimension(dm, &dsfield->height));
1169: dsfield->height++;
1170: PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc));
1171: dsfield->disc[0] = disc;
1172: PetscCall(PetscObjectReference((PetscObject)vec));
1173: dsfield->vec = vec;
1174: if (dmDG) {
1175: dsfield->dmDG = dmDG;
1176: PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG));
1177: dsfield->discDG[0] = discDG;
1178: PetscCall(PetscObjectReference((PetscObject)vecDG));
1179: dsfield->vecDG = vecDG;
1180: }
1181: *field = b;
1182: PetscFunctionReturn(PETSC_SUCCESS);
1183: }
1185: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1186: {
1187: PetscFunctionBegin;
1188: PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field));
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }