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: const PetscInt *coneK;
853: PetscInt eStart, minOrient, maxOrient, numOrient;
854: PetscInt *orients;
855: PetscReal **orientPoints;
856: PetscReal *cellPoints;
857: PetscReal *dummyWeights;
858: PetscQuadrature cellQuad = NULL;
860: PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc));
861: PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
862: PetscCall(PetscObjectGetClassId(faceDisc, &faceId));
863: PetscCall(PetscObjectGetClassId(cellDisc, &cellId));
864: PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported");
865: PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp));
866: PetscCall(PetscDualSpaceGetDM(dsp, &K));
867: PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL));
868: PetscCall(DMPlexGetCellType(K, eStart, &ct));
869: PetscCall(DMPlexGetConeSize(K, 0, &coneSize));
870: PetscCall(DMPlexGetCone(K, 0, &coneK));
871: PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts));
872: PetscCall(PetscMalloc1(dE * Nq, &cellPoints));
873: PetscCall(PetscMalloc1(Nq, &dummyWeights));
874: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
875: PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights));
876: minOrient = PETSC_MAX_INT;
877: maxOrient = PETSC_MIN_INT;
878: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
879: PetscInt point = points[p];
880: PetscInt numSupp, numChildren;
881: const PetscInt *supp;
883: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
884: PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
885: PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
886: PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
887: PetscCall(DMPlexGetSupport(dm, point, &supp));
888: for (s = 0; s < numSupp; s++) {
889: PetscInt cell = supp[s];
890: PetscInt numCone;
891: const PetscInt *cone, *orient;
893: PetscCall(DMPlexGetConeSize(dm, cell, &numCone));
894: PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
895: PetscCall(DMPlexGetCone(dm, cell, &cone));
896: PetscCall(DMPlexGetConeOrientation(dm, cell, &orient));
897: for (f = 0; f < coneSize; f++) {
898: if (cone[f] == point) break;
899: }
900: co[p][s][0] = f;
901: co[p][s][1] = orient[f];
902: co[p][s][2] = cell;
903: minOrient = PetscMin(minOrient, orient[f]);
904: maxOrient = PetscMax(maxOrient, orient[f]);
905: }
906: for (; s < 2; s++) {
907: co[p][s][0] = -1;
908: co[p][s][1] = -1;
909: co[p][s][2] = -1;
910: }
911: }
912: numOrient = maxOrient + 1 - minOrient;
913: PetscCall(DMPlexGetCone(K, 0, &coneK));
914: /* count all (face,orientation) doubles that appear */
915: PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints));
916: for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f]));
917: for (p = 0; p < numFaces; p++) {
918: for (s = 0; s < 2; s++) {
919: if (co[p][s][0] >= 0) {
920: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
921: orients[co[p][s][1] - minOrient]++;
922: }
923: }
924: }
925: for (o = 0; o < numOrient; o++) {
926: if (orients[o]) {
927: PetscInt orient = o + minOrient;
928: PetscInt q;
930: PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o]));
931: /* rotate the quadrature points appropriately */
932: switch (ct) {
933: case DM_POLYTOPE_POINT:
934: break;
935: case DM_POLYTOPE_SEGMENT:
936: if (orient == -2 || orient == 1) {
937: for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
938: } else {
939: for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
940: }
941: break;
942: case DM_POLYTOPE_TRIANGLE:
943: for (q = 0; q < Nq; q++) {
944: PetscReal lambda[3];
945: PetscReal lambdao[3];
947: /* convert to barycentric */
948: lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
949: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
950: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
951: if (orient >= 0) {
952: for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
953: } else {
954: for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
955: }
956: /* convert to coordinates */
957: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
958: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
959: }
960: break;
961: case DM_POLYTOPE_QUADRILATERAL:
962: for (q = 0; q < Nq; q++) {
963: PetscReal xi[2], xio[2];
964: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
966: xi[0] = geom->xi[2 * q];
967: xi[1] = geom->xi[2 * q + 1];
968: switch (oabs) {
969: case 1:
970: xio[0] = xi[1];
971: xio[1] = -xi[0];
972: break;
973: case 2:
974: xio[0] = -xi[0];
975: xio[1] = -xi[1];
976: break;
977: case 3:
978: xio[0] = -xi[1];
979: xio[1] = xi[0];
980: break;
981: case 0:
982: default:
983: xio[0] = xi[0];
984: xio[1] = xi[1];
985: break;
986: }
987: if (orient < 0) xio[0] = -xio[0];
988: orientPoints[o][2 * q + 0] = xio[0];
989: orientPoints[o][2 * q + 1] = xio[1];
990: }
991: break;
992: default:
993: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
994: }
995: }
996: }
997: for (f = 0; f < coneSize; f++) {
998: PetscInt face = coneK[f];
999: PetscReal v0[3];
1000: PetscReal J[9], detJ;
1001: PetscInt numCells, offset;
1002: PetscInt *cells;
1003: IS suppIS;
1005: PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ));
1006: for (o = 0; o <= numOrient; o++) {
1007: PetscFEGeom *cellGeom;
1009: if (!counts[f][o]) continue;
1010: /* If this (face,orientation) double appears,
1011: * convert the face quadrature points into volume quadrature points */
1012: for (q = 0; q < Nq; q++) {
1013: PetscReal xi0[3] = {-1., -1., -1.};
1015: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1016: }
1017: for (p = 0, numCells = 0; p < numFaces; p++) {
1018: for (s = 0; s < 2; s++) {
1019: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1020: }
1021: }
1022: PetscCall(PetscMalloc1(numCells, &cells));
1023: for (p = 0, offset = 0; p < numFaces; p++) {
1024: for (s = 0; s < 2; s++) {
1025: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1026: }
1027: }
1028: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
1029: PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom));
1030: for (p = 0, offset = 0; p < numFaces; p++) {
1031: for (s = 0; s < 2; s++) {
1032: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1033: for (q = 0; q < Nq * dE * dE; q++) {
1034: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1035: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1036: }
1037: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1038: offset++;
1039: }
1040: }
1041: }
1042: PetscCall(PetscFEGeomDestroy(&cellGeom));
1043: PetscCall(ISDestroy(&suppIS));
1044: PetscCall(PetscFree(cells));
1045: }
1046: }
1047: for (o = 0; o < numOrient; o++) {
1048: if (orients[o]) PetscCall(PetscFree(orientPoints[o]));
1049: }
1050: PetscCall(PetscFree2(orients, orientPoints));
1051: PetscCall(PetscQuadratureDestroy(&cellQuad));
1052: for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f]));
1053: PetscCall(PetscFree2(co, counts));
1054: }
1055: PetscCall(ISRestoreIndices(pointIS, &points));
1056: PetscCall(ISDestroy(&cellIS));
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1061: {
1062: PetscFunctionBegin;
1063: field->ops->destroy = DMFieldDestroy_DS;
1064: field->ops->evaluate = DMFieldEvaluate_DS;
1065: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1066: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1067: field->ops->getDegree = DMFieldGetDegree_DS;
1068: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1069: field->ops->view = DMFieldView_DS;
1070: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1075: {
1076: DMField_DS *dsfield;
1078: PetscFunctionBegin;
1079: PetscCall(PetscNew(&dsfield));
1080: field->data = dsfield;
1081: PetscCall(DMFieldInitialize_DS(field));
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1086: {
1087: DMField b;
1088: DMField_DS *dsfield;
1089: PetscObject disc = NULL, discDG = NULL;
1090: PetscSection section;
1091: PetscBool isContainer = PETSC_FALSE;
1092: PetscClassId id = -1;
1093: PetscInt numComponents = -1, dsNumFields;
1095: PetscFunctionBegin;
1096: PetscCall(DMGetLocalSection(dm, §ion));
1097: PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents));
1098: PetscCall(DMGetNumFields(dm, &dsNumFields));
1099: if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc));
1100: if (dsNumFields && dmDG) {
1101: PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG));
1102: PetscCall(PetscObjectReference(discDG));
1103: }
1104: if (disc) {
1105: PetscCall(PetscObjectGetClassId(disc, &id));
1106: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1107: }
1108: if (!disc || isContainer) {
1109: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1110: PetscFE fe;
1111: DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1112: PetscInt dim, cStart, cEnd, cellHeight;
1114: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1115: PetscCall(DMGetDimension(dm, &dim));
1116: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1117: if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct));
1118: PetscCallMPI(MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm));
1119: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe));
1120: PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view"));
1121: disc = (PetscObject)fe;
1122: } else PetscCall(PetscObjectReference(disc));
1123: PetscCall(PetscObjectGetClassId(disc, &id));
1124: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents));
1125: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1126: PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b));
1127: PetscCall(DMFieldSetType(b, DMFIELDDS));
1128: dsfield = (DMField_DS *)b->data;
1129: dsfield->fieldNum = fieldNum;
1130: PetscCall(DMGetDimension(dm, &dsfield->height));
1131: dsfield->height++;
1132: PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc));
1133: dsfield->disc[0] = disc;
1134: PetscCall(PetscObjectReference((PetscObject)vec));
1135: dsfield->vec = vec;
1136: if (dmDG) {
1137: dsfield->dmDG = dmDG;
1138: PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG));
1139: dsfield->discDG[0] = discDG;
1140: PetscCall(PetscObjectReference((PetscObject)vecDG));
1141: dsfield->vecDG = vecDG;
1142: }
1143: *field = b;
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1148: {
1149: PetscFunctionBegin;
1150: PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }