Actual source code: fe.c

  1: /* Basis Jet Tabulation

3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6: as a prime basis.

8:   \psi_i = \sum_k \alpha_{ki} \phi_k

10: Our nodal basis is defined in terms of the dual basis $n_j$

12:   n_j \cdot \psi_i = \delta_{ji}

14: and we may act on the first equation to obtain

16:   n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17:        \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18:                  I = V \alpha

20: so the coefficients of the nodal basis in the prime basis are

22:    \alpha = V^{-1}

24: We will define the dual basis vectors $n_j$ using a quadrature rule.

26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28: be implemented exactly as in FIAT using functionals $L_j$.

30: I will have to count the degrees correctly for the Legendre product when we are on simplices.

32: We will have three objects:
33:  - Space, P: this just need point evaluation I think
34:  - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35:  - FEM: This keeps {P, P', Q}
36: */
37: #include <petsc/private/petscfeimpl.h>
38: #include <petscdmplex.h>

40: PetscBool FEcite = PETSC_FALSE;
41: const char FECitation[] = "@article{kirby2004,\n"
42:                           "  title   = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
43:                           "  journal = {ACM Transactions on Mathematical Software},\n"
44:                           "  author  = {Robert C. Kirby},\n"
45:                           "  volume  = {30},\n"
46:                           "  number  = {4},\n"
47:                           "  pages   = {502--516},\n"
48:                           "  doi     = {10.1145/1039813.1039820},\n"
49:                           "  year    = {2004}\n}\n";

51: PetscClassId PETSCFE_CLASSID = 0;

53: PetscLogEvent PETSCFE_SetUp;

55: PetscFunctionList PetscFEList              = NULL;
56: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

58: /*@C
59:   PetscFERegister - Adds a new PetscFE implementation

61:   Not Collective

63:   Input Parameters:
64: + name        - The name of a new user-defined creation routine
65: - create_func - The creation routine itself

67:   Notes:
68:   PetscFERegister() may be called multiple times to add several user-defined PetscFEs

70:   Sample usage:
71: .vb
72:     PetscFERegister("my_fe", MyPetscFECreate);
73: .ve

75:   Then, your PetscFE type can be chosen with the procedural interface via
76: .vb
77:     PetscFECreate(MPI_Comm, PetscFE *);
78:     PetscFESetType(PetscFE, "my_fe");
79: .ve
80:    or at runtime via the option
81: .vb
82:     -petscfe_type my_fe
83: .ve

87: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

89: @*/
90: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
91: {

96:   return(0);
97: }

99: /*@C
100:   PetscFESetType - Builds a particular PetscFE

102:   Collective on fem

104:   Input Parameters:
105: + fem  - The PetscFE object
106: - name - The kind of FEM space

108:   Options Database Key:
109: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types

111:   Level: intermediate

113: .seealso: PetscFEGetType(), PetscFECreate()
114: @*/
115: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
116: {
117:   PetscErrorCode (*r)(PetscFE);
118:   PetscBool      match;

123:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
124:   if (match) return(0);

126:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
127:   PetscFunctionListFind(PetscFEList, name, &r);
128:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);

130:   if (fem->ops->destroy) {
131:     (*fem->ops->destroy)(fem);
132:     fem->ops->destroy = NULL;
133:   }
134:   (*r)(fem);
135:   PetscObjectChangeTypeName((PetscObject) fem, name);
136:   return(0);
137: }

139: /*@C
140:   PetscFEGetType - Gets the PetscFE type name (as a string) from the object.

142:   Not Collective

144:   Input Parameter:
145: . fem  - The PetscFE

147:   Output Parameter:
148: . name - The PetscFE type name

150:   Level: intermediate

152: .seealso: PetscFESetType(), PetscFECreate()
153: @*/
154: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
155: {

161:   if (!PetscFERegisterAllCalled) {
162:     PetscFERegisterAll();
163:   }
164:   *name = ((PetscObject) fem)->type_name;
165:   return(0);
166: }

168: /*@C
169:    PetscFEViewFromOptions - View from Options

171:    Collective on PetscFE

173:    Input Parameters:
174: +  A - the PetscFE object
175: .  obj - Optional object
176: -  name - command line option

178:    Level: intermediate
179: .seealso:  PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
180: @*/
181: PetscErrorCode  PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
182: {

187:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
188:   return(0);
189: }

191: /*@C
192:   PetscFEView - Views a PetscFE

194:   Collective on fem

196:   Input Parameters:
197: + fem - the PetscFE object to view
198: - viewer   - the viewer

200:   Level: beginner

202: .seealso PetscFEDestroy()
203: @*/
204: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
205: {
206:   PetscBool      iascii;

212:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
213:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
214:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
215:   if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
216:   return(0);
217: }

219: /*@
220:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

222:   Collective on fem

224:   Input Parameter:
225: . fem - the PetscFE object to set options for

227:   Options Database:
228: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
229: - -petscfe_num_batches - the number of cell batches to integrate serially

231:   Level: intermediate

233: .seealso PetscFEView()
234: @*/
235: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
236: {
237:   const char    *defaultType;
238:   char           name[256];
239:   PetscBool      flg;

244:   if (!((PetscObject) fem)->type_name) {
245:     defaultType = PETSCFEBASIC;
246:   } else {
247:     defaultType = ((PetscObject) fem)->type_name;
248:   }
249:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

251:   PetscObjectOptionsBegin((PetscObject) fem);
252:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
253:   if (flg) {
254:     PetscFESetType(fem, name);
255:   } else if (!((PetscObject) fem)->type_name) {
256:     PetscFESetType(fem, defaultType);
257:   }
258:   PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
259:   PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
260:   if (fem->ops->setfromoptions) {
261:     (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
262:   }
264:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
265:   PetscOptionsEnd();
266:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
267:   return(0);
268: }

270: /*@C
271:   PetscFESetUp - Construct data structures for the PetscFE

273:   Collective on fem

275:   Input Parameter:
276: . fem - the PetscFE object to setup

278:   Level: intermediate

280: .seealso PetscFEView(), PetscFEDestroy()
281: @*/
282: PetscErrorCode PetscFESetUp(PetscFE fem)
283: {

288:   if (fem->setupcalled) return(0);
289:   PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
290:   fem->setupcalled = PETSC_TRUE;
291:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
292:   PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
293:   return(0);
294: }

296: /*@
297:   PetscFEDestroy - Destroys a PetscFE object

299:   Collective on fem

301:   Input Parameter:
302: . fem - the PetscFE object to destroy

304:   Level: beginner

306: .seealso PetscFEView()
307: @*/
308: PetscErrorCode PetscFEDestroy(PetscFE *fem)
309: {

313:   if (!*fem) return(0);

316:   if (--((PetscObject)(*fem))->refct > 0) {*fem = NULL; return(0);}
317:   ((PetscObject) (*fem))->refct = 0;

319:   if ((*fem)->subspaces) {
320:     PetscInt dim, d;

322:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
323:     for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
324:   }
325:   PetscFree((*fem)->subspaces);
326:   PetscFree((*fem)->invV);
327:   PetscTabulationDestroy(&(*fem)->T);
328:   PetscTabulationDestroy(&(*fem)->Tf);
329:   PetscTabulationDestroy(&(*fem)->Tc);
330:   PetscSpaceDestroy(&(*fem)->basisSpace);
331:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
334: #ifdef PETSC_HAVE_LIBCEED
335:   CeedBasisDestroy(&(*fem)->ceedBasis);
336:   CeedDestroy(&(*fem)->ceed);
337: #endif

339:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
341:   return(0);
342: }

344: /*@
345:   PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().

347:   Collective

349:   Input Parameter:
350: . comm - The communicator for the PetscFE object

352:   Output Parameter:
353: . fem - The PetscFE object

355:   Level: beginner

357: .seealso: PetscFESetType(), PETSCFEGALERKIN
358: @*/
359: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
360: {
361:   PetscFE        f;

366:   PetscCitationsRegister(FECitation,&FEcite);
367:   *fem = NULL;
368:   PetscFEInitializePackage();

370:   PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);

372:   f->basisSpace    = NULL;
373:   f->dualSpace     = NULL;
374:   f->numComponents = 1;
375:   f->subspaces     = NULL;
376:   f->invV          = NULL;
377:   f->T             = NULL;
378:   f->Tf            = NULL;
379:   f->Tc            = NULL;
382:   f->blockSize     = 0;
383:   f->numBlocks     = 1;
384:   f->batchSize     = 0;
385:   f->numBatches    = 1;

387:   *fem = f;
388:   return(0);
389: }

391: /*@
392:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

394:   Not collective

396:   Input Parameter:
397: . fem - The PetscFE object

399:   Output Parameter:
400: . dim - The spatial dimension

402:   Level: intermediate

404: .seealso: PetscFECreate()
405: @*/
406: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
407: {
408:   DM             dm;

414:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
415:   DMGetDimension(dm, dim);
416:   return(0);
417: }

419: /*@
420:   PetscFESetNumComponents - Sets the number of components in the element

422:   Not collective

424:   Input Parameters:
425: + fem - The PetscFE object
426: - comp - The number of field components

428:   Level: intermediate

430: .seealso: PetscFECreate()
431: @*/
432: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
433: {
436:   fem->numComponents = comp;
437:   return(0);
438: }

440: /*@
441:   PetscFEGetNumComponents - Returns the number of components in the element

443:   Not collective

445:   Input Parameter:
446: . fem - The PetscFE object

448:   Output Parameter:
449: . comp - The number of field components

451:   Level: intermediate

453: .seealso: PetscFECreate()
454: @*/
455: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
456: {
460:   *comp = fem->numComponents;
461:   return(0);
462: }

464: /*@
465:   PetscFESetTileSizes - Sets the tile sizes for evaluation

467:   Not collective

469:   Input Parameters:
470: + fem - The PetscFE object
471: . blockSize - The number of elements in a block
472: . numBlocks - The number of blocks in a batch
473: . batchSize - The number of elements in a batch
474: - numBatches - The number of batches in a chunk

476:   Level: intermediate

478: .seealso: PetscFECreate()
479: @*/
480: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
481: {
484:   fem->blockSize  = blockSize;
485:   fem->numBlocks  = numBlocks;
486:   fem->batchSize  = batchSize;
487:   fem->numBatches = numBatches;
488:   return(0);
489: }

491: /*@
492:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

494:   Not collective

496:   Input Parameter:
497: . fem - The PetscFE object

499:   Output Parameters:
500: + blockSize - The number of elements in a block
501: . numBlocks - The number of blocks in a batch
502: . batchSize - The number of elements in a batch
503: - numBatches - The number of batches in a chunk

505:   Level: intermediate

507: .seealso: PetscFECreate()
508: @*/
509: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
510: {
517:   if (blockSize)  *blockSize  = fem->blockSize;
518:   if (numBlocks)  *numBlocks  = fem->numBlocks;
519:   if (batchSize)  *batchSize  = fem->batchSize;
520:   if (numBatches) *numBatches = fem->numBatches;
521:   return(0);
522: }

524: /*@
525:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

527:   Not collective

529:   Input Parameter:
530: . fem - The PetscFE object

532:   Output Parameter:
533: . sp - The PetscSpace object

535:   Level: intermediate

537: .seealso: PetscFECreate()
538: @*/
539: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
540: {
544:   *sp = fem->basisSpace;
545:   return(0);
546: }

548: /*@
549:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

551:   Not collective

553:   Input Parameters:
554: + fem - The PetscFE object
555: - sp - The PetscSpace object

557:   Level: intermediate

559: .seealso: PetscFECreate()
560: @*/
561: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
562: {

568:   PetscSpaceDestroy(&fem->basisSpace);
569:   fem->basisSpace = sp;
570:   PetscObjectReference((PetscObject) fem->basisSpace);
571:   return(0);
572: }

574: /*@
575:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

577:   Not collective

579:   Input Parameter:
580: . fem - The PetscFE object

582:   Output Parameter:
583: . sp - The PetscDualSpace object

585:   Level: intermediate

587: .seealso: PetscFECreate()
588: @*/
589: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
590: {
594:   *sp = fem->dualSpace;
595:   return(0);
596: }

598: /*@
599:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

601:   Not collective

603:   Input Parameters:
604: + fem - The PetscFE object
605: - sp - The PetscDualSpace object

607:   Level: intermediate

609: .seealso: PetscFECreate()
610: @*/
611: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
612: {

618:   PetscDualSpaceDestroy(&fem->dualSpace);
619:   fem->dualSpace = sp;
620:   PetscObjectReference((PetscObject) fem->dualSpace);
621:   return(0);
622: }

624: /*@

627:   Not collective

629:   Input Parameter:
630: . fem - The PetscFE object

632:   Output Parameter:
633: . q - The PetscQuadrature object

635:   Level: intermediate

637: .seealso: PetscFECreate()
638: @*/
640: {
645:   return(0);
646: }

648: /*@

651:   Not collective

653:   Input Parameters:
654: + fem - The PetscFE object
655: - q - The PetscQuadrature object

657:   Level: intermediate

659: .seealso: PetscFECreate()
660: @*/
662: {
663:   PetscInt       Nc, qNc;

668:   if (q == fem->quadrature) return(0);
669:   PetscFEGetNumComponents(fem, &Nc);
671:   if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
672:   PetscTabulationDestroy(&fem->T);
673:   PetscTabulationDestroy(&fem->Tc);
674:   PetscObjectReference((PetscObject) q);
677:   return(0);
678: }

680: /*@
681:   PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces

683:   Not collective

685:   Input Parameter:
686: . fem - The PetscFE object

688:   Output Parameter:
689: . q - The PetscQuadrature object

691:   Level: intermediate

693: .seealso: PetscFECreate()
694: @*/
696: {
701:   return(0);
702: }

704: /*@
705:   PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces

707:   Not collective

709:   Input Parameters:
710: + fem - The PetscFE object
711: - q - The PetscQuadrature object

713:   Level: intermediate

715: .seealso: PetscFECreate()
716: @*/
718: {
719:   PetscInt       Nc, qNc;

724:   PetscFEGetNumComponents(fem, &Nc);
726:   if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
727:   PetscTabulationDestroy(&fem->Tf);
730:   PetscObjectReference((PetscObject) q);
731:   return(0);
732: }

734: /*@

737:   Not collective

739:   Input Parameters:
740: + sfe - The PetscFE source for the quadratures
741: - tfe - The PetscFE target for the quadratures

743:   Level: intermediate

746: @*/
747: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
748: {
750:   PetscErrorCode  ierr;

759:   return(0);
760: }

762: /*@C
763:   PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension

765:   Not collective

767:   Input Parameter:
768: . fem - The PetscFE object

770:   Output Parameter:
771: . numDof - Array with the number of dofs per dimension

773:   Level: intermediate

775: .seealso: PetscFECreate()
776: @*/
777: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
778: {

784:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
785:   return(0);
786: }

788: /*@C
789:   PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell

791:   Not collective

793:   Input Parameters:
794: + fem - The PetscFE object
795: - k   - The highest derivative we need to tabulate, very often 1

797:   Output Parameter:
798: . T - The basis function values and derivatives at quadrature points

800:   Note:
801: $T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 802:$ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
803: $T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e 805: Level: intermediate 807: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy() 808: @*/ 809: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T) 810: { 811: PetscInt npoints; 812: const PetscReal *points; 813: PetscErrorCode ierr; 818: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL); 819: if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);} 820: if (fem->T && k > fem->T->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->T->K); 821: *T = fem->T; 822: return(0); 823: } 825: /*@C 826: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell 828: Not collective 830: Input Parameters: 831: + fem - The PetscFE object 832: - k - The highest derivative we need to tabulate, very often 1 834: Output Parameters: 835: . Tf - The basis function values and derivatives at face quadrature points 837: Note: 838:$ T->T[0] = Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
839: $T->T[1] = Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d 840:$ T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e

842:   Level: intermediate

844: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
845: @*/
846: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
847: {
848:   PetscErrorCode   ierr;

853:   if (!fem->Tf) {
854:     const PetscReal  xi0[3] = {-1., -1., -1.};
855:     PetscReal        v0[3], J[9], detJ;
857:     PetscDualSpace   sp;
858:     DM               dm;
859:     const PetscInt  *faces;
860:     PetscInt         dim, numFaces, f, npoints, q;
861:     const PetscReal *points;
862:     PetscReal       *facePoints;

864:     PetscFEGetDualSpace(fem, &sp);
865:     PetscDualSpaceGetDM(sp, &dm);
866:     DMGetDimension(dm, &dim);
867:     DMPlexGetConeSize(dm, 0, &numFaces);
868:     DMPlexGetCone(dm, 0, &faces);
870:     if (fq) {
871:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
872:       PetscMalloc1(numFaces*npoints*dim, &facePoints);
873:       for (f = 0; f < numFaces; ++f) {
874:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
875:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
876:       }
877:       PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
878:       PetscFree(facePoints);
879:     }
880:   }
881:   if (fem->Tf && k > fem->Tf->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->Tf->K);
882:   *Tf = fem->Tf;
883:   return(0);
884: }

886: /*@C
887:   PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points

889:   Not collective

891:   Input Parameter:
892: . fem - The PetscFE object

894:   Output Parameters:
895: . Tc - The basis function values at face centroid points

897:   Note:
898: $T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c 900: Level: intermediate 902: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() 903: @*/ 904: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc) 905: { 906: PetscErrorCode ierr; 911: if (!fem->Tc) { 912: PetscDualSpace sp; 913: DM dm; 914: const PetscInt *cone; 915: PetscReal *centroids; 916: PetscInt dim, numFaces, f; 918: PetscFEGetDualSpace(fem, &sp); 919: PetscDualSpaceGetDM(sp, &dm); 920: DMGetDimension(dm, &dim); 921: DMPlexGetConeSize(dm, 0, &numFaces); 922: DMPlexGetCone(dm, 0, &cone); 923: PetscMalloc1(numFaces*dim, &centroids); 924: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);} 925: PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc); 926: PetscFree(centroids); 927: } 928: *Tc = fem->Tc; 929: return(0); 930: } 932: /*@C 933: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. 935: Not collective 937: Input Parameters: 938: + fem - The PetscFE object 939: . nrepl - The number of replicas 940: . npoints - The number of tabulation points in a replica 941: . points - The tabulation point coordinates 942: - K - The number of derivatives calculated 944: Output Parameter: 945: . T - The basis function values and derivatives at tabulation points 947: Note: 948:$ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
949: $T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d 950:$ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

952:   Level: intermediate

954: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
955: @*/
956: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
957: {
958:   DM               dm;
959:   PetscDualSpace   Q;
960:   PetscInt         Nb;   /* Dimension of FE space P */
961:   PetscInt         Nc;   /* Field components */
962:   PetscInt         cdim; /* Reference coordinate dimension */
963:   PetscInt         k;
964:   PetscErrorCode   ierr;

967:   if (!npoints || !fem->dualSpace || K < 0) {
968:     *T = NULL;
969:     return(0);
970:   }
974:   PetscFEGetDualSpace(fem, &Q);
975:   PetscDualSpaceGetDM(Q, &dm);
976:   DMGetDimension(dm, &cdim);
977:   PetscDualSpaceGetDimension(Q, &Nb);
978:   PetscFEGetNumComponents(fem, &Nc);
979:   PetscMalloc1(1, T);
980:   (*T)->K    = !cdim ? 0 : K;
981:   (*T)->Nr   = nrepl;
982:   (*T)->Np   = npoints;
983:   (*T)->Nb   = Nb;
984:   (*T)->Nc   = Nc;
985:   (*T)->cdim = cdim;
986:   PetscMalloc1((*T)->K+1, &(*T)->T);
987:   for (k = 0; k <= (*T)->K; ++k) {
988:     PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
989:   }
990:   (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
991:   return(0);
992: }

994: /*@C
995:   PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

997:   Not collective

999:   Input Parameters:
1000: + fem     - The PetscFE object
1001: . npoints - The number of tabulation points
1002: . points  - The tabulation point coordinates
1003: . K       - The number of derivatives calculated
1004: - T       - An existing tabulation object with enough allocated space

1006:   Output Parameter:
1007: . T - The basis function values and derivatives at tabulation points

1009:   Note:
1010: $T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 1011:$ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1012: $T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e 1014: Level: intermediate 1016: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() 1017: @*/ 1018: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 1019: { 1023: if (!npoints || !fem->dualSpace || K < 0) return(0); 1027: if (PetscDefined(USE_DEBUG)) { 1028: DM dm; 1029: PetscDualSpace Q; 1030: PetscInt Nb; /* Dimension of FE space P */ 1031: PetscInt Nc; /* Field components */ 1032: PetscInt cdim; /* Reference coordinate dimension */ 1034: PetscFEGetDualSpace(fem, &Q); 1035: PetscDualSpaceGetDM(Q, &dm); 1036: DMGetDimension(dm, &cdim); 1037: PetscDualSpaceGetDimension(Q, &Nb); 1038: PetscFEGetNumComponents(fem, &Nc); 1039: if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K); 1040: if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb); 1041: if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc); 1042: if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim); 1043: } 1044: T->Nr = 1; 1045: T->Np = npoints; 1046: (*fem->ops->createtabulation)(fem, npoints, points, K, T); 1047: return(0); 1048: } 1050: /*@C 1051: PetscTabulationDestroy - Frees memory from the associated tabulation. 1053: Not collective 1055: Input Parameter: 1056: . T - The tabulation 1058: Level: intermediate 1060: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation() 1061: @*/ 1062: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T) 1063: { 1064: PetscInt k; 1069: if (!T || !(*T)) return(0); 1070: for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);} 1071: PetscFree((*T)->T); 1072: PetscFree(*T); 1073: *T = NULL; 1074: return(0); 1075: } 1077: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE) 1078: { 1079: PetscSpace bsp, bsubsp; 1080: PetscDualSpace dsp, dsubsp; 1081: PetscInt dim, depth, numComp, i, j, coneSize, order; 1082: PetscFEType type; 1083: DM dm; 1084: DMLabel label; 1085: PetscReal *xi, *v, *J, detJ; 1086: const char *name; 1087: PetscQuadrature origin, fullQuad, subQuad; 1093: PetscFEGetBasisSpace(fe,&bsp); 1094: PetscFEGetDualSpace(fe,&dsp); 1095: PetscDualSpaceGetDM(dsp,&dm); 1096: DMGetDimension(dm,&dim); 1097: DMPlexGetDepthLabel(dm,&label); 1098: DMLabelGetValue(label,refPoint,&depth); 1099: PetscCalloc1(depth,&xi); 1100: PetscMalloc1(dim,&v); 1101: PetscMalloc1(dim*dim,&J); 1102: for (i = 0; i < depth; i++) xi[i] = 0.; 1103: PetscQuadratureCreate(PETSC_COMM_SELF,&origin); 1104: PetscQuadratureSetData(origin,depth,0,1,xi,NULL); 1105: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ); 1106: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */ 1107: for (i = 1; i < dim; i++) { 1108: for (j = 0; j < depth; j++) { 1109: J[i * depth + j] = J[i * dim + j]; 1110: } 1111: } 1112: PetscQuadratureDestroy(&origin); 1113: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp); 1114: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp); 1115: PetscSpaceSetUp(bsubsp); 1116: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE); 1117: PetscFEGetType(fe,&type); 1118: PetscFESetType(*trFE,type); 1119: PetscFEGetNumComponents(fe,&numComp); 1120: PetscFESetNumComponents(*trFE,numComp); 1121: PetscFESetBasisSpace(*trFE,bsubsp); 1122: PetscFESetDualSpace(*trFE,dsubsp); 1123: PetscObjectGetName((PetscObject) fe, &name); 1124: if (name) {PetscFESetName(*trFE, name);} 1125: PetscFEGetQuadrature(fe,&fullQuad); 1126: PetscQuadratureGetOrder(fullQuad,&order); 1127: DMPlexGetConeSize(dm,refPoint,&coneSize); 1128: if (coneSize == 2 * depth) { 1129: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad); 1130: } else { 1131: PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad); 1132: } 1133: PetscFESetQuadrature(*trFE,subQuad); 1134: PetscFESetUp(*trFE); 1135: PetscQuadratureDestroy(&subQuad); 1136: PetscSpaceDestroy(&bsubsp); 1137: return(0); 1138: } 1140: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE) 1141: { 1142: PetscInt hStart, hEnd; 1143: PetscDualSpace dsp; 1144: DM dm; 1150: *trFE = NULL; 1151: PetscFEGetDualSpace(fe,&dsp); 1152: PetscDualSpaceGetDM(dsp,&dm); 1153: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd); 1154: if (hEnd <= hStart) return(0); 1155: PetscFECreatePointTrace(fe,hStart,trFE); 1156: return(0); 1157: } 1159: /*@ 1160: PetscFEGetDimension - Get the dimension of the finite element space on a cell 1162: Not collective 1164: Input Parameter: 1165: . fe - The PetscFE 1167: Output Parameter: 1168: . dim - The dimension 1170: Level: intermediate 1172: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 1173: @*/ 1174: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 1175: { 1181: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);} 1182: return(0); 1183: } 1185: /*@C 1186: PetscFEPushforward - Map the reference element function to real space 1188: Input Parameters: 1189: + fe - The PetscFE 1190: . fegeom - The cell geometry 1191: . Nv - The number of function values 1192: - vals - The function values 1194: Output Parameter: 1195: . vals - The transformed function values 1197: Level: advanced 1199: Note: This just forwards the call onto PetscDualSpacePushforward(). 1201: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1203: .seealso: PetscDualSpacePushforward() 1204: @*/ 1205: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1206: { 1210: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals); 1211: return(0); 1212: } 1214: /*@C 1215: PetscFEPushforwardGradient - Map the reference element function gradient to real space 1217: Input Parameters: 1218: + fe - The PetscFE 1219: . fegeom - The cell geometry 1220: . Nv - The number of function gradient values 1221: - vals - The function gradient values 1223: Output Parameter: 1224: . vals - The transformed function gradient values 1226: Level: advanced 1228: Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). 1230: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1232: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() 1233: @*/ 1234: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1235: { 1239: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals); 1240: return(0); 1241: } 1243: /*@C 1244: PetscFEPushforwardHessian - Map the reference element function Hessian to real space 1246: Input Parameters: 1247: + fe - The PetscFE 1248: . fegeom - The cell geometry 1249: . Nv - The number of function Hessian values 1250: - vals - The function Hessian values 1252: Output Parameter: 1253: . vals - The transformed function Hessian values 1255: Level: advanced 1257: Note: This just forwards the call onto PetscDualSpacePushforwardHessian(). 1259: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1261: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward() 1262: @*/ 1263: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1264: { 1268: PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals); 1269: return(0); 1270: } 1272: /* 1273: Purpose: Compute element vector for chunk of elements 1275: Input: 1276: Sizes: 1277: Ne: number of elements 1278: Nf: number of fields 1279: PetscFE 1280: dim: spatial dimension 1281: Nb: number of basis functions 1282: Nc: number of field components 1283: PetscQuadrature 1284: Nq: number of quadrature points 1286: Geometry: 1287: PetscFEGeom[Ne] possibly *Nq 1288: PetscReal v0s[dim] 1289: PetscReal n[dim] 1290: PetscReal jacobians[dim*dim] 1291: PetscReal jacobianInverses[dim*dim] 1292: PetscReal jacobianDeterminants 1293: FEM: 1294: PetscFE 1295: PetscQuadrature 1296: PetscReal quadPoints[Nq*dim] 1297: PetscReal quadWeights[Nq] 1298: PetscReal basis[Nq*Nb*Nc] 1299: PetscReal basisDer[Nq*Nb*Nc*dim] 1300: PetscScalar coefficients[Ne*Nb*Nc] 1301: PetscScalar elemVec[Ne*Nb*Nc] 1303: Problem: 1304: PetscInt f: the active field 1305: f0, f1 1307: Work Space: 1308: PetscFE 1309: PetscScalar f0[Nq*dim]; 1310: PetscScalar f1[Nq*dim*dim]; 1311: PetscScalar u[Nc]; 1312: PetscScalar gradU[Nc*dim]; 1313: PetscReal x[dim]; 1314: PetscScalar realSpaceDer[dim]; 1316: Purpose: Compute element vector for N_cb batches of elements 1318: Input: 1319: Sizes: 1320: N_cb: Number of serial cell batches 1322: Geometry: 1323: PetscReal v0s[Ne*dim] 1324: PetscReal jacobians[Ne*dim*dim] possibly *Nq 1325: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 1326: PetscReal jacobianDeterminants[Ne] possibly *Nq 1327: FEM: 1328: static PetscReal quadPoints[Nq*dim] 1329: static PetscReal quadWeights[Nq] 1330: static PetscReal basis[Nq*Nb*Nc] 1331: static PetscReal basisDer[Nq*Nb*Nc*dim] 1332: PetscScalar coefficients[Ne*Nb*Nc] 1333: PetscScalar elemVec[Ne*Nb*Nc] 1335: ex62.c: 1336: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 1337: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 1338: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 1339: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 1341: ex52.c: 1342: PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1343: PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1345: ex52_integrateElement.cu 1346: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 1348: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 1349: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1350: PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1352: ex52_integrateElementOpenCL.c: 1353: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 1354: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1355: PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1357: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 1358: */ 1360: /*@C 1361: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration 1363: Not collective 1365: Input Parameters: 1366: + prob - The PetscDS specifying the discretizations and continuum functions 1367: . field - The field being integrated 1368: . Ne - The number of elements in the chunk 1369: . cgeom - The cell geometry for each cell in the chunk 1370: . coefficients - The array of FEM basis coefficients for the elements 1371: . probAux - The PetscDS specifying the auxiliary discretizations 1372: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1374: Output Parameter: 1375: . integral - the integral for this field 1377: Level: intermediate 1379: .seealso: PetscFEIntegrateResidual() 1380: @*/ 1381: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1382: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1383: { 1384: PetscFE fe; 1389: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe); 1390: if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);} 1391: return(0); 1392: } 1394: /*@C 1395: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration 1397: Not collective 1399: Input Parameters: 1400: + prob - The PetscDS specifying the discretizations and continuum functions 1401: . field - The field being integrated 1402: . obj_func - The function to be integrated 1403: . Ne - The number of elements in the chunk 1404: . fgeom - The face geometry for each face in the chunk 1405: . coefficients - The array of FEM basis coefficients for the elements 1406: . probAux - The PetscDS specifying the auxiliary discretizations 1407: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1409: Output Parameter: 1410: . integral - the integral for this field 1412: Level: intermediate 1414: .seealso: PetscFEIntegrateResidual() 1415: @*/ 1416: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, 1417: void (*obj_func)(PetscInt, PetscInt, PetscInt, 1418: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1419: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1420: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1421: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1422: { 1423: PetscFE fe; 1428: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe); 1429: if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);} 1430: return(0); 1431: } 1433: /*@C 1434: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 1436: Not collective 1438: Input Parameters: 1439: + ds - The PetscDS specifying the discretizations and continuum functions 1440: . key - The (label+value, field) being integrated 1441: . Ne - The number of elements in the chunk 1442: . cgeom - The cell geometry for each cell in the chunk 1443: . coefficients - The array of FEM basis coefficients for the elements 1444: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1445: . probAux - The PetscDS specifying the auxiliary discretizations 1446: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1447: - t - The time 1449: Output Parameter: 1450: . elemVec - the element residual vectors from each element 1452: Note: 1453:$ Loop over batch of elements (e):
1454: $Loop over quadrature points (q): 1455:$     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1456: $Call f_0 and f_1 1457:$   Loop over element vector entries (f,fc --> i):
1458: $elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 1460: Level: intermediate 1462: .seealso: PetscFEIntegrateResidual() 1463: @*/ 1464: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1465: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1466: { 1467: PetscFE fe; 1472: PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe); 1473: if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);} 1474: return(0); 1475: } 1477: /*@C 1478: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary 1480: Not collective 1482: Input Parameters: 1483: + ds - The PetscDS specifying the discretizations and continuum functions 1484: . wf - The PetscWeakForm object holding the pointwise functions 1485: . key - The (label+value, field) being integrated 1486: . Ne - The number of elements in the chunk 1487: . fgeom - The face geometry for each cell in the chunk 1488: . coefficients - The array of FEM basis coefficients for the elements 1489: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1490: . probAux - The PetscDS specifying the auxiliary discretizations 1491: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1492: - t - The time 1494: Output Parameter: 1495: . elemVec - the element residual vectors from each element 1497: Level: intermediate 1499: .seealso: PetscFEIntegrateResidual() 1500: @*/ 1501: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1502: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1503: { 1504: PetscFE fe; 1509: PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe); 1510: if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);} 1511: return(0); 1512: } 1514: /*@C 1515: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1517: Not collective 1519: Input Parameters: 1520: + prob - The PetscDS specifying the discretizations and continuum functions 1521: . key - The (label+value, field) being integrated 1522: . s - The side of the cell being integrated, 0 for negative and 1 for positive 1523: . Ne - The number of elements in the chunk 1524: . fgeom - The face geometry for each cell in the chunk 1525: . coefficients - The array of FEM basis coefficients for the elements 1526: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1527: . probAux - The PetscDS specifying the auxiliary discretizations 1528: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1529: - t - The time 1531: Output Parameter 1532: . elemVec - the element residual vectors from each element 1534: Level: developer 1536: .seealso: PetscFEIntegrateResidual() 1537: @*/ 1538: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, 1539: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1540: { 1541: PetscFE fe; 1546: PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe); 1547: if (fe->ops->integratehybridresidual) {(*fe->ops->integratehybridresidual)(prob, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);} 1548: return(0); 1549: } 1551: /*@C 1552: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 1554: Not collective 1556: Input Parameters: 1557: + ds - The PetscDS specifying the discretizations and continuum functions 1558: . jtype - The type of matrix pointwise functions that should be used 1559: . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1560: . Ne - The number of elements in the chunk 1561: . cgeom - The cell geometry for each cell in the chunk 1562: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1563: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1564: . probAux - The PetscDS specifying the auxiliary discretizations 1565: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1566: . t - The time 1567: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1569: Output Parameter: 1570: . elemMat - the element matrices for the Jacobian from each element 1572: Note: 1573:$ Loop over batch of elements (e):
1574: $Loop over element matrix entries (f,fc,g,gc --> i,j): 1575:$     Loop over quadrature points (q):
1576: $Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1577:$         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1578: $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1579:$                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1580: $+ \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1581: Level: intermediate 1583: .seealso: PetscFEIntegrateResidual() 1584: @*/ 1585: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1586: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1587: { 1588: PetscFE fe; 1589: PetscInt Nf; 1594: PetscDSGetNumFields(ds, &Nf); 1595: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe); 1596: if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);} 1597: return(0); 1598: } 1600: /*@C 1601: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1603: Not collective 1605: Input Parameters: 1606: + ds - The PetscDS specifying the discretizations and continuum functions 1607: . wf - The PetscWeakForm holding the pointwise functions 1608: . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1609: . Ne - The number of elements in the chunk 1610: . fgeom - The face geometry for each cell in the chunk 1611: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1612: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1613: . probAux - The PetscDS specifying the auxiliary discretizations 1614: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1615: . t - The time 1616: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1618: Output Parameter: 1619: . elemMat - the element matrices for the Jacobian from each element 1621: Note: 1622:$ Loop over batch of elements (e):
1623: $Loop over element matrix entries (f,fc,g,gc --> i,j): 1624:$     Loop over quadrature points (q):
1625: $Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1626:$         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1627: $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1628:$                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1629: $+ \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1630: Level: intermediate 1632: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1633: @*/ 1634: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1635: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1636: { 1637: PetscFE fe; 1638: PetscInt Nf; 1643: PetscDSGetNumFields(ds, &Nf); 1644: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe); 1645: if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);} 1646: return(0); 1647: } 1649: /*@C 1650: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1652: Not collective 1654: Input Parameters: 1655: + ds - The PetscDS specifying the discretizations and continuum functions 1656: . jtype - The type of matrix pointwise functions that should be used 1657: . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1658: . Ne - The number of elements in the chunk 1659: . fgeom - The face geometry for each cell in the chunk 1660: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1661: . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1662: . probAux - The PetscDS specifying the auxiliary discretizations 1663: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1664: . t - The time 1665: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1667: Output Parameter 1668: . elemMat - the element matrices for the Jacobian from each element 1670: Note: 1671:$ Loop over batch of elements (e):
1672: $Loop over element matrix entries (f,fc,g,gc --> i,j): 1673:$     Loop over quadrature points (q):
1674: $Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1675:$         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1676: $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1677:$                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1678: \$                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1679:   Level: developer

1681: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1682: @*/
1683: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1684:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1685: {
1686:   PetscFE        fe;
1687:   PetscInt       Nf;

1692:   PetscDSGetNumFields(ds, &Nf);
1693:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1694:   if (fe->ops->integratehybridjacobian) {(*fe->ops->integratehybridjacobian)(ds, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1695:   return(0);
1696: }

1698: /*@
1699:   PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height

1701:   Input Parameters:
1702: + fe     - The finite element space
1703: - height - The height of the Plex point

1705:   Output Parameter:
1706: . subfe  - The subspace of this FE space

1708:   Note: For example, if we want the subspace of this space for a face, we would choose height = 1.

1712: .seealso: PetscFECreateDefault()
1713: @*/
1714: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1715: {
1716:   PetscSpace      P, subP;
1717:   PetscDualSpace  Q, subQ;
1719:   PetscFEType     fetype;
1720:   PetscInt        dim, Nc;
1721:   PetscErrorCode  ierr;

1726:   if (height == 0) {
1727:     *subfe = fe;
1728:     return(0);
1729:   }
1730:   PetscFEGetBasisSpace(fe, &P);
1731:   PetscFEGetDualSpace(fe, &Q);
1732:   PetscFEGetNumComponents(fe, &Nc);
1734:   PetscDualSpaceGetDimension(Q, &dim);
1735:   if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
1736:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1737:   if (height <= dim) {
1738:     if (!fe->subspaces[height-1]) {
1739:       PetscFE     sub = NULL;
1740:       const char *name;

1742:       PetscSpaceGetHeightSubspace(P, height, &subP);
1743:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1744:       if (subQ) {
1745:         PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1746:         PetscObjectGetName((PetscObject) fe,  &name);
1747:         PetscObjectSetName((PetscObject) sub,  name);
1748:         PetscFEGetType(fe, &fetype);
1749:         PetscFESetType(sub, fetype);
1750:         PetscFESetBasisSpace(sub, subP);
1751:         PetscFESetDualSpace(sub, subQ);
1752:         PetscFESetNumComponents(sub, Nc);
1753:         PetscFESetUp(sub);
1755:       }
1756:       fe->subspaces[height-1] = sub;
1757:     }
1758:     *subfe = fe->subspaces[height-1];
1759:   } else {
1760:     *subfe = NULL;
1761:   }
1762:   return(0);
1763: }

1765: /*@
1766:   PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1767:   to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1768:   sparsity). It is also used to create an interpolation between regularly refined meshes.

1770:   Collective on fem

1772:   Input Parameter:
1773: . fe - The initial PetscFE

1775:   Output Parameter:
1776: . feRef - The refined PetscFE

1780: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1781: @*/
1782: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1783: {
1784:   PetscSpace       P, Pref;
1785:   PetscDualSpace   Q, Qref;
1786:   DM               K, Kref;
1788:   const PetscReal *v0, *jac;
1789:   PetscInt         numComp, numSubelements;
1790:   PetscInt         cStart, cEnd, c;
1791:   PetscDualSpace  *cellSpaces;
1792:   PetscErrorCode   ierr;

1795:   PetscFEGetBasisSpace(fe, &P);
1796:   PetscFEGetDualSpace(fe, &Q);
1798:   PetscDualSpaceGetDM(Q, &K);
1799:   /* Create space */
1800:   PetscObjectReference((PetscObject) P);
1801:   Pref = P;
1802:   /* Create dual space */
1803:   PetscDualSpaceDuplicate(Q, &Qref);
1804:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1805:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1806:   PetscDualSpaceSetDM(Qref, Kref);
1807:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1808:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1809:   /* TODO: fix for non-uniform refinement */
1810:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1811:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1812:   PetscFree(cellSpaces);
1813:   DMDestroy(&Kref);
1814:   PetscDualSpaceSetUp(Qref);
1815:   /* Create element */
1816:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1817:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1818:   PetscFESetBasisSpace(*feRef, Pref);
1819:   PetscFESetDualSpace(*feRef, Qref);
1820:   PetscFEGetNumComponents(fe,    &numComp);
1821:   PetscFESetNumComponents(*feRef, numComp);
1822:   PetscFESetUp(*feRef);
1823:   PetscSpaceDestroy(&Pref);
1824:   PetscDualSpaceDestroy(&Qref);
1826:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1827:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1830:   return(0);
1831: }

1833: /*@C
1834:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1836:   Collective

1838:   Input Parameters:
1839: + comm      - The MPI comm
1840: . dim       - The spatial dimension
1841: . Nc        - The number of components
1842: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1843: . prefix    - The options prefix, or NULL
1844: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1846:   Output Parameter:
1847: . fem - The PetscFE object

1849:   Note:
1850:   Each subobject is SetFromOption() during creation, so that the object may be customized from the command line, using the prefix specified above. See the links below for the particular options available.

1852:   Level: beginner

1854: .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1855: @*/
1856: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1857: {
1859:   DM              K;
1860:   PetscSpace      P;
1861:   PetscDualSpace  Q;
1863:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1864:   PetscErrorCode  ierr;

1867:   /* Create space */
1868:   PetscSpaceCreate(comm, &P);
1869:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1870:   PetscSpacePolynomialSetTensor(P, tensor);
1871:   PetscSpaceSetNumComponents(P, Nc);
1872:   PetscSpaceSetNumVariables(P, dim);
1873:   PetscSpaceSetFromOptions(P);
1874:   PetscSpaceSetUp(P);
1875:   PetscSpaceGetDegree(P, &order, NULL);
1876:   PetscSpacePolynomialGetTensor(P, &tensor);
1877:   /* Create dual space */
1878:   PetscDualSpaceCreate(comm, &Q);
1879:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1880:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1881:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1882:   PetscDualSpaceSetDM(Q, K);
1883:   DMDestroy(&K);
1884:   PetscDualSpaceSetNumComponents(Q, Nc);
1885:   PetscDualSpaceSetOrder(Q, order);
1886:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1887:   PetscDualSpaceSetFromOptions(Q);
1888:   PetscDualSpaceSetUp(Q);
1889:   /* Create element */
1890:   PetscFECreate(comm, fem);
1891:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1892:   PetscFESetBasisSpace(*fem, P);
1893:   PetscFESetDualSpace(*fem, Q);
1894:   PetscFESetNumComponents(*fem, Nc);
1895:   PetscFESetFromOptions(*fem);
1896:   PetscFESetUp(*fem);
1897:   PetscSpaceDestroy(&P);
1898:   PetscDualSpaceDestroy(&Q);
1899:   /* Create quadrature (with specified order if given) */
1900:   qorder = qorder >= 0 ? qorder : order;
1901:   PetscObjectOptionsBegin((PetscObject)*fem);
1903:   PetscOptionsEnd();
1904:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1905:   if (isSimplex) {
1908:   } else {
1911:   }
1916:   return(0);
1917: }

1919: /*@
1920:   PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k

1922:   Collective

1924:   Input Parameters:
1925: + comm      - The MPI comm
1926: . dim       - The spatial dimension
1927: . Nc        - The number of components
1928: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1929: . k         - The degree k of the space
1930: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1932:   Output Parameter:
1933: . fem       - The PetscFE object

1935:   Level: beginner

1937:   Notes:
1938:   For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k.

1940: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1941: @*/
1942: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1943: {
1945:   DM              K;
1946:   PetscSpace      P;
1947:   PetscDualSpace  Q;
1949:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1950:   char            name[64];
1951:   PetscErrorCode  ierr;

1954:   /* Create space */
1955:   PetscSpaceCreate(comm, &P);
1956:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1957:   PetscSpacePolynomialSetTensor(P, tensor);
1958:   PetscSpaceSetNumComponents(P, Nc);
1959:   PetscSpaceSetNumVariables(P, dim);
1960:   PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1961:   PetscSpaceSetUp(P);
1962:   /* Create dual space */
1963:   PetscDualSpaceCreate(comm, &Q);
1964:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1965:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1966:   PetscDualSpaceSetDM(Q, K);
1967:   DMDestroy(&K);
1968:   PetscDualSpaceSetNumComponents(Q, Nc);
1969:   PetscDualSpaceSetOrder(Q, k);
1970:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1971:   PetscDualSpaceSetUp(Q);
1972:   /* Create finite element */
1973:   PetscFECreate(comm, fem);
1974:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1975:   PetscObjectSetName((PetscObject) *fem, name);
1976:   PetscFESetType(*fem, PETSCFEBASIC);
1977:   PetscFESetBasisSpace(*fem, P);
1978:   PetscFESetDualSpace(*fem, Q);
1979:   PetscFESetNumComponents(*fem, Nc);
1980:   PetscFESetUp(*fem);
1981:   PetscSpaceDestroy(&P);
1982:   PetscDualSpaceDestroy(&Q);
1983:   /* Create quadrature (with specified order if given) */
1984:   qorder = qorder >= 0 ? qorder : k;
1985:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1986:   if (isSimplex) {
1989:   } else {
1992:   }
1997:   /* Set finite element name */
1998:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1999:   PetscFESetName(*fem, name);
2000:   return(0);
2001: }

2003: /*@C
2004:   PetscFESetName - Names the FE and its subobjects

2006:   Not collective

2008:   Input Parameters:
2009: + fe   - The PetscFE
2010: - name - The name

2012:   Level: intermediate

2014: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
2015: @*/
2016: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2017: {
2018:   PetscSpace     P;
2019:   PetscDualSpace Q;

2023:   PetscFEGetBasisSpace(fe, &P);
2024:   PetscFEGetDualSpace(fe, &Q);
2025:   PetscObjectSetName((PetscObject) fe, name);
2026:   PetscObjectSetName((PetscObject) P,  name);
2027:   PetscObjectSetName((PetscObject) Q,  name);
2028:   return(0);
2029: }

2031: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2032: {
2033:   PetscInt       dOffset = 0, fOffset = 0, f, g;

2036:   for (f = 0; f < Nf; ++f) {
2037:     PetscFE          fe;
2038:     const PetscInt   k    = ds->jetDegree[f];
2039:     const PetscInt   cdim = T[f]->cdim;
2040:     const PetscInt   Nq   = T[f]->Np;
2041:     const PetscInt   Nbf  = T[f]->Nb;
2042:     const PetscInt   Ncf  = T[f]->Nc;
2043:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2044:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2045:     const PetscReal *Hq   = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL;
2046:     PetscInt         hOffset = 0, b, c, d;

2048:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2049:     for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2050:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2051:     for (b = 0; b < Nbf; ++b) {
2052:       for (c = 0; c < Ncf; ++c) {
2053:         const PetscInt cidx = b*Ncf+c;

2055:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2056:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2057:       }
2058:     }
2059:     if (k > 1) {
2060:       for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim;
2061:       for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0;
2062:       for (b = 0; b < Nbf; ++b) {
2063:         for (c = 0; c < Ncf; ++c) {
2064:           const PetscInt cidx = b*Ncf+c;

2066:           for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b];
2067:         }
2068:       }
2069:       PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);
2070:     }
2071:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2073:     if (u_t) {
2074:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2075:       for (b = 0; b < Nbf; ++b) {
2076:         for (c = 0; c < Ncf; ++c) {
2077:           const PetscInt cidx = b*Ncf+c;

2079:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2080:         }
2081:       }
2082:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2083:     }
2084:     fOffset += Ncf;
2085:     dOffset += Nbf;
2086:   }
2087:   return 0;
2088: }

