Actual source code: dmfieldds.c

  1: #include <petsc/private/dmfieldimpl.h>
  2: #include <petsc/private/petscfeimpl.h>
  3: #include <petscfe.h>
  4: #include <petscdmplex.h>
  5: #include <petscds.h>

  7: typedef struct _n_DMField_DS {
  8:   PetscBool    multifieldVec;
  9:   PetscInt     height;   /* Point height at which we want values and number of discretizations */
 10:   PetscInt     fieldNum; /* Number in DS of field which we evaluate */
 11:   PetscObject *disc;     /* Discretizations of this field at each height */
 12:   Vec          vec;      /* Field values */
 13:   DM           dmDG;     /* DM for the DG values */
 14:   PetscObject *discDG;   /* DG Discretizations of this field at each height */
 15:   Vec          vecDG;    /* DG Field values */
 16: } DMField_DS;

 18: static PetscErrorCode DMFieldDestroy_DS(DMField field)
 19: {
 20:   DMField_DS *dsfield = (DMField_DS *)field->data;
 21:   PetscInt    i;

 23:   PetscFunctionBegin;
 24:   PetscCall(VecDestroy(&dsfield->vec));
 25:   for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->disc[i]));
 26:   PetscCall(PetscFree(dsfield->disc));
 27:   PetscCall(VecDestroy(&dsfield->vecDG));
 28:   if (dsfield->discDG)
 29:     for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->discDG[i]));
 30:   PetscCall(PetscFree(dsfield->discDG));
 31:   PetscCall(PetscFree(dsfield));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: static PetscErrorCode DMFieldView_DS(DMField field, PetscViewer viewer)
 36: {
 37:   DMField_DS *dsfield = (DMField_DS *)field->data;
 38:   PetscObject disc;
 39:   PetscBool   iascii;

 41:   PetscFunctionBegin;
 42:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 43:   disc = dsfield->disc[0];
 44:   if (iascii) {
 45:     PetscCall(PetscViewerASCIIPrintf(viewer, "PetscDS field %" PetscInt_FMT "\n", dsfield->fieldNum));
 46:     PetscCall(PetscViewerASCIIPushTab(viewer));
 47:     PetscCall(PetscObjectView(disc, viewer));
 48:     PetscCall(PetscViewerASCIIPopTab(viewer));
 49:   }
 50:   PetscCall(PetscViewerASCIIPushTab(viewer));
 51:   PetscCheck(!dsfield->multifieldVec, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "View of subfield not implemented yet");
 52:   PetscCall(VecView(dsfield->vec, viewer));
 53:   PetscCall(PetscViewerASCIIPopTab(viewer));
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 55: }

 57: static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject discList[], PetscObject *disc)
 58: {
 59:   PetscFunctionBegin;
 60:   if (!discList[height]) {
 61:     PetscClassId id;

 63:     PetscCall(PetscObjectGetClassId(discList[0], &id));
 64:     if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height]));
 65:   }
 66:   *disc = discList[height];
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: /* y[m,c] = A[m,n,c] . b[n] */
 71: #define DMFieldDSdot(y, A, b, m, n, c, cast) \
 72:   do { \
 73:     PetscInt _i, _j, _k; \
 74:     for (_i = 0; _i < (m); _i++) { \
 75:       for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] = 0.; \
 76:       for (_j = 0; _j < (n); _j++) { \
 77:         for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \
 78:       } \
 79:     } \
 80:   } while (0)

 82: /*
 83:   Since this is used for coordinates, we need to allow for the possibility that values come from multiple sections/Vecs, so that we can have DG version of the coordinates for periodicity. This reproduces DMPlexGetCellCoordinates_Internal().
 84: */
 85: static PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
 86: {
 87:   DMField_DS        *dsfield = (DMField_DS *)field->data;
 88:   DM                 fdm     = dsfield->dmDG;
 89:   PetscSection       s       = NULL;
 90:   const PetscScalar *cvalues;
 91:   PetscInt           pStart, pEnd;

 93:   PetscFunctionBeginHot;
 94:   *isDG   = PETSC_FALSE;
 95:   *Nc     = 0;
 96:   *array  = NULL;
 97:   *values = NULL;
 98:   /* Check for cellwise section */
 99:   if (fdm) PetscCall(DMGetLocalSection(fdm, &s));
100:   if (!s) goto cg;
101:   /* Check that the cell exists in the cellwise section */
102:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
103:   if (cell < pStart || cell >= pEnd) goto cg;
104:   /* Check for cellwise coordinates for this cell */
105:   PetscCall(PetscSectionGetDof(s, cell, Nc));
106:   if (!*Nc) goto cg;
107:   /* Check for cellwise coordinates */
108:   if (!dsfield->vecDG) goto cg;
109:   /* Get cellwise coordinates */
110:   PetscCall(VecGetArrayRead(dsfield->vecDG, array));
111:   PetscCall(DMPlexPointLocalRead(fdm, cell, *array, &cvalues));
112:   PetscCall(DMGetWorkArray(fdm, *Nc, MPIU_SCALAR, values));
113:   PetscCall(PetscArraycpy(*values, cvalues, *Nc));
114:   PetscCall(VecRestoreArrayRead(dsfield->vecDG, array));
115:   *isDG = PETSC_TRUE;
116:   PetscFunctionReturn(PETSC_SUCCESS);
117: cg:
118:   /* Use continuous values */
119:   PetscCall(DMFieldGetDM(field, &fdm));
120:   PetscCall(DMGetLocalSection(fdm, &s));
121:   PetscCall(PetscSectionGetField(s, dsfield->fieldNum, &s));
122:   PetscCall(DMPlexVecGetClosure(fdm, s, dsfield->vec, cell, Nc, values));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: static PetscErrorCode DMFieldRestoreClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
127: {
128:   DMField_DS  *dsfield = (DMField_DS *)field->data;
129:   DM           fdm;
130:   PetscSection s;

132:   PetscFunctionBeginHot;
133:   if (*isDG) {
134:     PetscCall(DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values));
135:   } else {
136:     PetscCall(DMFieldGetDM(field, &fdm));
137:     PetscCall(DMGetLocalSection(fdm, &s));
138:     PetscCall(DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, (PetscScalar **)values));
139:   }
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
144: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
145: {
146:   DMField_DS      *dsfield = (DMField_DS *)field->data;
147:   DM               dm;
148:   PetscObject      disc;
149:   PetscClassId     classid;
150:   PetscInt         nq, nc, dim, meshDim, numCells;
151:   PetscSection     section;
152:   const PetscReal *qpoints;
153:   PetscBool        isStride;
154:   const PetscInt  *points = NULL;
155:   PetscInt         sfirst = -1, stride = -1;

157:   PetscFunctionBeginHot;
158:   dm = field->dm;
159:   nc = field->numComponents;
160:   PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL));
161:   PetscCall(DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc));
162:   PetscCall(DMGetDimension(dm, &meshDim));
163:   PetscCall(DMGetLocalSection(dm, &section));
164:   PetscCall(PetscSectionGetField(section, dsfield->fieldNum, &section));
165:   PetscCall(PetscObjectGetClassId(disc, &classid));
166:   /* TODO: batch */
167:   PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
168:   PetscCall(ISGetLocalSize(pointIS, &numCells));
169:   if (isStride) PetscCall(ISStrideGetInfo(pointIS, &sfirst, &stride));
170:   else PetscCall(ISGetIndices(pointIS, &points));
171:   if (classid == PETSCFE_CLASSID) {
172:     PetscFE         fe = (PetscFE)disc;
173:     PetscInt        feDim, i;
174:     PetscInt        K = H ? 2 : (D ? 1 : (B ? 0 : -1));
175:     PetscTabulation T;

177:     PetscCall(PetscFEGetDimension(fe, &feDim));
178:     PetscCall(PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T));
179:     for (i = 0; i < numCells; i++) {
180:       PetscInt           c = isStride ? (sfirst + i * stride) : points[i];
181:       PetscInt           closureSize;
182:       const PetscScalar *array;
183:       PetscScalar       *elem = NULL;
184:       PetscBool          isDG;

186:       PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
187:       if (B) {
188:         /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
189:         if (type == PETSC_SCALAR) {
190:           PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i];

192:           DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
193:         } else {
194:           PetscReal *cB = &((PetscReal *)B)[nc * nq * i];

196:           DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
197:         }
198:       }
199:       if (D) {
200:         if (type == PETSC_SCALAR) {
201:           PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i];

203:           DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
204:         } else {
205:           PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i];

207:           DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
208:         }
209:       }
210:       if (H) {
211:         if (type == PETSC_SCALAR) {
212:           PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i];

214:           DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
215:         } else {
216:           PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i];

218:           DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
219:         }
220:       }
221:       PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
222:     }
223:     PetscCall(PetscTabulationDestroy(&T));
224:   } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented");
225:   if (!isStride) PetscCall(ISRestoreIndices(pointIS, &points));
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

229: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
230: {
231:   DMField_DS        *dsfield = (DMField_DS *)field->data;
232:   PetscSF            cellSF  = NULL;
233:   const PetscSFNode *cells;
234:   PetscInt           c, nFound, numCells, feDim, nc;
235:   const PetscInt    *cellDegrees;
236:   const PetscScalar *pointsArray;
237:   PetscScalar       *cellPoints;
238:   PetscInt           gatherSize, gatherMax;
239:   PetscInt           dim, dimR, offset;
240:   MPI_Datatype       pointType;
241:   PetscObject        cellDisc;
242:   PetscFE            cellFE;
243:   PetscClassId       discID;
244:   PetscReal         *coordsReal, *coordsRef;
245:   PetscSection       section;
246:   PetscScalar       *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
247:   PetscReal         *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
248:   PetscReal         *v, *J, *invJ, *detJ;

250:   PetscFunctionBegin;
251:   nc = field->numComponents;
252:   PetscCall(DMGetLocalSection(field->dm, &section));
253:   PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
254:   PetscCall(PetscObjectGetClassId(cellDisc, &discID));
255:   PetscCheck(discID == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization type not supported");
256:   cellFE = (PetscFE)cellDisc;
257:   PetscCall(PetscFEGetDimension(cellFE, &feDim));
258:   PetscCall(DMGetCoordinateDim(field->dm, &dim));
259:   PetscCall(DMGetDimension(field->dm, &dimR));
260:   PetscCall(DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF));
261:   PetscCall(PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells));
262:   for (c = 0; c < nFound; c++) PetscCheck(cells[c].index >= 0, PetscObjectComm((PetscObject)points), PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " could not be located", c);
263:   PetscCall(PetscSFComputeDegreeBegin(cellSF, &cellDegrees));
264:   PetscCall(PetscSFComputeDegreeEnd(cellSF, &cellDegrees));
265:   for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
266:     gatherMax = PetscMax(gatherMax, cellDegrees[c]);
267:     gatherSize += cellDegrees[c];
268:   }
269:   PetscCall(PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef));
270:   PetscCall(PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ));
271:   if (datatype == PETSC_SCALAR) PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs));
272:   else PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr));

274:   PetscCallMPI(MPI_Type_contiguous(dim, MPIU_SCALAR, &pointType));
275:   PetscCallMPI(MPI_Type_commit(&pointType));
276:   PetscCall(VecGetArrayRead(points, &pointsArray));
277:   PetscCall(PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints));
278:   PetscCall(PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints));
279:   PetscCall(VecRestoreArrayRead(points, &pointsArray));
280:   for (c = 0, offset = 0; c < numCells; c++) {
281:     PetscInt nq = cellDegrees[c], p;

283:     if (nq) {
284:       PetscInt           K = H ? 2 : (D ? 1 : (B ? 0 : -1));
285:       PetscTabulation    T;
286:       PetscQuadrature    quad;
287:       const PetscScalar *array;
288:       PetscScalar       *elem = NULL;
289:       PetscReal         *quadPoints;
290:       PetscBool          isDG;
291:       PetscInt           closureSize, d, e, f, g;

293:       for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
294:       PetscCall(DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef));
295:       PetscCall(PetscFECreateTabulation(cellFE, 1, nq, coordsRef, K, &T));
296:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
297:       PetscCall(PetscMalloc1(dimR * nq, &quadPoints));
298:       for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
299:       PetscCall(PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL));
300:       PetscCall(DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ));
301:       PetscCall(PetscQuadratureDestroy(&quad));
302:       PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
303:       if (B) {
304:         if (datatype == PETSC_SCALAR) {
305:           PetscScalar *cB = &cellBs[nc * offset];

307:           DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
308:         } else {
309:           PetscReal *cB = &cellBr[nc * offset];

311:           DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
312:         }
313:       }
314:       if (D) {
315:         if (datatype == PETSC_SCALAR) {
316:           PetscScalar *cD = &cellDs[nc * dim * offset];

318:           DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
319:           for (p = 0; p < nq; p++) {
320:             for (g = 0; g < nc; g++) {
321:               PetscScalar vs[3];

323:               for (d = 0; d < dimR; d++) {
324:                 vs[d] = 0.;
325:                 for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
326:               }
327:               for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = vs[d];
328:             }
329:           }
330:         } else {
331:           PetscReal *cD = &cellDr[nc * dim * offset];

333:           DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
334:           for (p = 0; p < nq; p++) {
335:             for (g = 0; g < nc; g++) {
336:               for (d = 0; d < dimR; d++) {
337:                 v[d] = 0.;
338:                 for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
339:               }
340:               for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = v[d];
341:             }
342:           }
343:         }
344:       }
345:       if (H) {
346:         if (datatype == PETSC_SCALAR) {
347:           PetscScalar *cH = &cellHs[nc * dim * dim * offset];

349:           DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
350:           for (p = 0; p < nq; p++) {
351:             for (g = 0; g < nc * dimR; g++) {
352:               PetscScalar vs[3];

354:               for (d = 0; d < dimR; d++) {
355:                 vs[d] = 0.;
356:                 for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
357:               }
358:               for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = vs[d];
359:             }
360:             for (g = 0; g < nc; g++) {
361:               for (f = 0; f < dimR; f++) {
362:                 PetscScalar vs[3];

364:                 for (d = 0; d < dimR; d++) {
365:                   vs[d] = 0.;
366:                   for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
367:                 }
368:                 for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
369:               }
370:             }
371:           }
372:         } else {
373:           PetscReal *cH = &cellHr[nc * dim * dim * offset];

375:           DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
376:           for (p = 0; p < nq; p++) {
377:             for (g = 0; g < nc * dimR; g++) {
378:               for (d = 0; d < dimR; d++) {
379:                 v[d] = 0.;
380:                 for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
381:               }
382:               for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = v[d];
383:             }
384:             for (g = 0; g < nc; g++) {
385:               for (f = 0; f < dimR; f++) {
386:                 for (d = 0; d < dimR; d++) {
387:                   v[d] = 0.;
388:                   for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
389:                 }
390:                 for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
391:               }
392:             }
393:           }
394:         }
395:       }
396:       PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem));
397:       PetscCall(PetscTabulationDestroy(&T));
398:     }
399:     offset += nq;
400:   }
401:   {
402:     MPI_Datatype origtype;

404:     if (datatype == PETSC_SCALAR) {
405:       origtype = MPIU_SCALAR;
406:     } else {
407:       origtype = MPIU_REAL;
408:     }
409:     if (B) {
410:       MPI_Datatype Btype;

412:       PetscCallMPI(MPI_Type_contiguous(nc, origtype, &Btype));
413:       PetscCallMPI(MPI_Type_commit(&Btype));
414:       PetscCall(PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
415:       PetscCall(PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B));
416:       PetscCallMPI(MPI_Type_free(&Btype));
417:     }
418:     if (D) {
419:       MPI_Datatype Dtype;

421:       PetscCallMPI(MPI_Type_contiguous(nc * dim, origtype, &Dtype));
422:       PetscCallMPI(MPI_Type_commit(&Dtype));
423:       PetscCall(PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
424:       PetscCall(PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D));
425:       PetscCallMPI(MPI_Type_free(&Dtype));
426:     }
427:     if (H) {
428:       MPI_Datatype Htype;

430:       PetscCallMPI(MPI_Type_contiguous(nc * dim * dim, origtype, &Htype));
431:       PetscCallMPI(MPI_Type_commit(&Htype));
432:       PetscCall(PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
433:       PetscCall(PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H));
434:       PetscCallMPI(MPI_Type_free(&Htype));
435:     }
436:   }
437:   PetscCall(PetscFree4(v, J, invJ, detJ));
438:   PetscCall(PetscFree3(cellBr, cellDr, cellHr));
439:   PetscCall(PetscFree3(cellBs, cellDs, cellHs));
440:   PetscCall(PetscFree3(cellPoints, coordsReal, coordsRef));
441:   PetscCallMPI(MPI_Type_free(&pointType));
442:   PetscCall(PetscSFDestroy(&cellSF));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
447: {
448:   DMField_DS      *dsfield = (DMField_DS *)field->data;
449:   PetscInt         h, imin;
450:   PetscInt         dim;
451:   PetscClassId     id;
452:   PetscQuadrature  quad = NULL;
453:   PetscInt         maxDegree;
454:   PetscFEGeom     *geom;
455:   PetscInt         Nq, Nc, dimC, qNc, N;
456:   PetscInt         numPoints;
457:   void            *qB = NULL, *qD = NULL, *qH = NULL;
458:   const PetscReal *weights;
459:   MPI_Datatype     mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
460:   PetscObject      disc;
461:   DMField          coordField;

463:   PetscFunctionBegin;
464:   Nc = field->numComponents;
465:   PetscCall(DMGetCoordinateDim(field->dm, &dimC));
466:   PetscCall(DMGetDimension(field->dm, &dim));
467:   PetscCall(ISGetLocalSize(pointIS, &numPoints));
468:   PetscCall(ISGetMinMax(pointIS, &imin, NULL));
469:   for (h = 0; h < dsfield->height; h++) {
470:     PetscInt hEnd;

472:     PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
473:     if (imin < hEnd) break;
474:   }
475:   dim -= h;
476:   PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
477:   PetscCall(PetscObjectGetClassId(disc, &id));
478:   PetscCheck(id == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported");
479:   PetscCall(DMGetCoordinateField(field->dm, &coordField));
480:   PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree));
481:   if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad));
482:   if (!quad) PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, &quad));
483:   PetscCheck(quad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages");
484:   PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom));
485:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights));
486:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components");
487:   N = numPoints * Nq * Nc;
488:   if (B) PetscCall(DMGetWorkArray(field->dm, N, mpitype, &qB));
489:   if (D) PetscCall(DMGetWorkArray(field->dm, N * dimC, mpitype, &qD));
490:   if (H) PetscCall(DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
491:   PetscCall(DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH));
492:   if (B) {
493:     PetscInt i, j, k;

495:     if (type == PETSC_SCALAR) {
496:       PetscScalar *sB  = (PetscScalar *)B;
497:       PetscScalar *sqB = (PetscScalar *)qB;

499:       for (i = 0; i < numPoints; i++) {
500:         PetscReal vol = 0.;

502:         for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.;
503:         for (k = 0; k < Nq; k++) {
504:           vol += geom->detJ[i * Nq + k] * weights[k];
505:           for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j];
506:         }
507:         for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
508:       }
509:     } else {
510:       PetscReal *rB  = (PetscReal *)B;
511:       PetscReal *rqB = (PetscReal *)qB;

513:       for (i = 0; i < numPoints; i++) {
514:         PetscReal vol = 0.;

516:         for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.;
517:         for (k = 0; k < Nq; k++) {
518:           vol += geom->detJ[i * Nq + k] * weights[k];
519:           for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j];
520:         }
521:         for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
522:       }
523:     }
524:   }
525:   if (D) {
526:     PetscInt i, j, k, l, m;

528:     if (type == PETSC_SCALAR) {
529:       PetscScalar *sD  = (PetscScalar *)D;
530:       PetscScalar *sqD = (PetscScalar *)qD;

532:       for (i = 0; i < numPoints; i++) {
533:         PetscReal vol = 0.;

535:         for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
536:         for (k = 0; k < Nq; k++) {
537:           vol += geom->detJ[i * Nq + k] * weights[k];
538:           for (j = 0; j < Nc; j++) {
539:             PetscScalar pD[3] = {0., 0., 0.};

541:             for (l = 0; l < dimC; l++) {
542:               for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m];
543:             }
544:             for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
545:           }
546:         }
547:         for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
548:       }
549:     } else {
550:       PetscReal *rD  = (PetscReal *)D;
551:       PetscReal *rqD = (PetscReal *)qD;

553:       for (i = 0; i < numPoints; i++) {
554:         PetscReal vol = 0.;

556:         for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
557:         for (k = 0; k < Nq; k++) {
558:           vol += geom->detJ[i * Nq + k] * weights[k];
559:           for (j = 0; j < Nc; j++) {
560:             PetscReal pD[3] = {0., 0., 0.};

562:             for (l = 0; l < dimC; l++) {
563:               for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m];
564:             }
565:             for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
566:           }
567:         }
568:         for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
569:       }
570:     }
571:   }
572:   if (H) {
573:     PetscInt i, j, k, l, m, q, r;

575:     if (type == PETSC_SCALAR) {
576:       PetscScalar *sH  = (PetscScalar *)H;
577:       PetscScalar *sqH = (PetscScalar *)qH;

579:       for (i = 0; i < numPoints; i++) {
580:         PetscReal vol = 0.;

582:         for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.;
583:         for (k = 0; k < Nq; k++) {
584:           const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];

586:           vol += geom->detJ[i * Nq + k] * weights[k];
587:           for (j = 0; j < Nc; j++) {
588:             PetscScalar pH[3][3] = {
589:               {0., 0., 0.},
590:               {0., 0., 0.},
591:               {0., 0., 0.}
592:             };
593:             const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];

595:             for (l = 0; l < dimC; l++) {
596:               for (m = 0; m < dimC; m++) {
597:                 for (q = 0; q < dim; q++) {
598:                   for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
599:                 }
600:               }
601:             }
602:             for (l = 0; l < dimC; l++) {
603:               for (m = 0; m < dimC; m++) sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
604:             }
605:           }
606:         }
607:         for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
608:       }
609:     } else {
610:       PetscReal *rH  = (PetscReal *)H;
611:       PetscReal *rqH = (PetscReal *)qH;

613:       for (i = 0; i < numPoints; i++) {
614:         PetscReal vol = 0.;

616:         for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.;
617:         for (k = 0; k < Nq; k++) {
618:           const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];

620:           vol += geom->detJ[i * Nq + k] * weights[k];
621:           for (j = 0; j < Nc; j++) {
622:             PetscReal pH[3][3] = {
623:               {0., 0., 0.},
624:               {0., 0., 0.},
625:               {0., 0., 0.}
626:             };
627:             const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];

629:             for (l = 0; l < dimC; l++) {
630:               for (m = 0; m < dimC; m++) {
631:                 for (q = 0; q < dim; q++) {
632:                   for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
633:                 }
634:               }
635:             }
636:             for (l = 0; l < dimC; l++) {
637:               for (m = 0; m < dimC; m++) rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
638:             }
639:           }
640:         }
641:         for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
642:       }
643:     }
644:   }
645:   if (B) PetscCall(DMRestoreWorkArray(field->dm, N, mpitype, &qB));
646:   if (D) PetscCall(DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD));
647:   if (H) PetscCall(DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH));
648:   PetscCall(PetscFEGeomDestroy(&geom));
649:   PetscCall(PetscQuadratureDestroy(&quad));
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
654: {
655:   DMField_DS  *dsfield;
656:   PetscObject  disc;
657:   PetscInt     h, imin, imax;
658:   PetscClassId id;

660:   PetscFunctionBegin;
661:   dsfield = (DMField_DS *)field->data;
662:   PetscCall(ISGetMinMax(pointIS, &imin, &imax));
663:   if (imin >= imax) {
664:     h = 0;
665:   } else {
666:     for (h = 0; h < dsfield->height; h++) {
667:       PetscInt hEnd;

669:       PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd));
670:       if (imin < hEnd) break;
671:     }
672:   }
673:   PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
674:   PetscCall(PetscObjectGetClassId(disc, &id));
675:   if (id == PETSCFE_CLASSID) {
676:     PetscFE    fe = (PetscFE)disc;
677:     PetscSpace sp;

679:     PetscCall(PetscFEGetBasisSpace(fe, &sp));
680:     PetscCall(PetscSpaceGetDegree(sp, minDegree, maxDegree));
681:   }
682:   PetscFunctionReturn(PETSC_SUCCESS);
683: }

685: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
686: {
687:   DM              dm = field->dm;
688:   const PetscInt *points;
689:   DMPolytopeType  ct;
690:   PetscInt        dim, n;
691:   PetscBool       isplex;

693:   PetscFunctionBegin;
694:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
695:   PetscCall(ISGetLocalSize(pointIS, &n));
696:   if (isplex && n) {
697:     PetscCall(DMGetDimension(dm, &dim));
698:     PetscCall(ISGetIndices(pointIS, &points));
699:     PetscCall(DMPlexGetCellType(dm, points[0], &ct));
700:     switch (ct) {
701:     case DM_POLYTOPE_TRIANGLE:
702:     case DM_POLYTOPE_TETRAHEDRON:
703:       PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad));
704:       break;
705:     default:
706:       PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
707:     }
708:     PetscCall(ISRestoreIndices(pointIS, &points));
709:   } else PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, quad));
710:   PetscFunctionReturn(PETSC_SUCCESS);
711: }

713: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
714: {
715:   PetscInt     h, dim, imax, imin, cellHeight;
716:   DM           dm;
717:   DMField_DS  *dsfield;
718:   PetscObject  disc;
719:   PetscFE      fe;
720:   PetscClassId id;

722:   PetscFunctionBegin;
723:   dm      = field->dm;
724:   dsfield = (DMField_DS *)field->data;
725:   PetscCall(ISGetMinMax(pointIS, &imin, &imax));
726:   PetscCall(DMGetDimension(dm, &dim));
727:   for (h = 0; h <= dim; h++) {
728:     PetscInt hStart, hEnd;

730:     PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
731:     if (imax >= hStart && imin < hEnd) break;
732:   }
733:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
734:   h -= cellHeight;
735:   *quad = NULL;
736:   if (h < dsfield->height) {
737:     PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc));
738:     PetscCall(PetscObjectGetClassId(disc, &id));
739:     if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
740:     fe = (PetscFE)disc;
741:     PetscCall(PetscFEGetQuadrature(fe, quad));
742:     PetscCall(PetscObjectReference((PetscObject)*quad));
743:   }
744:   PetscFunctionReturn(PETSC_SUCCESS);
745: }

747: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
748: {
749:   const PetscInt *points;
750:   PetscInt        p, dim, dE, numFaces, Nq;
751:   PetscInt        maxDegree;
752:   DMLabel         depthLabel;
753:   IS              cellIS;
754:   DM              dm = field->dm;

756:   PetscFunctionBegin;
757:   dim = geom->dim;
758:   dE  = geom->dimEmbed;
759:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
760:   PetscCall(DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS));
761:   PetscCall(DMFieldGetDegree(field, cellIS, NULL, &maxDegree));
762:   PetscCall(ISGetIndices(pointIS, &points));
763:   numFaces = geom->numCells;
764:   Nq       = geom->numPoints;
765:   /* First, set local faces and flip normals so that they are outward for the first supporting cell */
766:   for (p = 0; p < numFaces; p++) {
767:     PetscInt        point = points[p];
768:     PetscInt        suppSize, s, coneSize, c, numChildren;
769:     const PetscInt *supp;

771:     PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
772:     PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
773:     PetscCall(DMPlexGetSupportSize(dm, point, &suppSize));
774:     PetscCheck(suppSize <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, suppSize);
775:     if (!suppSize) continue;
776:     PetscCall(DMPlexGetSupport(dm, point, &supp));
777:     for (s = 0; s < suppSize; ++s) {
778:       const PetscInt *cone, *ornt;

780:       PetscCall(DMPlexGetConeSize(dm, supp[s], &coneSize));
781:       PetscCall(DMPlexGetOrientedCone(dm, supp[s], &cone, &ornt));
782:       for (c = 0; c < coneSize; ++c)
783:         if (cone[c] == point) break;
784:       geom->face[p][s * 2 + 0] = c;
785:       geom->face[p][s * 2 + 1] = ornt[c];
786:       PetscCall(DMPlexRestoreOrientedCone(dm, supp[s], &cone, &ornt));
787:       PetscCheck(c != coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid connectivity: point %" PetscInt_FMT " not found in cone of support point %" PetscInt_FMT, point, supp[s]);
788:     }
789:     if (geom->face[p][1] < 0) {
790:       PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;

792:       for (q = 0; q < Np; ++q)
793:         for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
794:     }
795:   }
796:   if (maxDegree <= 1) {
797:     PetscInt     numCells, offset, *cells;
798:     PetscFEGeom *cellGeom;
799:     IS           suppIS;

801:     for (p = 0, numCells = 0; p < numFaces; p++) {
802:       PetscInt point = points[p];
803:       PetscInt numSupp, numChildren;

805:       PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
806:       PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
807:       PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
808:       PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
809:       numCells += numSupp;
810:     }
811:     PetscCall(PetscMalloc1(numCells, &cells));
812:     for (p = 0, offset = 0; p < numFaces; p++) {
813:       PetscInt        point = points[p];
814:       PetscInt        numSupp, s;
815:       const PetscInt *supp;

817:       PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
818:       PetscCall(DMPlexGetSupport(dm, point, &supp));
819:       for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
820:     }
821:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
822:     PetscCall(DMFieldCreateFEGeom(field, suppIS, quad, PETSC_FALSE, &cellGeom));
823:     for (p = 0, offset = 0; p < numFaces; p++) {
824:       PetscInt        point = points[p];
825:       PetscInt        numSupp, s, q;
826:       const PetscInt *supp;

828:       PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
829:       PetscCall(DMPlexGetSupport(dm, point, &supp));
830:       for (s = 0; s < numSupp; s++, offset++) {
831:         for (q = 0; q < Nq * dE * dE; q++) {
832:           geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
833:           geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
834:         }
835:         for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
836:       }
837:     }
838:     PetscCall(PetscFEGeomDestroy(&cellGeom));
839:     PetscCall(ISDestroy(&suppIS));
840:     PetscCall(PetscFree(cells));
841:   } else {
842:     DMField_DS    *dsfield = (DMField_DS *)field->data;
843:     PetscObject    faceDisc, cellDisc;
844:     PetscClassId   faceId, cellId;
845:     PetscDualSpace dsp;
846:     DM             K;
847:     DMPolytopeType ct;
848:     PetscInt(*co)[2][3];
849:     PetscInt        coneSize;
850:     PetscInt      **counts;
851:     PetscInt        f, i, o, q, s;
852:     PetscBool       found = PETSC_FALSE;
853:     const PetscInt *coneK;
854:     PetscInt        eStart, minOrient, maxOrient, numOrient;
855:     PetscInt       *orients;
856:     PetscReal     **orientPoints;
857:     PetscReal      *cellPoints;
858:     PetscReal      *dummyWeights;
859:     PetscQuadrature cellQuad = NULL;

861:     PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc));
862:     PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc));
863:     PetscCall(PetscObjectGetClassId(faceDisc, &faceId));
864:     PetscCall(PetscObjectGetClassId(cellDisc, &cellId));
865:     PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported");
866:     PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp));
867:     PetscCall(PetscDualSpaceGetDM(dsp, &K));
868:     PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL));
869:     PetscCall(DMPlexGetCellType(K, eStart, &ct));
870:     PetscCall(DMPlexGetConeSize(K, 0, &coneSize));
871:     PetscCall(DMPlexGetCone(K, 0, &coneK));
872:     PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts));
873:     PetscCall(PetscMalloc1(dE * Nq, &cellPoints));
874:     PetscCall(PetscMalloc1(Nq, &dummyWeights));
875:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad));
876:     PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights));
877:     minOrient = PETSC_MAX_INT;
878:     maxOrient = PETSC_MIN_INT;
879:     for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
880:       PetscInt        point = points[p];
881:       PetscInt        numSupp, numChildren;
882:       const PetscInt *supp;

