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 `PetscFEType`

 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:   Sample usage:
 68: .vb
 69:     PetscFERegister("my_fe", MyPetscFECreate);
 70: .ve

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

 82:   Level: advanced

 84:   Note:
 85:   `PetscFERegister()` may be called multiple times to add several user-defined `PetscFE`s

 87: .seealso: `PetscFE`, `PetscFEType`, `PetscFERegisterAll()`, `PetscFERegisterDestroy()`
 88: @*/
 89: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
 90: {
 91:   PetscFunctionListAdd(&PetscFEList, sname, function);
 92:   return 0;
 93: }

 95: /*@C
 96:   PetscFESetType - Builds a particular `PetscFE`

 98:   Collective on fem

100:   Input Parameters:
101: + fem  - The `PetscFE` object
102: - name - The kind of FEM space

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

107:   Level: intermediate

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

117:   PetscObjectTypeCompare((PetscObject)fem, name, &match);
118:   if (match) return 0;

120:   if (!PetscFERegisterAllCalled) PetscFERegisterAll();
121:   PetscFunctionListFind(PetscFEList, name, &r);

124:   PetscTryTypeMethod(fem, destroy);
125:   fem->ops->destroy = NULL;

127:   (*r)(fem);
128:   PetscObjectChangeTypeName((PetscObject)fem, name);
129:   return 0;
130: }

132: /*@C
133:   PetscFEGetType - Gets the `PetscFEType` (as a string) from the `PetscFE` object.

135:   Not Collective

137:   Input Parameter:
138: . fem  - The `PetscFE`

140:   Output Parameter:
141: . name - The `PetscFEType` name

143:   Level: intermediate

145: .seealso: `PetscFEType`, `PetscFE`, `PetscFESetType()`, `PetscFECreate()`
146: @*/
147: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
148: {
151:   if (!PetscFERegisterAllCalled) PetscFERegisterAll();
152:   *name = ((PetscObject)fem)->type_name;
153:   return 0;
154: }

156: /*@C
157:    PetscFEViewFromOptions - View from a `PetscFE` based on values in the options database

159:    Collective on A

161:    Input Parameters:
162: +  A - the `PetscFE` object
163: .  obj - Optional object that provides the options prefix
164: -  name - command line option name

166:    Level: intermediate

168: .seealso: `PetscFE`, `PetscFEView()`, `PetscObjectViewFromOptions()`, `PetscFECreate()`
169: @*/
170: PetscErrorCode PetscFEViewFromOptions(PetscFE A, PetscObject obj, const char name[])
171: {
173:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
174:   return 0;
175: }

177: /*@C
178:   PetscFEView - Views a `PetscFE`

180:   Collective on fem

182:   Input Parameters:
183: + fem - the `PetscFE` object to view
184: - viewer   - the viewer

186:   Level: beginner

188: .seealso: `PetscFE`, `PetscViewer`, `PetscFEDestroy()`, `PetscFEViewFromOptions()`
189: @*/
190: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
191: {
192:   PetscBool iascii;

196:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fem), &viewer);
197:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
198:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
199:   PetscTryTypeMethod(fem, view, viewer);
200:   return 0;
201: }

203: /*@
204:   PetscFESetFromOptions - sets parameters in a `PetscFE` from the options database

206:   Collective on fem

208:   Input Parameter:
209: . fem - the `PetscFE` object to set options for

211:   Options Database Keys:
212: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
213: - -petscfe_num_batches - the number of cell batches to integrate serially

215:   Level: intermediate

217: .seealso: `PetscFEV`, `PetscFEView()`
218: @*/
219: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
220: {
221:   const char *defaultType;
222:   char        name[256];
223:   PetscBool   flg;

226:   if (!((PetscObject)fem)->type_name) {
227:     defaultType = PETSCFEBASIC;
228:   } else {
229:     defaultType = ((PetscObject)fem)->type_name;
230:   }
231:   if (!PetscFERegisterAllCalled) PetscFERegisterAll();

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

250: /*@C
251:   PetscFESetUp - Construct data structures for the `PetscFE` after the `PetscFEType` has been set

253:   Collective on fem

255:   Input Parameter:
256: . fem - the `PetscFE` object to setup

258:   Level: intermediate

260: .seealso: `PetscFE`, `PetscFEView()`, `PetscFEDestroy()`
261: @*/
262: PetscErrorCode PetscFESetUp(PetscFE fem)
263: {
265:   if (fem->setupcalled) return 0;
266:   PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
267:   fem->setupcalled = PETSC_TRUE;
268:   PetscTryTypeMethod(fem, setup);
269:   PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
270:   return 0;
271: }

273: /*@
274:   PetscFEDestroy - Destroys a `PetscFE` object

276:   Collective on fem

278:   Input Parameter:
279: . fem - the `PetscFE` object to destroy

281:   Level: beginner

283: .seealso: `PetscFE`, `PetscFEView()`
284: @*/
285: PetscErrorCode PetscFEDestroy(PetscFE *fem)
286: {
287:   if (!*fem) return 0;

290:   if (--((PetscObject)(*fem))->refct > 0) {
291:     *fem = NULL;
292:     return 0;
293:   }
294:   ((PetscObject)(*fem))->refct = 0;

296:   if ((*fem)->subspaces) {
297:     PetscInt dim, d;

299:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
300:     for (d = 0; d < dim; ++d) PetscFEDestroy(&(*fem)->subspaces[d]);
301:   }
302:   PetscFree((*fem)->subspaces);
303:   PetscFree((*fem)->invV);
304:   PetscTabulationDestroy(&(*fem)->T);
305:   PetscTabulationDestroy(&(*fem)->Tf);
306:   PetscTabulationDestroy(&(*fem)->Tc);
307:   PetscSpaceDestroy(&(*fem)->basisSpace);
308:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
309:   PetscQuadratureDestroy(&(*fem)->quadrature);
310:   PetscQuadratureDestroy(&(*fem)->faceQuadrature);
311: #ifdef PETSC_HAVE_LIBCEED
312:   CeedBasisDestroy(&(*fem)->ceedBasis);
313:   CeedDestroy(&(*fem)->ceed);
314: #endif

316:   PetscTryTypeMethod((*fem), destroy);
317:   PetscHeaderDestroy(fem);
318:   return 0;
319: }

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

324:   Collective

326:   Input Parameter:
327: . comm - The communicator for the `PetscFE` object

329:   Output Parameter:
330: . fem - The `PetscFE` object

332:   Level: beginner

334: .seealso: `PetscFE`, `PetscFEType`, `PetscFESetType()`, `PetscFECreateDefault()`, `PETSCFEGALERKIN`
335: @*/
336: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
337: {
338:   PetscFE f;

341:   PetscCitationsRegister(FECitation, &FEcite);
342:   *fem = NULL;
343:   PetscFEInitializePackage();

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

347:   f->basisSpace    = NULL;
348:   f->dualSpace     = NULL;
349:   f->numComponents = 1;
350:   f->subspaces     = NULL;
351:   f->invV          = NULL;
352:   f->T             = NULL;
353:   f->Tf            = NULL;
354:   f->Tc            = NULL;
355:   PetscArrayzero(&f->quadrature, 1);
356:   PetscArrayzero(&f->faceQuadrature, 1);
357:   f->blockSize  = 0;
358:   f->numBlocks  = 1;
359:   f->batchSize  = 0;
360:   f->numBatches = 1;

362:   *fem = f;
363:   return 0;
364: }

366: /*@
367:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

369:   Not collective

371:   Input Parameter:
372: . fem - The `PetscFE` object

374:   Output Parameter:
375: . dim - The spatial dimension

377:   Level: intermediate

379: .seealso: `PetscFE`, `PetscFECreate()`
380: @*/
381: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
382: {
383:   DM dm;

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

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

395:   Not collective

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

401:   Level: intermediate

403: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`, `PetscFEGetNumComponents()`
404: @*/
405: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
406: {
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: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`, `PetscFEGetNumComponents()`
426: @*/
427: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
428: {
431:   *comp = fem->numComponents;
432:   return 0;
433: }

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

438:   Not collective

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

447:   Level: intermediate

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

461: /*@
462:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

464:   Not collective

466:   Input Parameter:
467: . fem - The `PetscFE` object

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

475:   Level: intermediate

477: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFESetTileSizes()`
478: @*/
479: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
480: {
486:   if (blockSize) *blockSize = fem->blockSize;
487:   if (numBlocks) *numBlocks = fem->numBlocks;
488:   if (batchSize) *batchSize = fem->batchSize;
489:   if (numBatches) *numBatches = fem->numBatches;
490:   return 0;
491: }

493: /*@
494:   PetscFEGetBasisSpace - Returns the `PetscSpace` used for the approximation of the solution for the `PetscFE`

496:   Not collective

498:   Input Parameter:
499: . fem - The `PetscFE` object

501:   Output Parameter:
502: . sp - The `PetscSpace` object

504:   Level: intermediate

506: .seealso: `PetscFE`, `PetscSpace`, `PetscFECreate()`
507: @*/
508: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
509: {
512:   *sp = fem->basisSpace;
513:   return 0;
514: }

516: /*@
517:   PetscFESetBasisSpace - Sets the `PetscSpace` used for the approximation of the solution

519:   Not collective

521:   Input Parameters:
522: + fem - The `PetscFE` object
523: - sp - The `PetscSpace` object

525:   Level: intermediate

527:   Developer Note:
528:   There is `PetscFESetBasisSpace()` but the `PetscFESetDualSpace()`, likely the Basis is unneeded in the function name

530: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetDualSpace()`
531: @*/
532: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
533: {
536:   PetscSpaceDestroy(&fem->basisSpace);
537:   fem->basisSpace = sp;
538:   PetscObjectReference((PetscObject)fem->basisSpace);
539:   return 0;
540: }

542: /*@
543:   PetscFEGetDualSpace - Returns the `PetscDualSpace` used to define the inner product for a `PetscFE`

545:   Not collective

547:   Input Parameter:
548: . fem - The `PetscFE` object

550:   Output Parameter:
551: . sp - The `PetscDualSpace` object

553:   Level: intermediate

555: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
556: @*/
557: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
558: {
561:   *sp = fem->dualSpace;
562:   return 0;
563: }

565: /*@
566:   PetscFESetDualSpace - Sets the `PetscDualSpace` used to define the inner product

568:   Not collective

570:   Input Parameters:
571: + fem - The `PetscFE` object
572: - sp - The `PetscDualSpace` object

574:   Level: intermediate

576: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetBasisSpace()`
577: @*/
578: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
579: {
582:   PetscDualSpaceDestroy(&fem->dualSpace);
583:   fem->dualSpace = sp;
584:   PetscObjectReference((PetscObject)fem->dualSpace);
585:   return 0;
586: }

588: /*@
589:   PetscFEGetQuadrature - Returns the `PetscQuadrature` used to calculate inner products

591:   Not collective

593:   Input Parameter:
594: . fem - The `PetscFE` object

596:   Output Parameter:
597: . q - The `PetscQuadrature` object

599:   Level: intermediate

601: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`
602: @*/
603: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
604: {
607:   *q = fem->quadrature;
608:   return 0;
609: }

611: /*@
612:   PetscFESetQuadrature - Sets the `PetscQuadrature` used to calculate inner products

614:   Not collective

616:   Input Parameters:
617: + fem - The `PetscFE` object
618: - q - The `PetscQuadrature` object

620:   Level: intermediate

622: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFEGetFaceQuadrature()`
623: @*/
624: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
625: {
626:   PetscInt Nc, qNc;

629:   if (q == fem->quadrature) return 0;
630:   PetscFEGetNumComponents(fem, &Nc);
631:   PetscQuadratureGetNumComponents(q, &qNc);
633:   PetscTabulationDestroy(&fem->T);
634:   PetscTabulationDestroy(&fem->Tc);
635:   PetscObjectReference((PetscObject)q);
636:   PetscQuadratureDestroy(&fem->quadrature);
637:   fem->quadrature = q;
638:   return 0;
639: }

641: /*@
642:   PetscFEGetFaceQuadrature - Returns the `PetscQuadrature` used to calculate inner products on faces

644:   Not collective

646:   Input Parameter:
647: . fem - The `PetscFE` object

649:   Output Parameter:
650: . q - The `PetscQuadrature` object

652:   Level: intermediate

654:   Developer Note:
655:   There is a special face quadrature but not edge, likely this API would benefit from a refactorization

657: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
658: @*/
659: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
660: {
663:   *q = fem->faceQuadrature;
664:   return 0;
665: }

667: /*@
668:   PetscFESetFaceQuadrature - Sets the `PetscQuadrature` used to calculate inner products on faces

670:   Not collective

672:   Input Parameters:
673: + fem - The `PetscFE` object
674: - q - The `PetscQuadrature` object

676:   Level: intermediate

678: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
679: @*/
680: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
681: {
682:   PetscInt Nc, qNc;

685:   PetscFEGetNumComponents(fem, &Nc);
686:   PetscQuadratureGetNumComponents(q, &qNc);
688:   PetscTabulationDestroy(&fem->Tf);
689:   PetscQuadratureDestroy(&fem->faceQuadrature);
690:   fem->faceQuadrature = q;
691:   PetscObjectReference((PetscObject)q);
692:   return 0;
693: }

695: /*@
696:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature to a new `PetscFE`

698:   Not collective

700:   Input Parameters:
701: + sfe - The `PetscFE` source for the quadratures
702: - tfe - The `PetscFE` target for the quadratures

704:   Level: intermediate

706: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
707: @*/
708: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
709: {
710:   PetscQuadrature q;

714:   PetscFEGetQuadrature(sfe, &q);
715:   PetscFESetQuadrature(tfe, q);
716:   PetscFEGetFaceQuadrature(sfe, &q);
717:   PetscFESetFaceQuadrature(tfe, q);
718:   return 0;
719: }

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

724:   Not collective

726:   Input Parameter:
727: . fem - The `PetscFE` object

729:   Output Parameter:
730: . numDof - Array with the number of dofs per dimension

732:   Level: intermediate

734: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
735: @*/
736: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
737: {
740:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
741:   return 0;
742: }

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

747:   Not collective

749:   Input Parameters:
750: + fem - The `PetscFE` object
751: - k   - The highest derivative we need to tabulate, very often 1

753:   Output Parameter:
754: . T - The basis function values and derivatives at quadrature points

756:   Level: intermediate

758:   Note:
759: .vb
760:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
761:   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
762:   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
763: .ve

765: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
766: @*/
767: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
768: {
769:   PetscInt         npoints;
770:   const PetscReal *points;

774:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
775:   if (!fem->T) PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);
777:   *T = fem->T;
778:   return 0;
779: }

781: /*@C
782:   PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell

784:   Not collective

786:   Input Parameters:
787: + fem - The `PetscFE` object
788: - k   - The highest derivative we need to tabulate, very often 1

790:   Output Parameters:
791: . Tf - The basis function values and derivatives at face quadrature points

793:   Level: intermediate

795:   Note:
796: .vb
797:   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
798:   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
799:   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
800: .ve

802: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
803: @*/
804: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
805: {
808:   if (!fem->Tf) {
809:     const PetscReal  xi0[3] = {-1., -1., -1.};
810:     PetscReal        v0[3], J[9], detJ;
811:     PetscQuadrature  fq;
812:     PetscDualSpace   sp;
813:     DM               dm;
814:     const PetscInt  *faces;
815:     PetscInt         dim, numFaces, f, npoints, q;
816:     const PetscReal *points;
817:     PetscReal       *facePoints;

819:     PetscFEGetDualSpace(fem, &sp);
820:     PetscDualSpaceGetDM(sp, &dm);
821:     DMGetDimension(dm, &dim);
822:     DMPlexGetConeSize(dm, 0, &numFaces);
823:     DMPlexGetCone(dm, 0, &faces);
824:     PetscFEGetFaceQuadrature(fem, &fq);
825:     if (fq) {
826:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
827:       PetscMalloc1(numFaces * npoints * dim, &facePoints);
828:       for (f = 0; f < numFaces; ++f) {
829:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
830:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * npoints + q) * dim]);
831:       }
832:       PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
833:       PetscFree(facePoints);
834:     }
835:   }
837:   *Tf = fem->Tf;
838:   return 0;
839: }

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

