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