884:       PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL));
885:       PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
886:       PetscCall(DMPlexGetSupportSize(dm, point, &numSupp));
887:       PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp);
888:       PetscCall(DMPlexGetSupport(dm, point, &supp));
889:       for (s = 0; s < numSupp; s++) {
890:         PetscInt        cell = supp[s];
891:         PetscInt        numCone;
892:         const PetscInt *cone, *orient;

894:         PetscCall(DMPlexGetConeSize(dm, cell, &numCone));
895:         // When we extract submeshes, we hang cells from the side that are not fully realized. We ignore these
896:         if (numCone == 1) {
897:           co[p][s][0] = -1;
898:           co[p][s][1] = -1;
899:           co[p][s][2] = -1;
900:           continue;
901:         }
902:         PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
903:         PetscCall(DMPlexGetCone(dm, cell, &cone));
904:         PetscCall(DMPlexGetConeOrientation(dm, cell, &orient));
905:         for (f = 0; f < coneSize; f++) {
906:           if (cone[f] == point) break;
907:         }
908:         co[p][s][0] = f;
909:         co[p][s][1] = orient[f];
910:         co[p][s][2] = cell;
911:         minOrient   = PetscMin(minOrient, orient[f]);
912:         maxOrient   = PetscMax(maxOrient, orient[f]);
913:         found       = PETSC_TRUE;
914:       }
915:       for (; s < 2; s++) {
916:         co[p][s][0] = -1;
917:         co[p][s][1] = -1;
918:         co[p][s][2] = -1;
919:       }
920:     }
921:     numOrient = found ? maxOrient + 1 - minOrient : 0;
922:     PetscCall(DMPlexGetCone(K, 0, &coneK));
923:     /* count all (face,orientation) doubles that appear */
924:     PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints));
925:     for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f]));
926:     for (p = 0; p < numFaces; p++) {
927:       for (s = 0; s < 2; s++) {
928:         if (co[p][s][0] >= 0) {
929:           counts[co[p][s][0]][co[p][s][1] - minOrient]++;
930:           orients[co[p][s][1] - minOrient]++;
931:         }
932:       }
933:     }
934:     for (o = 0; o < numOrient; o++) {
935:       if (orients[o]) {
936:         PetscInt orient = o + minOrient;
937:         PetscInt q;

939:         PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o]));
940:         /* rotate the quadrature points appropriately */
941:         switch (ct) {
942:         case DM_POLYTOPE_POINT:
943:           break;
944:         case DM_POLYTOPE_SEGMENT:
945:           if (orient == -2 || orient == 1) {
946:             for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
947:           } else {
948:             for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
949:           }
950:           break;
951:         case DM_POLYTOPE_TRIANGLE:
952:           for (q = 0; q < Nq; q++) {
953:             PetscReal lambda[3];
954:             PetscReal lambdao[3];

956:             /* convert to barycentric */
957:             lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
958:             lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
959:             lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
960:             if (orient >= 0) {
961:               for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
962:             } else {
963:               for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
964:             }
965:             /* convert to coordinates */
966:             orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
967:             orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
968:           }
969:           break;
970:         case DM_POLYTOPE_QUADRILATERAL:
971:           for (q = 0; q < Nq; q++) {
972:             PetscReal xi[2], xio[2];
973:             PetscInt  oabs = (orient >= 0) ? orient : -(orient + 1);

975:             xi[0] = geom->xi[2 * q];
976:             xi[1] = geom->xi[2 * q + 1];
977:             switch (oabs) {
978:             case 1:
979:               xio[0] = xi[1];
980:               xio[1] = -xi[0];
981:               break;
982:             case 2:
983:               xio[0] = -xi[0];
984:               xio[1] = -xi[1];
985:               break;
986:             case 3:
987:               xio[0] = -xi[1];
988:               xio[1] = xi[0];
989:               break;
990:             case 0:
991:             default:
992:               xio[0] = xi[0];
993:               xio[1] = xi[1];
994:               break;
995:             }
996:             if (orient < 0) xio[0] = -xio[0];
997:             orientPoints[o][2 * q + 0] = xio[0];
998:             orientPoints[o][2 * q + 1] = xio[1];
999:           }
1000:           break;
1001:         default:
1002:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
1003:         }
1004:       }
1005:     }
1006:     for (f = 0; f < coneSize; f++) {
1007:       PetscInt  face = coneK[f];
1008:       PetscReal v0[3];
1009:       PetscReal J[9], detJ;
1010:       PetscInt  numCells, offset;
1011:       PetscInt *cells;
1012:       IS        suppIS;

1014:       PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ));
1015:       for (o = 0; o <= numOrient; o++) {
1016:         PetscFEGeom *cellGeom;

1018:         if (!counts[f][o]) continue;
1019:         /* If this (face,orientation) double appears,
1020:          * convert the face quadrature points into volume quadrature points */
1021:         for (q = 0; q < Nq; q++) {
1022:           PetscReal xi0[3] = {-1., -1., -1.};

1024:           CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1025:         }
1026:         for (p = 0, numCells = 0; p < numFaces; p++) {
1027:           for (s = 0; s < 2; s++) {
1028:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1029:           }
1030:         }
1031:         PetscCall(PetscMalloc1(numCells, &cells));
1032:         for (p = 0, offset = 0; p < numFaces; p++) {
1033:           for (s = 0; s < 2; s++) {
1034:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1035:           }
1036:         }
1037:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS));
1038:         PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom));
1039:         for (p = 0, offset = 0; p < numFaces; p++) {
1040:           for (s = 0; s < 2; s++) {
1041:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1042:               for (q = 0; q < Nq * dE * dE; q++) {
1043:                 geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
1044:                 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1045:               }
1046:               for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1047:               offset++;
1048:             }
1049:           }
1050:         }
1051:         PetscCall(PetscFEGeomDestroy(&cellGeom));
1052:         PetscCall(ISDestroy(&suppIS));
1053:         PetscCall(PetscFree(cells));
1054:       }
1055:     }
1056:     for (o = 0; o < numOrient; o++) {
1057:       if (orients[o]) PetscCall(PetscFree(orientPoints[o]));
1058:     }
1059:     PetscCall(PetscFree2(orients, orientPoints));
1060:     PetscCall(PetscQuadratureDestroy(&cellQuad));
1061:     for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f]));
1062:     PetscCall(PetscFree2(co, counts));
1063:   }
1064:   PetscCall(ISRestoreIndices(pointIS, &points));
1065:   PetscCall(ISDestroy(&cellIS));
1066:   PetscFunctionReturn(PETSC_SUCCESS);
1067: }