844:   Not collective

846:   Input Parameter:
847: . fem - The `PetscFE` object

849:   Output Parameters:
850: . Tc - The basis function values at face centroid points

852:   Level: intermediate

854:   Note:
855: .vb
856:   T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
857: .ve

859: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
860: @*/
861: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
862: {
865:   if (!fem->Tc) {
866:     PetscDualSpace  sp;
867:     DM              dm;
868:     const PetscInt *cone;
869:     PetscReal      *centroids;
870:     PetscInt        dim, numFaces, f;

872:     PetscFEGetDualSpace(fem, &sp);
873:     PetscDualSpaceGetDM(sp, &dm);
874:     DMGetDimension(dm, &dim);
875:     DMPlexGetConeSize(dm, 0, &numFaces);
876:     DMPlexGetCone(dm, 0, &cone);
877:     PetscMalloc1(numFaces * dim, &centroids);
878:     for (f = 0; f < numFaces; ++f) DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f * dim], NULL);
879:     PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
880:     PetscFree(centroids);
881:   }
882:   *Tc = fem->Tc;
883:   return 0;
884: }

886: /*@C
887:   PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

889:   Not collective

891:   Input Parameters:
892: + fem     - The `PetscFE` object
893: . nrepl   - The number of replicas
894: . npoints - The number of tabulation points in a replica
895: . points  - The tabulation point coordinates
896: - K       - The number of derivatives calculated

898:   Output Parameter:
899: . T - The basis function values and derivatives at tabulation points

901:   Level: intermediate

903:   Note:
904: .vb
905:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
906:   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
907:   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

909: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
910: @*/
911: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
912: {
913:   DM             dm;
914:   PetscDualSpace Q;
915:   PetscInt       Nb;   /* Dimension of FE space P */
916:   PetscInt       Nc;   /* Field components */
917:   PetscInt       cdim; /* Reference coordinate dimension */
918:   PetscInt       k;

920:   if (!npoints || !fem->dualSpace || K < 0) {
921:     *T = NULL;
922:     return 0;
923:   }
927:   PetscFEGetDualSpace(fem, &Q);
928:   PetscDualSpaceGetDM(Q, &dm);
929:   DMGetDimension(dm, &cdim);
930:   PetscDualSpaceGetDimension(Q, &Nb);
931:   PetscFEGetNumComponents(fem, &Nc);
932:   PetscMalloc1(1, T);
933:   (*T)->K    = !cdim ? 0 : K;
934:   (*T)->Nr   = nrepl;
935:   (*T)->Np   = npoints;
936:   (*T)->Nb   = Nb;
937:   (*T)->Nc   = Nc;
938:   (*T)->cdim = cdim;
939:   PetscMalloc1((*T)->K + 1, &(*T)->T);
940:   for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
941:   PetscUseTypeMethod(fem, createtabulation, nrepl * npoints, points, K, *T);
942:   return 0;
943: }

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

948:   Not collective

950:   Input Parameters:
951: + fem     - The `PetscFE` object
952: . npoints - The number of tabulation points
953: . points  - The tabulation point coordinates
954: . K       - The number of derivatives calculated
955: - T       - An existing tabulation object with enough allocated space

957:   Output Parameter:
958: . T - The basis function values and derivatives at tabulation points

960:   Level: intermediate

962:   Note:
963: .vb
964:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
965:   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
966:   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
967: .ve

969: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
970: @*/
971: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
972: {
974:   if (!npoints || !fem->dualSpace || K < 0) return 0;
978:   if (PetscDefined(USE_DEBUG)) {
979:     DM             dm;
980:     PetscDualSpace Q;
981:     PetscInt       Nb;   /* Dimension of FE space P */
982:     PetscInt       Nc;   /* Field components */
983:     PetscInt       cdim; /* Reference coordinate dimension */

985:     PetscFEGetDualSpace(fem, &Q);
986:     PetscDualSpaceGetDM(Q, &dm);
987:     DMGetDimension(dm, &cdim);
988:     PetscDualSpaceGetDimension(Q, &Nb);
989:     PetscFEGetNumComponents(fem, &Nc);
994:   }
995:   T->Nr = 1;
996:   T->Np = npoints;
997:   PetscUseTypeMethod(fem, createtabulation, npoints, points, K, T);
998:   return 0;
999: }

1001: /*@C
1002:   PetscTabulationDestroy - Frees memory from the associated tabulation.

1004:   Not collective

1006:   Input Parameter:
1007: . T - The tabulation

1009:   Level: intermediate

1011: .seealso: `PetscTabulation`, `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
1012: @*/
1013: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1014: {
1015:   PetscInt k;

1018:   if (!T || !(*T)) return 0;
1019:   for (k = 0; k <= (*T)->K; ++k) PetscFree((*T)->T[k]);
1020:   PetscFree((*T)->T);
1021:   PetscFree(*T);
1022:   *T = NULL;
1023:   return 0;
1024: }

1026: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1027: {
1028:   PetscSpace      bsp, bsubsp;
1029:   PetscDualSpace  dsp, dsubsp;
1030:   PetscInt        dim, depth, numComp, i, j, coneSize, order;
1031:   PetscFEType     type;
1032:   DM              dm;
1033:   DMLabel         label;
1034:   PetscReal      *xi, *v, *J, detJ;
1035:   const char     *name;
1036:   PetscQuadrature origin, fullQuad, subQuad;

1040:   PetscFEGetBasisSpace(fe, &bsp);
1041:   PetscFEGetDualSpace(fe, &dsp);
1042:   PetscDualSpaceGetDM(dsp, &dm);
1043:   DMGetDimension(dm, &dim);
1044:   DMPlexGetDepthLabel(dm, &label);
1045:   DMLabelGetValue(label, refPoint, &depth);
1046:   PetscCalloc1(depth, &xi);
1047:   PetscMalloc1(dim, &v);
1048:   PetscMalloc1(dim * dim, &J);
1049:   for (i = 0; i < depth; i++) xi[i] = 0.;
1050:   PetscQuadratureCreate(PETSC_COMM_SELF, &origin);
1051:   PetscQuadratureSetData(origin, depth, 0, 1, xi, NULL);
1052:   DMPlexComputeCellGeometryFEM(dm, refPoint, origin, v, J, NULL, &detJ);
1053:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1054:   for (i = 1; i < dim; i++) {
1055:     for (j = 0; j < depth; j++) J[i * depth + j] = J[i * dim + j];
1056:   }
1057:   PetscQuadratureDestroy(&origin);
1058:   PetscDualSpaceGetPointSubspace(dsp, refPoint, &dsubsp);
1059:   PetscSpaceCreateSubspace(bsp, dsubsp, v, J, NULL, NULL, PETSC_OWN_POINTER, &bsubsp);
1060:   PetscSpaceSetUp(bsubsp);
1061:   PetscFECreate(PetscObjectComm((PetscObject)fe), trFE);
1062:   PetscFEGetType(fe, &type);
1063:   PetscFESetType(*trFE, type);
1064:   PetscFEGetNumComponents(fe, &numComp);
1065:   PetscFESetNumComponents(*trFE, numComp);
1066:   PetscFESetBasisSpace(*trFE, bsubsp);
1067:   PetscFESetDualSpace(*trFE, dsubsp);
1068:   PetscObjectGetName((PetscObject)fe, &name);
1069:   if (name) PetscFESetName(*trFE, name);
1070:   PetscFEGetQuadrature(fe, &fullQuad);
1071:   PetscQuadratureGetOrder(fullQuad, &order);
1072:   DMPlexGetConeSize(dm, refPoint, &coneSize);
1073:   if (coneSize == 2 * depth) PetscDTGaussTensorQuadrature(depth, 1, (order + 1) / 2, -1., 1., &subQuad);
1074:   else PetscDTStroudConicalQuadrature(depth, 1, (order + 1) / 2, -1., 1., &subQuad);
1075:   PetscFESetQuadrature(*trFE, subQuad);
1076:   PetscFESetUp(*trFE);
1077:   PetscQuadratureDestroy(&subQuad);
1078:   PetscSpaceDestroy(&bsubsp);
1079:   return 0;
1080: }

1082: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1083: {
1084:   PetscInt       hStart, hEnd;
1085:   PetscDualSpace dsp;
1086:   DM             dm;

1090:   *trFE = NULL;
1091:   PetscFEGetDualSpace(fe, &dsp);
1092:   PetscDualSpaceGetDM(dsp, &dm);
1093:   DMPlexGetHeightStratum(dm, height, &hStart, &hEnd);
1094:   if (hEnd <= hStart) return 0;
1095:   PetscFECreatePointTrace(fe, hStart, trFE);
1096:   return 0;
1097: }

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

1102:   Not collective

1104:   Input Parameter:
1105: . fe - The `PetscFE`

1107:   Output Parameter:
1108: . dim - The dimension

1110:   Level: intermediate

1112: .seealso: `PetscFE`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1113: @*/
1114: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1115: {
1118:   PetscTryTypeMethod(fem, getdimension, dim);
1119:   return 0;
1120: }

1122: /*@C
1123:   PetscFEPushforward - Map the reference element function to real space

1125:   Input Parameters:
1126: + fe     - The `PetscFE`
1127: . fegeom - The cell geometry
1128: . Nv     - The number of function values
1129: - vals   - The function values

1131:   Output Parameter:
1132: . vals   - The transformed function values

1134:   Level: advanced

1136:   Notes:
1137:   This just forwards the call onto `PetscDualSpacePushforward()`.

1139:   It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1141: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscDualSpacePushforward()`
1142: @*/
1143: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1144: {
1146:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1147:   return 0;
1148: }

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

1153:   Input Parameters:
1154: + fe     - The `PetscFE`
1155: . fegeom - The cell geometry
1156: . Nv     - The number of function gradient values
1157: - vals   - The function gradient values

1159:   Output Parameter:
1160: . vals   - The transformed function gradient values

1162:   Level: advanced

1164:   Notes:
1165:   This just forwards the call onto `PetscDualSpacePushforwardGradient()`.

1167:   It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1169: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1170: @*/
1171: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1172: {
1174:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1175:   return 0;
1176: }

1178: /*@C
1179:   PetscFEPushforwardHessian - Map the reference element function Hessian to real space

1181:   Input Parameters:
1182: + fe     - The `PetscFE`
1183: . fegeom - The cell geometry
1184: . Nv     - The number of function Hessian values
1185: - vals   - The function Hessian values

1187:   Output Parameter:
1188: . vals   - The transformed function Hessian values

1190:   Level: advanced

1192:   Notes:
1193:   This just forwards the call onto `PetscDualSpacePushforwardHessian()`.

1195:   It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1197:   Developer Note:
1198:   It is unclear why all these one line convenience routines are desirable

1200: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1201: @*/
1202: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1203: {
1205:   PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1206:   return 0;
1207: }

1209: /*
1210: Purpose: Compute element vector for chunk of elements

1212: Input:
1213:   Sizes:
1214:      Ne:  number of elements
1215:      Nf:  number of fields
1216:      PetscFE
1217:        dim: spatial dimension
1218:        Nb:  number of basis functions
1219:        Nc:  number of field components
1220:        PetscQuadrature
1221:          Nq:  number of quadrature points

1223:   Geometry:
1224:      PetscFEGeom[Ne] possibly *Nq
1225:        PetscReal v0s[dim]
1226:        PetscReal n[dim]
1227:        PetscReal jacobians[dim*dim]
1228:        PetscReal jacobianInverses[dim*dim]
1229:        PetscReal jacobianDeterminants
1230:   FEM:
1231:      PetscFE
1232:        PetscQuadrature
1233:          PetscReal   quadPoints[Nq*dim]
1234:          PetscReal   quadWeights[Nq]
1235:        PetscReal   basis[Nq*Nb*Nc]
1236:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1237:      PetscScalar coefficients[Ne*Nb*Nc]
1238:      PetscScalar elemVec[Ne*Nb*Nc]

1240:   Problem:
1241:      PetscInt f: the active field
1242:      f0, f1

1244:   Work Space:
1245:      PetscFE
1246:        PetscScalar f0[Nq*dim];
1247:        PetscScalar f1[Nq*dim*dim];
1248:        PetscScalar u[Nc];
1249:        PetscScalar gradU[Nc*dim];
1250:        PetscReal   x[dim];
1251:        PetscScalar realSpaceDer[dim];

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

1255: Input:
1256:   Sizes:
1257:      N_cb: Number of serial cell batches

1259:   Geometry:
1260:      PetscReal v0s[Ne*dim]
1261:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1262:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1263:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1264:   FEM:
1265:      static PetscReal   quadPoints[Nq*dim]
1266:      static PetscReal   quadWeights[Nq]
1267:      static PetscReal   basis[Nq*Nb*Nc]
1268:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1269:      PetscScalar coefficients[Ne*Nb*Nc]
1270:      PetscScalar elemVec[Ne*Nb*Nc]

1272: ex62.c:
1273:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1274:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1275:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1276:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1278: ex52.c:
1279:   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)
1280:   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)

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

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

1289: ex52_integrateElementOpenCL.c:
1290: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1291:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1292:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1300:   Not collective

1302:   Input Parameters:
1303: + prob         - The `PetscDS` specifying the discretizations and continuum functions
1304: . field        - The field being integrated
1305: . Ne           - The number of elements in the chunk
1306: . cgeom        - The cell geometry for each cell in the chunk
1307: . coefficients - The array of FEM basis coefficients for the elements
1308: . probAux      - The `PetscDS` specifying the auxiliary discretizations
1309: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1311:   Output Parameter:
1312: . integral     - the integral for this field

1314:   Level: intermediate

1316:   Developer Note:
1317:   The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.

1319: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrateBd()`
1320: @*/
1321: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1322: {
1323:   PetscFE fe;

1326:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1327:   if (fe->ops->integrate) (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);
1328:   return 0;
1329: }

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

1334:   Not collective

1336:   Input Parameters:
1337: + prob         - The `PetscDS` specifying the discretizations and continuum functions
1338: . field        - The field being integrated
1339: . obj_func     - The function to be integrated
1340: . Ne           - The number of elements in the chunk
1341: . fgeom        - The face geometry for each face in the chunk
1342: . coefficients - The array of FEM basis coefficients for the elements
1343: . probAux      - The `PetscDS` specifying the auxiliary discretizations
1344: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1346:   Output Parameter:
1347: . integral     - the integral for this field

1349:   Level: intermediate

1351:   Developer Note:
1352:   The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.