2090: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2091: {
2092:   PetscInt       dOffset = 0, fOffset = 0, g;

2095:   for (g = 0; g < 2*Nf-1; ++g) {
2096:     if (!T[g/2]) continue;
2097:     {
2098:     PetscFE          fe;
2099:     const PetscInt   f   = g/2;
2100:     const PetscInt   cdim = T[f]->cdim;
2101:     const PetscInt   Nq   = T[f]->Np;
2102:     const PetscInt   Nbf  = T[f]->Nb;
2103:     const PetscInt   Ncf  = T[f]->Nc;
2104:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2105:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2106:     PetscInt         b, c, d;

2108:     fe = (PetscFE) ds->disc[f];
2109:     for (c = 0; c < Ncf; ++c)      u[fOffset+c] = 0.0;
2110:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2111:     for (b = 0; b < Nbf; ++b) {
2112:       for (c = 0; c < Ncf; ++c) {
2113:         const PetscInt cidx = b*Ncf+c;

2115:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2116:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2117:       }
2118:     }
2119:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2121:     if (u_t) {
2122:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2123:       for (b = 0; b < Nbf; ++b) {
2124:         for (c = 0; c < Ncf; ++c) {
2125:           const PetscInt cidx = b*Ncf+c;

2127:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2128:         }
2129:       }
2130:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2131:     }
2132:     fOffset += Ncf;
2133:     dOffset += Nbf;
2134:   }
2135:   }
2136:   return 0;
2137: }

2139: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2140: {
2141:   PetscFE         fe;
2142:   PetscTabulation Tc;
2143:   PetscInt        b, c;
2144:   PetscErrorCode  ierr;

2146:   if (!prob) return 0;
2147:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2148:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2149:   {
2150:     const PetscReal *faceBasis = Tc->T[0];
2151:     const PetscInt   Nb        = Tc->Nb;
2152:     const PetscInt   Nc        = Tc->Nc;

2154:     for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2155:     for (b = 0; b < Nb; ++b) {
2156:       for (c = 0; c < Nc; ++c) {
2157:         u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2158:       }
2159:     }
2160:   }
2161:   return 0;
2162: }

2164: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2165: {
2166:   PetscFEGeom      pgeom;
2167:   const PetscInt   dEt      = T->cdim;
2168:   const PetscInt   dE       = fegeom->dimEmbed;
2169:   const PetscInt   Nq       = T->Np;
2170:   const PetscInt   Nb       = T->Nb;
2171:   const PetscInt   Nc       = T->Nc;
2172:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2173:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dEt];
2174:   PetscInt         q, b, c, d;
2175:   PetscErrorCode   ierr;

2177:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2178:   for (q = 0; q < Nq; ++q) {
2179:     for (b = 0; b < Nb; ++b) {
2180:       for (c = 0; c < Nc; ++c) {
2181:         const PetscInt bcidx = b*Nc+c;

2183:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2184:         for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dEt+bcidx*dEt+d];
2185:       }
2186:     }
2187:     PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);
2188:     PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);
2190:     for (b = 0; b < Nb; ++b) {
2191:       for (c = 0; c < Nc; ++c) {
2192:         const PetscInt bcidx = b*Nc+c;
2193:         const PetscInt qcidx = q*Nc+c;

2195:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2196:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2197:       }
2198:     }
2199:   }
2200:   return(0);
2201: }

