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 isascii;
41: PetscFunctionBegin;
42: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
43: disc = dsfield->disc[0];
44: if (isascii) {
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, feDim, i, K = H ? 2 : (D ? 1 : (B ? 0 : -1));
152: PetscSection section;
153: const PetscReal *qpoints;
154: PetscBool isStride;
155: const PetscInt *points = NULL;
156: PetscInt sfirst = -1, stride = -1;
157: PetscFE fe;
158: PetscTabulation T;
160: PetscFunctionBeginHot;
161: dm = field->dm;
162: nc = field->numComponents;
163: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL));
164: PetscCall(DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc));
165: PetscCall(DMGetDimension(dm, &meshDim));
166: PetscCall(DMGetLocalSection(dm, §ion));
167: PetscCall(PetscSectionGetField(section, dsfield->fieldNum, §ion));
168: PetscCall(PetscObjectGetClassId(disc, &classid));
169: /* TODO: batch */
170: PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
171: PetscCall(ISGetLocalSize(pointIS, &numCells));
172: if (isStride) PetscCall(ISStrideGetInfo(pointIS, &sfirst, &stride));
173: else PetscCall(ISGetIndices(pointIS, &points));
174: PetscCheck(classid == PETSCFE_CLASSID, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented");
175: fe = (PetscFE)disc;
176: PetscCall(PetscFEGetDimension(fe, &feDim));
177: PetscCall(PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T));
178: for (i = 0; i < numCells; i++) {
179: PetscInt c = isStride ? (sfirst + i * stride) : points[i];
180: PetscInt closureSize;
181: const PetscScalar *array;
182: PetscScalar *elem = NULL;
183: PetscBool isDG;
185: PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
186: if (B) {
187: /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
188: if (type == PETSC_SCALAR) {
189: PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i];
191: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
192: } else {
193: PetscReal *cB = &((PetscReal *)B)[nc * nq * i];
195: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
196: }
197: }
198: if (D) {
199: if (type == PETSC_SCALAR) {
200: PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i];
202: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
203: } else {
204: PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i];
206: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
207: }
208: }
209: if (H) {
210: if (type == PETSC_SCALAR) {
211: PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i];
213: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
214: } else {
215: PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i];
217: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
218: }
219: }
220: PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
221: }
222: PetscCall(PetscTabulationDestroy(&T));
223: if (!isStride) PetscCall(ISRestoreIndices(pointIS, &points));
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
228: {
229: DMField_DS *dsfield = (DMField_DS *)field->data;
230: PetscSF cellSF = NULL;
231: const PetscSFNode *cells;
232: PetscInt c, nFound, numCells, feDim, nc;
233: const PetscInt *cellDegrees;
234: const PetscScalar *pointsArray;
235: PetscScalar *cellPoints;
236: PetscInt gatherSize, gatherMax;
237: PetscInt dim, dimR, offset;
238: MPI_Datatype pointType;
239: PetscObject cellDisc;
240: PetscFE cellFE;
241: PetscClassId discID;
242: PetscReal *coordsReal, *coordsRef;
243: PetscSection section;
244: PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
245: PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
246: PetscReal *v, *J, *invJ, *detJ;
247: PetscMPIInt idim;
249: PetscFunctionBegin;
250: nc = field->numComponents;
251: PetscCall(DMGetLocalSection(field->dm, §ion));
252: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
253: PetscCall(PetscObjectGetClassId(cellDisc, &discID));
254: PetscCheck(discID == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization type not supported");
255: cellFE = (PetscFE)cellDisc;
256: PetscCall(PetscFEGetDimension(cellFE, &feDim));
257: PetscCall(DMGetCoordinateDim(field->dm, &dim));
258: PetscCall(DMGetDimension(field->dm, &dimR));
259: PetscCall(DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF));
260: PetscCall(PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells));
261: 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);
262: PetscCall(PetscSFComputeDegreeBegin(cellSF, &cellDegrees));
263: PetscCall(PetscSFComputeDegreeEnd(cellSF, &cellDegrees));
264: for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
265: gatherMax = PetscMax(gatherMax, cellDegrees[c]);
266: gatherSize += cellDegrees[c];
267: }
268: PetscCall(PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef));
269: PetscCall(PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ));
270: 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));
271: 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));
273: PetscCall(PetscMPIIntCast(dim, &idim));
274: PetscCallMPI(MPI_Type_contiguous(idim, MPIU_SCALAR, &pointType));
275: PetscCallMPI(MPI_Type_commit(&pointType));
276: PetscCall(VecGetArrayRead(points, &pointsArray));
277: PetscCall(PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints));
278: PetscCall(PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints));
279: PetscCall(VecRestoreArrayRead(points, &pointsArray));
280: for (c = 0, offset = 0; c < numCells; c++) {
281: PetscInt nq = cellDegrees[c], p;
283: if (nq) {
284: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
285: PetscTabulation T;
286: PetscQuadrature quad;
287: const PetscScalar *array;
288: PetscScalar *elem = NULL;
289: PetscReal *quadPoints;
290: PetscBool isDG;
291: PetscInt closureSize, d, e, f, g;
293: for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
294: PetscCall(DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef));
295: PetscCall(PetscFECreateTabulation(cellFE, 1, nq, coordsRef, K, &T));
296: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
297: PetscCall(PetscMalloc1(dimR * nq, &quadPoints));
298: for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
299: PetscCall(PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL));
300: PetscCall(DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ));
301: PetscCall(PetscQuadratureDestroy(&quad));
302: PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
303: if (B) {
304: if (datatype == PETSC_SCALAR) {
305: PetscScalar *cB = &cellBs[nc * offset];
307: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
308: } else {
309: PetscReal *cB = &cellBr[nc * offset];
311: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
312: }
313: }
314: if (D) {
315: if (datatype == PETSC_SCALAR) {
316: PetscScalar *cD = &cellDs[nc * dim * offset];
318: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
319: for (p = 0; p < nq; p++) {
320: for (g = 0; g < nc; g++) {
321: PetscScalar vs[3];
323: for (d = 0; d < dimR; d++) {
324: vs[d] = 0.;
325: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
326: }
327: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = vs[d];
328: }
329: }
330: } else {
331: PetscReal *cD = &cellDr[nc * dim * offset];
333: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
334: for (p = 0; p < nq; p++) {
335: for (g = 0; g < nc; g++) {
336: for (d = 0; d < dimR; d++) {
337: v[d] = 0.;
338: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
339: }
340: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = v[d];
341: }
342: }
343: }
344: }
345: if (H) {
346: if (datatype == PETSC_SCALAR) {
347: PetscScalar *cH = &cellHs[nc * dim * dim * offset];
349: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
350: for (p = 0; p < nq; p++) {
351: for (g = 0; g < nc * dimR; g++) {
352: PetscScalar vs[3];
354: for (d = 0; d < dimR; d++) {
355: vs[d] = 0.;
356: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
357: }
358: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = vs[d];
359: }
360: for (g = 0; g < nc; g++) {
361: for (f = 0; f < dimR; f++) {
362: PetscScalar vs[3];
364: for (d = 0; d < dimR; d++) {
365: vs[d] = 0.;
366: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
367: }
368: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
369: }
370: }
371: }
372: } else {
373: PetscReal *cH = &cellHr[nc * dim * dim * offset];
375: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
376: for (p = 0; p < nq; p++) {
377: for (g = 0; g < nc * dimR; g++) {
378: for (d = 0; d < dimR; d++) {
379: v[d] = 0.;
380: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
381: }
382: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = v[d];
383: }
384: for (g = 0; g < nc; g++) {
385: for (f = 0; f < dimR; f++) {
386: for (d = 0; d < dimR; d++) {
387: v[d] = 0.;
388: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
389: }
390: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
391: }
392: }
393: }
394: }
395: }
396: PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
397: PetscCall(PetscTabulationDestroy(&T));
398: }
399: offset += nq;
400: }
401: {
402: MPI_Datatype origtype;
404: if (datatype == PETSC_SCALAR) {
405: origtype = MPIU_SCALAR;
406: } else {
407: origtype = MPIU_REAL;
408: }
409: if (B) {
410: MPI_Datatype Btype;
412: PetscCall(PetscMPIIntCast(nc, &idim));
413: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Btype));
414: PetscCallMPI(MPI_Type_commit(&Btype));
415: PetscCall(PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
416: PetscCall(PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
417: PetscCallMPI(MPI_Type_free(&Btype));
418: }
419: if (D) {
420: MPI_Datatype Dtype;
422: PetscCall(PetscMPIIntCast(nc * dim, &idim));
423: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Dtype));
424: PetscCallMPI(MPI_Type_commit(&Dtype));
425: PetscCall(PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
426: PetscCall(PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
427: PetscCallMPI(MPI_Type_free(&Dtype));
428: }
429: if (H) {
430: MPI_Datatype Htype;
432: PetscCall(PetscMPIIntCast(nc * dim * dim, &idim));
433: PetscCallMPI(MPI_Type_contiguous(idim, origtype, &Htype));
434: PetscCallMPI(MPI_Type_commit(&Htype));
435: PetscCall(PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
436: PetscCall(PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
437: PetscCallMPI(MPI_Type_free(&Htype));
438: }
439: }
440: PetscCall(PetscFree4(v, J, invJ, detJ));
441: PetscCall(PetscFree3(cellBr, cellDr, cellHr));
442: PetscCall(PetscFree3(cellBs, cellDs, cellHs));
443: PetscCall(PetscFree3(cellPoints, coordsReal, coordsRef));
444: PetscCallMPI(MPI_Type_free(&pointType));
445: PetscCall(PetscSFDestroy(&cellSF));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
450: {
451: DMField_DS *dsfield = (DMField_DS *)field->data;
452: PetscInt h, imin;
453: PetscInt dim;
454: PetscClassId id;
455: PetscQuadrature quad = NULL;
456: PetscInt maxDegree;
457: PetscFEGeom *geom;
458: PetscInt Nq, Nc, dimC, qNc, N;
459: PetscInt numPoints;
460: void *qB = NULL, *qD = NULL, *qH = NULL;
461: const PetscReal *weights;
462: MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
463: PetscObject disc;
464: DMField coordField;
466: PetscFunctionBegin;
467: Nc = field->numComponents;
468: PetscCall(DMGetCoordinateDim(field->dm, &dimC));
469: PetscCall(DMGetDimension(field->dm, &dim));
470: PetscCall(ISGetLocalSize(pointIS, &numPoints));
471: PetscCall(ISGetMinMax(pointIS, &imin, NULL));
472: for (h = 0; h < dsfield->height; h++) {
473: PetscInt hEnd;
475: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
476: if (imin < hEnd) break;
477: }
478: dim -= h;
479: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
480: PetscCall(PetscObjectGetClassId(disc, &id));
481: PetscCheck(id == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported");
482: PetscCall(DMGetCoordinateField(field->dm, &coordField));
483: PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
484: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
485: if (!quad) PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, &quad));
486: PetscCheck(quad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages");
487: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FEGEOM_BASIC, &geom));
488: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights));
489: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components");
490: N = numPoints * Nq * Nc;
491: if (B) PetscCall(DMGetWorkArray(field->dm, N, mpitype, &qB));
492: if (D) PetscCall(DMGetWorkArray(field->dm, N * dimC, mpitype, &qD));
493: if (H) PetscCall(DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
494: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH));
495: if (B) {
496: PetscInt i, j, k;
498: if (type == PETSC_SCALAR) {
499: PetscScalar *sB = (PetscScalar *)B;
500: PetscScalar *sqB = (PetscScalar *)qB;
502: for (i = 0; i < numPoints; i++) {
503: PetscReal vol = 0.;
505: for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.;
506: for (k = 0; k < Nq; k++) {
507: vol += geom->detJ[i * Nq + k] * weights[k];
508: for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j];
509: }
510: for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
511: }
512: } else {
513: PetscReal *rB = (PetscReal *)B;
514: PetscReal *rqB = (PetscReal *)qB;
516: for (i = 0; i < numPoints; i++) {
517: PetscReal vol = 0.;
519: for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.;
520: for (k = 0; k < Nq; k++) {
521: vol += geom->detJ[i * Nq + k] * weights[k];
522: for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j];
523: }
524: for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
525: }
526: }
527: }
528: if (D) {
529: PetscInt i, j, k, l, m;
531: if (type == PETSC_SCALAR) {
532: PetscScalar *sD = (PetscScalar *)D;
533: PetscScalar *sqD = (PetscScalar *)qD;
535: for (i = 0; i < numPoints; i++) {
536: PetscReal vol = 0.;
538: for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
539: for (k = 0; k < Nq; k++) {
540: vol += geom->detJ[i * Nq + k] * weights[k];
541: for (j = 0; j < Nc; j++) {
542: PetscScalar pD[3] = {0., 0., 0.};
544: for (l = 0; l < dimC; l++) {
545: 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];
546: }
547: for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
548: }
549: }
550: for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
551: }
552: } else {
553: PetscReal *rD = (PetscReal *)D;
554: PetscReal *rqD = (PetscReal *)qD;
556: for (i = 0; i < numPoints; i++) {
557: PetscReal vol = 0.;
559: for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
560: for (k = 0; k < Nq; k++) {
561: vol += geom->detJ[i * Nq + k] * weights[k];
562: for (j = 0; j < Nc; j++) {
563: PetscReal pD[3] = {0., 0., 0.};
565: for (l = 0; l < dimC; l++) {
566: 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];
567: }
568: for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
569: }
570: }
571: for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
572: }
573: }
574: }
575: if (H) {
576: PetscInt i, j, k, l, m, q, r;
578: if (type == PETSC_SCALAR) {
579: PetscScalar *sH = (PetscScalar *)H;
580: PetscScalar *sqH = (PetscScalar *)qH;
582: for (i = 0; i < numPoints; i++) {
583: PetscReal vol = 0.;
585: for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.;
586: for (k = 0; k < Nq; k++) {
587: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
589: vol += geom->detJ[i * Nq + k] * weights[k];
590: for (j = 0; j < Nc; j++) {
591: PetscScalar pH[3][3] = {
592: {0., 0., 0.},
593: {0., 0., 0.},
594: {0., 0., 0.}
595: };
596: const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];
598: for (l = 0; l < dimC; l++) {
599: for (m = 0; m < dimC; m++) {
600: for (q = 0; q < dim; q++) {
601: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
602: }
603: }
604: }
605: for (l = 0; l < dimC; l++) {
606: 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];
607: }
608: }
609: }
610: for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
611: }
612: } else {
613: PetscReal *rH = (PetscReal *)H;
614: PetscReal *rqH = (PetscReal *)qH;
616: for (i = 0; i < numPoints; i++) {
617: PetscReal vol = 0.;
619: for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.;
620: for (k = 0; k < Nq; k++) {
621: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
623: vol += geom->detJ[i * Nq + k] * weights[k];
624: for (j = 0; j < Nc; j++) {
625: PetscReal pH[3][3] = {
626: {0., 0., 0.},
627: {0., 0., 0.},
628: {0., 0., 0.}
629: };
630: const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];
632: for (l = 0; l < dimC; l++) {
633: for (m = 0; m < dimC; m++) {
634: for (q = 0; q < dim; q++) {
635: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
636: }
637: }
638: }
639: for (l = 0; l < dimC; l++) {
640: 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];
641: }
642: }
643: }
644: for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
645: }
646: }
647: }
648: if (B) PetscCall(DMRestoreWorkArray(field->dm, N, mpitype, &qB));
649: if (D) PetscCall(DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD));
650: if (H) PetscCall(DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
651: PetscCall(PetscFEGeomDestroy(&geom));
652: PetscCall(PetscQuadratureDestroy(&quad));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
657: {
658: DMField_DS *dsfield;
659: PetscObject disc;
660: PetscInt h, imin, imax;
661: PetscClassId id;
663: PetscFunctionBegin;
664: dsfield = (DMField_DS *)field->data;
665: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
666: if (imin >= imax) {
667: h = 0;
668: } else {
669: for (h = 0; h < dsfield->height; h++) {
670: PetscInt hEnd;
672: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
673: if (imin < hEnd) break;
674: }
675: }
676: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
677: PetscCall(PetscObjectGetClassId(disc, &id));
678: if (id == PETSCFE_CLASSID) {
679: PetscFE fe = (PetscFE)disc;
680: PetscSpace sp;
682: PetscCall(PetscFEGetBasisSpace(fe, &sp));
683: PetscCall(PetscSpaceGetDegree(sp, minDegree, maxDegree));
684: }
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
689: {
690: DM dm = field->dm;
691: const PetscInt *points;
692: DMPolytopeType ct;
693: PetscInt dim, n;
694: PetscBool isplex;
696: PetscFunctionBegin;
697: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
698: PetscCall(ISGetLocalSize(pointIS, &n));
699: if (isplex && n) {
700: PetscCall(DMGetDimension(dm, &dim));
701: PetscCall(ISGetIndices(pointIS, &points));
702: PetscCall(DMPlexGetCellType(dm, points[0], &ct));
703: switch (ct) {
704: case DM_POLYTOPE_TRIANGLE:
705: case DM_POLYTOPE_TETRAHEDRON:
706: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad));
707: break;
708: default:
709: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
710: }
711: PetscCall(ISRestoreIndices(pointIS, &points));
712: } else PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, quad));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
717: {
718: PetscInt h, dim, imax, imin, cellHeight;
719: DM dm;
720: DMField_DS *dsfield;
721: PetscObject disc;
722: PetscFE fe;
723: PetscClassId id;
725: PetscFunctionBegin;
726: dm = field->dm;
727: dsfield = (DMField_DS *)field->data;
728: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
729: PetscCall(DMGetDimension(dm, &dim));
730: for (h = 0; h <= dim; h++) {
731: PetscInt hStart, hEnd;
733: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
734: if (imax >= hStart && imin < hEnd) break;
735: }
736: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
737: h -= cellHeight;
738: *quad = NULL;
739: if (h < dsfield->height) {
740: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
741: PetscCall(PetscObjectGetClassId(disc, &id));
742: if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
743: fe = (PetscFE)disc;
744: PetscCall(PetscFEGetQuadrature(fe, quad));
745: PetscCall(PetscObjectReference((PetscObject)*quad));
746: }
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: static PetscErrorCode DMFieldCreateDefaultFaceQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
751: {
752: PetscInt h, dim, imax, imin, cellHeight;
753: DM dm;
754: DMField_DS *dsfield;
755: PetscObject disc;
756: PetscFE fe;
757: PetscClassId id;
759: PetscFunctionBegin;
760: dm = field->dm;
761: dsfield = (DMField_DS *)field->data;
762: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
763: PetscCall(DMGetDimension(dm, &dim));
764: for (h = 0; h <= dim; h++) {
765: PetscInt hStart, hEnd;
767: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
768: if (imax >= hStart && imin < hEnd) break;
769: }
770: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
771: h -= cellHeight;
772: *quad = NULL;
773: if (h < dsfield->height) {
774: PetscQuadrature fq;
776: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
777: PetscCall(PetscObjectGetClassId(disc, &id));
778: if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
779: fe = (PetscFE)disc;
780: PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
781: PetscCall(PetscFEExpandFaceQuadrature(fe, fq, quad));
782: }
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
787: {
788: const PetscInt *points;
789: PetscInt p, dim, dE, numFaces, Nq;
790: PetscInt maxDegree;
791: DMLabel depthLabel;
792: IS cellIS;
793: DM dm = field->dm;
795: PetscFunctionBegin;
796: dim = geom->dim;
797: dE = geom->dimEmbed;
798: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
799: PetscCall(DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS));
800: PetscCall(DMFieldGetDegree(field, cellIS, NULL, &maxDegree));
801: PetscCall(ISGetIndices(pointIS, &points));
802: numFaces = geom->numCells;
803: Nq = geom->numPoints;
804: /* First, set local faces and flip normals so that they are outward for the first supporting cell */
805: for (p = 0; p < numFaces; p++) {
806: PetscInt point = points[p];
807: PetscInt suppSize, s, coneSize, c, numChildren;
808: const PetscInt *supp;
810: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
811: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
812: PetscCall(DMPlexGetSupportSize(dm, point, &suppSize));
813: PetscCheck(suppSize <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, suppSize);
814: if (!suppSize) continue;
815: PetscCall(DMPlexGetSupport(dm, point, &supp));
816: for (s = 0; s < suppSize; ++s) {
817: const PetscInt *cone, *ornt;
819: PetscCall(DMPlexGetConeSize(dm, supp[s], &coneSize));
820: PetscCall(DMPlexGetOrientedCone(dm, supp[s], &cone, &ornt));
821: for (c = 0; c < coneSize; ++c)
822: if (cone[c] == point) break;
823: geom->face[p][s * 2 + 0] = c;
824: geom->face[p][s * 2 + 1] = ornt[c];
825: PetscCall(DMPlexRestoreOrientedCone(dm, supp[s], &cone, &ornt));
826: 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]);
827: }
828: if (geom->face[p][1] < 0) {
829: PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;
831: for (q = 0; q < Np; ++q)
832: for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
833: }
834: }
835: if (maxDegree <= 1) {
836: PetscQuadrature cellQuad = NULL;
837: PetscInt numCells, offset, *cells;
838: PetscFEGeom *cellGeom;
839: IS suppIS;
841: if (quad) {
842: DM dm;
843: PetscReal *points, *weights;
844: PetscInt tdim, Nc, Np;
846: PetscCall(DMFieldGetDM(field, &dm));
847: PetscCall(DMGetDimension(dm, &tdim));
848: if (tdim > dim) {
849: // Make a compatible cell quadrature (points don't matter since its affine)
850: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
851: PetscCall(PetscQuadratureGetData(quad, NULL, &Nc, &Np, NULL, NULL));
852: PetscCall(PetscCalloc1((dim + 1) * Np, &points));
853: PetscCall(PetscCalloc1(Nc * Np, &weights));
854: PetscCall(PetscQuadratureSetData(cellQuad, dim + 1, Nc, Np, points, weights));
855: } else {
856: // TODO J will be wrong here, but other things need to be fixed
857: // This path comes from calling DMProjectBdFieldLabelLocal() in Plex ex5
858: PetscCall(PetscObjectReference((PetscObject)quad));
859: cellQuad = quad;
860: }
861: }
862: for (p = 0, numCells = 0; p < numFaces; p++) {
863: PetscInt point = points[p];
864: PetscInt numSupp, numChildren;
866: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
867: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
868: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
869: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
870: numCells += numSupp;
871: }
872: PetscCall(PetscMalloc1(numCells, &cells));
873: for (p = 0, offset = 0; p < numFaces; p++) {
874: PetscInt point = points[p];
875: PetscInt numSupp, s;
876: const PetscInt *supp;
878: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
879: PetscCall(DMPlexGetSupport(dm, point, &supp));
880: for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
881: }
882: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
883: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FEGEOM_BASIC, &cellGeom));
884: for (p = 0, offset = 0; p < numFaces; p++) {
885: PetscInt point = points[p];
886: PetscInt numSupp, s, q;
887: const PetscInt *supp;
889: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
890: PetscCall(DMPlexGetSupport(dm, point, &supp));
891: for (s = 0; s < numSupp; s++, offset++) {
892: for (q = 0; q < Nq * dE * dE; q++) {
893: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
894: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
895: }
896: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
897: }
898: }
899: PetscCall(PetscFEGeomDestroy(&cellGeom));
900: PetscCall(PetscQuadratureDestroy(&cellQuad));
901: PetscCall(ISDestroy(&suppIS));
902: PetscCall(PetscFree(cells));
903: } else {
904: DMField_DS *dsfield = (DMField_DS *)field->data;
905: PetscObject faceDisc, cellDisc;
906: PetscClassId faceId, cellId;
907: PetscDualSpace dsp;
908: DM K;
909: DMPolytopeType ct;
910: PetscInt (*co)[2][3];
911: PetscInt coneSize;
912: PetscInt **counts;
913: PetscInt f, i, o, q, s;
914: PetscBool found = PETSC_FALSE;
915: const PetscInt *coneK;
916: PetscInt eStart, minOrient, maxOrient, numOrient;
917: PetscInt *orients;
918: PetscReal **orientPoints;
919: PetscReal *cellPoints;
920: PetscReal *dummyWeights;
921: PetscQuadrature cellQuad = NULL;
923: PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc));
924: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
925: PetscCall(PetscObjectGetClassId(faceDisc, &faceId));
926: PetscCall(PetscObjectGetClassId(cellDisc, &cellId));
927: PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported");
928: PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp));
929: PetscCall(PetscDualSpaceGetDM(dsp, &K));
930: PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL));
931: PetscCall(DMPlexGetCellType(K, eStart, &ct));
932: PetscCall(DMPlexGetConeSize(K, 0, &coneSize));
933: PetscCall(DMPlexGetCone(K, 0, &coneK));
934: PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts));
935: PetscCall(PetscMalloc1(dE * Nq, &cellPoints));
936: PetscCall(PetscMalloc1(Nq, &dummyWeights));
937: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
938: PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights));
939: minOrient = PETSC_INT_MAX;
940: maxOrient = PETSC_INT_MIN;
941: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
942: PetscInt point = points[p];
943: PetscInt numSupp, numChildren;
944: const PetscInt *supp;
946: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
947: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
948: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
949: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
950: PetscCall(DMPlexGetSupport(dm, point, &supp));
951: for (s = 0; s < numSupp; s++) {
952: PetscInt cell = supp[s];
953: PetscInt numCone;
954: const PetscInt *cone, *orient;
956: PetscCall(DMPlexGetConeSize(dm, cell, &numCone));
957: // When we extract submeshes, we hang cells from the side that are not fully realized. We ignore these
958: if (numCone == 1) {
959: co[p][s][0] = -1;
960: co[p][s][1] = -1;
961: co[p][s][2] = -1;
962: continue;
963: }
964: PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
965: PetscCall(DMPlexGetCone(dm, cell, &cone));
966: PetscCall(DMPlexGetConeOrientation(dm, cell, &orient));
967: for (f = 0; f < coneSize; f++) {
968: if (cone[f] == point) break;
969: }
970: co[p][s][0] = f;
971: co[p][s][1] = orient[f];
972: co[p][s][2] = cell;
973: minOrient = PetscMin(minOrient, orient[f]);
974: maxOrient = PetscMax(maxOrient, orient[f]);
975: found = PETSC_TRUE;
976: }
977: for (; s < 2; s++) {
978: co[p][s][0] = -1;
979: co[p][s][1] = -1;
980: co[p][s][2] = -1;
981: }
982: }
983: numOrient = found ? maxOrient + 1 - minOrient : 0;
984: PetscCall(DMPlexGetCone(K, 0, &coneK));
985: /* count all (face,orientation) doubles that appear */
986: PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints));
987: for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f]));
988: for (p = 0; p < numFaces; p++) {
989: for (s = 0; s < 2; s++) {
990: if (co[p][s][0] >= 0) {
991: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
992: orients[co[p][s][1] - minOrient]++;
993: }
994: }
995: }
996: for (o = 0; o < numOrient; o++) {
997: if (orients[o]) {
998: PetscInt orient = o + minOrient;
999: PetscInt q;
1001: PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o]));
1002: /* rotate the quadrature points appropriately */
1003: switch (ct) {
1004: case DM_POLYTOPE_POINT:
1005: break;
1006: case DM_POLYTOPE_SEGMENT:
1007: if (orient == -2 || orient == 1) {
1008: for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
1009: } else {
1010: for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
1011: }
1012: break;
1013: case DM_POLYTOPE_TRIANGLE:
1014: for (q = 0; q < Nq; q++) {
1015: PetscReal lambda[3];
1016: PetscReal lambdao[3];
1018: /* convert to barycentric */
1019: lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
1020: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
1021: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
1022: if (orient >= 0) {
1023: for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
1024: } else {
1025: for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
1026: }
1027: /* convert to coordinates */
1028: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
1029: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
1030: }
1031: break;
1032: case DM_POLYTOPE_QUADRILATERAL:
1033: for (q = 0; q < Nq; q++) {
1034: PetscReal xi[2], xio[2];
1035: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
1037: xi[0] = geom->xi[2 * q];
1038: xi[1] = geom->xi[2 * q + 1];
1039: switch (oabs) {
1040: case 1:
1041: xio[0] = xi[1];
1042: xio[1] = -xi[0];
1043: break;
1044: case 2:
1045: xio[0] = -xi[0];
1046: xio[1] = -xi[1];
1047: break;
1048: case 3:
1049: xio[0] = -xi[1];
1050: xio[1] = xi[0];
1051: break;
1052: case 0:
1053: default:
1054: xio[0] = xi[0];
1055: xio[1] = xi[1];
1056: break;
1057: }
1058: if (orient < 0) xio[0] = -xio[0];
1059: orientPoints[o][2 * q + 0] = xio[0];
1060: orientPoints[o][2 * q + 1] = xio[1];
1061: }
1062: break;
1063: default:
1064: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
1065: }
1066: }
1067: }
1068: for (f = 0; f < coneSize; f++) {
1069: PetscInt face = coneK[f];
1070: PetscReal v0[3];
1071: PetscReal J[9], detJ;
1072: PetscInt numCells, offset;
1073: PetscInt *cells;
1074: IS suppIS;
1076: PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ));
1077: for (o = 0; o <= numOrient; o++) {
1078: PetscFEGeom *cellGeom;
1080: if (!counts[f][o]) continue;
1081: /* If this (face,orientation) double appears,
1082: * convert the face quadrature points into volume quadrature points */
1083: for (q = 0; q < Nq; q++) {
1084: PetscReal xi0[3] = {-1., -1., -1.};
1086: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1087: }
1088: for (p = 0, numCells = 0; p < numFaces; p++) {
1089: for (s = 0; s < 2; s++) {
1090: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1091: }
1092: }
1093: PetscCall(PetscMalloc1(numCells, &cells));
1094: for (p = 0, offset = 0; p < numFaces; p++) {
1095: for (s = 0; s < 2; s++) {
1096: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1097: }
1098: }
1099: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
1100: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FEGEOM_BASIC, &cellGeom));
1101: for (p = 0, offset = 0; p < numFaces; p++) {
1102: for (s = 0; s < 2; s++) {
1103: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1104: for (q = 0; q < Nq * dE * dE; q++) {
1105: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1106: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1107: }
1108: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1109: offset++;
1110: }
1111: }
1112: }
1113: PetscCall(PetscFEGeomDestroy(&cellGeom));
1114: PetscCall(ISDestroy(&suppIS));
1115: PetscCall(PetscFree(cells));
1116: }
1117: }
1118: for (o = 0; o < numOrient; o++) {
1119: if (orients[o]) PetscCall(PetscFree(orientPoints[o]));
1120: }
1121: PetscCall(PetscFree2(orients, orientPoints));
1122: PetscCall(PetscQuadratureDestroy(&cellQuad));
1123: for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f]));
1124: PetscCall(PetscFree2(co, counts));
1125: }
1126: PetscCall(ISRestoreIndices(pointIS, &points));
1127: PetscCall(ISDestroy(&cellIS));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1132: {
1133: PetscFunctionBegin;
1134: field->ops->destroy = DMFieldDestroy_DS;
1135: field->ops->evaluate = DMFieldEvaluate_DS;
1136: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1137: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1138: field->ops->getDegree = DMFieldGetDegree_DS;
1139: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1140: field->ops->createDefaultFaceQuadrature = DMFieldCreateDefaultFaceQuadrature_DS;
1141: field->ops->view = DMFieldView_DS;
1142: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1143: PetscFunctionReturn(PETSC_SUCCESS);
1144: }
1146: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1147: {
1148: DMField_DS *dsfield;
1150: PetscFunctionBegin;
1151: PetscCall(PetscNew(&dsfield));
1152: field->data = dsfield;
1153: PetscCall(DMFieldInitialize_DS(field));
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1158: {
1159: DMField b;
1160: DMField_DS *dsfield;
1161: PetscObject disc = NULL, discDG = NULL;
1162: PetscSection section;
1163: PetscBool isContainer = PETSC_FALSE;
1164: PetscClassId id = -1;
1165: PetscInt numComponents = -1, dsNumFields;
1167: PetscFunctionBegin;
1168: PetscCall(DMGetLocalSection(dm, §ion));
1169: PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents));
1170: PetscCall(DMGetNumFields(dm, &dsNumFields));
1171: if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc));
1172: if (dsNumFields && dmDG) {
1173: PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG));
1174: PetscCall(PetscObjectReference(discDG));
1175: }
1176: if (disc) {
1177: PetscCall(PetscObjectGetClassId(disc, &id));
1178: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1179: }
1180: if (!disc || isContainer) {
1181: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1182: PetscFE fe;
1183: DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1184: PetscInt dim, cStart, cEnd, cellHeight;
1186: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1187: PetscCall(DMGetDimension(dm, &dim));
1188: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1189: if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct));
1190: PetscCallMPI(MPIU_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm));
1191: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe));
1192: PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view"));
1193: disc = (PetscObject)fe;
1194: } else PetscCall(PetscObjectReference(disc));
1195: PetscCall(PetscObjectGetClassId(disc, &id));
1196: PetscCheck(id == PETSCFE_CLASSID, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1197: PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents));
1198: PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b));
1199: PetscCall(DMFieldSetType(b, DMFIELDDS));
1200: dsfield = (DMField_DS *)b->data;
1201: dsfield->fieldNum = fieldNum;
1202: PetscCall(DMGetDimension(dm, &dsfield->height));
1203: dsfield->height++;
1204: PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc));
1205: dsfield->disc[0] = disc;
1206: PetscCall(PetscObjectReference((PetscObject)vec));
1207: dsfield->vec = vec;
1208: if (dmDG) {
1209: dsfield->dmDG = dmDG;
1210: PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG));
1211: dsfield->discDG[0] = discDG;
1212: PetscCall(PetscObjectReference((PetscObject)vecDG));
1213: dsfield->vecDG = vecDG;
1214: }
1215: *field = b;
1216: PetscFunctionReturn(PETSC_SUCCESS);
1217: }
1219: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1220: {
1221: PetscFunctionBegin;
1222: PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field));
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }