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: if (!discList[height]) {
61: PetscClassId id;
63: PetscCall(PetscObjectGetClassId(discList[0], &id));
64: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height]));
65: }
66: *disc = discList[height];
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: /* y[m,c] = A[m,n,c] . b[n] */
71: #define DMFieldDSdot(y, A, b, m, n, c, cast) \
72: do { \
73: PetscInt _i, _j, _k; \
74: for (_i = 0; _i < (m); _i++) { \
75: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] = 0.; \
76: for (_j = 0; _j < (n); _j++) { \
77: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \
78: } \
79: } \
80: } while (0)
82: /*
83: 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().
84: */
85: static PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
86: {
87: DMField_DS *dsfield = (DMField_DS *)field->data;
88: DM fdm = dsfield->dmDG;
89: PetscSection s = NULL;
90: const PetscScalar *cvalues;
91: PetscInt pStart, pEnd;
93: PetscFunctionBeginHot;
94: *isDG = PETSC_FALSE;
95: *Nc = 0;
96: *array = NULL;
97: *values = NULL;
98: /* Check for cellwise section */
99: if (fdm) PetscCall(DMGetLocalSection(fdm, &s));
100: if (!s) goto cg;
101: /* Check that the cell exists in the cellwise section */
102: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
103: if (cell < pStart || cell >= pEnd) goto cg;
104: /* Check for cellwise coordinates for this cell */
105: PetscCall(PetscSectionGetDof(s, cell, Nc));
106: if (!*Nc) goto cg;
107: /* Check for cellwise coordinates */
108: if (!dsfield->vecDG) goto cg;
109: /* Get cellwise coordinates */
110: PetscCall(VecGetArrayRead(dsfield->vecDG, array));
111: PetscCall(DMPlexPointLocalRead(fdm, cell, *array, &cvalues));
112: PetscCall(DMGetWorkArray(fdm, *Nc, MPIU_SCALAR, values));
113: PetscCall(PetscArraycpy(*values, cvalues, *Nc));
114: PetscCall(VecRestoreArrayRead(dsfield->vecDG, array));
115: *isDG = PETSC_TRUE;
116: PetscFunctionReturn(PETSC_SUCCESS);
117: cg:
118: /* Use continuous values */
119: PetscCall(DMFieldGetDM(field, &fdm));
120: PetscCall(DMGetLocalSection(fdm, &s));
121: PetscCall(PetscSectionGetField(s, dsfield->fieldNum, &s));
122: PetscCall(DMPlexVecGetClosure(fdm, s, dsfield->vec, cell, Nc, values));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode DMFieldRestoreClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
127: {
128: DMField_DS *dsfield = (DMField_DS *)field->data;
129: DM fdm;
130: PetscSection s;
132: PetscFunctionBeginHot;
133: if (*isDG) {
134: PetscCall(DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values));
135: } else {
136: PetscCall(DMFieldGetDM(field, &fdm));
137: PetscCall(DMGetLocalSection(fdm, &s));
138: PetscCall(DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, (PetscScalar **)values));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
144: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
145: {
146: DMField_DS *dsfield = (DMField_DS *)field->data;
147: DM dm;
148: PetscObject disc;
149: PetscClassId classid;
150: PetscInt nq, nc, dim, meshDim, numCells;
151: PetscSection section;
152: const PetscReal *qpoints;
153: PetscBool isStride;
154: const PetscInt *points = NULL;
155: PetscInt sfirst = -1, stride = -1;
157: PetscFunctionBeginHot;
158: dm = field->dm;
159: nc = field->numComponents;
160: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL));
161: PetscCall(DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc));
162: PetscCall(DMGetDimension(dm, &meshDim));
163: PetscCall(DMGetLocalSection(dm, §ion));
164: PetscCall(PetscSectionGetField(section, dsfield->fieldNum, §ion));
165: PetscCall(PetscObjectGetClassId(disc, &classid));
166: /* TODO: batch */
167: PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
168: PetscCall(ISGetLocalSize(pointIS, &numCells));
169: if (isStride) PetscCall(ISStrideGetInfo(pointIS, &sfirst, &stride));
170: else PetscCall(ISGetIndices(pointIS, &points));
171: if (classid == PETSCFE_CLASSID) {
172: PetscFE fe = (PetscFE)disc;
173: PetscInt feDim, i;
174: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
175: PetscTabulation T;
177: PetscCall(PetscFEGetDimension(fe, &feDim));
178: PetscCall(PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T));
179: for (i = 0; i < numCells; i++) {
180: PetscInt c = isStride ? (sfirst + i * stride) : points[i];
181: PetscInt closureSize;
182: const PetscScalar *array;
183: PetscScalar *elem = NULL;
184: PetscBool isDG;
186: PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
187: if (B) {
188: /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
189: if (type == PETSC_SCALAR) {
190: PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i];
192: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
193: } else {
194: PetscReal *cB = &((PetscReal *)B)[nc * nq * i];
196: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
197: }
198: }
199: if (D) {
200: if (type == PETSC_SCALAR) {
201: PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i];
203: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
204: } else {
205: PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i];
207: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
208: }
209: }
210: if (H) {
211: if (type == PETSC_SCALAR) {
212: PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i];
214: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
215: } else {
216: PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i];
218: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
219: }
220: }
221: PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
222: }
223: PetscCall(PetscTabulationDestroy(&T));
224: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented");
225: if (!isStride) PetscCall(ISRestoreIndices(pointIS, &points));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
230: {
231: DMField_DS *dsfield = (DMField_DS *)field->data;
232: PetscSF cellSF = NULL;
233: const PetscSFNode *cells;
234: PetscInt c, nFound, numCells, feDim, nc;
235: const PetscInt *cellDegrees;
236: const PetscScalar *pointsArray;
237: PetscScalar *cellPoints;
238: PetscInt gatherSize, gatherMax;
239: PetscInt dim, dimR, offset;
240: MPI_Datatype pointType;
241: PetscObject cellDisc;
242: PetscFE cellFE;
243: PetscClassId discID;
244: PetscReal *coordsReal, *coordsRef;
245: PetscSection section;
246: PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
247: PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
248: PetscReal *v, *J, *invJ, *detJ;
250: PetscFunctionBegin;
251: nc = field->numComponents;
252: PetscCall(DMGetLocalSection(field->dm, §ion));
253: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
254: PetscCall(PetscObjectGetClassId(cellDisc, &discID));
255: PetscCheck(discID == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization type not supported");
256: cellFE = (PetscFE)cellDisc;
257: PetscCall(PetscFEGetDimension(cellFE, &feDim));
258: PetscCall(DMGetCoordinateDim(field->dm, &dim));
259: PetscCall(DMGetDimension(field->dm, &dimR));
260: PetscCall(DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF));
261: PetscCall(PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells));
262: 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);
263: PetscCall(PetscSFComputeDegreeBegin(cellSF, &cellDegrees));
264: PetscCall(PetscSFComputeDegreeEnd(cellSF, &cellDegrees));
265: for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
266: gatherMax = PetscMax(gatherMax, cellDegrees[c]);
267: gatherSize += cellDegrees[c];
268: }
269: PetscCall(PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef));
270: PetscCall(PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ));
271: if (datatype == PETSC_SCALAR) PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs));
272: else PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr));
274: PetscCallMPI(MPI_Type_contiguous(dim, 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: PetscCallMPI(MPI_Type_contiguous(nc, origtype, &Btype));
413: PetscCallMPI(MPI_Type_commit(&Btype));
414: PetscCall(PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
415: PetscCall(PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
416: PetscCallMPI(MPI_Type_free(&Btype));
417: }
418: if (D) {
419: MPI_Datatype Dtype;
421: PetscCallMPI(MPI_Type_contiguous(nc * dim, origtype, &Dtype));
422: PetscCallMPI(MPI_Type_commit(&Dtype));
423: PetscCall(PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
424: PetscCall(PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
425: PetscCallMPI(MPI_Type_free(&Dtype));
426: }
427: if (H) {
428: MPI_Datatype Htype;
430: PetscCallMPI(MPI_Type_contiguous(nc * dim * dim, origtype, &Htype));
431: PetscCallMPI(MPI_Type_commit(&Htype));
432: PetscCall(PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
433: PetscCall(PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
434: PetscCallMPI(MPI_Type_free(&Htype));
435: }
436: }
437: PetscCall(PetscFree4(v, J, invJ, detJ));
438: PetscCall(PetscFree3(cellBr, cellDr, cellHr));
439: PetscCall(PetscFree3(cellBs, cellDs, cellHs));
440: PetscCall(PetscFree3(cellPoints, coordsReal, coordsRef));
441: PetscCallMPI(MPI_Type_free(&pointType));
442: PetscCall(PetscSFDestroy(&cellSF));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
447: {
448: DMField_DS *dsfield = (DMField_DS *)field->data;
449: PetscInt h, imin;
450: PetscInt dim;
451: PetscClassId id;
452: PetscQuadrature quad = NULL;
453: PetscInt maxDegree;
454: PetscFEGeom *geom;
455: PetscInt Nq, Nc, dimC, qNc, N;
456: PetscInt numPoints;
457: void *qB = NULL, *qD = NULL, *qH = NULL;
458: const PetscReal *weights;
459: MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
460: PetscObject disc;
461: DMField coordField;
463: PetscFunctionBegin;
464: Nc = field->numComponents;
465: PetscCall(DMGetCoordinateDim(field->dm, &dimC));
466: PetscCall(DMGetDimension(field->dm, &dim));
467: PetscCall(ISGetLocalSize(pointIS, &numPoints));
468: PetscCall(ISGetMinMax(pointIS, &imin, NULL));
469: for (h = 0; h < dsfield->height; h++) {
470: PetscInt hEnd;
472: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
473: if (imin < hEnd) break;
474: }
475: dim -= h;
476: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
477: PetscCall(PetscObjectGetClassId(disc, &id));
478: PetscCheck(id == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported");
479: PetscCall(DMGetCoordinateField(field->dm, &coordField));
480: PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
481: if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
482: if (!quad) PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, &quad));
483: PetscCheck(quad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages");
484: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom));
485: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights));
486: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components");
487: N = numPoints * Nq * Nc;
488: if (B) PetscCall(DMGetWorkArray(field->dm, N, mpitype, &qB));
489: if (D) PetscCall(DMGetWorkArray(field->dm, N * dimC, mpitype, &qD));
490: if (H) PetscCall(DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
491: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH));
492: if (B) {
493: PetscInt i, j, k;
495: if (type == PETSC_SCALAR) {
496: PetscScalar *sB = (PetscScalar *)B;
497: PetscScalar *sqB = (PetscScalar *)qB;
499: for (i = 0; i < numPoints; i++) {
500: PetscReal vol = 0.;
502: for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.;
503: for (k = 0; k < Nq; k++) {
504: vol += geom->detJ[i * Nq + k] * weights[k];
505: for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j];
506: }
507: for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
508: }
509: } else {
510: PetscReal *rB = (PetscReal *)B;
511: PetscReal *rqB = (PetscReal *)qB;
513: for (i = 0; i < numPoints; i++) {
514: PetscReal vol = 0.;
516: for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.;
517: for (k = 0; k < Nq; k++) {
518: vol += geom->detJ[i * Nq + k] * weights[k];
519: for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j];
520: }
521: for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
522: }
523: }
524: }
525: if (D) {
526: PetscInt i, j, k, l, m;
528: if (type == PETSC_SCALAR) {
529: PetscScalar *sD = (PetscScalar *)D;
530: PetscScalar *sqD = (PetscScalar *)qD;
532: for (i = 0; i < numPoints; i++) {
533: PetscReal vol = 0.;
535: for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
536: for (k = 0; k < Nq; k++) {
537: vol += geom->detJ[i * Nq + k] * weights[k];
538: for (j = 0; j < Nc; j++) {
539: PetscScalar pD[3] = {0., 0., 0.};
541: for (l = 0; l < dimC; l++) {
542: 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];
543: }
544: for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
545: }
546: }
547: for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
548: }
549: } else {
550: PetscReal *rD = (PetscReal *)D;
551: PetscReal *rqD = (PetscReal *)qD;
553: for (i = 0; i < numPoints; i++) {
554: PetscReal vol = 0.;
556: for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
557: for (k = 0; k < Nq; k++) {
558: vol += geom->detJ[i * Nq + k] * weights[k];
559: for (j = 0; j < Nc; j++) {
560: PetscReal pD[3] = {0., 0., 0.};
562: for (l = 0; l < dimC; l++) {
563: 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];
564: }
565: for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
566: }
567: }
568: for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
569: }
570: }
571: }
572: if (H) {
573: PetscInt i, j, k, l, m, q, r;
575: if (type == PETSC_SCALAR) {
576: PetscScalar *sH = (PetscScalar *)H;
577: PetscScalar *sqH = (PetscScalar *)qH;
579: for (i = 0; i < numPoints; i++) {
580: PetscReal vol = 0.;
582: for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.;
583: for (k = 0; k < Nq; k++) {
584: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
586: vol += geom->detJ[i * Nq + k] * weights[k];
587: for (j = 0; j < Nc; j++) {
588: PetscScalar pH[3][3] = {
589: {0., 0., 0.},
590: {0., 0., 0.},
591: {0., 0., 0.}
592: };
593: const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];
595: for (l = 0; l < dimC; l++) {
596: for (m = 0; m < dimC; m++) {
597: for (q = 0; q < dim; q++) {
598: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
599: }
600: }
601: }
602: for (l = 0; l < dimC; l++) {
603: 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];
604: }
605: }
606: }
607: for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
608: }
609: } else {
610: PetscReal *rH = (PetscReal *)H;
611: PetscReal *rqH = (PetscReal *)qH;
613: for (i = 0; i < numPoints; i++) {
614: PetscReal vol = 0.;
616: for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.;
617: for (k = 0; k < Nq; k++) {
618: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
620: vol += geom->detJ[i * Nq + k] * weights[k];
621: for (j = 0; j < Nc; j++) {
622: PetscReal pH[3][3] = {
623: {0., 0., 0.},
624: {0., 0., 0.},
625: {0., 0., 0.}
626: };
627: const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];
629: for (l = 0; l < dimC; l++) {
630: for (m = 0; m < dimC; m++) {
631: for (q = 0; q < dim; q++) {
632: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
633: }
634: }
635: }
636: for (l = 0; l < dimC; l++) {
637: 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];
638: }
639: }
640: }
641: for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
642: }
643: }
644: }
645: if (B) PetscCall(DMRestoreWorkArray(field->dm, N, mpitype, &qB));
646: if (D) PetscCall(DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD));
647: if (H) PetscCall(DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
648: PetscCall(PetscFEGeomDestroy(&geom));
649: PetscCall(PetscQuadratureDestroy(&quad));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
654: {
655: DMField_DS *dsfield;
656: PetscObject disc;
657: PetscInt h, imin, imax;
658: PetscClassId id;
660: PetscFunctionBegin;
661: dsfield = (DMField_DS *)field->data;
662: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
663: if (imin >= imax) {
664: h = 0;
665: } else {
666: for (h = 0; h < dsfield->height; h++) {
667: PetscInt hEnd;
669: PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
670: if (imin < hEnd) break;
671: }
672: }
673: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
674: PetscCall(PetscObjectGetClassId(disc, &id));
675: if (id == PETSCFE_CLASSID) {
676: PetscFE fe = (PetscFE)disc;
677: PetscSpace sp;
679: PetscCall(PetscFEGetBasisSpace(fe, &sp));
680: PetscCall(PetscSpaceGetDegree(sp, minDegree, maxDegree));
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
686: {
687: DM dm = field->dm;
688: const PetscInt *points;
689: DMPolytopeType ct;
690: PetscInt dim, n;
691: PetscBool isplex;
693: PetscFunctionBegin;
694: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
695: PetscCall(ISGetLocalSize(pointIS, &n));
696: if (isplex && n) {
697: PetscCall(DMGetDimension(dm, &dim));
698: PetscCall(ISGetIndices(pointIS, &points));
699: PetscCall(DMPlexGetCellType(dm, points[0], &ct));
700: switch (ct) {
701: case DM_POLYTOPE_TRIANGLE:
702: case DM_POLYTOPE_TETRAHEDRON:
703: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad));
704: break;
705: default:
706: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
707: }
708: PetscCall(ISRestoreIndices(pointIS, &points));
709: } else PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, quad));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
714: {
715: PetscInt h, dim, imax, imin, cellHeight;
716: DM dm;
717: DMField_DS *dsfield;
718: PetscObject disc;
719: PetscFE fe;
720: PetscClassId id;
722: PetscFunctionBegin;
723: dm = field->dm;
724: dsfield = (DMField_DS *)field->data;
725: PetscCall(ISGetMinMax(pointIS, &imin, &imax));
726: PetscCall(DMGetDimension(dm, &dim));
727: for (h = 0; h <= dim; h++) {
728: PetscInt hStart, hEnd;
730: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
731: if (imax >= hStart && imin < hEnd) break;
732: }
733: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
734: h -= cellHeight;
735: *quad = NULL;
736: if (h < dsfield->height) {
737: PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
738: PetscCall(PetscObjectGetClassId(disc, &id));
739: if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
740: fe = (PetscFE)disc;
741: PetscCall(PetscFEGetQuadrature(fe, quad));
742: PetscCall(PetscObjectReference((PetscObject)*quad));
743: }
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
748: {
749: const PetscInt *points;
750: PetscInt p, dim, dE, numFaces, Nq;
751: PetscInt maxDegree;
752: DMLabel depthLabel;
753: IS cellIS;
754: DM dm = field->dm;
756: PetscFunctionBegin;
757: dim = geom->dim;
758: dE = geom->dimEmbed;
759: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
760: PetscCall(DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS));
761: PetscCall(DMFieldGetDegree(field, cellIS, NULL, &maxDegree));
762: PetscCall(ISGetIndices(pointIS, &points));
763: numFaces = geom->numCells;
764: Nq = geom->numPoints;
765: /* First, set local faces and flip normals so that they are outward for the first supporting cell */
766: for (p = 0; p < numFaces; p++) {
767: PetscInt point = points[p];
768: PetscInt suppSize, s, coneSize, c, numChildren;
769: const PetscInt *supp;
771: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
772: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
773: PetscCall(DMPlexGetSupportSize(dm, point, &suppSize));
774: PetscCheck(suppSize <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, suppSize);
775: if (!suppSize) continue;
776: PetscCall(DMPlexGetSupport(dm, point, &supp));
777: for (s = 0; s < suppSize; ++s) {
778: const PetscInt *cone, *ornt;
780: PetscCall(DMPlexGetConeSize(dm, supp[s], &coneSize));
781: PetscCall(DMPlexGetOrientedCone(dm, supp[s], &cone, &ornt));
782: for (c = 0; c < coneSize; ++c)
783: if (cone[c] == point) break;
784: geom->face[p][s * 2 + 0] = c;
785: geom->face[p][s * 2 + 1] = ornt[c];
786: PetscCall(DMPlexRestoreOrientedCone(dm, supp[s], &cone, &ornt));
787: 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]);
788: }
789: if (geom->face[p][1] < 0) {
790: PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;
792: for (q = 0; q < Np; ++q)
793: for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
794: }
795: }
796: if (maxDegree <= 1) {
797: PetscInt numCells, offset, *cells;
798: PetscFEGeom *cellGeom;
799: IS suppIS;
801: for (p = 0, numCells = 0; p < numFaces; p++) {
802: PetscInt point = points[p];
803: PetscInt numSupp, numChildren;
805: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
806: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
807: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
808: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
809: numCells += numSupp;
810: }
811: PetscCall(PetscMalloc1(numCells, &cells));
812: for (p = 0, offset = 0; p < numFaces; p++) {
813: PetscInt point = points[p];
814: PetscInt numSupp, s;
815: const PetscInt *supp;
817: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
818: PetscCall(DMPlexGetSupport(dm, point, &supp));
819: for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
820: }
821: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
822: PetscCall(DMFieldCreateFEGeom(field, suppIS, quad, PETSC_FALSE, &cellGeom));
823: for (p = 0, offset = 0; p < numFaces; p++) {
824: PetscInt point = points[p];
825: PetscInt numSupp, s, q;
826: const PetscInt *supp;
828: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
829: PetscCall(DMPlexGetSupport(dm, point, &supp));
830: for (s = 0; s < numSupp; s++, offset++) {
831: for (q = 0; q < Nq * dE * dE; q++) {
832: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
833: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
834: }
835: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
836: }
837: }
838: PetscCall(PetscFEGeomDestroy(&cellGeom));
839: PetscCall(ISDestroy(&suppIS));
840: PetscCall(PetscFree(cells));
841: } else {
842: DMField_DS *dsfield = (DMField_DS *)field->data;
843: PetscObject faceDisc, cellDisc;
844: PetscClassId faceId, cellId;
845: PetscDualSpace dsp;
846: DM K;
847: DMPolytopeType ct;
848: PetscInt(*co)[2][3];
849: PetscInt coneSize;
850: PetscInt **counts;
851: PetscInt f, i, o, q, s;
852: PetscBool found = PETSC_FALSE;
853: const PetscInt *coneK;
854: PetscInt eStart, minOrient, maxOrient, numOrient;
855: PetscInt *orients;
856: PetscReal **orientPoints;
857: PetscReal *cellPoints;
858: PetscReal *dummyWeights;
859: PetscQuadrature cellQuad = NULL;
861: PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc));
862: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
863: PetscCall(PetscObjectGetClassId(faceDisc, &faceId));
864: PetscCall(PetscObjectGetClassId(cellDisc, &cellId));
865: PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported");
866: PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp));
867: PetscCall(PetscDualSpaceGetDM(dsp, &K));
868: PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL));
869: PetscCall(DMPlexGetCellType(K, eStart, &ct));
870: PetscCall(DMPlexGetConeSize(K, 0, &coneSize));
871: PetscCall(DMPlexGetCone(K, 0, &coneK));
872: PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts));
873: PetscCall(PetscMalloc1(dE * Nq, &cellPoints));
874: PetscCall(PetscMalloc1(Nq, &dummyWeights));
875: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
876: PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights));
877: minOrient = PETSC_MAX_INT;
878: maxOrient = PETSC_MIN_INT;
879: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
880: PetscInt point = points[p];
881: PetscInt numSupp, numChildren;
882: const PetscInt *supp;
884: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
885: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
886: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
887: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
888: PetscCall(DMPlexGetSupport(dm, point, &supp));
889: for (s = 0; s < numSupp; s++) {
890: PetscInt cell = supp[s];
891: PetscInt numCone;
892: const PetscInt *cone, *orient;
894: PetscCall(DMPlexGetConeSize(dm, cell, &numCone));
895: // When we extract submeshes, we hang cells from the side that are not fully realized. We ignore these
896: if (numCone == 1) {
897: co[p][s][0] = -1;
898: co[p][s][1] = -1;
899: co[p][s][2] = -1;
900: continue;
901: }
902: PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
903: PetscCall(DMPlexGetCone(dm, cell, &cone));
904: PetscCall(DMPlexGetConeOrientation(dm, cell, &orient));
905: for (f = 0; f < coneSize; f++) {
906: if (cone[f] == point) break;
907: }
908: co[p][s][0] = f;
909: co[p][s][1] = orient[f];
910: co[p][s][2] = cell;
911: minOrient = PetscMin(minOrient, orient[f]);
912: maxOrient = PetscMax(maxOrient, orient[f]);
913: found = PETSC_TRUE;
914: }
915: for (; s < 2; s++) {
916: co[p][s][0] = -1;
917: co[p][s][1] = -1;
918: co[p][s][2] = -1;
919: }
920: }
921: numOrient = found ? maxOrient + 1 - minOrient : 0;
922: PetscCall(DMPlexGetCone(K, 0, &coneK));
923: /* count all (face,orientation) doubles that appear */
924: PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints));
925: for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f]));
926: for (p = 0; p < numFaces; p++) {
927: for (s = 0; s < 2; s++) {
928: if (co[p][s][0] >= 0) {
929: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
930: orients[co[p][s][1] - minOrient]++;
931: }
932: }
933: }
934: for (o = 0; o < numOrient; o++) {
935: if (orients[o]) {
936: PetscInt orient = o + minOrient;
937: PetscInt q;
939: PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o]));
940: /* rotate the quadrature points appropriately */
941: switch (ct) {
942: case DM_POLYTOPE_POINT:
943: break;
944: case DM_POLYTOPE_SEGMENT:
945: if (orient == -2 || orient == 1) {
946: for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
947: } else {
948: for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
949: }
950: break;
951: case DM_POLYTOPE_TRIANGLE:
952: for (q = 0; q < Nq; q++) {
953: PetscReal lambda[3];
954: PetscReal lambdao[3];
956: /* convert to barycentric */
957: lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
958: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
959: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
960: if (orient >= 0) {
961: for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
962: } else {
963: for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
964: }
965: /* convert to coordinates */
966: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
967: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
968: }
969: break;
970: case DM_POLYTOPE_QUADRILATERAL:
971: for (q = 0; q < Nq; q++) {
972: PetscReal xi[2], xio[2];
973: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
975: xi[0] = geom->xi[2 * q];
976: xi[1] = geom->xi[2 * q + 1];
977: switch (oabs) {
978: case 1:
979: xio[0] = xi[1];
980: xio[1] = -xi[0];
981: break;
982: case 2:
983: xio[0] = -xi[0];
984: xio[1] = -xi[1];
985: break;
986: case 3:
987: xio[0] = -xi[1];
988: xio[1] = xi[0];
989: break;
990: case 0:
991: default:
992: xio[0] = xi[0];
993: xio[1] = xi[1];
994: break;
995: }
996: if (orient < 0) xio[0] = -xio[0];
997: orientPoints[o][2 * q + 0] = xio[0];
998: orientPoints[o][2 * q + 1] = xio[1];
999: }
1000: break;
1001: default:
1002: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
1003: }
1004: }
1005: }
1006: for (f = 0; f < coneSize; f++) {
1007: PetscInt face = coneK[f];
1008: PetscReal v0[3];
1009: PetscReal J[9], detJ;
1010: PetscInt numCells, offset;
1011: PetscInt *cells;
1012: IS suppIS;
1014: PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ));
1015: for (o = 0; o <= numOrient; o++) {
1016: PetscFEGeom *cellGeom;
1018: if (!counts[f][o]) continue;
1019: /* If this (face,orientation) double appears,
1020: * convert the face quadrature points into volume quadrature points */
1021: for (q = 0; q < Nq; q++) {
1022: PetscReal xi0[3] = {-1., -1., -1.};
1024: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1025: }
1026: for (p = 0, numCells = 0; p < numFaces; p++) {
1027: for (s = 0; s < 2; s++) {
1028: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1029: }
1030: }
1031: PetscCall(PetscMalloc1(numCells, &cells));
1032: for (p = 0, offset = 0; p < numFaces; p++) {
1033: for (s = 0; s < 2; s++) {
1034: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1035: }
1036: }
1037: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
1038: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom));
1039: for (p = 0, offset = 0; p < numFaces; p++) {
1040: for (s = 0; s < 2; s++) {
1041: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1042: for (q = 0; q < Nq * dE * dE; q++) {
1043: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1044: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1045: }
1046: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1047: offset++;
1048: }
1049: }
1050: }
1051: PetscCall(PetscFEGeomDestroy(&cellGeom));
1052: PetscCall(ISDestroy(&suppIS));
1053: PetscCall(PetscFree(cells));
1054: }
1055: }
1056: for (o = 0; o < numOrient; o++) {
1057: if (orients[o]) PetscCall(PetscFree(orientPoints[o]));
1058: }
1059: PetscCall(PetscFree2(orients, orientPoints));
1060: PetscCall(PetscQuadratureDestroy(&cellQuad));
1061: for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f]));
1062: PetscCall(PetscFree2(co, counts));
1063: }
1064: PetscCall(ISRestoreIndices(pointIS, &points));
1065: PetscCall(ISDestroy(&cellIS));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1070: {
1071: PetscFunctionBegin;
1072: field->ops->destroy = DMFieldDestroy_DS;
1073: field->ops->evaluate = DMFieldEvaluate_DS;
1074: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1075: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1076: field->ops->getDegree = DMFieldGetDegree_DS;
1077: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1078: field->ops->view = DMFieldView_DS;
1079: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1084: {
1085: DMField_DS *dsfield;
1087: PetscFunctionBegin;
1088: PetscCall(PetscNew(&dsfield));
1089: field->data = dsfield;
1090: PetscCall(DMFieldInitialize_DS(field));
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1095: {
1096: DMField b;
1097: DMField_DS *dsfield;
1098: PetscObject disc = NULL, discDG = NULL;
1099: PetscSection section;
1100: PetscBool isContainer = PETSC_FALSE;
1101: PetscClassId id = -1;
1102: PetscInt numComponents = -1, dsNumFields;
1104: PetscFunctionBegin;
1105: PetscCall(DMGetLocalSection(dm, §ion));
1106: PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents));
1107: PetscCall(DMGetNumFields(dm, &dsNumFields));
1108: if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc));
1109: if (dsNumFields && dmDG) {
1110: PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG));
1111: PetscCall(PetscObjectReference(discDG));
1112: }
1113: if (disc) {
1114: PetscCall(PetscObjectGetClassId(disc, &id));
1115: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1116: }
1117: if (!disc || isContainer) {
1118: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1119: PetscFE fe;
1120: DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1121: PetscInt dim, cStart, cEnd, cellHeight;
1123: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1124: PetscCall(DMGetDimension(dm, &dim));
1125: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1126: if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct));
1127: PetscCallMPI(MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm));
1128: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe));
1129: PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view"));
1130: disc = (PetscObject)fe;
1131: } else PetscCall(PetscObjectReference(disc));
1132: PetscCall(PetscObjectGetClassId(disc, &id));
1133: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents));
1134: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1135: PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b));
1136: PetscCall(DMFieldSetType(b, DMFIELDDS));
1137: dsfield = (DMField_DS *)b->data;
1138: dsfield->fieldNum = fieldNum;
1139: PetscCall(DMGetDimension(dm, &dsfield->height));
1140: dsfield->height++;
1141: PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc));
1142: dsfield->disc[0] = disc;
1143: PetscCall(PetscObjectReference((PetscObject)vec));
1144: dsfield->vec = vec;
1145: if (dmDG) {
1146: dsfield->dmDG = dmDG;
1147: PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG));
1148: dsfield->discDG[0] = discDG;
1149: PetscCall(PetscObjectReference((PetscObject)vecDG));
1150: dsfield->vecDG = vecDG;
1151: }
1152: *field = b;
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1157: {
1158: PetscFunctionBegin;
1159: PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }