Actual source code: fe.c

petsc-master 2019-10-14
Report Typos and Errors
  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: PetscFunctionList PetscFEList              = NULL;
 54: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

 56: /*@C
 57:   PetscFERegister - Adds a new PetscFE implementation

 59:   Not Collective

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

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

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

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

 83:   Level: advanced

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

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

 93:   PetscFunctionListAdd(&PetscFEList, sname, function);
 94:   return(0);
 95: }

 97: /*@C
 98:   PetscFESetType - Builds a particular PetscFE

100:   Collective on fem

102:   Input Parameters:
103: + fem  - The PetscFE object
104: - name - The kind of FEM space

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

109:   Level: intermediate

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

121:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
122:   if (match) return(0);

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

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

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

140:   Not Collective

142:   Input Parameter:
143: . fem  - The PetscFE

145:   Output Parameter:
146: . name - The PetscFE type name

148:   Level: intermediate

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

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

166: /*@C
167:   PetscFEView - Views a PetscFE

169:   Collective on fem

171:   Input Parameter:
172: + fem - the PetscFE object to view
173: - viewer   - the viewer

175:   Level: beginner

177: .seealso PetscFEDestroy()
178: @*/
179: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
180: {
181:   PetscBool      iascii;

187:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
188:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
189:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
190:   if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
191:   return(0);
192: }

194: /*@
195:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

197:   Collective on fem

199:   Input Parameter:
200: . fem - the PetscFE object to set options for

202:   Options Database:
203: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
204: - -petscfe_num_batches - the number of cell batches to integrate serially

206:   Level: intermediate

208: .seealso PetscFEView()
209: @*/
210: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
211: {
212:   const char    *defaultType;
213:   char           name[256];
214:   PetscBool      flg;

219:   if (!((PetscObject) fem)->type_name) {
220:     defaultType = PETSCFEBASIC;
221:   } else {
222:     defaultType = ((PetscObject) fem)->type_name;
223:   }
224:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

226:   PetscObjectOptionsBegin((PetscObject) fem);
227:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
228:   if (flg) {
229:     PetscFESetType(fem, name);
230:   } else if (!((PetscObject) fem)->type_name) {
231:     PetscFESetType(fem, defaultType);
232:   }
233:   PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
234:   PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
235:   if (fem->ops->setfromoptions) {
236:     (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
237:   }
238:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
239:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
240:   PetscOptionsEnd();
241:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
242:   return(0);
243: }

245: /*@C
246:   PetscFESetUp - Construct data structures for the PetscFE

248:   Collective on fem

250:   Input Parameter:
251: . fem - the PetscFE object to setup

253:   Level: intermediate

255: .seealso PetscFEView(), PetscFEDestroy()
256: @*/
257: PetscErrorCode PetscFESetUp(PetscFE fem)
258: {

263:   if (fem->setupcalled) return(0);
264:   fem->setupcalled = PETSC_TRUE;
265:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
266:   return(0);
267: }

269: /*@
270:   PetscFEDestroy - Destroys a PetscFE object

272:   Collective on fem

274:   Input Parameter:
275: . fem - the PetscFE object to destroy

277:   Level: beginner

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

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

289:   if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
290:   ((PetscObject) (*fem))->refct = 0;

292:   if ((*fem)->subspaces) {
293:     PetscInt dim, d;

295:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
296:     for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
297:   }
298:   PetscFree((*fem)->subspaces);
299:   PetscFree((*fem)->invV);
300:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
301:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
302:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
303:   PetscSpaceDestroy(&(*fem)->basisSpace);
304:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
305:   PetscQuadratureDestroy(&(*fem)->quadrature);
306:   PetscQuadratureDestroy(&(*fem)->faceQuadrature);

308:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
309:   PetscHeaderDestroy(fem);
310:   return(0);
311: }

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

316:   Collective

318:   Input Parameter:
319: . comm - The communicator for the PetscFE object

321:   Output Parameter:
322: . fem - The PetscFE object

324:   Level: beginner

326: .seealso: PetscFESetType(), PETSCFEGALERKIN
327: @*/
328: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
329: {
330:   PetscFE        f;

335:   PetscCitationsRegister(FECitation,&FEcite);
336:   *fem = NULL;
337:   PetscFEInitializePackage();

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

341:   f->basisSpace    = NULL;
342:   f->dualSpace     = NULL;
343:   f->numComponents = 1;
344:   f->subspaces     = NULL;
345:   f->invV          = NULL;
346:   f->B             = NULL;
347:   f->D             = NULL;
348:   f->H             = NULL;
349:   f->Bf            = NULL;
350:   f->Df            = NULL;
351:   f->Hf            = NULL;
352:   PetscArrayzero(&f->quadrature, 1);
353:   PetscArrayzero(&f->faceQuadrature, 1);
354:   f->blockSize     = 0;
355:   f->numBlocks     = 1;
356:   f->batchSize     = 0;
357:   f->numBatches    = 1;

359:   *fem = f;
360:   return(0);
361: }

363: /*@
364:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

366:   Not collective

368:   Input Parameter:
369: . fem - The PetscFE object

371:   Output Parameter:
372: . dim - The spatial dimension

374:   Level: intermediate

376: .seealso: PetscFECreate()
377: @*/
378: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
379: {
380:   DM             dm;

386:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
387:   DMGetDimension(dm, dim);
388:   return(0);
389: }

391: /*@
392:   PetscFESetNumComponents - Sets the number of components in the element

394:   Not collective

396:   Input Parameters:
397: + fem - The PetscFE object
398: - comp - The number of field components

400:   Level: intermediate

402: .seealso: PetscFECreate()
403: @*/
404: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
405: {
408:   fem->numComponents = comp;
409:   return(0);
410: }

412: /*@
413:   PetscFEGetNumComponents - Returns the number of components in the element

415:   Not collective

417:   Input Parameter:
418: . fem - The PetscFE object

420:   Output Parameter:
421: . comp - The number of field components

423:   Level: intermediate

425: .seealso: PetscFECreate()
426: @*/
427: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
428: {
432:   *comp = fem->numComponents;
433:   return(0);
434: }

436: /*@
437:   PetscFESetTileSizes - Sets the tile sizes for evaluation

439:   Not collective

441:   Input Parameters:
442: + fem - The PetscFE object
443: . blockSize - The number of elements in a block
444: . numBlocks - The number of blocks in a batch
445: . batchSize - The number of elements in a batch
446: - numBatches - The number of batches in a chunk

448:   Level: intermediate

450: .seealso: PetscFECreate()
451: @*/
452: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
453: {
456:   fem->blockSize  = blockSize;
457:   fem->numBlocks  = numBlocks;
458:   fem->batchSize  = batchSize;
459:   fem->numBatches = numBatches;
460:   return(0);
461: }

463: /*@
464:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

466:   Not collective

468:   Input Parameter:
469: . fem - The PetscFE object

471:   Output Parameters:
472: + blockSize - The number of elements in a block
473: . numBlocks - The number of blocks in a batch
474: . batchSize - The number of elements in a batch
475: - numBatches - The number of batches in a chunk

477:   Level: intermediate

479: .seealso: PetscFECreate()
480: @*/
481: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
482: {
489:   if (blockSize)  *blockSize  = fem->blockSize;
490:   if (numBlocks)  *numBlocks  = fem->numBlocks;
491:   if (batchSize)  *batchSize  = fem->batchSize;
492:   if (numBatches) *numBatches = fem->numBatches;
493:   return(0);
494: }

496: /*@
497:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

499:   Not collective

501:   Input Parameter:
502: . fem - The PetscFE object

504:   Output Parameter:
505: . sp - The PetscSpace object

507:   Level: intermediate

509: .seealso: PetscFECreate()
510: @*/
511: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
512: {
516:   *sp = fem->basisSpace;
517:   return(0);
518: }

520: /*@
521:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

523:   Not collective

525:   Input Parameters:
526: + fem - The PetscFE object
527: - sp - The PetscSpace object

529:   Level: intermediate

531: .seealso: PetscFECreate()
532: @*/
533: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
534: {

540:   PetscSpaceDestroy(&fem->basisSpace);
541:   fem->basisSpace = sp;
542:   PetscObjectReference((PetscObject) fem->basisSpace);
543:   return(0);
544: }

546: /*@
547:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

549:   Not collective

551:   Input Parameter:
552: . fem - The PetscFE object

554:   Output Parameter:
555: . sp - The PetscDualSpace object

557:   Level: intermediate

559: .seealso: PetscFECreate()
560: @*/
561: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
562: {
566:   *sp = fem->dualSpace;
567:   return(0);
568: }

570: /*@
571:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

573:   Not collective

575:   Input Parameters:
576: + fem - The PetscFE object
577: - sp - The PetscDualSpace object

579:   Level: intermediate

581: .seealso: PetscFECreate()
582: @*/
583: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
584: {

590:   PetscDualSpaceDestroy(&fem->dualSpace);
591:   fem->dualSpace = sp;
592:   PetscObjectReference((PetscObject) fem->dualSpace);
593:   return(0);
594: }

596: /*@
597:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

599:   Not collective

601:   Input Parameter:
602: . fem - The PetscFE object

604:   Output Parameter:
605: . q - The PetscQuadrature object

607:   Level: intermediate

609: .seealso: PetscFECreate()
610: @*/
611: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
612: {
616:   *q = fem->quadrature;
617:   return(0);
618: }

620: /*@
621:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

623:   Not collective

625:   Input Parameters:
626: + fem - The PetscFE object
627: - q - The PetscQuadrature object

629:   Level: intermediate

631: .seealso: PetscFECreate()
632: @*/
633: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
634: {
635:   PetscInt       Nc, qNc;

640:   PetscFEGetNumComponents(fem, &Nc);
641:   PetscQuadratureGetNumComponents(q, &qNc);
642:   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);
643:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
644:   PetscQuadratureDestroy(&fem->quadrature);
645:   fem->quadrature = q;
646:   PetscObjectReference((PetscObject) q);
647:   return(0);
648: }

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

653:   Not collective

655:   Input Parameter:
656: . fem - The PetscFE object

658:   Output Parameter:
659: . q - The PetscQuadrature object

661:   Level: intermediate

663: .seealso: PetscFECreate()
664: @*/
665: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
666: {
670:   *q = fem->faceQuadrature;
671:   return(0);
672: }

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

677:   Not collective

679:   Input Parameters:
680: + fem - The PetscFE object
681: - q - The PetscQuadrature object

683:   Level: intermediate

685: .seealso: PetscFECreate()
686: @*/
687: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
688: {

693:   PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
694:   PetscQuadratureDestroy(&fem->faceQuadrature);
695:   fem->faceQuadrature = q;
696:   PetscObjectReference((PetscObject) q);
697:   return(0);
698: }

700: /*@
701:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

703:   Not collective

705:   Input Parameters:
706: + sfe - The PetscFE source for the quadratures
707: - tfe - The PetscFE target for the quadratures

709:   Level: intermediate

711: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
712: @*/
713: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
714: {
715:   PetscQuadrature q;
716:   PetscErrorCode  ierr;

721:   PetscFEGetQuadrature(sfe, &q);
722:   PetscFESetQuadrature(tfe,  q);
723:   PetscFEGetFaceQuadrature(sfe, &q);
724:   PetscFESetFaceQuadrature(tfe,  q);
725:   return(0);
726: }

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

731:   Not collective

733:   Input Parameter:
734: . fem - The PetscFE object

736:   Output Parameter:
737: . numDof - Array with the number of dofs per dimension

739:   Level: intermediate

741: .seealso: PetscFECreate()
742: @*/
743: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
744: {

750:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
751:   return(0);
752: }

754: /*@C
755:   PetscFEGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points

757:   Not collective

759:   Input Parameter:
760: . fem - The PetscFE object

762:   Output Parameters:
763: + B - The basis function values at quadrature points
764: . D - The basis function derivatives at quadrature points
765: - H - The basis function second derivatives at quadrature points

767:   Note:
768: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
769: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
770: $ 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

772:   Level: intermediate

774: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
775: @*/
776: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
777: {
778:   PetscInt         npoints;
779:   const PetscReal *points;
780:   PetscErrorCode   ierr;

787:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
788:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
789:   if (B) *B = fem->B;
790:   if (D) *D = fem->D;
791:   if (H) *H = fem->H;
792:   return(0);
793: }

795: /*@C
796:   PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points

798:   Not collective

800:   Input Parameter:
801: . fem - The PetscFE object

803:   Output Parameters:
804: + B - The basis function values at face quadrature points
805: . D - The basis function derivatives at face quadrature points
806: - H - The basis function second derivatives at face quadrature points

808:   Note:
809: $ Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
810: $ 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
811: $ 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

813:   Level: intermediate

815: .seealso: PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
816: @*/
817: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
818: {
819:   PetscErrorCode   ierr;

826:   if (!fem->Bf) {
827:     const PetscReal  xi0[3] = {-1., -1., -1.};
828:     PetscReal        v0[3], J[9], detJ;
829:     PetscQuadrature  fq;
830:     PetscDualSpace   sp;
831:     DM               dm;
832:     const PetscInt  *faces;
833:     PetscInt         dim, numFaces, f, npoints, q;
834:     const PetscReal *points;
835:     PetscReal       *facePoints;

837:     PetscFEGetDualSpace(fem, &sp);
838:     PetscDualSpaceGetDM(sp, &dm);
839:     DMGetDimension(dm, &dim);
840:     DMPlexGetConeSize(dm, 0, &numFaces);
841:     DMPlexGetCone(dm, 0, &faces);
842:     PetscFEGetFaceQuadrature(fem, &fq);
843:     if (fq) {
844:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
845:       PetscMalloc1(numFaces*npoints*dim, &facePoints);
846:       for (f = 0; f < numFaces; ++f) {
847:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
848:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
849:       }
850:       PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
851:       PetscFree(facePoints);
852:     }
853:   }
854:   if (Bf) *Bf = fem->Bf;
855:   if (Df) *Df = fem->Df;
856:   if (Hf) *Hf = fem->Hf;
857:   return(0);
858: }

860: /*@C
861:   PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face centroid points

863:   Not collective

865:   Input Parameter:
866: . fem - The PetscFE object

868:   Output Parameters:
869: + B - The basis function values at face centroid points
870: . D - The basis function derivatives at face centroid points
871: - H - The basis function second derivatives at face centroid points

873:   Note:
874: $ Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
875: $ Df[((f*pdim + i)*Nc + c)*dim + d] is the derivative value at point f for basis function i, component c, in direction d
876: $ Hf[(((f*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f for basis function i, component c, in directions d and e

878:   Level: intermediate

880: .seealso: PetscFEGetFaceTabulation(), PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
881: @*/
882: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
883: {
884:   PetscErrorCode   ierr;

889:   if (!fem->F) {
890:     PetscDualSpace  sp;
891:     DM              dm;
892:     const PetscInt *cone;
893:     PetscReal      *centroids;
894:     PetscInt        dim, numFaces, f;

896:     PetscFEGetDualSpace(fem, &sp);
897:     PetscDualSpaceGetDM(sp, &dm);
898:     DMGetDimension(dm, &dim);
899:     DMPlexGetConeSize(dm, 0, &numFaces);
900:     DMPlexGetCone(dm, 0, &cone);
901:     PetscMalloc1(numFaces*dim, &centroids);
902:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
903:     PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
904:     PetscFree(centroids);
905:   }
906:   *F = fem->F;
907:   return(0);
908: }

910: /*@C
911:   PetscFEGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

913:   Not collective

915:   Input Parameters:
916: + fem     - The PetscFE object
917: . npoints - The number of tabulation points
918: - points  - The tabulation point coordinates

920:   Output Parameters:
921: + B - The basis function values at tabulation points
922: . D - The basis function derivatives at tabulation points
923: - H - The basis function second derivatives at tabulation points

925:   Note:
926: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
927: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
928: $ 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

930:   Level: intermediate

932: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
933: @*/
934: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
935: {
936:   DM               dm;
937:   PetscInt         pdim; /* Dimension of FE space P */
938:   PetscInt         dim;  /* Spatial dimension */
939:   PetscInt         comp; /* Field components */
940:   PetscErrorCode   ierr;

943:   if (!npoints) {
944:     if (B) *B = NULL;
945:     if (D) *D = NULL;
946:     if (H) *H = NULL;
947:     return(0);
948:   }
954:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
955:   DMGetDimension(dm, &dim);
956:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
957:   PetscFEGetNumComponents(fem, &comp);
958:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
959:   if (!dim) {
960:     if (D) *D = NULL;
961:     if (H) *H = NULL;
962:   } else {
963:     if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
964:     if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
965:   }
966:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
967:   return(0);
968: }

970: /*@C
971:   PetscFERestoreTabulation - Frees memory from the associated tabulation.

973:   Not collective

975:   Input Parameters:
976: + fem     - The PetscFE object
977: . npoints - The number of tabulation points
978: . points  - The tabulation point coordinates
979: . B - The basis function values at tabulation points
980: . D - The basis function derivatives at tabulation points
981: - H - The basis function second derivatives at tabulation points

983:   Note:
984: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
985: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
986: $ 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

988:   Level: intermediate

990: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation()
991: @*/
992: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
993: {
994:   DM             dm;

999:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
1000:   if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
1001:   if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
1002:   if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
1003:   return(0);
1004: }

1006: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1007: {
1008:   PetscSpace     bsp, bsubsp;
1009:   PetscDualSpace dsp, dsubsp;
1010:   PetscInt       dim, depth, numComp, i, j, coneSize, order;
1011:   PetscFEType    type;
1012:   DM             dm;
1013:   DMLabel        label;
1014:   PetscReal      *xi, *v, *J, detJ;
1015:   const char     *name;
1016:   PetscQuadrature origin, fullQuad, subQuad;

1022:   PetscFEGetBasisSpace(fe,&bsp);
1023:   PetscFEGetDualSpace(fe,&dsp);
1024:   PetscDualSpaceGetDM(dsp,&dm);
1025:   DMGetDimension(dm,&dim);
1026:   DMPlexGetDepthLabel(dm,&label);
1027:   DMLabelGetValue(label,refPoint,&depth);
1028:   PetscCalloc1(depth,&xi);
1029:   PetscMalloc1(dim,&v);
1030:   PetscMalloc1(dim*dim,&J);
1031:   for (i = 0; i < depth; i++) xi[i] = 0.;
1032:   PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1033:   PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1034:   DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1035:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1036:   for (i = 1; i < dim; i++) {
1037:     for (j = 0; j < depth; j++) {
1038:       J[i * depth + j] = J[i * dim + j];
1039:     }
1040:   }
1041:   PetscQuadratureDestroy(&origin);
1042:   PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1043:   PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1044:   PetscSpaceSetUp(bsubsp);
1045:   PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1046:   PetscFEGetType(fe,&type);
1047:   PetscFESetType(*trFE,type);
1048:   PetscFEGetNumComponents(fe,&numComp);
1049:   PetscFESetNumComponents(*trFE,numComp);
1050:   PetscFESetBasisSpace(*trFE,bsubsp);
1051:   PetscFESetDualSpace(*trFE,dsubsp);
1052:   PetscObjectGetName((PetscObject) fe, &name);
1053:   if (name) {PetscFESetName(*trFE, name);}
1054:   PetscFEGetQuadrature(fe,&fullQuad);
1055:   PetscQuadratureGetOrder(fullQuad,&order);
1056:   DMPlexGetConeSize(dm,refPoint,&coneSize);
1057:   if (coneSize == 2 * depth) {
1058:     PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1059:   } else {
1060:     PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1061:   }
1062:   PetscFESetQuadrature(*trFE,subQuad);
1063:   PetscFESetUp(*trFE);
1064:   PetscQuadratureDestroy(&subQuad);
1065:   PetscSpaceDestroy(&bsubsp);
1066:   return(0);
1067: }

1069: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1070: {
1071:   PetscInt       hStart, hEnd;
1072:   PetscDualSpace dsp;
1073:   DM             dm;

1079:   *trFE = NULL;
1080:   PetscFEGetDualSpace(fe,&dsp);
1081:   PetscDualSpaceGetDM(dsp,&dm);
1082:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1083:   if (hEnd <= hStart) return(0);
1084:   PetscFECreatePointTrace(fe,hStart,trFE);
1085:   return(0);
1086: }


1089: /*@
1090:   PetscFEGetDimension - Get the dimension of the finite element space on a cell

1092:   Not collective

1094:   Input Parameter:
1095: . fe - The PetscFE

1097:   Output Parameter:
1098: . dim - The dimension

1100:   Level: intermediate

1102: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1103: @*/
1104: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1105: {

1111:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1112:   return(0);
1113: }

1115: /*@C
1116:   PetscFEPushforward - Map the reference element function to real space

1118:   Input Parameters:
1119: + fe     - The PetscFE
1120: . fegeom - The cell geometry
1121: . Nv     - The number of function values
1122: - vals   - The function values

1124:   Output Parameter:
1125: . vals   - The transformed function values

1127:   Level: advanced

1129:   Note: This just forwards the call onto PetscDualSpacePushforward().

1131: .seealso: PetscDualSpacePushforward()
1132: @*/
1133: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1134: {

1138:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1139:   return(0);
1140: }

1142: /*@C
1143:   PetscFEPushforwardGradient - Map the reference element function gradient to real space

1145:   Input Parameters:
1146: + fe     - The PetscFE
1147: . fegeom - The cell geometry
1148: . Nv     - The number of function gradient values
1149: - vals   - The function gradient values

1151:   Output Parameter:
1152: . vals   - The transformed function gradient values

1154:   Level: advanced

1156:   Note: This just forwards the call onto PetscDualSpacePushforwardGradient().

1158: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1159: @*/
1160: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1161: {

1165:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1166:   return(0);
1167: }

1169: /*
1170: Purpose: Compute element vector for chunk of elements

1172: Input:
1173:   Sizes:
1174:      Ne:  number of elements
1175:      Nf:  number of fields
1176:      PetscFE
1177:        dim: spatial dimension
1178:        Nb:  number of basis functions
1179:        Nc:  number of field components
1180:        PetscQuadrature
1181:          Nq:  number of quadrature points

1183:   Geometry:
1184:      PetscFEGeom[Ne] possibly *Nq
1185:        PetscReal v0s[dim]
1186:        PetscReal n[dim]
1187:        PetscReal jacobians[dim*dim]
1188:        PetscReal jacobianInverses[dim*dim]
1189:        PetscReal jacobianDeterminants
1190:   FEM:
1191:      PetscFE
1192:        PetscQuadrature
1193:          PetscReal   quadPoints[Nq*dim]
1194:          PetscReal   quadWeights[Nq]
1195:        PetscReal   basis[Nq*Nb*Nc]
1196:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1197:      PetscScalar coefficients[Ne*Nb*Nc]
1198:      PetscScalar elemVec[Ne*Nb*Nc]

1200:   Problem:
1201:      PetscInt f: the active field
1202:      f0, f1

1204:   Work Space:
1205:      PetscFE
1206:        PetscScalar f0[Nq*dim];
1207:        PetscScalar f1[Nq*dim*dim];
1208:        PetscScalar u[Nc];
1209:        PetscScalar gradU[Nc*dim];
1210:        PetscReal   x[dim];
1211:        PetscScalar realSpaceDer[dim];

1213: Purpose: Compute element vector for N_cb batches of elements

1215: Input:
1216:   Sizes:
1217:      N_cb: Number of serial cell batches

1219:   Geometry:
1220:      PetscReal v0s[Ne*dim]
1221:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1222:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1223:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1224:   FEM:
1225:      static PetscReal   quadPoints[Nq*dim]
1226:      static PetscReal   quadWeights[Nq]
1227:      static PetscReal   basis[Nq*Nb*Nc]
1228:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1229:      PetscScalar coefficients[Ne*Nb*Nc]
1230:      PetscScalar elemVec[Ne*Nb*Nc]

1232: ex62.c:
1233:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1234:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1235:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1236:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1238: ex52.c:
1239:   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)
1240:   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)

1242: ex52_integrateElement.cu
1243: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)

1245: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1246:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1247:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

1249: ex52_integrateElementOpenCL.c:
1250: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1251:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1252:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

1254: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1255: */

1257: /*@C
1258:   PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration

1260:   Not collective

1262:   Input Parameters:
1263: + fem          - The PetscFE object for the field being integrated
1264: . prob         - The PetscDS specifying the discretizations and continuum functions
1265: . field        - The field being integrated
1266: . Ne           - The number of elements in the chunk
1267: . cgeom        - The cell geometry for each cell in the chunk
1268: . coefficients - The array of FEM basis coefficients for the elements
1269: . probAux      - The PetscDS specifying the auxiliary discretizations
1270: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1272:   Output Parameter
1273: . integral     - the integral for this field

1275:   Level: intermediate

1277: .seealso: PetscFEIntegrateResidual()
1278: @*/
1279: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1280:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1281: {
1282:   PetscFE        fe;

1287:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1288:   if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1289:   return(0);
1290: }

1292: /*@C
1293:   PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration

1295:   Not collective

1297:   Input Parameters:
1298: + fem          - The PetscFE object for the field being integrated
1299: . prob         - The PetscDS specifying the discretizations and continuum functions
1300: . field        - The field being integrated
1301: . obj_func     - The function to be integrated
1302: . Ne           - The number of elements in the chunk
1303: . fgeom        - The face geometry for each face in the chunk
1304: . coefficients - The array of FEM basis coefficients for the elements
1305: . probAux      - The PetscDS specifying the auxiliary discretizations
1306: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1308:   Output Parameter
1309: . integral     - the integral for this field

1311:   Level: intermediate

1313: .seealso: PetscFEIntegrateResidual()
1314: @*/
1315: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1316:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1317:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1318:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1319:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1320:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1321: {
1322:   PetscFE        fe;

1327:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1328:   if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1329:   return(0);
1330: }

1332: /*@C
1333:   PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration

1335:   Not collective

1337:   Input Parameters:
1338: + fem          - The PetscFE object for the field being integrated
1339: . prob         - The PetscDS specifying the discretizations and continuum functions
1340: . field        - The field being integrated
1341: . Ne           - The number of elements in the chunk
1342: . cgeom        - The cell geometry for each cell in the chunk
1343: . coefficients - The array of FEM basis coefficients for the elements
1344: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1345: . probAux      - The PetscDS specifying the auxiliary discretizations
1346: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1347: - t            - The time

1349:   Output Parameter
1350: . elemVec      - the element residual vectors from each element

1352:   Note:
1353: $ Loop over batch of elements (e):
1354: $   Loop over quadrature points (q):
1355: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1356: $     Call f_0 and f_1
1357: $   Loop over element vector entries (f,fc --> i):
1358: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1360:   Level: intermediate

1362: .seealso: PetscFEIntegrateResidual()
1363: @*/
1364: PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1365:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1366: {
1367:   PetscFE        fe;

1372:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1373:   if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1374:   return(0);
1375: }

1377: /*@C
1378:   PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary

1380:   Not collective

1382:   Input Parameters:
1383: + fem          - The PetscFE object for the field being integrated
1384: . prob         - The PetscDS specifying the discretizations and continuum functions
1385: . field        - The field being integrated
1386: . Ne           - The number of elements in the chunk
1387: . fgeom        - The face geometry for each cell in the chunk
1388: . coefficients - The array of FEM basis coefficients for the elements
1389: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1390: . probAux      - The PetscDS specifying the auxiliary discretizations
1391: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1392: - t            - The time

1394:   Output Parameter
1395: . elemVec      - the element residual vectors from each element

1397:   Level: intermediate

1399: .seealso: PetscFEIntegrateResidual()
1400: @*/
1401: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1402:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1403: {
1404:   PetscFE        fe;

1409:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1410:   if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1411:   return(0);
1412: }

1414: /*@C
1415:   PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration

1417:   Not collective

1419:   Input Parameters:
1420: + fem          - The PetscFE object for the field being integrated
1421: . prob         - The PetscDS specifying the discretizations and continuum functions
1422: . jtype        - The type of matrix pointwise functions that should be used
1423: . fieldI       - The test field being integrated
1424: . fieldJ       - The basis field being integrated
1425: . Ne           - The number of elements in the chunk
1426: . cgeom        - The cell geometry for each cell in the chunk
1427: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1428: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1429: . probAux      - The PetscDS specifying the auxiliary discretizations
1430: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1431: . t            - The time
1432: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1434:   Output Parameter
1435: . elemMat      - the element matrices for the Jacobian from each element

1437:   Note:
1438: $ Loop over batch of elements (e):
1439: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1440: $     Loop over quadrature points (q):
1441: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1442: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1443: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1444: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1445: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1446:   Level: intermediate

1448: .seealso: PetscFEIntegrateResidual()
1449: @*/
1450: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1451:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1452: {
1453:   PetscFE        fe;

1458:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1459:   if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1460:   return(0);
1461: }

1463: /*@C
1464:   PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration

1466:   Not collective

1468:   Input Parameters:
1469: . prob         - The PetscDS specifying the discretizations and continuum functions
1470: . fieldI       - The test field being integrated
1471: . fieldJ       - The basis field being integrated
1472: . Ne           - The number of elements in the chunk
1473: . fgeom        - The face geometry for each cell in the chunk
1474: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1475: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1476: . probAux      - The PetscDS specifying the auxiliary discretizations
1477: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1478: . t            - The time
1479: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1481:   Output Parameter
1482: . elemMat              - the element matrices for the Jacobian from each element

1484:   Note:
1485: $ Loop over batch of elements (e):
1486: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1487: $     Loop over quadrature points (q):
1488: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1489: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1490: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1491: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1492: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1493:   Level: intermediate

1495: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1496: @*/
1497: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1498:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1499: {
1500:   PetscFE        fe;

1505:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1506:   if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1507:   return(0);
1508: }

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

1513:   Input Parameters:
1514: + fe     - The finite element space
1515: - height - The height of the Plex point

1517:   Output Parameter:
1518: . subfe  - The subspace of this FE space

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

1522:   Level: advanced

1524: .seealso: PetscFECreateDefault()
1525: @*/
1526: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1527: {
1528:   PetscSpace      P, subP;
1529:   PetscDualSpace  Q, subQ;
1530:   PetscQuadrature subq;
1531:   PetscFEType     fetype;
1532:   PetscInt        dim, Nc;
1533:   PetscErrorCode  ierr;

1538:   if (height == 0) {
1539:     *subfe = fe;
1540:     return(0);
1541:   }
1542:   PetscFEGetBasisSpace(fe, &P);
1543:   PetscFEGetDualSpace(fe, &Q);
1544:   PetscFEGetNumComponents(fe, &Nc);
1545:   PetscFEGetFaceQuadrature(fe, &subq);
1546:   PetscDualSpaceGetDimension(Q, &dim);
1547:   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);}
1548:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1549:   if (height <= dim) {
1550:     if (!fe->subspaces[height-1]) {
1551:       PetscFE     sub;
1552:       const char *name;

1554:       PetscSpaceGetHeightSubspace(P, height, &subP);
1555:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1556:       PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1557:       PetscObjectGetName((PetscObject) fe,  &name);
1558:       PetscObjectSetName((PetscObject) sub,  name);
1559:       PetscFEGetType(fe, &fetype);
1560:       PetscFESetType(sub, fetype);
1561:       PetscFESetBasisSpace(sub, subP);
1562:       PetscFESetDualSpace(sub, subQ);
1563:       PetscFESetNumComponents(sub, Nc);
1564:       PetscFESetUp(sub);
1565:       PetscFESetQuadrature(sub, subq);
1566:       fe->subspaces[height-1] = sub;
1567:     }
1568:     *subfe = fe->subspaces[height-1];
1569:   } else {
1570:     *subfe = NULL;
1571:   }
1572:   return(0);
1573: }

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

1580:   Collective on fem

1582:   Input Parameter:
1583: . fe - The initial PetscFE

1585:   Output Parameter:
1586: . feRef - The refined PetscFE

1588:   Level: advanced

1590: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1591: @*/
1592: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1593: {
1594:   PetscSpace       P, Pref;
1595:   PetscDualSpace   Q, Qref;
1596:   DM               K, Kref;
1597:   PetscQuadrature  q, qref;
1598:   const PetscReal *v0, *jac;
1599:   PetscInt         numComp, numSubelements;
1600:   PetscErrorCode   ierr;

1603:   PetscFEGetBasisSpace(fe, &P);
1604:   PetscFEGetDualSpace(fe, &Q);
1605:   PetscFEGetQuadrature(fe, &q);
1606:   PetscDualSpaceGetDM(Q, &K);
1607:   /* Create space */
1608:   PetscObjectReference((PetscObject) P);
1609:   Pref = P;
1610:   /* Create dual space */
1611:   PetscDualSpaceDuplicate(Q, &Qref);
1612:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1613:   PetscDualSpaceSetDM(Qref, Kref);
1614:   DMDestroy(&Kref);
1615:   PetscDualSpaceSetUp(Qref);
1616:   /* Create element */
1617:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1618:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1619:   PetscFESetBasisSpace(*feRef, Pref);
1620:   PetscFESetDualSpace(*feRef, Qref);
1621:   PetscFEGetNumComponents(fe,    &numComp);
1622:   PetscFESetNumComponents(*feRef, numComp);
1623:   PetscFESetUp(*feRef);
1624:   PetscSpaceDestroy(&Pref);
1625:   PetscDualSpaceDestroy(&Qref);
1626:   /* Create quadrature */
1627:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1628:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1629:   PetscFESetQuadrature(*feRef, qref);
1630:   PetscQuadratureDestroy(&qref);
1631:   return(0);
1632: }

1634: /*@C
1635:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1637:   Collective

1639:   Input Parameters:
1640: + comm      - The MPI comm
1641: . dim       - The spatial dimension
1642: . Nc        - The number of components
1643: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1644: . prefix    - The options prefix, or NULL
1645: - qorder    - The quadrature order

1647:   Output Parameter:
1648: . fem - The PetscFE object

1650:   Level: beginner

1652: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1653: @*/
1654: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1655: {
1656:   PetscQuadrature q, fq;
1657:   DM              K;
1658:   PetscSpace      P;
1659:   PetscDualSpace  Q;
1660:   PetscInt        order, quadPointsPerEdge;
1661:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1662:   PetscErrorCode  ierr;

1665:   /* Create space */
1666:   PetscSpaceCreate(comm, &P);
1667:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1668:   PetscSpacePolynomialSetTensor(P, tensor);
1669:   PetscSpaceSetNumComponents(P, Nc);
1670:   PetscSpaceSetNumVariables(P, dim);
1671:   PetscSpaceSetFromOptions(P);
1672:   PetscSpaceSetUp(P);
1673:   PetscSpaceGetDegree(P, &order, NULL);
1674:   PetscSpacePolynomialGetTensor(P, &tensor);
1675:   /* Create dual space */
1676:   PetscDualSpaceCreate(comm, &Q);
1677:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1678:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1679:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1680:   PetscDualSpaceSetDM(Q, K);
1681:   DMDestroy(&K);
1682:   PetscDualSpaceSetNumComponents(Q, Nc);
1683:   PetscDualSpaceSetOrder(Q, order);
1684:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1685:   PetscDualSpaceSetFromOptions(Q);
1686:   PetscDualSpaceSetUp(Q);
1687:   /* Create element */
1688:   PetscFECreate(comm, fem);
1689:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1690:   PetscFESetBasisSpace(*fem, P);
1691:   PetscFESetDualSpace(*fem, Q);
1692:   PetscFESetNumComponents(*fem, Nc);
1693:   PetscFESetFromOptions(*fem);
1694:   PetscFESetUp(*fem);
1695:   PetscSpaceDestroy(&P);
1696:   PetscDualSpaceDestroy(&Q);
1697:   /* Create quadrature (with specified order if given) */
1698:   qorder = qorder >= 0 ? qorder : order;
1699:   PetscObjectOptionsBegin((PetscObject)*fem);
1700:   PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1701:   PetscOptionsEnd();
1702:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1703:   if (isSimplex) {
1704:     PetscDTGaussJacobiQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1705:     PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1706:   } else {
1707:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1708:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1709:   }
1710:   PetscFESetQuadrature(*fem, q);
1711:   PetscFESetFaceQuadrature(*fem, fq);
1712:   PetscQuadratureDestroy(&q);
1713:   PetscQuadratureDestroy(&fq);
1714:   return(0);
1715: }

1717: /*@C
1718:   PetscFESetName - Names the FE and its subobjects

1720:   Not collective

1722:   Input Parameters:
1723: + fe   - The PetscFE
1724: - name - The name

1726:   Level: intermediate

1728: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1729: @*/
1730: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1731: {
1732:   PetscSpace     P;
1733:   PetscDualSpace Q;

1737:   PetscFEGetBasisSpace(fe, &P);
1738:   PetscFEGetDualSpace(fe, &Q);
1739:   PetscObjectSetName((PetscObject) fe, name);
1740:   PetscObjectSetName((PetscObject) P,  name);
1741:   PetscObjectSetName((PetscObject) Q,  name);
1742:   return(0);
1743: }

1745: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1746: {
1747:   PetscInt       dOffset = 0, fOffset = 0, f;

1750:   for (f = 0; f < Nf; ++f) {
1751:     PetscFE          fe;
1752:     const PetscInt   Nbf = Nb[f], Ncf = Nc[f];
1753:     const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
1754:     const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
1755:     PetscInt         b, c, d;

1757:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
1758:     for (c = 0; c < Ncf; ++c)     u[fOffset+c] = 0.0;
1759:     for (d = 0; d < dim*Ncf; ++d) u_x[fOffset*dim+d] = 0.0;
1760:     for (b = 0; b < Nbf; ++b) {
1761:       for (c = 0; c < Ncf; ++c) {
1762:         const PetscInt cidx = b*Ncf+c;

1764:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
1765:         for (d = 0; d < dim; ++d) u_x[(fOffset+c)*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
1766:       }
1767:     }
1768:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
1769:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*dim]);
1770:     if (u_t) {
1771:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
1772:       for (b = 0; b < Nbf; ++b) {
1773:         for (c = 0; c < Ncf; ++c) {
1774:           const PetscInt cidx = b*Ncf+c;

1776:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
1777:         }
1778:       }
1779:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
1780:     }
1781:     fOffset += Ncf;
1782:     dOffset += Nbf;
1783:   }
1784:   return 0;
1785: }