1354: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrate()`
1355: @*/
1356: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, void (*obj_func)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1357: {
1358:   PetscFE fe;

1361:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1362:   if (fe->ops->integratebd) (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);
1363:   return 0;
1364: }

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

1369:   Not collective

1371:   Input Parameters:
1372: + ds           - The PetscDS specifying the discretizations and continuum functions
1373: . key          - The (label+value, field) being integrated
1374: . Ne           - The number of elements in the chunk
1375: . cgeom        - The cell geometry for each cell in the chunk
1376: . coefficients - The array of FEM basis coefficients for the elements
1377: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1378: . probAux      - The PetscDS specifying the auxiliary discretizations
1379: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1380: - t            - The time

1382:   Output Parameter:
1383: . elemVec      - the element residual vectors from each element

1385:   Level: intermediate

1387:   Note:
1388: .vb
1389:   Loop over batch of elements (e):
1390:     Loop over quadrature points (q):
1391:       Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1392:       Call f_0 and f_1
1393:     Loop over element vector entries (f,fc --> i):
1394:       elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1395: .ve

1397: .seealso: `PetscFEIntegrateResidual()`
1398: @*/
1399: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1400: {
1401:   PetscFE fe;

1405:   PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1406:   if (fe->ops->integrateresidual) (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1407:   return 0;
1408: }

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

1413:   Not collective

1415:   Input Parameters:
1416: + ds           - The PetscDS specifying the discretizations and continuum functions
1417: . wf           - The PetscWeakForm object holding the pointwise functions
1418: . key          - The (label+value, field) being integrated
1419: . Ne           - The number of elements in the chunk
1420: . fgeom        - The face geometry for each cell in the chunk
1421: . coefficients - The array of FEM basis coefficients for the elements
1422: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1423: . probAux      - The PetscDS specifying the auxiliary discretizations
1424: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1425: - t            - The time

1427:   Output Parameter:
1428: . elemVec      - the element residual vectors from each element

1430:   Level: intermediate

1432: .seealso: `PetscFEIntegrateResidual()`
1433: @*/
1434: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1435: {
1436:   PetscFE fe;

1439:   PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1440:   if (fe->ops->integratebdresidual) (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1441:   return 0;
1442: }

1444: /*@C
1445:   PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration

1447:   Not collective

1449:   Input Parameters:
1450: + prob         - The PetscDS specifying the discretizations and continuum functions
1451: . key          - The (label+value, field) being integrated
1452: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1453: . Ne           - The number of elements in the chunk
1454: . fgeom        - The face geometry for each cell in the chunk
1455: . coefficients - The array of FEM basis coefficients for the elements
1456: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1457: . probAux      - The PetscDS specifying the auxiliary discretizations
1458: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1459: - t            - The time

1461:   Output Parameter
1462: . elemVec      - the element residual vectors from each element

1464:   Level: developer

1466: .seealso: `PetscFEIntegrateResidual()`
1467: @*/
1468: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1469: {
1470:   PetscFE fe;

1473:   PetscDSGetDiscretization(prob, key.field, (PetscObject *)&fe);
1474:   if (fe->ops->integratehybridresidual) (*fe->ops->integratehybridresidual)(prob, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1475:   return 0;
1476: }

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

1481:   Not collective

1483:   Input Parameters:
1484: + ds           - The PetscDS specifying the discretizations and continuum functions
1485: . jtype        - The type of matrix pointwise functions that should be used
1486: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1487: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1488: . Ne           - The number of elements in the chunk
1489: . cgeom        - The cell geometry for each cell in the chunk
1490: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1491: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1492: . probAux      - The PetscDS specifying the auxiliary discretizations
1493: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1494: . t            - The time
1495: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1497:   Output Parameter:
1498: . elemMat      - the element matrices for the Jacobian from each element

1500:   Level: intermediate

1502:   Note:
1503: .vb
1504:   Loop over batch of elements (e):
1505:     Loop over element matrix entries (f,fc,g,gc --> i,j):
1506:       Loop over quadrature points (q):
1507:         Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1508:           elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1509:                        + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1510:                        + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1511:                        + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1512: .ve

1514: .seealso: `PetscFEIntegrateResidual()`
1515: @*/
1516: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1517: {
1518:   PetscFE  fe;
1519:   PetscInt Nf;

1522:   PetscDSGetNumFields(ds, &Nf);
1523:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1524:   if (fe->ops->integratejacobian) (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1525:   return 0;
1526: }

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

1531:   Not collective

1533:   Input Parameters:
1534: + ds           - The PetscDS specifying the discretizations and continuum functions
1535: . wf           - The PetscWeakForm holding the pointwise functions
1536: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1537: . Ne           - The number of elements in the chunk
1538: . fgeom        - The face geometry for each cell in the chunk
1539: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1540: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1541: . probAux      - The PetscDS specifying the auxiliary discretizations
1542: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1543: . t            - The time
1544: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1546:   Output Parameter:
1547: . elemMat              - the element matrices for the Jacobian from each element

1549:   Level: intermediate

1551:   Note:
1552: .vb
1553:   Loop over batch of elements (e):
1554:     Loop over element matrix entries (f,fc,g,gc --> i,j):
1555:       Loop over quadrature points (q):
1556:         Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1557:           elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1558:                        + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1559:                        + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1560:                        + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1561: .ve

1563: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1564: @*/
1565: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1566: {
1567:   PetscFE  fe;
1568:   PetscInt Nf;

1571:   PetscDSGetNumFields(ds, &Nf);
1572:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1573:   if (fe->ops->integratebdjacobian) (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1574:   return 0;
1575: }

1577: /*@C
1578:   PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration

1580:   Not collective

1582:   Input Parameters:
1583: + ds           - The PetscDS specifying the discretizations and continuum functions
1584: . jtype        - The type of matrix pointwise functions that should be used
1585: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1586: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1587: . Ne           - The number of elements in the chunk
1588: . fgeom        - The face geometry for each cell in the chunk
1589: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1590: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1591: . probAux      - The PetscDS specifying the auxiliary discretizations
1592: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1593: . t            - The time
1594: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1596:   Output Parameter
1597: . elemMat              - the element matrices for the Jacobian from each element

1599:   Level: developer

1601:   Note:
1602: .vb
1603:   Loop over batch of elements (e):
1604:     Loop over element matrix entries (f,fc,g,gc --> i,j):
1605:       Loop over quadrature points (q):
1606:         Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1607:           elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1608:                        + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1609:                        + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1610:                        + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1611: .ve

1613: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1614: @*/
1615: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1616: {
1617:   PetscFE  fe;
1618:   PetscInt Nf;

1621:   PetscDSGetNumFields(ds, &Nf);
1622:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1623:   if (fe->ops->integratehybridjacobian) (*fe->ops->integratehybridjacobian)(ds, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1624:   return 0;
1625: }

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

1630:   Input Parameters:
1631: + fe     - The finite element space
1632: - height - The height of the Plex point

1634:   Output Parameter:
1635: . subfe  - The subspace of this FE space

1637:   Level: advanced

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

1642: .seealso: `PetscFECreateDefault()`
1643: @*/
1644: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1645: {
1646:   PetscSpace      P, subP;
1647:   PetscDualSpace  Q, subQ;
1648:   PetscQuadrature subq;
1649:   PetscFEType     fetype;
1650:   PetscInt        dim, Nc;

1654:   if (height == 0) {
1655:     *subfe = fe;
1656:     return 0;
1657:   }
1658:   PetscFEGetBasisSpace(fe, &P);
1659:   PetscFEGetDualSpace(fe, &Q);
1660:   PetscFEGetNumComponents(fe, &Nc);
1661:   PetscFEGetFaceQuadrature(fe, &subq);
1662:   PetscDualSpaceGetDimension(Q, &dim);
1664:   if (!fe->subspaces) PetscCalloc1(dim, &fe->subspaces);
1665:   if (height <= dim) {
1666:     if (!fe->subspaces[height - 1]) {
1667:       PetscFE     sub = NULL;
1668:       const char *name;

1670:       PetscSpaceGetHeightSubspace(P, height, &subP);
1671:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1672:       if (subQ) {
1673:         PetscFECreate(PetscObjectComm((PetscObject)fe), &sub);
1674:         PetscObjectGetName((PetscObject)fe, &name);
1675:         PetscObjectSetName((PetscObject)sub, name);
1676:         PetscFEGetType(fe, &fetype);
1677:         PetscFESetType(sub, fetype);
1678:         PetscFESetBasisSpace(sub, subP);
1679:         PetscFESetDualSpace(sub, subQ);
1680:         PetscFESetNumComponents(sub, Nc);
1681:         PetscFESetUp(sub);
1682:         PetscFESetQuadrature(sub, subq);
1683:       }
1684:       fe->subspaces[height - 1] = sub;
1685:     }
1686:     *subfe = fe->subspaces[height - 1];
1687:   } else {
1688:     *subfe = NULL;
1689:   }
1690:   return 0;
1691: }

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

1698:   Collective on fem

1700:   Input Parameter:
1701: . fe - The initial PetscFE

1703:   Output Parameter:
1704: . feRef - The refined PetscFE

1706:   Level: advanced

1708: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1709: @*/
1710: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1711: {
1712:   PetscSpace       P, Pref;
1713:   PetscDualSpace   Q, Qref;
1714:   DM               K, Kref;
1715:   PetscQuadrature  q, qref;
1716:   const PetscReal *v0, *jac;
1717:   PetscInt         numComp, numSubelements;
1718:   PetscInt         cStart, cEnd, c;
1719:   PetscDualSpace  *cellSpaces;

1721:   PetscFEGetBasisSpace(fe, &P);
1722:   PetscFEGetDualSpace(fe, &Q);
1723:   PetscFEGetQuadrature(fe, &q);
1724:   PetscDualSpaceGetDM(Q, &K);
1725:   /* Create space */
1726:   PetscObjectReference((PetscObject)P);
1727:   Pref = P;
1728:   /* Create dual space */
1729:   PetscDualSpaceDuplicate(Q, &Qref);
1730:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1731:   DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref);
1732:   PetscDualSpaceSetDM(Qref, Kref);
1733:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1734:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1735:   /* TODO: fix for non-uniform refinement */
1736:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1737:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1738:   PetscFree(cellSpaces);
1739:   DMDestroy(&Kref);
1740:   PetscDualSpaceSetUp(Qref);
1741:   /* Create element */
1742:   PetscFECreate(PetscObjectComm((PetscObject)fe), feRef);
1743:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1744:   PetscFESetBasisSpace(*feRef, Pref);
1745:   PetscFESetDualSpace(*feRef, Qref);
1746:   PetscFEGetNumComponents(fe, &numComp);
1747:   PetscFESetNumComponents(*feRef, numComp);
1748:   PetscFESetUp(*feRef);
1749:   PetscSpaceDestroy(&Pref);
1750:   PetscDualSpaceDestroy(&Qref);
1751:   /* Create quadrature */
1752:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1753:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1754:   PetscFESetQuadrature(*feRef, qref);
1755:   PetscQuadratureDestroy(&qref);
1756:   return 0;
1757: }

1759: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1760: {
1761:   PetscSpace     P;
1762:   PetscDualSpace Q;
1763:   DM             K;
1764:   DMPolytopeType ct;
1765:   PetscInt       degree;
1766:   char           name[64];

1768:   PetscFEGetBasisSpace(fe, &P);
1769:   PetscSpaceGetDegree(P, &degree, NULL);
1770:   PetscFEGetDualSpace(fe, &Q);
1771:   PetscDualSpaceGetDM(Q, &K);
1772:   DMPlexGetCellType(K, 0, &ct);
1773:   switch (ct) {
1774:   case DM_POLYTOPE_SEGMENT:
1775:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1776:   case DM_POLYTOPE_QUADRILATERAL:
1777:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1778:   case DM_POLYTOPE_HEXAHEDRON:
1779:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1780:     PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree);
1781:     break;
1782:   case DM_POLYTOPE_TRIANGLE:
1783:   case DM_POLYTOPE_TETRAHEDRON:
1784:     PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree);
1785:     break;
1786:   case DM_POLYTOPE_TRI_PRISM:
1787:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1788:     PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree);
1789:     break;
1790:   default:
1791:     PetscSNPrintf(name, sizeof(name), "FE");
1792:   }
1793:   PetscFESetName(fe, name);
1794:   return 0;
1795: }

1797: static PetscErrorCode PetscFECreateDefaultQuadrature_Private(PetscInt dim, DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq)
1798: {
1799:   const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1);

1801:   switch (ct) {
1802:   case DM_POLYTOPE_SEGMENT:
1803:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1804:   case DM_POLYTOPE_QUADRILATERAL:
1805:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1806:   case DM_POLYTOPE_HEXAHEDRON:
1807:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1808:     PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1809:     PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1810:     break;
1811:   case DM_POLYTOPE_TRIANGLE:
1812:   case DM_POLYTOPE_TETRAHEDRON:
1813:     PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1814:     PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1815:     break;
1816:   case DM_POLYTOPE_TRI_PRISM:
1817:   case DM_POLYTOPE_TRI_PRISM_TENSOR: {
1818:     PetscQuadrature q1, q2;

1820:     PetscDTStroudConicalQuadrature(2, 1, quadPointsPerEdge, -1.0, 1.0, &q1);
1821:     PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2);
1822:     PetscDTTensorQuadratureCreate(q1, q2, q);
1823:     PetscQuadratureDestroy(&q1);
1824:     PetscQuadratureDestroy(&q2);
1825:   }
1826:     PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1827:     /* TODO Need separate quadratures for each face */
1828:     break;
1829:   default:
1830:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
1831:   }
1832:   return 0;
1833: }

1835: /*@
1836:   PetscFECreateFromSpaces - Create a `PetscFE` from the basis and dual spaces

1838:   Collective

1840:   Input Parameters:
1841: + P  - The basis space
1842: . Q  - The dual space
1843: . q  - The cell quadrature
1844: - fq - The face quadrature

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

1849:   Level: beginner

1851:   Note:
1852:   The `PetscFE` takes ownership of these spaces by calling destroy on each. They should not be used after this call, and for borrowed references from `PetscFEGetSpace()` and the like, the caller must use `PetscObjectReference` before this call.

1854: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`,
1855:           `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1856: @*/
1857: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1858: {
1859:   PetscInt    Nc;
1860:   const char *prefix;

1862:   PetscFECreate(PetscObjectComm((PetscObject)P), fem);
1863:   PetscObjectGetOptionsPrefix((PetscObject)P, &prefix);
1864:   PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix);
1865:   PetscFESetType(*fem, PETSCFEBASIC);
1866:   PetscFESetBasisSpace(*fem, P);
1867:   PetscFESetDualSpace(*fem, Q);
1868:   PetscSpaceGetNumComponents(P, &Nc);
1869:   PetscFESetNumComponents(*fem, Nc);
1870:   PetscFESetUp(*fem);
1871:   PetscSpaceDestroy(&P);
1872:   PetscDualSpaceDestroy(&Q);
1873:   PetscFESetQuadrature(*fem, q);
1874:   PetscFESetFaceQuadrature(*fem, fq);
1875:   PetscQuadratureDestroy(&q);
1876:   PetscQuadratureDestroy(&fq);
1877:   PetscFESetDefaultName_Private(*fem);
1878:   return 0;
1879: }

1881: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1882: {
1883:   DM              K;
1884:   PetscSpace      P;
1885:   PetscDualSpace  Q;
1886:   PetscQuadrature q, fq;
1887:   PetscBool       tensor;

1891:   switch (ct) {
1892:   case DM_POLYTOPE_SEGMENT:
1893:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1894:   case DM_POLYTOPE_QUADRILATERAL:
1895:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1896:   case DM_POLYTOPE_HEXAHEDRON:
1897:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1898:     tensor = PETSC_TRUE;
1899:     break;
1900:   default:
1901:     tensor = PETSC_FALSE;
1902:   }
1903:   /* Create space */
1904:   PetscSpaceCreate(comm, &P);
1905:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1906:   PetscObjectSetOptionsPrefix((PetscObject)P, prefix);
1907:   PetscSpacePolynomialSetTensor(P, tensor);
1908:   PetscSpaceSetNumComponents(P, Nc);
1909:   PetscSpaceSetNumVariables(P, dim);
1910:   if (degree >= 0) {
1911:     PetscSpaceSetDegree(P, degree, PETSC_DETERMINE);
1912:     if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1913:       PetscSpace Pend, Pside;

1915:       PetscSpaceCreate(comm, &Pend);
1916:       PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL);
1917:       PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE);
1918:       PetscSpaceSetNumComponents(Pend, Nc);
1919:       PetscSpaceSetNumVariables(Pend, dim - 1);
1920:       PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE);
1921:       PetscSpaceCreate(comm, &Pside);
1922:       PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL);
1923:       PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE);
1924:       PetscSpaceSetNumComponents(Pside, 1);
1925:       PetscSpaceSetNumVariables(Pside, 1);
1926:       PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE);
1927:       PetscSpaceSetType(P, PETSCSPACETENSOR);
1928:       PetscSpaceTensorSetNumSubspaces(P, 2);
1929:       PetscSpaceTensorSetSubspace(P, 0, Pend);
1930:       PetscSpaceTensorSetSubspace(P, 1, Pside);
1931:       PetscSpaceDestroy(&Pend);
1932:       PetscSpaceDestroy(&Pside);
1933:     }
1934:   }
1935:   if (setFromOptions) PetscSpaceSetFromOptions(P);
1936:   PetscSpaceSetUp(P);
1937:   PetscSpaceGetDegree(P, &degree, NULL);
1938:   PetscSpacePolynomialGetTensor(P, &tensor);
1939:   PetscSpaceGetNumComponents(P, &Nc);
1940:   /* Create dual space */
1941:   PetscDualSpaceCreate(comm, &Q);
1942:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1943:   PetscObjectSetOptionsPrefix((PetscObject)Q, prefix);
1944:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1945:   PetscDualSpaceSetDM(Q, K);
1946:   DMDestroy(&K);
1947:   PetscDualSpaceSetNumComponents(Q, Nc);
1948:   PetscDualSpaceSetOrder(Q, degree);
1949:   /* TODO For some reason, we need a tensor dualspace with wedges */
1950:   PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE);
1951:   if (setFromOptions) PetscDualSpaceSetFromOptions(Q);
1952:   PetscDualSpaceSetUp(Q);
1953:   /* Create quadrature */
1954:   qorder = qorder >= 0 ? qorder : degree;
1955:   if (setFromOptions) {
1956:     PetscObjectOptionsBegin((PetscObject)P);
1957:     PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0);
1958:     PetscOptionsEnd();
1959:   }
1960:   PetscFECreateDefaultQuadrature_Private(dim, ct, qorder, &q, &fq);
1961:   /* Create finite element */
1962:   PetscFECreateFromSpaces(P, Q, q, fq, fem);
1963:   if (setFromOptions) PetscFESetFromOptions(*fem);
1964:   return 0;
1965: }

1967: /*@C
1968:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1970:   Collective

1972:   Input Parameters:
1973: + comm      - The MPI comm
1974: . dim       - The spatial dimension
1975: . Nc        - The number of components
1976: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1977: . prefix    - The options prefix, or NULL
1978: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1980:   Output Parameter:
1981: . fem - The PetscFE object

1983:   Level: beginner

1985:   Note:
1986:   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.

1988: .seealso: `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1989: @*/
1990: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1991: {
1992:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1993:   return 0;
1994: }

