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:   VecDestroy(&dsfield->vec);
 24:   for (i = 0; i < dsfield->height; i++) PetscObjectDereference(dsfield->disc[i]);
 25:   PetscFree(dsfield->disc);
 26:   VecDestroy(&dsfield->vecDG);
 27:   if (dsfield->discDG)
 28:     for (i = 0; i < dsfield->height; i++) PetscObjectDereference(dsfield->discDG[i]);
 29:   PetscFree(dsfield->discDG);
 30:   PetscFree(dsfield);
 31:   return 0;
 32: }

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

 40:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 41:   disc = dsfield->disc[0];
 42:   if (iascii) {
 43:     PetscViewerASCIIPrintf(viewer, "PetscDS field %" PetscInt_FMT "\n", dsfield->fieldNum);
 44:     PetscViewerASCIIPushTab(viewer);
 45:     PetscObjectView(disc, viewer);
 46:     PetscViewerASCIIPopTab(viewer);
 47:   }
 48:   PetscViewerASCIIPushTab(viewer);
 50:   VecView(dsfield->vec, viewer);
 51:   PetscViewerASCIIPopTab(viewer);
 52:   return 0;
 53: }

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

 60:     PetscObjectGetClassId(discList[0], &id);
 61:     if (id == PETSCFE_CLASSID) PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height]);
 62:   }
 63:   *disc = discList[height];
 64:   return 0;
 65: }

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

 79: /*
 80:   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().
 81: */
 82: PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
 83: {
 84:   DMField_DS        *dsfield = (DMField_DS *)field->data;
 85:   DM                 fdm     = dsfield->dmDG;
 86:   PetscSection       s       = NULL;
 87:   const PetscScalar *cvalues;
 88:   PetscInt           pStart, pEnd;

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

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

130:   if (*isDG) {
131:     DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values);
132:   } else {
133:     DMFieldGetDM(field, &fdm);
134:     DMGetLocalSection(fdm, &s);
135:     DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, (PetscScalar **)values);
136:   }
137:   return 0;
138: }

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

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

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

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

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

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

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

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

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

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

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

247:   nc = field->numComponents;
248:   DMGetLocalSection(field->dm, &section);
249:   DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc);
250:   PetscObjectGetClassId(cellDisc, &discID);
252:   cellFE = (PetscFE)cellDisc;
253:   PetscFEGetDimension(cellFE, &feDim);
254:   DMGetCoordinateDim(field->dm, &dim);
255:   DMGetDimension(field->dm, &dimR);
256:   DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF);
257:   PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells);
259:   PetscSFComputeDegreeBegin(cellSF, &cellDegrees);
260:   PetscSFComputeDegreeEnd(cellSF, &cellDegrees);
261:   for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
262:     gatherMax = PetscMax(gatherMax, cellDegrees[c]);
263:     gatherSize += cellDegrees[c];
264:   }
265:   PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef);
266:   PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ);
267:   if (datatype == PETSC_SCALAR) PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs);
268:   else PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr);

270:   MPI_Type_contiguous(dim, MPIU_SCALAR, &pointType);
271:   MPI_Type_commit(&pointType);
272:   VecGetArrayRead(points, &pointsArray);
273:   PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints);
274:   PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints);
275:   VecRestoreArrayRead(points, &pointsArray);
276:   for (c = 0, offset = 0; c < numCells; c++) {
277:     PetscInt nq = cellDegrees[c], p;

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

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

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

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

314:           DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
315:           for (p = 0; p < nq; p++) {
316:             for (g = 0; g < nc; g++) {
317:               PetscScalar vs[3];

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

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

345:           DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
346:           for (p = 0; p < nq; p++) {
347:             for (g = 0; g < nc * dimR; g++) {
348:               PetscScalar vs[3];

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

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

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

400:     if (datatype == PETSC_SCALAR) {
401:       origtype = MPIU_SCALAR;
402:     } else {
403:       origtype = MPIU_REAL;
404:     }
405:     if (B) {
406:       MPI_Datatype Btype;

408:       MPI_Type_contiguous(nc, origtype, &Btype);
409:       MPI_Type_commit(&Btype);
410:       PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B);
411:       PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B);
412:       MPI_Type_free(&Btype);
413:     }
414:     if (D) {
415:       MPI_Datatype Dtype;

417:       MPI_Type_contiguous(nc * dim, origtype, &Dtype);
418:       MPI_Type_commit(&Dtype);
419:       PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D);
420:       PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D);
421:       MPI_Type_free(&Dtype);
422:     }
423:     if (H) {
424:       MPI_Datatype Htype;

426:       MPI_Type_contiguous(nc * dim * dim, origtype, &Htype);
427:       MPI_Type_commit(&Htype);
428:       PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H);
429:       PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H);
430:       MPI_Type_free(&Htype);
431:     }
432:   }
433:   PetscFree4(v, J, invJ, detJ);
434:   PetscFree3(cellBr, cellDr, cellHr);
435:   PetscFree3(cellBs, cellDs, cellHs);
436:   PetscFree3(cellPoints, coordsReal, coordsRef);
437:   MPI_Type_free(&pointType);
438:   PetscSFDestroy(&cellSF);
439:   return 0;
440: }

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

459:   Nc = field->numComponents;
460:   DMGetCoordinateDim(field->dm, &dimC);
461:   DMGetDimension(field->dm, &dim);
462:   ISGetLocalSize(pointIS, &numPoints);
463:   ISGetMinMax(pointIS, &imin, NULL);
464:   for (h = 0; h < dsfield->height; h++) {
465:     PetscInt hEnd;

467:     DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd);
468:     if (imin < hEnd) break;
469:   }
470:   dim -= h;
471:   DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
472:   PetscObjectGetClassId(disc, &id);
474:   DMGetCoordinateField(field->dm, &coordField);
475:   DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
476:   if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
477:   if (!quad) DMFieldCreateDefaultQuadrature(field, pointIS, &quad);
479:   DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
480:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights);
482:   N = numPoints * Nq * Nc;
483:   if (B) DMGetWorkArray(field->dm, N, mpitype, &qB);
484:   if (D) DMGetWorkArray(field->dm, N * dimC, mpitype, &qD);
485:   if (H) DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
486:   DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH);
487:   if (B) {
488:     PetscInt i, j, k;

490:     if (type == PETSC_SCALAR) {
491:       PetscScalar *sB  = (PetscScalar *)B;
492:       PetscScalar *sqB = (PetscScalar *)qB;

494:       for (i = 0; i < numPoints; i++) {
495:         PetscReal vol = 0.;

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

508:       for (i = 0; i < numPoints; i++) {
509:         PetscReal vol = 0.;

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

523:     if (type == PETSC_SCALAR) {
524:       PetscScalar *sD  = (PetscScalar *)D;
525:       PetscScalar *sqD = (PetscScalar *)qD;

527:       for (i = 0; i < numPoints; i++) {
528:         PetscReal vol = 0.;

530:         for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
531:         for (k = 0; k < Nq; k++) {
532:           vol += geom->detJ[i * Nq + k] * weights[k];
533:           for (j = 0; j < Nc; j++) {
534:             PetscScalar pD[3] = {0., 0., 0.};

536:             for (l = 0; l < dimC; l++) {
537:               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];
538:             }
539:             for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
540:           }
541:         }
542:         for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
543:       }
544:     } else {
545:       PetscReal *rD  = (PetscReal *)D;
546:       PetscReal *rqD = (PetscReal *)qD;

548:       for (i = 0; i < numPoints; i++) {
549:         PetscReal vol = 0.;

551:         for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
552:         for (k = 0; k < Nq; k++) {
553:           vol += geom->detJ[i * Nq + k] * weights[k];
554:           for (j = 0; j < Nc; j++) {
555:             PetscReal pD[3] = {0., 0., 0.};

557:             for (l = 0; l < dimC; l++) {
558:               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];
559:             }
560:             for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
561:           }
562:         }
563:         for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
564:       }
565:     }
566:   }
567:   if (H) {
568:     PetscInt i, j, k, l, m, q, r;

570:     if (type == PETSC_SCALAR) {
571:       PetscScalar *sH  = (PetscScalar *)H;
572:       PetscScalar *sqH = (PetscScalar *)qH;

574:       for (i = 0; i < numPoints; i++) {
575:         PetscReal vol = 0.;

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

581:           vol += geom->detJ[i * Nq + k] * weights[k];
582:           for (j = 0; j < Nc; j++) {
583:             PetscScalar pH[3][3] = {
584:               {0., 0., 0.},
585:               {0., 0., 0.},
586:               {0., 0., 0.}
587:             };
588:             const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];

590:             for (l = 0; l < dimC; l++) {
591:               for (m = 0; m < dimC; m++) {
592:                 for (q = 0; q < dim; q++) {
593:                   for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
594:                 }
595:               }
596:             }
597:             for (l = 0; l < dimC; l++) {
598:               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];
599:             }
600:           }
601:         }
602:         for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
603:       }
604:     } else {
605:       PetscReal *rH  = (PetscReal *)H;
606:       PetscReal *rqH = (PetscReal *)qH;

608:       for (i = 0; i < numPoints; i++) {
609:         PetscReal vol = 0.;

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

615:           vol += geom->detJ[i * Nq + k] * weights[k];
616:           for (j = 0; j < Nc; j++) {
617:             PetscReal pH[3][3] = {
618:               {0., 0., 0.},
619:               {0., 0., 0.},
620:               {0., 0., 0.}
621:             };
622:             const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];

624:             for (l = 0; l < dimC; l++) {
625:               for (m = 0; m < dimC; m++) {
626:                 for (q = 0; q < dim; q++) {
627:                   for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
628:                 }
629:               }
630:             }
631:             for (l = 0; l < dimC; l++) {
632:               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];
633:             }
634:           }
635:         }
636:         for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
637:       }
638:     }
639:   }
640:   if (B) DMRestoreWorkArray(field->dm, N, mpitype, &qB);
641:   if (D) DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD);
642:   if (H) DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
643:   PetscFEGeomDestroy(&geom);
644:   PetscQuadratureDestroy(&quad);
645:   return 0;
646: }

648: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
649: {
650:   DMField_DS  *dsfield;
651:   PetscObject  disc;
652:   PetscInt     h, imin, imax;
653:   PetscClassId id;

655:   dsfield = (DMField_DS *)field->data;
656:   ISGetMinMax(pointIS, &imin, &imax);
657:   if (imin >= imax) {
658:     h = 0;
659:   } else {
660:     for (h = 0; h < dsfield->height; h++) {
661:       PetscInt hEnd;

663:       DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd);
664:       if (imin < hEnd) break;
665:     }
666:   }
667:   DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
668:   PetscObjectGetClassId(disc, &id);
669:   if (id == PETSCFE_CLASSID) {
670:     PetscFE    fe = (PetscFE)disc;
671:     PetscSpace sp;

673:     PetscFEGetBasisSpace(fe, &sp);
674:     PetscSpaceGetDegree(sp, minDegree, maxDegree);
675:   }
676:   return 0;
677: }

679: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
680: {
681:   DM              dm = field->dm;
682:   const PetscInt *points;
683:   DMPolytopeType  ct;
684:   PetscInt        dim, n;
685:   PetscBool       isplex;

687:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
688:   ISGetLocalSize(pointIS, &n);
689:   if (isplex && n) {
690:     DMGetDimension(dm, &dim);
691:     ISGetIndices(pointIS, &points);
692:     DMPlexGetCellType(dm, points[0], &ct);
693:     switch (ct) {
694:     case DM_POLYTOPE_TRIANGLE:
695:     case DM_POLYTOPE_TETRAHEDRON:
696:       PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad);
697:       break;
698:     default:
699:       PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);
700:     }
701:     ISRestoreIndices(pointIS, &points);
702:   } else DMFieldCreateDefaultQuadrature(field, pointIS, quad);
703:   return 0;
704: }

706: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
707: {
708:   PetscInt     h, dim, imax, imin, cellHeight;
709:   DM           dm;
710:   DMField_DS  *dsfield;
711:   PetscObject  disc;
712:   PetscFE      fe;
713:   PetscClassId id;

715:   dm      = field->dm;
716:   dsfield = (DMField_DS *)field->data;
717:   ISGetMinMax(pointIS, &imin, &imax);
718:   DMGetDimension(dm, &dim);
719:   for (h = 0; h <= dim; h++) {
720:     PetscInt hStart, hEnd;

722:     DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
723:     if (imax >= hStart && imin < hEnd) break;
724:   }
725:   DMPlexGetVTKCellHeight(dm, &cellHeight);
726:   h -= cellHeight;
727:   *quad = NULL;
728:   if (h < dsfield->height) {
729:     DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
730:     PetscObjectGetClassId(disc, &id);
731:     if (id != PETSCFE_CLASSID) return 0;
732:     fe = (PetscFE)disc;
733:     PetscFEGetQuadrature(fe, quad);
734:     PetscObjectReference((PetscObject)*quad);
735:   }
736:   return 0;
737: }

739: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
740: {
741:   const PetscInt *points;
742:   PetscInt        p, dim, dE, numFaces, Nq;
743:   PetscInt        maxDegree;
744:   DMLabel         depthLabel;
745:   IS              cellIS;
746:   DM              dm = field->dm;

748:   dim = geom->dim;
749:   dE  = geom->dimEmbed;
750:   DMPlexGetDepthLabel(dm, &depthLabel);
751:   DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);
752:   DMFieldGetDegree(field, cellIS, NULL, &maxDegree);
753:   ISGetIndices(pointIS, &points);
754:   numFaces = geom->numCells;
755:   Nq       = geom->numPoints;
756:   /* First, set local faces and flip normals so that they are outward for the first supporting cell */
757:   for (p = 0; p < numFaces; p++) {
758:     PetscInt        point = points[p];
759:     PetscInt        suppSize, s, coneSize, c, numChildren;
760:     const PetscInt *supp, *cone, *ornt;

762:     DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
764:     DMPlexGetSupportSize(dm, point, &suppSize);
766:     if (!suppSize) continue;
767:     DMPlexGetSupport(dm, point, &supp);
768:     for (s = 0; s < suppSize; ++s) {
769:       DMPlexGetConeSize(dm, supp[s], &coneSize);
770:       DMPlexGetCone(dm, supp[s], &cone);
771:       for (c = 0; c < coneSize; ++c)
772:         if (cone[c] == point) break;
774:       geom->face[p][s] = c;
775:     }
776:     DMPlexGetConeOrientation(dm, supp[0], &ornt);
777:     if (ornt[geom->face[p][0]] < 0) {
778:       PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;

780:       for (q = 0; q < Np; ++q)
781:         for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
782:     }
783:   }
784:   if (maxDegree <= 1) {
785:     PetscInt     numCells, offset, *cells;
786:     PetscFEGeom *cellGeom;
787:     IS           suppIS;

789:     for (p = 0, numCells = 0; p < numFaces; p++) {
790:       PetscInt point = points[p];
791:       PetscInt numSupp, numChildren;

793:       DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
795:       DMPlexGetSupportSize(dm, point, &numSupp);
797:       numCells += numSupp;
798:     }
799:     PetscMalloc1(numCells, &cells);
800:     for (p = 0, offset = 0; p < numFaces; p++) {
801:       PetscInt        point = points[p];
802:       PetscInt        numSupp, s;
803:       const PetscInt *supp;

805:       DMPlexGetSupportSize(dm, point, &numSupp);
806:       DMPlexGetSupport(dm, point, &supp);
807:       for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
808:     }
809:     ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS);
810:     DMFieldCreateFEGeom(field, suppIS, quad, PETSC_FALSE, &cellGeom);
811:     for (p = 0, offset = 0; p < numFaces; p++) {
812:       PetscInt        point = points[p];
813:       PetscInt        numSupp, s, q;
814:       const PetscInt *supp;

816:       DMPlexGetSupportSize(dm, point, &numSupp);
817:       DMPlexGetSupport(dm, point, &supp);
818:       for (s = 0; s < numSupp; s++, offset++) {
819:         for (q = 0; q < Nq * dE * dE; q++) {
820:           geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
821:           geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
822:         }
823:         for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
824:       }
825:     }
826:     PetscFEGeomDestroy(&cellGeom);
827:     ISDestroy(&suppIS);
828:     PetscFree(cells);
829:   } else {
830:     DMField_DS    *dsfield = (DMField_DS *)field->data;
831:     PetscObject    faceDisc, cellDisc;
832:     PetscClassId   faceId, cellId;
833:     PetscDualSpace dsp;
834:     DM             K;
835:     DMPolytopeType ct;
836:     PetscInt(*co)[2][3];
837:     PetscInt        coneSize;
838:     PetscInt      **counts;
839:     PetscInt        f, i, o, q, s;
840:     const PetscInt *coneK;
841:     PetscInt        eStart, minOrient, maxOrient, numOrient;
842:     PetscInt       *orients;
843:     PetscReal     **orientPoints;
844:     PetscReal      *cellPoints;
845:     PetscReal      *dummyWeights;
846:     PetscQuadrature cellQuad = NULL;

848:     DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc);
849:     DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc);
850:     PetscObjectGetClassId(faceDisc, &faceId);
851:     PetscObjectGetClassId(cellDisc, &cellId);
853:     PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);
854:     PetscDualSpaceGetDM(dsp, &K);
855:     DMPlexGetHeightStratum(K, 1, &eStart, NULL);
856:     DMPlexGetCellType(K, eStart, &ct);
857:     DMPlexGetConeSize(K, 0, &coneSize);
858:     DMPlexGetCone(K, 0, &coneK);
859:     PetscMalloc2(numFaces, &co, coneSize, &counts);
860:     PetscMalloc1(dE * Nq, &cellPoints);
861:     PetscMalloc1(Nq, &dummyWeights);
862:     PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad);
863:     PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);
864:     minOrient = PETSC_MAX_INT;
865:     maxOrient = PETSC_MIN_INT;
866:     for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
867:       PetscInt        point = points[p];
868:       PetscInt        numSupp, numChildren;
869:       const PetscInt *supp;

