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