1996: /*@C
1997:   PetscFECreateByCell - Create a PetscFE for basic FEM computation

1999:   Collective

2001:   Input Parameters:
2002: + comm   - The MPI comm
2003: . dim    - The spatial dimension
2004: . Nc     - The number of components
2005: . ct     - The celltype of the reference cell
2006: . prefix - The options prefix, or NULL
2007: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

2009:   Output Parameter:
2010: . fem - The PetscFE object

2012:   Level: beginner

2014:   Note:
2015:   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.

2017: .seealso: `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2018: @*/
2019: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
2020: {
2021:   PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
2022:   return 0;
2023: }

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

2028:   Collective

2030:   Input Parameters:
2031: + comm      - The MPI comm
2032: . dim       - The spatial dimension
2033: . Nc        - The number of components
2034: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2035: . k         - The degree k of the space
2036: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

2038:   Output Parameter:
2039: . fem       - The PetscFE object

2041:   Level: beginner

2043:   Note:
2044:   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.

2046: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2047: @*/
2048: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2049: {
2050:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem);
2051:   return 0;
2052: }

2054: /*@
2055:   PetscFECreateLagrangeByCell - Create a PetscFE for the basic Lagrange space of degree k

2057:   Collective

2059:   Input Parameters:
2060: + comm      - The MPI comm
2061: . dim       - The spatial dimension
2062: . Nc        - The number of components
2063: . ct        - The celltype of the reference cell
2064: . k         - The degree k of the space
2065: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

2067:   Output Parameter:
2068: . fem       - The PetscFE object

2070:   Level: beginner

2072:   Note:
2073:   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.

2075: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2076: @*/
2077: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2078: {
2079:   PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem);
2080:   return 0;
2081: }