2203: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2204: {
2205:   const PetscInt   dE       = T->cdim;
2206:   const PetscInt   Nq       = T->Np;
2207:   const PetscInt   Nb       = T->Nb;
2208:   const PetscInt   Nc       = T->Nc;
2209:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2210:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2211:   PetscInt         q, b, c, d;
2212:   PetscErrorCode   ierr;

2214:   for (b = 0; b < Nb; ++b) elemVec[Nb*s+b] = 0.0;
2215:   for (q = 0; q < Nq; ++q) {
2216:     for (b = 0; b < Nb; ++b) {
2217:       for (c = 0; c < Nc; ++c) {
2218:         const PetscInt bcidx = b*Nc+c;

2220:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2221:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2222:       }
2223:     }
2224:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2226:     for (b = 0; b < Nb; ++b) {
2227:       for (c = 0; c < Nc; ++c) {
2228:         const PetscInt bcidx = b*Nc+c;
2229:         const PetscInt qcidx = q*Nc+c;

2231:         elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2232:         for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2233:       }
2234:     }
2235:   }
2236:   return(0);
2237: }

2239: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2240: {
2241:   const PetscInt   dE        = TI->cdim;
2242:   const PetscInt   NqI       = TI->Np;
2243:   const PetscInt   NbI       = TI->Nb;
2244:   const PetscInt   NcI       = TI->Nc;
2245:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2246:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2247:   const PetscInt   NqJ       = TJ->Np;
2248:   const PetscInt   NbJ       = TJ->Nb;
2249:   const PetscInt   NcJ       = TJ->Nc;
2250:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2251:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2252:   PetscInt         f, fc, g, gc, df, dg;
2253:   PetscErrorCode   ierr;

2255:   for (f = 0; f < NbI; ++f) {
2256:     for (fc = 0; fc < NcI; ++fc) {
2257:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2259:       tmpBasisI[fidx] = basisI[fidx];
2260:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2261:     }
2262:   }
2263:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2265:   for (g = 0; g < NbJ; ++g) {
2266:     for (gc = 0; gc < NcJ; ++gc) {
2267:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2269:       tmpBasisJ[gidx] = basisJ[gidx];
2270:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2271:     }
2272:   }
2273:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2275:   for (f = 0; f < NbI; ++f) {
2276:     for (fc = 0; fc < NcI; ++fc) {
2277:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2278:       const PetscInt i    = offsetI+f; /* Element matrix row */
2279:       for (g = 0; g < NbJ; ++g) {
2280:         for (gc = 0; gc < NcJ; ++gc) {
2281:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2282:           const PetscInt j    = offsetJ+g; /* Element matrix column */
2283:           const PetscInt fOff = eOffset+i*totDim+j;

2285:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2286:           for (df = 0; df < dE; ++df) {
2287:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2288:             elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2289:             for (dg = 0; dg < dE; ++dg) {
2290:               elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2291:             }
2292:           }
2293:         }
2294:       }
2295:     }
2296:   }
2297:   return(0);
2298: }