871:       DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
873:       DMPlexGetSupportSize(dm, point, &numSupp);
875:       DMPlexGetSupport(dm, point, &supp);
876:       for (s = 0; s < numSupp; s++) {
877:         PetscInt        cell = supp[s];
878:         PetscInt        numCone;
879:         const PetscInt *cone, *orient;

881:         DMPlexGetConeSize(dm, cell, &numCone);
883:         DMPlexGetCone(dm, cell, &cone);
884:         DMPlexGetConeOrientation(dm, cell, &orient);
885:         for (f = 0; f < coneSize; f++) {
886:           if (cone[f] == point) break;
887:         }
888:         co[p][s][0] = f;
889:         co[p][s][1] = orient[f];
890:         co[p][s][2] = cell;
891:         minOrient   = PetscMin(minOrient, orient[f]);
892:         maxOrient   = PetscMax(maxOrient, orient[f]);
893:       }
894:       for (; s < 2; s++) {
895:         co[p][s][0] = -1;
896:         co[p][s][1] = -1;
897:         co[p][s][2] = -1;
898:       }
899:     }
900:     numOrient = maxOrient + 1 - minOrient;
901:     DMPlexGetCone(K, 0, &coneK);
902:     /* count all (face,orientation) doubles that appear */
903:     PetscCalloc2(numOrient, &orients, numOrient, &orientPoints);
904:     for (f = 0; f < coneSize; f++) PetscCalloc1(numOrient + 1, &counts[f]);
905:     for (p = 0; p < numFaces; p++) {
906:       for (s = 0; s < 2; s++) {
907:         if (co[p][s][0] >= 0) {
908:           counts[co[p][s][0]][co[p][s][1] - minOrient]++;
909:           orients[co[p][s][1] - minOrient]++;
910:         }
911:       }
912:     }
913:     for (o = 0; o < numOrient; o++) {
914:       if (orients[o]) {
915:         PetscInt orient = o + minOrient;
916:         PetscInt q;

918:         PetscMalloc1(Nq * dim, &orientPoints[o]);
919:         /* rotate the quadrature points appropriately */
920:         switch (ct) {
921:         case DM_POLYTOPE_POINT:
922:           break;
923:         case DM_POLYTOPE_SEGMENT:
924:           if (orient == -2 || orient == 1) {
925:             for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
926:           } else {
927:             for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
928:           }
929:           break;
930:         case DM_POLYTOPE_TRIANGLE:
931:           for (q = 0; q < Nq; q++) {
932:             PetscReal lambda[3];
933:             PetscReal lambdao[3];

935:             /* convert to barycentric */
936:             lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
937:             lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
938:             lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
939:             if (orient >= 0) {
940:               for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
941:             } else {
942:               for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
943:             }
944:             /* convert to coordinates */
945:             orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
946:             orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
947:           }
948:           break;
949:         case DM_POLYTOPE_QUADRILATERAL:
950:           for (q = 0; q < Nq; q++) {
951:             PetscReal xi[2], xio[2];
952:             PetscInt  oabs = (orient >= 0) ? orient : -(orient + 1);

954:             xi[0] = geom->xi[2 * q];
955:             xi[1] = geom->xi[2 * q + 1];
956:             switch (oabs) {
957:             case 1:
958:               xio[0] = xi[1];
959:               xio[1] = -xi[0];
960:               break;
961:             case 2:
962:               xio[0] = -xi[0];
963:               xio[1] = -xi[1];
964:               break;
965:             case 3:
966:               xio[0] = -xi[1];
967:               xio[1] = xi[0];
968:               break;
969:             case 0:
970:             default:
971:               xio[0] = xi[0];
972:               xio[1] = xi[1];
973:               break;
974:             }
975:             if (orient < 0) xio[0] = -xio[0];
976:             orientPoints[o][2 * q + 0] = xio[0];
977:             orientPoints[o][2 * q + 1] = xio[1];
978:           }
979:           break;
980:         default:
981:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
982:         }
983:       }
984:     }
985:     for (f = 0; f < coneSize; f++) {
986:       PetscInt  face = coneK[f];
987:       PetscReal v0[3];
988:       PetscReal J[9], detJ;
989:       PetscInt  numCells, offset;
990:       PetscInt *cells;
991:       IS        suppIS;

993:       DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ);
994:       for (o = 0; o <= numOrient; o++) {
995:         PetscFEGeom *cellGeom;

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

1003:           CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1004:         }
1005:         for (p = 0, numCells = 0; p < numFaces; p++) {
1006:           for (s = 0; s < 2; s++) {
1007:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1008:           }
1009:         }
1010:         PetscMalloc1(numCells, &cells);
1011:         for (p = 0, offset = 0; p < numFaces; p++) {
1012:           for (s = 0; s < 2; s++) {
1013:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1014:           }
1015:         }
1016:         ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS);
1017:         DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom);
1018:         for (p = 0, offset = 0; p < numFaces; p++) {
1019:           for (s = 0; s < 2; s++) {
1020:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1021:               for (q = 0; q < Nq * dE * dE; q++) {
1022:                 geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
1023:                 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1024:               }
1025:               for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1026:               offset++;
1027:             }
1028:           }
1029:         }
1030:         PetscFEGeomDestroy(&cellGeom);
1031:         ISDestroy(&suppIS);
1032:         PetscFree(cells);
1033:       }
1034:     }
1035:     for (o = 0; o < numOrient; o++) {
1036:       if (orients[o]) PetscFree(orientPoints[o]);
1037:     }
1038:     PetscFree2(orients, orientPoints);
1039:     PetscQuadratureDestroy(&cellQuad);
1040:     for (f = 0; f < coneSize; f++) PetscFree(counts[f]);
1041:     PetscFree2(co, counts);
1042:   }
1043:   ISRestoreIndices(pointIS, &points);
1044:   ISDestroy(&cellIS);
1045:   return 0;
1046: }