2083: /*@C
2084:   PetscFESetName - Names the FE and its subobjects

2086:   Not collective

2088:   Input Parameters:
2089: + fe   - The PetscFE
2090: - name - The name

2092:   Level: intermediate

2094: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2095: @*/
2096: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2097: {
2098:   PetscSpace     P;
2099:   PetscDualSpace Q;

2101:   PetscFEGetBasisSpace(fe, &P);
2102:   PetscFEGetDualSpace(fe, &Q);
2103:   PetscObjectSetName((PetscObject)fe, name);
2104:   PetscObjectSetName((PetscObject)P, name);
2105:   PetscObjectSetName((PetscObject)Q, name);
2106:   return 0;
2107: }

2109: 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[])
2110: {
2111:   PetscInt dOffset = 0, fOffset = 0, f, g;

2113:   for (f = 0; f < Nf; ++f) {
2114:     PetscFE          fe;
2115:     const PetscInt   k       = ds->jetDegree[f];
2116:     const PetscInt   cdim    = T[f]->cdim;
2117:     const PetscInt   Nq      = T[f]->Np;
2118:     const PetscInt   Nbf     = T[f]->Nb;
2119:     const PetscInt   Ncf     = T[f]->Nc;
2120:     const PetscReal *Bq      = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2121:     const PetscReal *Dq      = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2122:     const PetscReal *Hq      = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2123:     PetscInt         hOffset = 0, b, c, d;

2125:     PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
2126:     for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2127:     for (d = 0; d < cdim * Ncf; ++d) u_x[fOffset * cdim + d] = 0.0;
2128:     for (b = 0; b < Nbf; ++b) {
2129:       for (c = 0; c < Ncf; ++c) {
2130:         const PetscInt cidx = b * Ncf + c;

2132:         u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2133:         for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * cdim + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2134:       }
2135:     }
2136:     if (k > 1) {
2137:       for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * cdim;
2138:       for (d = 0; d < cdim * cdim * Ncf; ++d) u_x[hOffset + fOffset * cdim * cdim + d] = 0.0;
2139:       for (b = 0; b < Nbf; ++b) {
2140:         for (c = 0; c < Ncf; ++c) {
2141:           const PetscInt cidx = b * Ncf + c;

2143:           for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * cdim * cdim + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2144:         }
2145:       }
2146:       PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * cdim * cdim]);
2147:     }
2148:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2149:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * cdim]);
2150:     if (u_t) {
2151:       for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2152:       for (b = 0; b < Nbf; ++b) {
2153:         for (c = 0; c < Ncf; ++c) {
2154:           const PetscInt cidx = b * Ncf + c;

2156:           u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2157:         }
2158:       }
2159:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2160:     }
2161:     fOffset += Ncf;
2162:     dOffset += Nbf;
2163:   }
2164:   return 0;
2165: }

2167: 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[])
2168: {
2169:   PetscInt dOffset = 0, fOffset = 0, f, g;

2171:   /* f is the field number in the DS, g is the field number in u[] */
2172:   for (f = 0, g = 0; f < Nf; ++f) {
2173:     PetscFE          fe  = (PetscFE)ds->disc[f];
2174:     const PetscInt   dEt = T[f]->cdim;
2175:     const PetscInt   dE  = fegeom->dimEmbed;
2176:     const PetscInt   Nq  = T[f]->Np;
2177:     const PetscInt   Nbf = T[f]->Nb;
2178:     const PetscInt   Ncf = T[f]->Nc;
2179:     const PetscReal *Bq  = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2180:     const PetscReal *Dq  = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2181:     PetscBool        isCohesive;
2182:     PetscInt         Ns, s;

2184:     if (!T[f]) continue;
2185:     PetscDSGetCohesive(ds, f, &isCohesive);
2186:     Ns = isCohesive ? 1 : 2;
2187:     for (s = 0; s < Ns; ++s, ++g) {
2188:       PetscInt b, c, d;

2190:       for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2191:       for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2192:       for (b = 0; b < Nbf; ++b) {
2193:         for (c = 0; c < Ncf; ++c) {
2194:           const PetscInt cidx = b * Ncf + c;

2196:           u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2197:           for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2198:         }
2199:       }
2200:       PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2201:       PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]);
2202:       if (u_t) {
2203:         for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2204:         for (b = 0; b < Nbf; ++b) {
2205:           for (c = 0; c < Ncf; ++c) {
2206:             const PetscInt cidx = b * Ncf + c;

2208:             u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2209:           }
2210:         }
2211:         PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2212:       }
2213:       fOffset += Ncf;
2214:       dOffset += Nbf;
2215:     }
2216:   }
2217:   return 0;
2218: }