2300: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2301: {
2302:   const PetscInt   dE        = TI->cdim;
2303:   const PetscInt   NqI       = TI->Np;
2304:   const PetscInt   NbI       = TI->Nb;
2305:   const PetscInt   NcI       = TI->Nc;
2306:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2307:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2308:   const PetscInt   NqJ       = TJ->Np;
2309:   const PetscInt   NbJ       = TJ->Nb;
2310:   const PetscInt   NcJ       = TJ->Nc;
2311:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2312:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2313:   const PetscInt   Ns = isHybridI ? 1 : 2;
2314:   const PetscInt   Nt = isHybridJ ? 1 : 2;
2315:   PetscInt         f, fc, g, gc, df, dg, s, t;
2316:   PetscErrorCode   ierr;

2318:   for (f = 0; f < NbI; ++f) {
2319:     for (fc = 0; fc < NcI; ++fc) {
2320:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2322:       tmpBasisI[fidx] = basisI[fidx];
2323:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2324:     }
2325:   }
2326:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2328:   for (g = 0; g < NbJ; ++g) {
2329:     for (gc = 0; gc < NcJ; ++gc) {
2330:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2332:       tmpBasisJ[gidx] = basisJ[gidx];
2333:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2334:     }
2335:   }
2336:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2338:   for (s = 0; s < Ns; ++s) {
2339:     for (f = 0; f < NbI; ++f) {
2340:       for (fc = 0; fc < NcI; ++fc) {
2341:         const PetscInt sc   = NcI*s+fc;  /* components from each side of the surface */
2342:         const PetscInt fidx = f*NcI+fc;  /* Test function basis index */
2343:         const PetscInt i    = offsetI+NbI*s+f; /* Element matrix row */
2344:         for (t = 0; t < Nt; ++t) {
2345:           for (g = 0; g < NbJ; ++g) {
2346:             for (gc = 0; gc < NcJ; ++gc) {
2347:               const PetscInt tc   = NcJ*t+gc;  /* components from each side of the surface */
2348:               const PetscInt gidx = g*NcJ+gc;  /* Trial function basis index */
2349:               const PetscInt j    = offsetJ+NbJ*t+g; /* Element matrix column */
2350:               const PetscInt fOff = eOffset+i*totDim+j;

2352:               elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2353:               for (df = 0; df < dE; ++df) {
2354:                 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2355:                 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2356:                 for (dg = 0; dg < dE; ++dg) {
2357:                   elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2358:                 }
2359:               }
2360:             }
2361:           }
2362:         }
2363:       }
2364:     }
2365:   }
2366:   return(0);
2367: }

2370: {
2371:   PetscDualSpace  dsp;
2372:   DM              dm;
2374:   PetscInt        dim, cdim, Nq;
2375:   PetscErrorCode  ierr;

2378:   PetscFEGetDualSpace(fe, &dsp);
2379:   PetscDualSpaceGetDM(dsp, &dm);
2380:   DMGetDimension(dm, &dim);
2381:   DMGetCoordinateDim(dm, &cdim);
2385:   PetscMalloc1(Nq*cdim,      &cgeom->v);
2386:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2387:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2388:   PetscMalloc1(Nq,           &cgeom->detJ);
2389:   cgeom->dim       = dim;
2390:   cgeom->dimEmbed  = cdim;
2391:   cgeom->numCells  = 1;
2392:   cgeom->numPoints = Nq;
2393:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2394:   return(0);
2395: }

2397: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2398: {

2402:   PetscFree(cgeom->v);
2403:   PetscFree(cgeom->J);
2404:   PetscFree(cgeom->invJ);
2405:   PetscFree(cgeom->detJ);
2406:   return(0);
2407: }