1069: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1070: {
1071:   PetscFunctionBegin;
1072:   field->ops->destroy                 = DMFieldDestroy_DS;
1073:   field->ops->evaluate                = DMFieldEvaluate_DS;
1074:   field->ops->evaluateFE              = DMFieldEvaluateFE_DS;
1075:   field->ops->evaluateFV              = DMFieldEvaluateFV_DS;
1076:   field->ops->getDegree               = DMFieldGetDegree_DS;
1077:   field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1078:   field->ops->view                    = DMFieldView_DS;
1079:   field->ops->computeFaceData         = DMFieldComputeFaceData_DS;
1080:   PetscFunctionReturn(PETSC_SUCCESS);
1081: }

1083: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1084: {
1085:   DMField_DS *dsfield;

1087:   PetscFunctionBegin;
1088:   PetscCall(PetscNew(&dsfield));
1089:   field->data = dsfield;
1090:   PetscCall(DMFieldInitialize_DS(field));
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

1094: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1095: {
1096:   DMField      b;
1097:   DMField_DS  *dsfield;
1098:   PetscObject  disc = NULL, discDG = NULL;
1099:   PetscSection section;
1100:   PetscBool    isContainer   = PETSC_FALSE;
1101:   PetscClassId id            = -1;
1102:   PetscInt     numComponents = -1, dsNumFields;

1104:   PetscFunctionBegin;
1105:   PetscCall(DMGetLocalSection(dm, &section));
1106:   PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents));
1107:   PetscCall(DMGetNumFields(dm, &dsNumFields));
1108:   if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc));
1109:   if (dsNumFields && dmDG) {
1110:     PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG));
1111:     PetscCall(PetscObjectReference(discDG));
1112:   }
1113:   if (disc) {
1114:     PetscCall(PetscObjectGetClassId(disc, &id));
1115:     isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1116:   }
1117:   if (!disc || isContainer) {
1118:     MPI_Comm       comm = PetscObjectComm((PetscObject)dm);
1119:     PetscFE        fe;
1120:     DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1121:     PetscInt       dim, cStart, cEnd, cellHeight;

1123:     PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1124:     PetscCall(DMGetDimension(dm, &dim));
1125:     PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1126:     if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct));
1127:     PetscCallMPI(MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm));
1128:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe));
1129:     PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view"));
1130:     disc = (PetscObject)fe;
1131:   } else PetscCall(PetscObjectReference(disc));
1132:   PetscCall(PetscObjectGetClassId(disc, &id));
1133:   if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents));
1134:   else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1135:   PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b));
1136:   PetscCall(DMFieldSetType(b, DMFIELDDS));
1137:   dsfield           = (DMField_DS *)b->data;
1138:   dsfield->fieldNum = fieldNum;
1139:   PetscCall(DMGetDimension(dm, &dsfield->height));
1140:   dsfield->height++;
1141:   PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc));
1142:   dsfield->disc[0] = disc;
1143:   PetscCall(PetscObjectReference((PetscObject)vec));
1144:   dsfield->vec = vec;
1145:   if (dmDG) {
1146:     dsfield->dmDG = dmDG;
1147:     PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG));
1148:     dsfield->discDG[0] = discDG;
1149:     PetscCall(PetscObjectReference((PetscObject)vecDG));
1150:     dsfield->vecDG = vecDG;
1151:   }
1152:   *field = b;
1153:   PetscFunctionReturn(PETSC_SUCCESS);
1154: }

1156: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1157: {
1158:   PetscFunctionBegin;
1159:   PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field));
1160:   PetscFunctionReturn(PETSC_SUCCESS);
1161: }