2220: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2221: {
2222:   PetscFE         fe;
2223:   PetscTabulation Tc;
2224:   PetscInt        b, c;

2226:   if (!prob) return 0;
2227:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
2228:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2229:   {
2230:     const PetscReal *faceBasis = Tc->T[0];
2231:     const PetscInt   Nb        = Tc->Nb;
2232:     const PetscInt   Nc        = Tc->Nc;

2234:     for (c = 0; c < Nc; ++c) u[c] = 0.0;
2235:     for (b = 0; b < Nb; ++b) {
2236:       for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2237:     }
2238:   }
2239:   return 0;
2240: }

2242: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2243: {
2244:   PetscFEGeom      pgeom;
2245:   const PetscInt   dEt      = T->cdim;
2246:   const PetscInt   dE       = fegeom->dimEmbed;
2247:   const PetscInt   Nq       = T->Np;
2248:   const PetscInt   Nb       = T->Nb;
2249:   const PetscInt   Nc       = T->Nc;
2250:   const PetscReal *basis    = &T->T[0][r * Nq * Nb * Nc];
2251:   const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt];
2252:   PetscInt         q, b, c, d;

2254:   for (q = 0; q < Nq; ++q) {
2255:     for (b = 0; b < Nb; ++b) {
2256:       for (c = 0; c < Nc; ++c) {
2257:         const PetscInt bcidx = b * Nc + c;

2259:         tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2260:         for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d];
2261:         for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0;
2262:       }
2263:     }
2264:     PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);
2265:     PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);
2266:     PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer);
2267:     for (b = 0; b < Nb; ++b) {
2268:       for (c = 0; c < Nc; ++c) {
2269:         const PetscInt bcidx = b * Nc + c;
2270:         const PetscInt qcidx = q * Nc + c;

2272:         elemVec[b] += tmpBasis[bcidx] * f0[qcidx];
2273:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2274:       }
2275:     }
2276:   }
2277:   return (0);
2278: }

2280: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2281: {
2282:   const PetscInt   dE       = T->cdim;
2283:   const PetscInt   Nq       = T->Np;
2284:   const PetscInt   Nb       = T->Nb;
2285:   const PetscInt   Nc       = T->Nc;
2286:   const PetscReal *basis    = &T->T[0][r * Nq * Nb * Nc];
2287:   const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];
2288:   PetscInt         q, b, c, d;

2290:   for (q = 0; q < Nq; ++q) {
2291:     for (b = 0; b < Nb; ++b) {
2292:       for (c = 0; c < Nc; ++c) {
2293:         const PetscInt bcidx = b * Nc + c;

2295:         tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2296:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2297:       }
2298:     }
2299:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2300:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2301:     for (b = 0; b < Nb; ++b) {
2302:       for (c = 0; c < Nc; ++c) {
2303:         const PetscInt bcidx = b * Nc + c;
2304:         const PetscInt qcidx = q * Nc + c;

2306:         elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2307:         for (d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2308:       }
2309:     }
2310:   }
2311:   return (0);
2312: }

2314: 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[])
2315: {
2316:   const PetscInt   dE        = TI->cdim;
2317:   const PetscInt   NqI       = TI->Np;
2318:   const PetscInt   NbI       = TI->Nb;
2319:   const PetscInt   NcI       = TI->Nc;
2320:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2321:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2322:   const PetscInt   NqJ       = TJ->Np;
2323:   const PetscInt   NbJ       = TJ->Nb;
2324:   const PetscInt   NcJ       = TJ->Nc;
2325:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2326:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2327:   PetscInt         f, fc, g, gc, df, dg;

2329:   for (f = 0; f < NbI; ++f) {
2330:     for (fc = 0; fc < NcI; ++fc) {
2331:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2333:       tmpBasisI[fidx] = basisI[fidx];
2334:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2335:     }
2336:   }
2337:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2338:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2339:   for (g = 0; g < NbJ; ++g) {
2340:     for (gc = 0; gc < NcJ; ++gc) {
2341:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2343:       tmpBasisJ[gidx] = basisJ[gidx];
2344:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2345:     }
2346:   }
2347:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2348:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2349:   for (f = 0; f < NbI; ++f) {
2350:     for (fc = 0; fc < NcI; ++fc) {
2351:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2352:       const PetscInt i    = offsetI + f;  /* Element matrix row */
2353:       for (g = 0; g < NbJ; ++g) {
2354:         for (gc = 0; gc < NcJ; ++gc) {
2355:           const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2356:           const PetscInt j    = offsetJ + g;  /* Element matrix column */
2357:           const PetscInt fOff = eOffset + i * totDim + j;

2359:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2360:           for (df = 0; df < dE; ++df) {
2361:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2362:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2363:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2364:           }
2365:         }
2366:       }
2367:     }
2368:   }
2369:   return (0);
2370: }

2372: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt s, 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[])
2373: {
2374:   const PetscInt   dE        = TI->cdim;
2375:   const PetscInt   NqI       = TI->Np;
2376:   const PetscInt   NbI       = TI->Nb;
2377:   const PetscInt   NcI       = TI->Nc;
2378:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2379:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2380:   const PetscInt   NqJ       = TJ->Np;
2381:   const PetscInt   NbJ       = TJ->Nb;
2382:   const PetscInt   NcJ       = TJ->Nc;
2383:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2384:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2385:   const PetscInt   so        = isHybridI ? 0 : s;
2386:   const PetscInt   to        = isHybridJ ? 0 : s;
2387:   PetscInt         f, fc, g, gc, df, dg;

2389:   for (f = 0; f < NbI; ++f) {
2390:     for (fc = 0; fc < NcI; ++fc) {
2391:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2393:       tmpBasisI[fidx] = basisI[fidx];
2394:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2395:     }
2396:   }
2397:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2398:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2399:   for (g = 0; g < NbJ; ++g) {
2400:     for (gc = 0; gc < NcJ; ++gc) {
2401:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2403:       tmpBasisJ[gidx] = basisJ[gidx];
2404:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2405:     }
2406:   }
2407:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2408:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2409:   for (f = 0; f < NbI; ++f) {
2410:     for (fc = 0; fc < NcI; ++fc) {
2411:       const PetscInt fidx = f * NcI + fc;           /* Test function basis index */
2412:       const PetscInt i    = offsetI + NbI * so + f; /* Element matrix row */
2413:       for (g = 0; g < NbJ; ++g) {
2414:         for (gc = 0; gc < NcJ; ++gc) {
2415:           const PetscInt gidx = g * NcJ + gc;           /* Trial function basis index */
2416:           const PetscInt j    = offsetJ + NbJ * to + g; /* Element matrix column */
2417:           const PetscInt fOff = eOffset + i * totDim + j;

2419:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2420:           for (df = 0; df < dE; ++df) {
2421:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2422:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2423:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2424:           }
2425:         }
2426:       }
2427:     }
2428:   }
2429:   return (0);
2430: }

2432: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2433: {
2434:   PetscDualSpace  dsp;
2435:   DM              dm;
2436:   PetscQuadrature quadDef;
2437:   PetscInt        dim, cdim, Nq;

2439:   PetscFEGetDualSpace(fe, &dsp);
2440:   PetscDualSpaceGetDM(dsp, &dm);
2441:   DMGetDimension(dm, &dim);
2442:   DMGetCoordinateDim(dm, &cdim);
2443:   PetscFEGetQuadrature(fe, &quadDef);
2444:   quad = quad ? quad : quadDef;
2445:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2446:   PetscMalloc1(Nq * cdim, &cgeom->v);
2447:   PetscMalloc1(Nq * cdim * cdim, &cgeom->J);
2448:   PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ);
2449:   PetscMalloc1(Nq, &cgeom->detJ);
2450:   cgeom->dim       = dim;
2451:   cgeom->dimEmbed  = cdim;
2452:   cgeom->numCells  = 1;
2453:   cgeom->numPoints = Nq;
2454:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2455:   return 0;
2456: }

2458: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2459: {
2460:   PetscFree(cgeom->v);
2461:   PetscFree(cgeom->J);
2462:   PetscFree(cgeom->invJ);
2463:   PetscFree(cgeom->detJ);
2464:   return 0;
2465: }