1787: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
1788: {
1789:   PetscFE        fe;
1790:   PetscReal     *faceBasis;
1791:   PetscInt       Nb, Nc, b, c;

1794:   if (!prob) return 0;
1795:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1796:   PetscFEGetDimension(fe, &Nb);
1797:   PetscFEGetNumComponents(fe, &Nc);
1798:   PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
1799:   for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
1800:   for (b = 0; b < Nb; ++b) {
1801:     for (c = 0; c < Nc; ++c) {
1802:       const PetscInt cidx = b*Nc+c;

1804:       u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
1805:     }
1806:   }
1807:   return 0;
1808: }

1810: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
1811: {
1812:   PetscInt       q, b, c, d;

1815:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
1816:   for (q = 0; q < Nq; ++q) {
1817:     for (b = 0; b < Nb; ++b) {
1818:       for (c = 0; c < Nc; ++c) {
1819:         const PetscInt bcidx = b*Nc+c;

1821:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
1822:         for (d = 0; d < dim; ++d) tmpBasisDer[bcidx*dim+d] = basisDer[q*Nb*Nc*dim+bcidx*dim+d];
1823:       }
1824:     }
1825:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
1826:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
1827:     for (b = 0; b < Nb; ++b) {
1828:       for (c = 0; c < Nc; ++c) {
1829:         const PetscInt bcidx = b*Nc+c;
1830:         const PetscInt qcidx = q*Nc+c;

1832:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
1833:         for (d = 0; d < dim; ++d) elemVec[b] += tmpBasisDer[bcidx*dim+d]*f1[qcidx*dim+d];
1834:       }
1835:     }
1836:   }
1837:   return(0);
1838: }

1840: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt dim, PetscInt NbI, PetscInt NcI, const PetscReal basisI[], const PetscReal basisDerI[], PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscInt NbJ, PetscInt NcJ, const PetscReal basisJ[], const PetscReal basisDerJ[], 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[])
1841: {
1842:   PetscInt       f, fc, g, gc, df, dg;

1845:   for (f = 0; f < NbI; ++f) {
1846:     for (fc = 0; fc < NcI; ++fc) {
1847:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

1849:       tmpBasisI[fidx] = basisI[fidx];
1850:       for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dim+df] = basisDerI[fidx*dim+df];
1851:     }
1852:   }
1853:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
1854:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
1855:   for (g = 0; g < NbJ; ++g) {
1856:     for (gc = 0; gc < NcJ; ++gc) {
1857:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

1859:       tmpBasisJ[gidx] = basisJ[gidx];
1860:       for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dim+dg] = basisDerJ[gidx*dim+dg];
1861:     }
1862:   }
1863:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
1864:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
1865:   for (f = 0; f < NbI; ++f) {
1866:     for (fc = 0; fc < NcI; ++fc) {
1867:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
1868:       const PetscInt i    = offsetI+f; /* Element matrix row */
1869:       for (g = 0; g < NbJ; ++g) {
1870:         for (gc = 0; gc < NcJ; ++gc) {
1871:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
1872:           const PetscInt j    = offsetJ+g; /* Element matrix column */
1873:           const PetscInt fOff = eOffset+i*totDim+j;

1875:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
1876:           for (df = 0; df < dim; ++df) {
1877:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dim+df]*tmpBasisDerJ[gidx*dim+df];
1878:             elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g2[(fc*NcJ+gc)*dim+df]*tmpBasisJ[gidx];
1879:             for (dg = 0; dg < dim; ++dg) {
1880:               elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g3[((fc*NcJ+gc)*dim+df)*dim+dg]*tmpBasisDerJ[gidx*dim+dg];
1881:             }
1882:           }
1883:         }
1884:       }
1885:     }
1886:   }
1887:   return(0);
1888: }

1890: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
1891: {
1892:   PetscDualSpace  dsp;
1893:   DM              dm;
1894:   PetscQuadrature quadDef;
1895:   PetscInt        dim, cdim, Nq;
1896:   PetscErrorCode  ierr;

1899:   PetscFEGetDualSpace(fe, &dsp);
1900:   PetscDualSpaceGetDM(dsp, &dm);
1901:   DMGetDimension(dm, &dim);
1902:   DMGetCoordinateDim(dm, &cdim);
1903:   PetscFEGetQuadrature(fe, &quadDef);
1904:   quad = quad ? quad : quadDef;
1905:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
1906:   PetscMalloc1(Nq*cdim,      &cgeom->v);
1907:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
1908:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
1909:   PetscMalloc1(Nq,           &cgeom->detJ);
1910:   cgeom->dim       = dim;
1911:   cgeom->dimEmbed  = cdim;
1912:   cgeom->numCells  = 1;
1913:   cgeom->numPoints = Nq;
1914:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
1915:   return(0);
1916: }

1918: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
1919: {

1923:   PetscFree(cgeom->v);
1924:   PetscFree(cgeom->J);
1925:   PetscFree(cgeom->invJ);
1926:   PetscFree(cgeom->detJ);
1927:   return(0);
1928: }