1048: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1049: {
1050:   field->ops->destroy                 = DMFieldDestroy_DS;
1051:   field->ops->evaluate                = DMFieldEvaluate_DS;
1052:   field->ops->evaluateFE              = DMFieldEvaluateFE_DS;
1053:   field->ops->evaluateFV              = DMFieldEvaluateFV_DS;
1054:   field->ops->getDegree               = DMFieldGetDegree_DS;
1055:   field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1056:   field->ops->view                    = DMFieldView_DS;
1057:   field->ops->computeFaceData         = DMFieldComputeFaceData_DS;
1058:   return 0;
1059: }

1061: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1062: {
1063:   DMField_DS *dsfield;

1065:   PetscNew(&dsfield);
1066:   field->data = dsfield;
1067:   DMFieldInitialize_DS(field);
1068:   return 0;
1069: }

1071: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1072: {
1073:   DMField      b;
1074:   DMField_DS  *dsfield;
1075:   PetscObject  disc = NULL, discDG = NULL;
1076:   PetscSection section;
1077:   PetscBool    isContainer   = PETSC_FALSE;
1078:   PetscClassId id            = -1;
1079:   PetscInt     numComponents = -1, dsNumFields;

1081:   DMGetLocalSection(dm, &section);
1082:   PetscSectionGetFieldComponents(section, fieldNum, &numComponents);
1083:   DMGetNumFields(dm, &dsNumFields);
1084:   if (dsNumFields) DMGetField(dm, fieldNum, NULL, &disc);
1085:   if (dsNumFields && dmDG) {
1086:     DMGetField(dmDG, fieldNum, NULL, &discDG);
1087:     PetscObjectReference(discDG);
1088:   }
1089:   if (disc) {
1090:     PetscObjectGetClassId(disc, &id);
1091:     isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1092:   }
1093:   if (!disc || isContainer) {
1094:     MPI_Comm       comm = PetscObjectComm((PetscObject)dm);
1095:     PetscFE        fe;
1096:     DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1097:     PetscInt       dim, cStart, cEnd, cellHeight;

1099:     DMPlexGetVTKCellHeight(dm, &cellHeight);
1100:     DMGetDimension(dm, &dim);
1101:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1102:     if (cEnd > cStart) DMPlexGetCellType(dm, cStart, &locct);
1103:     MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm);
1104:     PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe);
1105:     PetscFEViewFromOptions(fe, NULL, "-field_fe_view");
1106:     disc = (PetscObject)fe;
1107:   } else PetscObjectReference(disc);
1108:   PetscObjectGetClassId(disc, &id);
1109:   if (id == PETSCFE_CLASSID) PetscFEGetNumComponents((PetscFE)disc, &numComponents);
1110:   else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1111:   DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b);
1112:   DMFieldSetType(b, DMFIELDDS);
1113:   dsfield           = (DMField_DS *)b->data;
1114:   dsfield->fieldNum = fieldNum;
1115:   DMGetDimension(dm, &dsfield->height);
1116:   dsfield->height++;
1117:   PetscCalloc1(dsfield->height, &dsfield->disc);
1118:   dsfield->disc[0] = disc;
1119:   PetscObjectReference((PetscObject)vec);
1120:   dsfield->vec = vec;
1121:   if (dmDG) {
1122:     dsfield->dmDG = dmDG;
1123:     PetscCalloc1(dsfield->height, &dsfield->discDG);
1124:     dsfield->discDG[0] = discDG;
1125:     PetscObjectReference((PetscObject)vecDG);
1126:     dsfield->vecDG = vecDG;
1127:   }
1128:   *field = b;
1129:   return 0;
1130: }

1132: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1133: {
1134:   DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field);
1135:   return 0;
1136: }