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: + sname    - The name of a new user-defined creation routine
 65: - function - The creation routine

 67:   Example 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:   PetscFunctionBegin;
 92:   PetscCall(PetscFunctionListAdd(&PetscFEList, sname, function));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

 99:   Collective

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

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

108:   Level: intermediate

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

117:   PetscFunctionBegin;
119:   PetscCall(PetscObjectTypeCompare((PetscObject)fem, name, &match));
120:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

122:   if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
123:   PetscCall(PetscFunctionListFind(PetscFEList, name, &r));
124:   PetscCheck(r, PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);

126:   PetscTryTypeMethod(fem, destroy);
127:   fem->ops->destroy = NULL;

129:   PetscCall((*r)(fem));
130:   PetscCall(PetscObjectChangeTypeName((PetscObject)fem, name));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

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

137:   Not Collective

139:   Input Parameter:
140: . fem - The `PetscFE`

142:   Output Parameter:
143: . name - The `PetscFEType` name

145:   Level: intermediate

147: .seealso: `PetscFEType`, `PetscFE`, `PetscFESetType()`, `PetscFECreate()`
148: @*/
149: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
150: {
151:   PetscFunctionBegin;
153:   PetscAssertPointer(name, 2);
154:   if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
155:   *name = ((PetscObject)fem)->type_name;
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

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

162:   Collective

164:   Input Parameters:
165: + A    - the `PetscFE` object
166: . obj  - Optional object that provides the options prefix
167: - name - command line option name

169:   Level: intermediate

171: .seealso: `PetscFE`, `PetscFEView()`, `PetscObjectViewFromOptions()`, `PetscFECreate()`
172: @*/
173: PetscErrorCode PetscFEViewFromOptions(PetscFE A, PetscObject obj, const char name[])
174: {
175:   PetscFunctionBegin;
177:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

181: /*@C
182:   PetscFEView - Views a `PetscFE`

184:   Collective

186:   Input Parameters:
187: + fem    - the `PetscFE` object to view
188: - viewer - the viewer

190:   Level: beginner

192: .seealso: `PetscFE`, `PetscViewer`, `PetscFEDestroy()`, `PetscFEViewFromOptions()`
193: @*/
194: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
195: {
196:   PetscBool iascii;

198:   PetscFunctionBegin;
201:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fem), &viewer));
202:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer));
203:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
204:   PetscTryTypeMethod(fem, view, viewer);
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

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

211:   Collective

213:   Input Parameter:
214: . fem - the `PetscFE` object to set options for

216:   Options Database Keys:
217: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
218: - -petscfe_num_batches - the number of cell batches to integrate serially

220:   Level: intermediate

222: .seealso: `PetscFEV`, `PetscFEView()`
223: @*/
224: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
225: {
226:   const char *defaultType;
227:   char        name[256];
228:   PetscBool   flg;

230:   PetscFunctionBegin;
232:   if (!((PetscObject)fem)->type_name) {
233:     defaultType = PETSCFEBASIC;
234:   } else {
235:     defaultType = ((PetscObject)fem)->type_name;
236:   }
237:   if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());

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

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

259:   Collective

261:   Input Parameter:
262: . fem - the `PetscFE` object to setup

264:   Level: intermediate

266: .seealso: `PetscFE`, `PetscFEView()`, `PetscFEDestroy()`
267: @*/
268: PetscErrorCode PetscFESetUp(PetscFE fem)
269: {
270:   PetscFunctionBegin;
272:   if (fem->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
273:   PetscCall(PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0));
274:   fem->setupcalled = PETSC_TRUE;
275:   PetscTryTypeMethod(fem, setup);
276:   PetscCall(PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0));
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: /*@
281:   PetscFEDestroy - Destroys a `PetscFE` object

283:   Collective

285:   Input Parameter:
286: . fem - the `PetscFE` object to destroy

288:   Level: beginner

290: .seealso: `PetscFE`, `PetscFEView()`
291: @*/
292: PetscErrorCode PetscFEDestroy(PetscFE *fem)
293: {
294:   PetscFunctionBegin;
295:   if (!*fem) PetscFunctionReturn(PETSC_SUCCESS);

298:   if (--((PetscObject)*fem)->refct > 0) {
299:     *fem = NULL;
300:     PetscFunctionReturn(PETSC_SUCCESS);
301:   }
302:   ((PetscObject)*fem)->refct = 0;

304:   if ((*fem)->subspaces) {
305:     PetscInt dim, d;

307:     PetscCall(PetscDualSpaceGetDimension((*fem)->dualSpace, &dim));
308:     for (d = 0; d < dim; ++d) PetscCall(PetscFEDestroy(&(*fem)->subspaces[d]));
309:   }
310:   PetscCall(PetscFree((*fem)->subspaces));
311:   PetscCall(PetscFree((*fem)->invV));
312:   PetscCall(PetscTabulationDestroy(&(*fem)->T));
313:   PetscCall(PetscTabulationDestroy(&(*fem)->Tf));
314:   PetscCall(PetscTabulationDestroy(&(*fem)->Tc));
315:   PetscCall(PetscSpaceDestroy(&(*fem)->basisSpace));
316:   PetscCall(PetscDualSpaceDestroy(&(*fem)->dualSpace));
317:   PetscCall(PetscQuadratureDestroy(&(*fem)->quadrature));
318:   PetscCall(PetscQuadratureDestroy(&(*fem)->faceQuadrature));
319: #ifdef PETSC_HAVE_LIBCEED
320:   PetscCallCEED(CeedBasisDestroy(&(*fem)->ceedBasis));
321:   PetscCallCEED(CeedDestroy(&(*fem)->ceed));
322: #endif

324:   PetscTryTypeMethod(*fem, destroy);
325:   PetscCall(PetscHeaderDestroy(fem));
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

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

332:   Collective

334:   Input Parameter:
335: . comm - The communicator for the `PetscFE` object

337:   Output Parameter:
338: . fem - The `PetscFE` object

340:   Level: beginner

342: .seealso: `PetscFE`, `PetscFEType`, `PetscFESetType()`, `PetscFECreateDefault()`, `PETSCFEGALERKIN`
343: @*/
344: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
345: {
346:   PetscFE f;

348:   PetscFunctionBegin;
349:   PetscAssertPointer(fem, 2);
350:   PetscCall(PetscCitationsRegister(FECitation, &FEcite));
351:   *fem = NULL;
352:   PetscCall(PetscFEInitializePackage());

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

356:   f->basisSpace    = NULL;
357:   f->dualSpace     = NULL;
358:   f->numComponents = 1;
359:   f->subspaces     = NULL;
360:   f->invV          = NULL;
361:   f->T             = NULL;
362:   f->Tf            = NULL;
363:   f->Tc            = NULL;
364:   PetscCall(PetscArrayzero(&f->quadrature, 1));
365:   PetscCall(PetscArrayzero(&f->faceQuadrature, 1));
366:   f->blockSize  = 0;
367:   f->numBlocks  = 1;
368:   f->batchSize  = 0;
369:   f->numBatches = 1;

371:   *fem = f;
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: /*@
376:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

378:   Not Collective

380:   Input Parameter:
381: . fem - The `PetscFE` object

383:   Output Parameter:
384: . dim - The spatial dimension

386:   Level: intermediate

388: .seealso: `PetscFE`, `PetscFECreate()`
389: @*/
390: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
391: {
392:   DM dm;

394:   PetscFunctionBegin;
396:   PetscAssertPointer(dim, 2);
397:   PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
398:   PetscCall(DMGetDimension(dm, dim));
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: /*@
403:   PetscFESetNumComponents - Sets the number of field components in the element

405:   Not Collective

407:   Input Parameters:
408: + fem  - The `PetscFE` object
409: - comp - The number of field components

411:   Level: intermediate

413: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`, `PetscFEGetNumComponents()`
414: @*/
415: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
416: {
417:   PetscFunctionBegin;
419:   fem->numComponents = comp;
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: /*@
424:   PetscFEGetNumComponents - Returns the number of components in the element

426:   Not Collective

428:   Input Parameter:
429: . fem - The `PetscFE` object

431:   Output Parameter:
432: . comp - The number of field components

434:   Level: intermediate

436: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`
437: @*/
438: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
439: {
440:   PetscFunctionBegin;
442:   PetscAssertPointer(comp, 2);
443:   *comp = fem->numComponents;
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: /*@
448:   PetscFESetTileSizes - Sets the tile sizes for evaluation

450:   Not Collective

452:   Input Parameters:
453: + fem        - The `PetscFE` object
454: . blockSize  - The number of elements in a block
455: . numBlocks  - The number of blocks in a batch
456: . batchSize  - The number of elements in a batch
457: - numBatches - The number of batches in a chunk

459:   Level: intermediate

461: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetTileSizes()`
462: @*/
463: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
464: {
465:   PetscFunctionBegin;
467:   fem->blockSize  = blockSize;
468:   fem->numBlocks  = numBlocks;
469:   fem->batchSize  = batchSize;
470:   fem->numBatches = numBatches;
471:   PetscFunctionReturn(PETSC_SUCCESS);
472: }

474: /*@
475:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

477:   Not Collective

479:   Input Parameter:
480: . fem - The `PetscFE` object

482:   Output Parameters:
483: + blockSize  - The number of elements in a block
484: . numBlocks  - The number of blocks in a batch
485: . batchSize  - The number of elements in a batch
486: - numBatches - The number of batches in a chunk

488:   Level: intermediate

490: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFESetTileSizes()`
491: @*/
492: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
493: {
494:   PetscFunctionBegin;
496:   if (blockSize) PetscAssertPointer(blockSize, 2);
497:   if (numBlocks) PetscAssertPointer(numBlocks, 3);
498:   if (batchSize) PetscAssertPointer(batchSize, 4);
499:   if (numBatches) PetscAssertPointer(numBatches, 5);
500:   if (blockSize) *blockSize = fem->blockSize;
501:   if (numBlocks) *numBlocks = fem->numBlocks;
502:   if (batchSize) *batchSize = fem->batchSize;
503:   if (numBatches) *numBatches = fem->numBatches;
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

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

510:   Not Collective

512:   Input Parameter:
513: . fem - The `PetscFE` object

515:   Output Parameter:
516: . sp - The `PetscSpace` object

518:   Level: intermediate

520: .seealso: `PetscFE`, `PetscSpace`, `PetscFECreate()`
521: @*/
522: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
523: {
524:   PetscFunctionBegin;
526:   PetscAssertPointer(sp, 2);
527:   *sp = fem->basisSpace;
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: /*@
532:   PetscFESetBasisSpace - Sets the `PetscSpace` used for the approximation of the solution

534:   Not Collective

536:   Input Parameters:
537: + fem - The `PetscFE` object
538: - sp  - The `PetscSpace` object

540:   Level: intermediate

542:   Developer Notes:
543:   There is `PetscFESetBasisSpace()` but the `PetscFESetDualSpace()`, likely the Basis is unneeded in the function name

545: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetDualSpace()`
546: @*/
547: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
548: {
549:   PetscFunctionBegin;
552:   PetscCall(PetscSpaceDestroy(&fem->basisSpace));
553:   fem->basisSpace = sp;
554:   PetscCall(PetscObjectReference((PetscObject)fem->basisSpace));
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }

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

561:   Not Collective

563:   Input Parameter:
564: . fem - The `PetscFE` object

566:   Output Parameter:
567: . sp - The `PetscDualSpace` object

569:   Level: intermediate

571: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
572: @*/
573: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
574: {
575:   PetscFunctionBegin;
577:   PetscAssertPointer(sp, 2);
578:   *sp = fem->dualSpace;
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: /*@
583:   PetscFESetDualSpace - Sets the `PetscDualSpace` used to define the inner product

585:   Not Collective

587:   Input Parameters:
588: + fem - The `PetscFE` object
589: - sp  - The `PetscDualSpace` object

591:   Level: intermediate

593: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetBasisSpace()`
594: @*/
595: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
596: {
597:   PetscFunctionBegin;
600:   PetscCall(PetscDualSpaceDestroy(&fem->dualSpace));
601:   fem->dualSpace = sp;
602:   PetscCall(PetscObjectReference((PetscObject)fem->dualSpace));
603:   PetscFunctionReturn(PETSC_SUCCESS);
604: }

606: /*@
607:   PetscFEGetQuadrature - Returns the `PetscQuadrature` used to calculate inner products

609:   Not Collective

611:   Input Parameter:
612: . fem - The `PetscFE` object

614:   Output Parameter:
615: . q - The `PetscQuadrature` object

617:   Level: intermediate

619: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`
620: @*/
621: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
622: {
623:   PetscFunctionBegin;
625:   PetscAssertPointer(q, 2);
626:   *q = fem->quadrature;
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: /*@
631:   PetscFESetQuadrature - Sets the `PetscQuadrature` used to calculate inner products

633:   Not Collective

635:   Input Parameters:
636: + fem - The `PetscFE` object
637: - q   - The `PetscQuadrature` object

639:   Level: intermediate

641: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFEGetFaceQuadrature()`
642: @*/
643: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
644: {
645:   PetscInt Nc, qNc;

647:   PetscFunctionBegin;
649:   if (q == fem->quadrature) PetscFunctionReturn(PETSC_SUCCESS);
650:   PetscCall(PetscFEGetNumComponents(fem, &Nc));
651:   PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
652:   PetscCheck(!(qNc != 1) || !(Nc != qNc), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_SIZ, "FE components %" PetscInt_FMT " != Quadrature components %" PetscInt_FMT " and non-scalar quadrature", Nc, qNc);
653:   PetscCall(PetscTabulationDestroy(&fem->T));
654:   PetscCall(PetscTabulationDestroy(&fem->Tc));
655:   PetscCall(PetscObjectReference((PetscObject)q));
656:   PetscCall(PetscQuadratureDestroy(&fem->quadrature));
657:   fem->quadrature = q;
658:   PetscFunctionReturn(PETSC_SUCCESS);
659: }

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

664:   Not Collective

666:   Input Parameter:
667: . fem - The `PetscFE` object

669:   Output Parameter:
670: . q - The `PetscQuadrature` object

672:   Level: intermediate

674:   Developer Notes:
675:   There is a special face quadrature but not edge, likely this API would benefit from a refactorization

677: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
678: @*/
679: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
680: {
681:   PetscFunctionBegin;
683:   PetscAssertPointer(q, 2);
684:   *q = fem->faceQuadrature;
685:   PetscFunctionReturn(PETSC_SUCCESS);
686: }

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

691:   Not Collective

693:   Input Parameters:
694: + fem - The `PetscFE` object
695: - q   - The `PetscQuadrature` object

697:   Level: intermediate

699: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`
700: @*/
701: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
702: {
703:   PetscInt Nc, qNc;

705:   PetscFunctionBegin;
707:   if (q == fem->faceQuadrature) PetscFunctionReturn(PETSC_SUCCESS);
708:   PetscCall(PetscFEGetNumComponents(fem, &Nc));
709:   PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
710:   PetscCheck(!(qNc != 1) || !(Nc != qNc), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_SIZ, "FE components %" PetscInt_FMT " != Quadrature components %" PetscInt_FMT " and non-scalar quadrature", Nc, qNc);
711:   PetscCall(PetscTabulationDestroy(&fem->Tf));
712:   PetscCall(PetscObjectReference((PetscObject)q));
713:   PetscCall(PetscQuadratureDestroy(&fem->faceQuadrature));
714:   fem->faceQuadrature = q;
715:   PetscFunctionReturn(PETSC_SUCCESS);
716: }

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

721:   Not Collective

723:   Input Parameters:
724: + sfe - The `PetscFE` source for the quadratures
725: - tfe - The `PetscFE` target for the quadratures

727:   Level: intermediate

729: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
730: @*/
731: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
732: {
733:   PetscQuadrature q;

735:   PetscFunctionBegin;
738:   PetscCall(PetscFEGetQuadrature(sfe, &q));
739:   PetscCall(PetscFESetQuadrature(tfe, q));
740:   PetscCall(PetscFEGetFaceQuadrature(sfe, &q));
741:   PetscCall(PetscFESetFaceQuadrature(tfe, q));
742:   PetscFunctionReturn(PETSC_SUCCESS);
743: }

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

748:   Not Collective

750:   Input Parameter:
751: . fem - The `PetscFE` object

753:   Output Parameter:
754: . numDof - Array of length `dim` with the number of dofs in each dimension

756:   Level: intermediate

758: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
759: @*/
760: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt *numDof[])
761: {
762:   PetscFunctionBegin;
764:   PetscAssertPointer(numDof, 2);
765:   PetscCall(PetscDualSpaceGetNumDof(fem->dualSpace, numDof));
766:   PetscFunctionReturn(PETSC_SUCCESS);
767: }

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

772:   Not Collective

774:   Input Parameters:
775: + fem - The `PetscFE` object
776: - k   - The highest derivative we need to tabulate, very often 1

778:   Output Parameter:
779: . T - The basis function values and derivatives at quadrature points

781:   Level: intermediate

783:   Note:
784: .vb
785:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
786:   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
787:   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
788: .ve

790: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
791: @*/
792: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
793: {
794:   PetscInt         npoints;
795:   const PetscReal *points;

797:   PetscFunctionBegin;
799:   PetscAssertPointer(T, 3);
800:   PetscCall(PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL));
801:   if (!fem->T) PetscCall(PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T));
802:   PetscCheck(!fem->T || k <= fem->T->K || (!fem->T->cdim && !fem->T->K), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %" PetscInt_FMT " derivatives, but only tabulated %" PetscInt_FMT, k, fem->T->K);
803:   *T = fem->T;
804:   PetscFunctionReturn(PETSC_SUCCESS);
805: }

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

810:   Not Collective

812:   Input Parameters:
813: + fem - The `PetscFE` object
814: - k   - The highest derivative we need to tabulate, very often 1

816:   Output Parameter:
817: . Tf - The basis function values and derivatives at face quadrature points

819:   Level: intermediate

821:   Note:
822: .vb
823:   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
824:   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
825:   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
826: .ve

828: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
829: @*/
830: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
831: {
832:   PetscFunctionBegin;
834:   PetscAssertPointer(Tf, 3);
835:   if (!fem->Tf) {
836:     const PetscReal  xi0[3] = {-1., -1., -1.};
837:     PetscReal        v0[3], J[9], detJ;
838:     PetscQuadrature  fq;
839:     PetscDualSpace   sp;
840:     DM               dm;
841:     const PetscInt  *faces;
842:     PetscInt         dim, numFaces, f, npoints, q;
843:     const PetscReal *points;
844:     PetscReal       *facePoints;

846:     PetscCall(PetscFEGetDualSpace(fem, &sp));
847:     PetscCall(PetscDualSpaceGetDM(sp, &dm));
848:     PetscCall(DMGetDimension(dm, &dim));
849:     PetscCall(DMPlexGetConeSize(dm, 0, &numFaces));
850:     PetscCall(DMPlexGetCone(dm, 0, &faces));
851:     PetscCall(PetscFEGetFaceQuadrature(fem, &fq));
852:     if (fq) {
853:       PetscCall(PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL));
854:       PetscCall(PetscMalloc1(numFaces * npoints * dim, &facePoints));
855:       for (f = 0; f < numFaces; ++f) {
856:         PetscCall(DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ));
857:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * npoints + q) * dim]);
858:       }
859:       PetscCall(PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf));
860:       PetscCall(PetscFree(facePoints));
861:     }
862:   }
863:   PetscCheck(!fem->Tf || k <= fem->Tf->K, PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %" PetscInt_FMT " derivatives, but only tabulated %" PetscInt_FMT, k, fem->Tf->K);
864:   *Tf = fem->Tf;
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

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

871:   Not Collective

873:   Input Parameter:
874: . fem - The `PetscFE` object

876:   Output Parameter:
877: . Tc - The basis function values at face centroid points

879:   Level: intermediate

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

886: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
887: @*/
888: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
889: {
890:   PetscFunctionBegin;
892:   PetscAssertPointer(Tc, 2);
893:   if (!fem->Tc) {
894:     PetscDualSpace  sp;
895:     DM              dm;
896:     const PetscInt *cone;
897:     PetscReal      *centroids;
898:     PetscInt        dim, numFaces, f;

900:     PetscCall(PetscFEGetDualSpace(fem, &sp));
901:     PetscCall(PetscDualSpaceGetDM(sp, &dm));
902:     PetscCall(DMGetDimension(dm, &dim));
903:     PetscCall(DMPlexGetConeSize(dm, 0, &numFaces));
904:     PetscCall(DMPlexGetCone(dm, 0, &cone));
905:     PetscCall(PetscMalloc1(numFaces * dim, &centroids));
906:     for (f = 0; f < numFaces; ++f) PetscCall(DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f * dim], NULL));
907:     PetscCall(PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc));
908:     PetscCall(PetscFree(centroids));
909:   }
910:   *Tc = fem->Tc;
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

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

917:   Not Collective

919:   Input Parameters:
920: + fem     - The `PetscFE` object
921: . nrepl   - The number of replicas
922: . npoints - The number of tabulation points in a replica
923: . points  - The tabulation point coordinates
924: - K       - The number of derivatives calculated

926:   Output Parameter:
927: . T - The basis function values and derivatives at tabulation points

929:   Level: intermediate

931:   Note:
932: .vb
933:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
934:   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
935:   T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis
936:   T->function i, component c, in directions d and e
937: .ve

939: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
940: @*/
941: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
942: {
943:   DM             dm;
944:   PetscDualSpace Q;
945:   PetscInt       Nb;   /* Dimension of FE space P */
946:   PetscInt       Nc;   /* Field components */
947:   PetscInt       cdim; /* Reference coordinate dimension */
948:   PetscInt       k;

950:   PetscFunctionBegin;
951:   if (!npoints || !fem->dualSpace || K < 0) {
952:     *T = NULL;
953:     PetscFunctionReturn(PETSC_SUCCESS);
954:   }
956:   PetscAssertPointer(points, 4);
957:   PetscAssertPointer(T, 6);
958:   PetscCall(PetscFEGetDualSpace(fem, &Q));
959:   PetscCall(PetscDualSpaceGetDM(Q, &dm));
960:   PetscCall(DMGetDimension(dm, &cdim));
961:   PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
962:   PetscCall(PetscFEGetNumComponents(fem, &Nc));
963:   PetscCall(PetscMalloc1(1, T));
964:   (*T)->K    = !cdim ? 0 : K;
965:   (*T)->Nr   = nrepl;
966:   (*T)->Np   = npoints;
967:   (*T)->Nb   = Nb;
968:   (*T)->Nc   = Nc;
969:   (*T)->cdim = cdim;
970:   PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
971:   for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscCalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
972:   PetscUseTypeMethod(fem, createtabulation, nrepl * npoints, points, K, *T);
973:   PetscFunctionReturn(PETSC_SUCCESS);
974: }

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

979:   Not Collective

981:   Input Parameters:
982: + fem     - The `PetscFE` object
983: . npoints - The number of tabulation points
984: . points  - The tabulation point coordinates
985: . K       - The number of derivatives calculated
986: - T       - An existing tabulation object with enough allocated space

988:   Output Parameter:
989: . T - The basis function values and derivatives at tabulation points

991:   Level: intermediate

993:   Note:
994: .vb
995:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
996:   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
997:   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
998: .ve

1000: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
1001: @*/
1002: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1003: {
1004:   PetscFunctionBeginHot;
1005:   if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(PETSC_SUCCESS);
1007:   PetscAssertPointer(points, 3);
1008:   PetscAssertPointer(T, 5);
1009:   if (PetscDefined(USE_DEBUG)) {
1010:     DM             dm;
1011:     PetscDualSpace Q;
1012:     PetscInt       Nb;   /* Dimension of FE space P */
1013:     PetscInt       Nc;   /* Field components */
1014:     PetscInt       cdim; /* Reference coordinate dimension */

1016:     PetscCall(PetscFEGetDualSpace(fem, &Q));
1017:     PetscCall(PetscDualSpaceGetDM(Q, &dm));
1018:     PetscCall(DMGetDimension(dm, &cdim));
1019:     PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
1020:     PetscCall(PetscFEGetNumComponents(fem, &Nc));
1021:     PetscCheck(T->K == (!cdim ? 0 : K), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %" PetscInt_FMT " must match requested K %" PetscInt_FMT, T->K, !cdim ? 0 : K);
1022:     PetscCheck(T->Nb == Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %" PetscInt_FMT " must match requested Nb %" PetscInt_FMT, T->Nb, Nb);
1023:     PetscCheck(T->Nc == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %" PetscInt_FMT " must match requested Nc %" PetscInt_FMT, T->Nc, Nc);
1024:     PetscCheck(T->cdim == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %" PetscInt_FMT " must match requested cdim %" PetscInt_FMT, T->cdim, cdim);
1025:   }
1026:   T->Nr = 1;
1027:   T->Np = npoints;
1028:   PetscUseTypeMethod(fem, createtabulation, npoints, points, K, T);
1029:   PetscFunctionReturn(PETSC_SUCCESS);
1030: }

1032: /*@C
1033:   PetscTabulationDestroy - Frees memory from the associated tabulation.

1035:   Not Collective

1037:   Input Parameter:
1038: . T - The tabulation

1040:   Level: intermediate

1042: .seealso: `PetscTabulation`, `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
1043: @*/
1044: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1045: {
1046:   PetscInt k;

1048:   PetscFunctionBegin;
1049:   PetscAssertPointer(T, 1);
1050:   if (!T || !(*T)) PetscFunctionReturn(PETSC_SUCCESS);
1051:   for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscFree((*T)->T[k]));
1052:   PetscCall(PetscFree((*T)->T));
1053:   PetscCall(PetscFree(*T));
1054:   *T = NULL;
1055:   PetscFunctionReturn(PETSC_SUCCESS);
1056: }

1058: static PetscErrorCode PetscFECreatePointTraceDefault_Internal(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1059: {
1060:   PetscSpace      bsp, bsubsp;
1061:   PetscDualSpace  dsp, dsubsp;
1062:   PetscInt        dim, depth, numComp, i, j, coneSize, order;
1063:   DM              dm;
1064:   DMLabel         label;
1065:   PetscReal      *xi, *v, *J, detJ;
1066:   const char     *name;
1067:   PetscQuadrature origin, fullQuad, subQuad;

1069:   PetscFunctionBegin;
1070:   PetscCall(PetscFEGetBasisSpace(fe, &bsp));
1071:   PetscCall(PetscFEGetDualSpace(fe, &dsp));
1072:   PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1073:   PetscCall(DMGetDimension(dm, &dim));
1074:   PetscCall(DMPlexGetDepthLabel(dm, &label));
1075:   PetscCall(DMLabelGetValue(label, refPoint, &depth));
1076:   PetscCall(PetscCalloc1(depth, &xi));
1077:   PetscCall(PetscMalloc1(dim, &v));
1078:   PetscCall(PetscMalloc1(dim * dim, &J));
1079:   for (i = 0; i < depth; i++) xi[i] = 0.;
1080:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &origin));
1081:   PetscCall(PetscQuadratureSetData(origin, depth, 0, 1, xi, NULL));
1082:   PetscCall(DMPlexComputeCellGeometryFEM(dm, refPoint, origin, v, J, NULL, &detJ));
1083:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1084:   for (i = 1; i < dim; i++) {
1085:     for (j = 0; j < depth; j++) J[i * depth + j] = J[i * dim + j];
1086:   }
1087:   PetscCall(PetscQuadratureDestroy(&origin));
1088:   PetscCall(PetscDualSpaceGetPointSubspace(dsp, refPoint, &dsubsp));
1089:   PetscCall(PetscSpaceCreateSubspace(bsp, dsubsp, v, J, NULL, NULL, PETSC_OWN_POINTER, &bsubsp));
1090:   PetscCall(PetscSpaceSetUp(bsubsp));
1091:   PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), trFE));
1092:   PetscCall(PetscFESetType(*trFE, PETSCFEBASIC));
1093:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
1094:   PetscCall(PetscFESetNumComponents(*trFE, numComp));
1095:   PetscCall(PetscFESetBasisSpace(*trFE, bsubsp));
1096:   PetscCall(PetscFESetDualSpace(*trFE, dsubsp));
1097:   PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1098:   if (name) PetscCall(PetscFESetName(*trFE, name));
1099:   PetscCall(PetscFEGetQuadrature(fe, &fullQuad));
1100:   PetscCall(PetscQuadratureGetOrder(fullQuad, &order));
1101:   PetscCall(DMPlexGetConeSize(dm, refPoint, &coneSize));
1102:   if (coneSize == 2 * depth) PetscCall(PetscDTGaussTensorQuadrature(depth, 1, (order + 2) / 2, -1., 1., &subQuad));
1103:   else PetscCall(PetscDTSimplexQuadrature(depth, order, PETSCDTSIMPLEXQUAD_DEFAULT, &subQuad));
1104:   PetscCall(PetscFESetQuadrature(*trFE, subQuad));
1105:   PetscCall(PetscFESetUp(*trFE));
1106:   PetscCall(PetscQuadratureDestroy(&subQuad));
1107:   PetscCall(PetscSpaceDestroy(&bsubsp));
1108:   PetscFunctionReturn(PETSC_SUCCESS);
1109: }

1111: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1112: {
1113:   PetscFunctionBegin;
1115:   PetscAssertPointer(trFE, 3);
1116:   if (fe->ops->createpointtrace) PetscUseTypeMethod(fe, createpointtrace, refPoint, trFE);
1117:   else PetscCall(PetscFECreatePointTraceDefault_Internal(fe, refPoint, trFE));
1118:   PetscFunctionReturn(PETSC_SUCCESS);
1119: }

1121: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1122: {
1123:   PetscInt       hStart, hEnd;
1124:   PetscDualSpace dsp;
1125:   DM             dm;

1127:   PetscFunctionBegin;
1129:   PetscAssertPointer(trFE, 3);
1130:   *trFE = NULL;
1131:   PetscCall(PetscFEGetDualSpace(fe, &dsp));
1132:   PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1133:   PetscCall(DMPlexGetHeightStratum(dm, height, &hStart, &hEnd));
1134:   if (hEnd <= hStart) PetscFunctionReturn(PETSC_SUCCESS);
1135:   PetscCall(PetscFECreatePointTrace(fe, hStart, trFE));
1136:   PetscFunctionReturn(PETSC_SUCCESS);
1137: }

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

1142:   Not Collective

1144:   Input Parameter:
1145: . fem - The `PetscFE`

1147:   Output Parameter:
1148: . dim - The dimension

1150:   Level: intermediate

1152: .seealso: `PetscFE`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1153: @*/
1154: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1155: {
1156:   PetscFunctionBegin;
1158:   PetscAssertPointer(dim, 2);
1159:   PetscTryTypeMethod(fem, getdimension, dim);
1160:   PetscFunctionReturn(PETSC_SUCCESS);
1161: }

1163: /*@C
1164:   PetscFEPushforward - Map the reference element function to real space

1166:   Input Parameters:
1167: + fe     - The `PetscFE`
1168: . fegeom - The cell geometry
1169: . Nv     - The number of function values
1170: - vals   - The function values

1172:   Output Parameter:
1173: . vals - The transformed function values

1175:   Level: advanced

1177:   Notes:
1178:   This just forwards the call onto `PetscDualSpacePushforward()`.

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

1182: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscDualSpacePushforward()`
1183: @*/
1184: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1185: {
1186:   PetscFunctionBeginHot;
1187:   PetscCall(PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1188:   PetscFunctionReturn(PETSC_SUCCESS);
1189: }

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

1194:   Input Parameters:
1195: + fe     - The `PetscFE`
1196: . fegeom - The cell geometry
1197: . Nv     - The number of function gradient values
1198: - vals   - The function gradient values

1200:   Output Parameter:
1201: . vals - The transformed function gradient values

1203:   Level: advanced

1205:   Notes:
1206:   This just forwards the call onto `PetscDualSpacePushforwardGradient()`.

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

1210: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1211: @*/
1212: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1213: {
1214:   PetscFunctionBeginHot;
1215:   PetscCall(PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

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

1222:   Input Parameters:
1223: + fe     - The `PetscFE`
1224: . fegeom - The cell geometry
1225: . Nv     - The number of function Hessian values
1226: - vals   - The function Hessian values

1228:   Output Parameter:
1229: . vals - The transformed function Hessian values

1231:   Level: advanced

1233:   Notes:
1234:   This just forwards the call onto `PetscDualSpacePushforwardHessian()`.

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

1238:   Developer Notes:
1239:   It is unclear why all these one line convenience routines are desirable

1241: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1242: @*/
1243: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1244: {
1245:   PetscFunctionBeginHot;
1246:   PetscCall(PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1247:   PetscFunctionReturn(PETSC_SUCCESS);
1248: }

1250: /*
1251: Purpose: Compute element vector for chunk of elements

1253: Input:
1254:   Sizes:
1255:      Ne:  number of elements
1256:      Nf:  number of fields
1257:      PetscFE
1258:        dim: spatial dimension
1259:        Nb:  number of basis functions
1260:        Nc:  number of field components
1261:        PetscQuadrature
1262:          Nq:  number of quadrature points

1264:   Geometry:
1265:      PetscFEGeom[Ne] possibly *Nq
1266:        PetscReal v0s[dim]
1267:        PetscReal n[dim]
1268:        PetscReal jacobians[dim*dim]
1269:        PetscReal jacobianInverses[dim*dim]
1270:        PetscReal jacobianDeterminants
1271:   FEM:
1272:      PetscFE
1273:        PetscQuadrature
1274:          PetscReal   quadPoints[Nq*dim]
1275:          PetscReal   quadWeights[Nq]
1276:        PetscReal   basis[Nq*Nb*Nc]
1277:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1278:      PetscScalar coefficients[Ne*Nb*Nc]
1279:      PetscScalar elemVec[Ne*Nb*Nc]

1281:   Problem:
1282:      PetscInt f: the active field
1283:      f0, f1

1285:   Work Space:
1286:      PetscFE
1287:        PetscScalar f0[Nq*dim];
1288:        PetscScalar f1[Nq*dim*dim];
1289:        PetscScalar u[Nc];
1290:        PetscScalar gradU[Nc*dim];
1291:        PetscReal   x[dim];
1292:        PetscScalar realSpaceDer[dim];

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

1296: Input:
1297:   Sizes:
1298:      N_cb: Number of serial cell batches

1300:   Geometry:
1301:      PetscReal v0s[Ne*dim]
1302:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1303:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1304:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1305:   FEM:
1306:      static PetscReal   quadPoints[Nq*dim]
1307:      static PetscReal   quadWeights[Nq]
1308:      static PetscReal   basis[Nq*Nb*Nc]
1309:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1310:      PetscScalar coefficients[Ne*Nb*Nc]
1311:      PetscScalar elemVec[Ne*Nb*Nc]

1313: ex62.c:
1314:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1315:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1316:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1317:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1319: ex52.c:
1320:   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)
1321:   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)

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

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

1330: ex52_integrateElementOpenCL.c:
1331: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1332:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1333:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1341:   Not Collective

1343:   Input Parameters:
1344: + prob            - The `PetscDS` specifying the discretizations and continuum functions
1345: . field           - The field being integrated
1346: . Ne              - The number of elements in the chunk
1347: . cgeom           - The cell geometry for each cell in the chunk
1348: . coefficients    - The array of FEM basis coefficients for the elements
1349: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1350: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1352:   Output Parameter:
1353: . integral - the integral for this field

1355:   Level: intermediate

1357:   Developer Notes:
1358:   The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.

1360: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrateBd()`
1361: @*/
1362: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1363: {
1364:   PetscFE fe;

1366:   PetscFunctionBegin;
1368:   PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1369:   if (fe->ops->integrate) PetscCall((*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral));
1370:   PetscFunctionReturn(PETSC_SUCCESS);
1371: }

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

1376:   Not Collective

1378:   Input Parameters:
1379: + prob            - The `PetscDS` specifying the discretizations and continuum functions
1380: . field           - The field being integrated
1381: . obj_func        - The function to be integrated
1382: . Ne              - The number of elements in the chunk
1383: . geom            - The face geometry for each face in the chunk
1384: . coefficients    - The array of FEM basis coefficients for the elements
1385: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1386: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1388:   Output Parameter:
1389: . integral - the integral for this field

1391:   Level: intermediate

1393:   Developer Notes:
1394:   The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.

1396: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrate()`
1397: @*/
1398: 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[])
1399: {
1400:   PetscFE fe;

1402:   PetscFunctionBegin;
1404:   PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1405:   if (fe->ops->integratebd) PetscCall((*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral));
1406:   PetscFunctionReturn(PETSC_SUCCESS);
1407: }

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

1412:   Not Collective

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

1425:   Output Parameter:
1426: . elemVec - the element residual vectors from each element

1428:   Level: intermediate

1430:   Note:
1431: .vb
1432:   Loop over batch of elements (e):
1433:     Loop over quadrature points (q):
1434:       Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1435:       Call f_0 and f_1
1436:     Loop over element vector entries (f,fc --> i):
1437:       elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1438: .ve

1440: .seealso: `PetscFEIntegrateBdResidual()`
1441: @*/
1442: 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[])
1443: {
1444:   PetscFE fe;

1446:   PetscFunctionBeginHot;
1448:   PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1449:   if (fe->ops->integrateresidual) PetscCall((*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1450:   PetscFunctionReturn(PETSC_SUCCESS);
1451: }

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

1456:   Not Collective

1458:   Input Parameters:
1459: + ds              - The `PetscDS` specifying the discretizations and continuum functions
1460: . wf              - The PetscWeakForm object holding the pointwise functions
1461: . key             - The (label+value, field) being integrated
1462: . Ne              - The number of elements in the chunk
1463: . fgeom           - The face geometry for each cell in the chunk
1464: . coefficients    - The array of FEM basis coefficients for the elements
1465: . coefficients_t  - The array of FEM basis time derivative coefficients for the elements
1466: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1467: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1468: - t               - The time

1470:   Output Parameter:
1471: . elemVec - the element residual vectors from each element

1473:   Level: intermediate

1475: .seealso: `PetscFEIntegrateResidual()`
1476: @*/
1477: 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[])
1478: {
1479:   PetscFE fe;

1481:   PetscFunctionBegin;
1483:   PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1484:   if (fe->ops->integratebdresidual) PetscCall((*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1485:   PetscFunctionReturn(PETSC_SUCCESS);
1486: }

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

1491:   Not Collective

1493:   Input Parameters:
1494: + ds              - The `PetscDS` specifying the discretizations and continuum functions
1495: . dsIn            - The `PetscDS` specifying the discretizations and continuum functions for input
1496: . key             - The (label+value, field) being integrated
1497: . s               - The side of the cell being integrated, 0 for negative and 1 for positive
1498: . Ne              - The number of elements in the chunk
1499: . fgeom           - The face geometry for each cell in the chunk
1500: . coefficients    - The array of FEM basis coefficients for the elements
1501: . coefficients_t  - The array of FEM basis time derivative coefficients for the elements
1502: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1503: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1504: - t               - The time

1506:   Output Parameter:
1507: . elemVec - the element residual vectors from each element

1509:   Level: developer

1511: .seealso: `PetscFEIntegrateResidual()`
1512: @*/
1513: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1514: {
1515:   PetscFE fe;

1517:   PetscFunctionBegin;
1520:   PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1521:   if (fe->ops->integratehybridresidual) PetscCall((*fe->ops->integratehybridresidual)(ds, dsIn, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

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

1528:   Not Collective

1530:   Input Parameters:
1531: + ds              - The `PetscDS` specifying the discretizations and continuum functions
1532: . jtype           - The type of matrix pointwise functions that should be used
1533: . key             - The (label+value, fieldI*Nf + fieldJ) being integrated
1534: . Ne              - The number of elements in the chunk
1535: . cgeom           - The cell geometry for each cell in the chunk
1536: . coefficients    - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1537: . coefficients_t  - The array of FEM basis time derivative coefficients for the elements
1538: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1539: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1540: . t               - The time
1541: - u_tshift        - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1543:   Output Parameter:
1544: . elemMat - the element matrices for the Jacobian from each element

1546:   Level: intermediate

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

1560: .seealso: `PetscFEIntegrateResidual()`
1561: @*/
1562: 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[])
1563: {
1564:   PetscFE  fe;
1565:   PetscInt Nf;

1567:   PetscFunctionBegin;
1569:   PetscCall(PetscDSGetNumFields(ds, &Nf));
1570:   PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1571:   if (fe->ops->integratejacobian) PetscCall((*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1572:   PetscFunctionReturn(PETSC_SUCCESS);
1573: }

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

1578:   Not Collective

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

1594:   Output Parameter:
1595: . elemMat - the element matrices for the Jacobian from each element

1597:   Level: intermediate

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

1611: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1612: @*/
1613: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFEJacobianType jtype, 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[])
1614: {
1615:   PetscFE  fe;
1616:   PetscInt Nf;

1618:   PetscFunctionBegin;
1620:   PetscCall(PetscDSGetNumFields(ds, &Nf));
1621:   PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1622:   if (fe->ops->integratebdjacobian) PetscCall((*fe->ops->integratebdjacobian)(ds, wf, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1623:   PetscFunctionReturn(PETSC_SUCCESS);
1624: }

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

1629:   Not Collective

1631:   Input Parameters:
1632: + ds              - The `PetscDS` specifying the discretizations and continuum functions for the output
1633: . dsIn            - The `PetscDS` specifying the discretizations and continuum functions for the input
1634: . jtype           - The type of matrix pointwise functions that should be used
1635: . key             - The (label+value, fieldI*Nf + fieldJ) being integrated
1636: . s               - The side of the cell being integrated, 0 for negative and 1 for positive
1637: . Ne              - The number of elements in the chunk
1638: . fgeom           - The face geometry for each cell in the chunk
1639: . coefficients    - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1640: . coefficients_t  - The array of FEM basis time derivative coefficients for the elements
1641: . probAux         - The `PetscDS` specifying the auxiliary discretizations
1642: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1643: . t               - The time
1644: - u_tshift        - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1646:   Output Parameter:
1647: . elemMat - the element matrices for the Jacobian from each element

1649:   Level: developer

1651:   Note:
1652: .vb
1653:   Loop over batch of elements (e):
1654:     Loop over element matrix entries (f,fc,g,gc --> i,j):
1655:       Loop over quadrature points (q):
1656:         Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1657:           elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1658:                        + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1659:                        + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1660:                        + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1661: .ve

1663: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1664: @*/
1665: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscDS dsIn, 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[])
1666: {
1667:   PetscFE  fe;
1668:   PetscInt Nf;

1670:   PetscFunctionBegin;
1672:   PetscCall(PetscDSGetNumFields(ds, &Nf));
1673:   PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1674:   if (fe->ops->integratehybridjacobian) PetscCall((*fe->ops->integratehybridjacobian)(ds, dsIn, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1675:   PetscFunctionReturn(PETSC_SUCCESS);
1676: }

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

1681:   Input Parameters:
1682: + fe     - The finite element space
1683: - height - The height of the `DMPLEX` point

1685:   Output Parameter:
1686: . subfe - The subspace of this `PetscFE` space

1688:   Level: advanced

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

1693: .seealso: `PetscFECreateDefault()`
1694: @*/
1695: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1696: {
1697:   PetscSpace      P, subP;
1698:   PetscDualSpace  Q, subQ;
1699:   PetscQuadrature subq;
1700:   PetscInt        dim, Nc;

1702:   PetscFunctionBegin;
1704:   PetscAssertPointer(subfe, 3);
1705:   if (height == 0) {
1706:     *subfe = fe;
1707:     PetscFunctionReturn(PETSC_SUCCESS);
1708:   }
1709:   PetscCall(PetscFEGetBasisSpace(fe, &P));
1710:   PetscCall(PetscFEGetDualSpace(fe, &Q));
1711:   PetscCall(PetscFEGetNumComponents(fe, &Nc));
1712:   PetscCall(PetscFEGetFaceQuadrature(fe, &subq));
1713:   PetscCall(PetscDualSpaceGetDimension(Q, &dim));
1714:   PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
1715:   if (!fe->subspaces) PetscCall(PetscCalloc1(dim, &fe->subspaces));
1716:   if (height <= dim) {
1717:     if (!fe->subspaces[height - 1]) {
1718:       PetscFE     sub = NULL;
1719:       const char *name;

1721:       PetscCall(PetscSpaceGetHeightSubspace(P, height, &subP));
1722:       PetscCall(PetscDualSpaceGetHeightSubspace(Q, height, &subQ));
1723:       if (subQ) {
1724:         PetscCall(PetscObjectReference((PetscObject)subP));
1725:         PetscCall(PetscObjectReference((PetscObject)subQ));
1726:         PetscCall(PetscObjectReference((PetscObject)subq));
1727:         PetscCall(PetscFECreateFromSpaces(subP, subQ, subq, NULL, &sub));
1728:       }
1729:       if (sub) {
1730:         PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1731:         if (name) PetscCall(PetscFESetName(sub, name));
1732:       }
1733:       fe->subspaces[height - 1] = sub;
1734:     }
1735:     *subfe = fe->subspaces[height - 1];
1736:   } else {
1737:     *subfe = NULL;
1738:   }
1739:   PetscFunctionReturn(PETSC_SUCCESS);
1740: }

1742: /*@
1743:   PetscFERefine - Create a "refined" `PetscFE` object that refines the reference cell into
1744:   smaller copies.

1746:   Collective

1748:   Input Parameter:
1749: . fe - The initial `PetscFE`

1751:   Output Parameter:
1752: . feRef - The refined `PetscFE`

1754:   Level: advanced

1756:   Notes:
1757:   This is typically used to generate a preconditioner for a higher order method from a lower order method on a
1758:   refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1759:   interpolation between regularly refined meshes.

1761: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1762: @*/
1763: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1764: {
1765:   PetscSpace       P, Pref;
1766:   PetscDualSpace   Q, Qref;
1767:   DM               K, Kref;
1768:   PetscQuadrature  q, qref;
1769:   const PetscReal *v0, *jac;
1770:   PetscInt         numComp, numSubelements;
1771:   PetscInt         cStart, cEnd, c;
1772:   PetscDualSpace  *cellSpaces;

1774:   PetscFunctionBegin;
1775:   PetscCall(PetscFEGetBasisSpace(fe, &P));
1776:   PetscCall(PetscFEGetDualSpace(fe, &Q));
1777:   PetscCall(PetscFEGetQuadrature(fe, &q));
1778:   PetscCall(PetscDualSpaceGetDM(Q, &K));
1779:   /* Create space */
1780:   PetscCall(PetscObjectReference((PetscObject)P));
1781:   Pref = P;
1782:   /* Create dual space */
1783:   PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1784:   PetscCall(PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED));
1785:   PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref));
1786:   PetscCall(DMGetCoordinatesLocalSetUp(Kref));
1787:   PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1788:   PetscCall(DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd));
1789:   PetscCall(PetscMalloc1(cEnd - cStart, &cellSpaces));
1790:   /* TODO: fix for non-uniform refinement */
1791:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1792:   PetscCall(PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces));
1793:   PetscCall(PetscFree(cellSpaces));
1794:   PetscCall(DMDestroy(&Kref));
1795:   PetscCall(PetscDualSpaceSetUp(Qref));
1796:   /* Create element */
1797:   PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), feRef));
1798:   PetscCall(PetscFESetType(*feRef, PETSCFECOMPOSITE));
1799:   PetscCall(PetscFESetBasisSpace(*feRef, Pref));
1800:   PetscCall(PetscFESetDualSpace(*feRef, Qref));
1801:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
1802:   PetscCall(PetscFESetNumComponents(*feRef, numComp));
1803:   PetscCall(PetscFESetUp(*feRef));
1804:   PetscCall(PetscSpaceDestroy(&Pref));
1805:   PetscCall(PetscDualSpaceDestroy(&Qref));
1806:   /* Create quadrature */
1807:   PetscCall(PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL));
1808:   PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1809:   PetscCall(PetscFESetQuadrature(*feRef, qref));
1810:   PetscCall(PetscQuadratureDestroy(&qref));
1811:   PetscFunctionReturn(PETSC_SUCCESS);
1812: }

1814: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1815: {
1816:   PetscSpace     P;
1817:   PetscDualSpace Q;
1818:   DM             K;
1819:   DMPolytopeType ct;
1820:   PetscInt       degree;
1821:   char           name[64];

1823:   PetscFunctionBegin;
1824:   PetscCall(PetscFEGetBasisSpace(fe, &P));
1825:   PetscCall(PetscSpaceGetDegree(P, &degree, NULL));
1826:   PetscCall(PetscFEGetDualSpace(fe, &Q));
1827:   PetscCall(PetscDualSpaceGetDM(Q, &K));
1828:   PetscCall(DMPlexGetCellType(K, 0, &ct));
1829:   switch (ct) {
1830:   case DM_POLYTOPE_SEGMENT:
1831:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1832:   case DM_POLYTOPE_QUADRILATERAL:
1833:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1834:   case DM_POLYTOPE_HEXAHEDRON:
1835:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1836:     PetscCall(PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree));
1837:     break;
1838:   case DM_POLYTOPE_TRIANGLE:
1839:   case DM_POLYTOPE_TETRAHEDRON:
1840:     PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree));
1841:     break;
1842:   case DM_POLYTOPE_TRI_PRISM:
1843:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1844:     PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree));
1845:     break;
1846:   default:
1847:     PetscCall(PetscSNPrintf(name, sizeof(name), "FE"));
1848:   }
1849:   PetscCall(PetscFESetName(fe, name));
1850:   PetscFunctionReturn(PETSC_SUCCESS);
1851: }

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

1856:   Collective

1858:   Input Parameters:
1859: + P  - The basis space
1860: . Q  - The dual space
1861: . q  - The cell quadrature
1862: - fq - The face quadrature

1864:   Output Parameter:
1865: . fem - The `PetscFE` object

1867:   Level: beginner

1869:   Note:
1870:   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.

1872: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`,
1873:           `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1874: @*/
1875: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1876: {
1877:   PetscInt    Nc;
1878:   PetscInt    p_Ns = -1, p_Nc = -1, q_Ns = -1, q_Nc = -1;
1879:   PetscBool   p_is_uniform_sum = PETSC_FALSE, p_interleave_basis = PETSC_FALSE, p_interleave_components = PETSC_FALSE;
1880:   PetscBool   q_is_uniform_sum = PETSC_FALSE, q_interleave_basis = PETSC_FALSE, q_interleave_components = PETSC_FALSE;
1881:   const char *prefix;

1883:   PetscFunctionBegin;
1884:   PetscCall(PetscObjectTypeCompare((PetscObject)P, PETSCSPACESUM, &p_is_uniform_sum));
1885:   if (p_is_uniform_sum) {
1886:     PetscSpace subsp_0 = NULL;
1887:     PetscCall(PetscSpaceSumGetNumSubspaces(P, &p_Ns));
1888:     PetscCall(PetscSpaceGetNumComponents(P, &p_Nc));
1889:     PetscCall(PetscSpaceSumGetConcatenate(P, &p_is_uniform_sum));
1890:     PetscCall(PetscSpaceSumGetInterleave(P, &p_interleave_basis, &p_interleave_components));
1891:     for (PetscInt s = 0; s < p_Ns; s++) {
1892:       PetscSpace subsp;

1894:       PetscCall(PetscSpaceSumGetSubspace(P, s, &subsp));
1895:       if (!s) {
1896:         subsp_0 = subsp;
1897:       } else if (subsp != subsp_0) {
1898:         p_is_uniform_sum = PETSC_FALSE;
1899:       }
1900:     }
1901:   }
1902:   PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &q_is_uniform_sum));
1903:   if (q_is_uniform_sum) {
1904:     PetscDualSpace subsp_0 = NULL;
1905:     PetscCall(PetscDualSpaceSumGetNumSubspaces(Q, &q_Ns));
1906:     PetscCall(PetscDualSpaceGetNumComponents(Q, &q_Nc));
1907:     PetscCall(PetscDualSpaceSumGetConcatenate(Q, &q_is_uniform_sum));
1908:     PetscCall(PetscDualSpaceSumGetInterleave(Q, &q_interleave_basis, &q_interleave_components));
1909:     for (PetscInt s = 0; s < q_Ns; s++) {
1910:       PetscDualSpace subsp;

1912:       PetscCall(PetscDualSpaceSumGetSubspace(Q, s, &subsp));
1913:       if (!s) {
1914:         subsp_0 = subsp;
1915:       } else if (subsp != subsp_0) {
1916:         q_is_uniform_sum = PETSC_FALSE;
1917:       }
1918:     }
1919:   }
1920:   if (p_is_uniform_sum && q_is_uniform_sum && (p_interleave_basis == q_interleave_basis) && (p_interleave_components == q_interleave_components) && (p_Ns == q_Ns) && (p_Nc == q_Nc)) {
1921:     PetscSpace     scalar_space;
1922:     PetscDualSpace scalar_dspace;
1923:     PetscFE        scalar_fe;

1925:     PetscCall(PetscSpaceSumGetSubspace(P, 0, &scalar_space));
1926:     PetscCall(PetscDualSpaceSumGetSubspace(Q, 0, &scalar_dspace));
1927:     PetscCall(PetscObjectReference((PetscObject)scalar_space));
1928:     PetscCall(PetscObjectReference((PetscObject)scalar_dspace));
1929:     PetscCall(PetscObjectReference((PetscObject)q));
1930:     PetscCall(PetscObjectReference((PetscObject)fq));
1931:     PetscCall(PetscFECreateFromSpaces(scalar_space, scalar_dspace, q, fq, &scalar_fe));
1932:     PetscCall(PetscFECreateVector(scalar_fe, p_Ns, p_interleave_basis, p_interleave_components, fem));
1933:     PetscCall(PetscFEDestroy(&scalar_fe));
1934:   } else {
1935:     PetscCall(PetscFECreate(PetscObjectComm((PetscObject)P), fem));
1936:     PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
1937:   }
1938:   PetscCall(PetscSpaceGetNumComponents(P, &Nc));
1939:   PetscCall(PetscFESetNumComponents(*fem, Nc));
1940:   PetscCall(PetscFESetBasisSpace(*fem, P));
1941:   PetscCall(PetscFESetDualSpace(*fem, Q));
1942:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)P, &prefix));
1943:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
1944:   PetscCall(PetscFESetUp(*fem));
1945:   PetscCall(PetscSpaceDestroy(&P));
1946:   PetscCall(PetscDualSpaceDestroy(&Q));
1947:   PetscCall(PetscFESetQuadrature(*fem, q));
1948:   PetscCall(PetscFESetFaceQuadrature(*fem, fq));
1949:   PetscCall(PetscQuadratureDestroy(&q));
1950:   PetscCall(PetscQuadratureDestroy(&fq));
1951:   PetscCall(PetscFESetDefaultName_Private(*fem));
1952:   PetscFunctionReturn(PETSC_SUCCESS);
1953: }

1955: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1956: {
1957:   DM              K;
1958:   PetscSpace      P;
1959:   PetscDualSpace  Q;
1960:   PetscQuadrature q, fq;
1961:   PetscBool       tensor;

1963:   PetscFunctionBegin;
1964:   if (prefix) PetscAssertPointer(prefix, 5);
1965:   PetscAssertPointer(fem, 9);
1966:   switch (ct) {
1967:   case DM_POLYTOPE_SEGMENT:
1968:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1969:   case DM_POLYTOPE_QUADRILATERAL:
1970:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1971:   case DM_POLYTOPE_HEXAHEDRON:
1972:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1973:     tensor = PETSC_TRUE;
1974:     break;
1975:   default:
1976:     tensor = PETSC_FALSE;
1977:   }
1978:   /* Create space */
1979:   PetscCall(PetscSpaceCreate(comm, &P));
1980:   PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
1981:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
1982:   PetscCall(PetscSpacePolynomialSetTensor(P, tensor));
1983:   PetscCall(PetscSpaceSetNumComponents(P, Nc));
1984:   PetscCall(PetscSpaceSetNumVariables(P, dim));
1985:   if (degree >= 0) {
1986:     PetscCall(PetscSpaceSetDegree(P, degree, PETSC_DETERMINE));
1987:     if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1988:       PetscSpace Pend, Pside;

1990:       PetscCall(PetscSpaceSetNumComponents(P, 1));
1991:       PetscCall(PetscSpaceCreate(comm, &Pend));
1992:       PetscCall(PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL));
1993:       PetscCall(PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE));
1994:       PetscCall(PetscSpaceSetNumComponents(Pend, 1));
1995:       PetscCall(PetscSpaceSetNumVariables(Pend, dim - 1));
1996:       PetscCall(PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE));
1997:       PetscCall(PetscSpaceCreate(comm, &Pside));
1998:       PetscCall(PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL));
1999:       PetscCall(PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE));
2000:       PetscCall(PetscSpaceSetNumComponents(Pside, 1));
2001:       PetscCall(PetscSpaceSetNumVariables(Pside, 1));
2002:       PetscCall(PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE));
2003:       PetscCall(PetscSpaceSetType(P, PETSCSPACETENSOR));
2004:       PetscCall(PetscSpaceTensorSetNumSubspaces(P, 2));
2005:       PetscCall(PetscSpaceTensorSetSubspace(P, 0, Pend));
2006:       PetscCall(PetscSpaceTensorSetSubspace(P, 1, Pside));
2007:       PetscCall(PetscSpaceDestroy(&Pend));
2008:       PetscCall(PetscSpaceDestroy(&Pside));

2010:       if (Nc > 1) {
2011:         PetscSpace scalar_P = P;

2013:         PetscCall(PetscSpaceCreate(comm, &P));
2014:         PetscCall(PetscSpaceSetNumVariables(P, dim));
2015:         PetscCall(PetscSpaceSetNumComponents(P, Nc));
2016:         PetscCall(PetscSpaceSetType(P, PETSCSPACESUM));
2017:         PetscCall(PetscSpaceSumSetNumSubspaces(P, Nc));
2018:         PetscCall(PetscSpaceSumSetConcatenate(P, PETSC_TRUE));
2019:         PetscCall(PetscSpaceSumSetInterleave(P, PETSC_TRUE, PETSC_FALSE));
2020:         for (PetscInt i = 0; i < Nc; i++) PetscCall(PetscSpaceSumSetSubspace(P, i, scalar_P));
2021:         PetscCall(PetscSpaceDestroy(&scalar_P));
2022:       }
2023:     }
2024:   }
2025:   if (setFromOptions) PetscCall(PetscSpaceSetFromOptions(P));
2026:   PetscCall(PetscSpaceSetUp(P));
2027:   PetscCall(PetscSpaceGetDegree(P, &degree, NULL));
2028:   PetscCall(PetscSpacePolynomialGetTensor(P, &tensor));
2029:   PetscCall(PetscSpaceGetNumComponents(P, &Nc));
2030:   /* Create dual space */
2031:   PetscCall(PetscDualSpaceCreate(comm, &Q));
2032:   PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
2033:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
2034:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
2035:   PetscCall(PetscDualSpaceSetDM(Q, K));
2036:   PetscCall(DMDestroy(&K));
2037:   PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
2038:   PetscCall(PetscDualSpaceSetOrder(Q, degree));
2039:   PetscCall(PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE));
2040:   if (setFromOptions) PetscCall(PetscDualSpaceSetFromOptions(Q));
2041:   PetscCall(PetscDualSpaceSetUp(Q));
2042:   /* Create quadrature */
2043:   qorder = qorder >= 0 ? qorder : degree;
2044:   if (setFromOptions) {
2045:     PetscObjectOptionsBegin((PetscObject)P);
2046:     PetscCall(PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0));
2047:     PetscOptionsEnd();
2048:   }
2049:   PetscCall(PetscDTCreateDefaultQuadrature(ct, qorder, &q, &fq));
2050:   /* Create finite element */
2051:   PetscCall(PetscFECreateFromSpaces(P, Q, q, fq, fem));
2052:   if (setFromOptions) PetscCall(PetscFESetFromOptions(*fem));
2053:   PetscFunctionReturn(PETSC_SUCCESS);
2054: }

2056: /*@C
2057:   PetscFECreateDefault - Create a `PetscFE` for basic FEM computation

2059:   Collective

2061:   Input Parameters:
2062: + comm      - The MPI comm
2063: . dim       - The spatial dimension
2064: . Nc        - The number of components
2065: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2066: . prefix    - The options prefix, or `NULL`
2067: - qorder    - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree

2069:   Output Parameter:
2070: . fem - The `PetscFE` object

2072:   Level: beginner

2074:   Note:
2075:   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.

2077: .seealso: `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2078: @*/
2079: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
2080: {
2081:   PetscFunctionBegin;
2082:   PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2083:   PetscFunctionReturn(PETSC_SUCCESS);
2084: }

2086: /*@C
2087:   PetscFECreateByCell - Create a `PetscFE` for basic FEM computation

2089:   Collective

2091:   Input Parameters:
2092: + comm   - The MPI comm
2093: . dim    - The spatial dimension
2094: . Nc     - The number of components
2095: . ct     - The celltype of the reference cell
2096: . prefix - The options prefix, or `NULL`
2097: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree

2099:   Output Parameter:
2100: . fem - The `PetscFE` object

2102:   Level: beginner

2104:   Note:
2105:   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.

2107: .seealso: `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2108: @*/
2109: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
2110: {
2111:   PetscFunctionBegin;
2112:   PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2113:   PetscFunctionReturn(PETSC_SUCCESS);
2114: }

2116: /*@
2117:   PetscFECreateLagrange - Create a `PetscFE` for the basic Lagrange space of degree k

2119:   Collective

2121:   Input Parameters:
2122: + comm      - The MPI comm
2123: . dim       - The spatial dimension
2124: . Nc        - The number of components
2125: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2126: . k         - The degree k of the space
2127: - qorder    - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree

2129:   Output Parameter:
2130: . fem - The `PetscFE` object

2132:   Level: beginner

2134:   Note:
2135:   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.

2137: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2138: @*/
2139: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2140: {
2141:   PetscFunctionBegin;
2142:   PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem));
2143:   PetscFunctionReturn(PETSC_SUCCESS);
2144: }

2146: /*@
2147:   PetscFECreateLagrangeByCell - Create a `PetscFE` for the basic Lagrange space of degree k

2149:   Collective

2151:   Input Parameters:
2152: + comm   - The MPI comm
2153: . dim    - The spatial dimension
2154: . Nc     - The number of components
2155: . ct     - The celltype of the reference cell
2156: . k      - The degree k of the space
2157: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree

2159:   Output Parameter:
2160: . fem - The `PetscFE` object

2162:   Level: beginner

2164:   Note:
2165:   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.

2167: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2168: @*/
2169: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2170: {
2171:   PetscFunctionBegin;
2172:   PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem));
2173:   PetscFunctionReturn(PETSC_SUCCESS);
2174: }

2176: /*@C
2177:   PetscFESetName - Names the `PetscFE` and its subobjects

2179:   Not Collective

2181:   Input Parameters:
2182: + fe   - The `PetscFE`
2183: - name - The name

2185:   Level: intermediate

2187: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2188: @*/
2189: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2190: {
2191:   PetscSpace     P;
2192:   PetscDualSpace Q;

2194:   PetscFunctionBegin;
2195:   PetscCall(PetscFEGetBasisSpace(fe, &P));
2196:   PetscCall(PetscFEGetDualSpace(fe, &Q));
2197:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
2198:   PetscCall(PetscObjectSetName((PetscObject)P, name));
2199:   PetscCall(PetscObjectSetName((PetscObject)Q, name));
2200:   PetscFunctionReturn(PETSC_SUCCESS);
2201: }

2203: 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[])
2204: {
2205:   PetscInt dOffset = 0, fOffset = 0, f, g;

2207:   for (f = 0; f < Nf; ++f) {
2208:     PetscCheck(r < T[f]->Nr, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", r, T[f]->Nr);
2209:     PetscCheck(q < T[f]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", q, T[f]->Np);
2210:     PetscFE          fe;
2211:     const PetscInt   k       = ds->jetDegree[f];
2212:     const PetscInt   cdim    = T[f]->cdim;
2213:     const PetscInt   dE      = fegeom->dimEmbed;
2214:     const PetscInt   Nq      = T[f]->Np;
2215:     const PetscInt   Nbf     = T[f]->Nb;
2216:     const PetscInt   Ncf     = T[f]->Nc;
2217:     const PetscReal *Bq      = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2218:     const PetscReal *Dq      = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2219:     const PetscReal *Hq      = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2220:     PetscInt         hOffset = 0, b, c, d;

2222:     PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
2223:     for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2224:     for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2225:     for (b = 0; b < Nbf; ++b) {
2226:       for (c = 0; c < Ncf; ++c) {
2227:         const PetscInt cidx = b * Ncf + c;

2229:         u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2230:         for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2231:       }
2232:     }
2233:     if (k > 1) {
2234:       for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * dE;
2235:       for (d = 0; d < dE * dE * Ncf; ++d) u_x[hOffset + fOffset * dE * dE + d] = 0.0;
2236:       for (b = 0; b < Nbf; ++b) {
2237:         for (c = 0; c < Ncf; ++c) {
2238:           const PetscInt cidx = b * Ncf + c;

2240:           for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * dE * dE + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2241:         }
2242:       }
2243:       PetscCall(PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * dE * dE]));
2244:     }
2245:     PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset]));
2246:     PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]));
2247:     if (u_t) {
2248:       for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2249:       for (b = 0; b < Nbf; ++b) {
2250:         for (c = 0; c < Ncf; ++c) {
2251:           const PetscInt cidx = b * Ncf + c;

2253:           u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2254:         }
2255:       }
2256:       PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2257:     }
2258:     fOffset += Ncf;
2259:     dOffset += Nbf;
2260:   }
2261:   return PETSC_SUCCESS;
2262: }

2264: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt rc, PetscInt qc, PetscTabulation Tab[], const PetscInt rf[], const PetscInt qf[], PetscTabulation Tabf[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2265: {
2266:   PetscInt dOffset = 0, fOffset = 0, f, g;

2268:   /* f is the field number in the DS, g is the field number in u[] */
2269:   for (f = 0, g = 0; f < Nf; ++f) {
2270:     PetscBool isCohesive;
2271:     PetscInt  Ns, s;

2273:     if (!Tab[f]) continue;
2274:     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
2275:     Ns = isCohesive ? 1 : 2;
2276:     {
2277:       PetscTabulation T   = isCohesive ? Tab[f] : Tabf[f];
2278:       PetscFE         fe  = (PetscFE)ds->disc[f];
2279:       const PetscInt  dEt = T->cdim;
2280:       const PetscInt  dE  = fegeom->dimEmbed;
2281:       const PetscInt  Nq  = T->Np;
2282:       const PetscInt  Nbf = T->Nb;
2283:       const PetscInt  Ncf = T->Nc;

2285:       for (s = 0; s < Ns; ++s, ++g) {
2286:         const PetscInt   r  = isCohesive ? rc : rf[s];
2287:         const PetscInt   q  = isCohesive ? qc : qf[s];
2288:         const PetscReal *Bq = &T->T[0][(r * Nq + q) * Nbf * Ncf];
2289:         const PetscReal *Dq = &T->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2290:         PetscInt         b, c, d;

2292:         PetscCheck(r < T->Nr, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " Side %" PetscInt_FMT " Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", f, s, r, T->Nr);
2293:         PetscCheck(q < T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " Side %" PetscInt_FMT " Point number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", f, s, q, T->Np);
2294:         for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2295:         for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2296:         for (b = 0; b < Nbf; ++b) {
2297:           for (c = 0; c < Ncf; ++c) {
2298:             const PetscInt cidx = b * Ncf + c;

2300:             u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2301:             for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2302:           }
2303:         }
2304:         PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset]));
2305:         PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]));
2306:         if (u_t) {
2307:           for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2308:           for (b = 0; b < Nbf; ++b) {
2309:             for (c = 0; c < Ncf; ++c) {
2310:               const PetscInt cidx = b * Ncf + c;

2312:               u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2313:             }
2314:           }
2315:           PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2316:         }
2317:         fOffset += Ncf;
2318:         dOffset += Nbf;
2319:       }
2320:     }
2321:   }
2322:   return PETSC_SUCCESS;
2323: }

2325: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2326: {
2327:   PetscFE         fe;
2328:   PetscTabulation Tc;
2329:   PetscInt        b, c;

2331:   if (!prob) return PETSC_SUCCESS;
2332:   PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
2333:   PetscCall(PetscFEGetFaceCentroidTabulation(fe, &Tc));
2334:   {
2335:     const PetscReal *faceBasis = Tc->T[0];
2336:     const PetscInt   Nb        = Tc->Nb;
2337:     const PetscInt   Nc        = Tc->Nc;

2339:     for (c = 0; c < Nc; ++c) u[c] = 0.0;
2340:     for (b = 0; b < Nb; ++b) {
2341:       for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2342:     }
2343:   }
2344:   return PETSC_SUCCESS;
2345: }

2347: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2348: {
2349:   PetscFEGeom      pgeom;
2350:   const PetscInt   dEt      = T->cdim;
2351:   const PetscInt   dE       = fegeom->dimEmbed;
2352:   const PetscInt   Nq       = T->Np;
2353:   const PetscInt   Nb       = T->Nb;
2354:   const PetscInt   Nc       = T->Nc;
2355:   const PetscReal *basis    = &T->T[0][r * Nq * Nb * Nc];
2356:   const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt];
2357:   PetscInt         q, b, c, d;

2359:   for (q = 0; q < Nq; ++q) {
2360:     for (b = 0; b < Nb; ++b) {
2361:       for (c = 0; c < Nc; ++c) {
2362:         const PetscInt bcidx = b * Nc + c;

2364:         tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2365:         for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d];
2366:         for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0;
2367:       }
2368:     }
2369:     PetscCall(PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom));
2370:     PetscCall(PetscFEPushforward(fe, &pgeom, Nb, tmpBasis));
2371:     PetscCall(PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer));
2372:     for (b = 0; b < Nb; ++b) {
2373:       for (c = 0; c < Nc; ++c) {
2374:         const PetscInt bcidx = b * Nc + c;
2375:         const PetscInt qcidx = q * Nc + c;

2377:         elemVec[b] += tmpBasis[bcidx] * f0[qcidx];
2378:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2379:       }
2380:     }
2381:   }
2382:   return PETSC_SUCCESS;
2383: }

2385: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt side, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2386: {
2387:   const PetscInt   dE       = T->cdim;
2388:   const PetscInt   Nq       = T->Np;
2389:   const PetscInt   Nb       = T->Nb;
2390:   const PetscInt   Nc       = T->Nc;
2391:   const PetscReal *basis    = &T->T[0][r * Nq * Nb * Nc];
2392:   const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];

2394:   for (PetscInt q = 0; q < Nq; ++q) {
2395:     for (PetscInt b = 0; b < Nb; ++b) {
2396:       for (PetscInt c = 0; c < Nc; ++c) {
2397:         const PetscInt bcidx = b * Nc + c;

2399:         tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2400:         for (PetscInt d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2401:       }
2402:     }
2403:     PetscCall(PetscFEPushforward(fe, fegeom, Nb, tmpBasis));
2404:     // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2405:     // PetscCall(PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer));
2406:     if (side == 2) {
2407:       // Integrating over whole cohesive cell, so insert for both sides
2408:       for (PetscInt s = 0; s < 2; ++s) {
2409:         for (PetscInt b = 0; b < Nb; ++b) {
2410:           for (PetscInt c = 0; c < Nc; ++c) {
2411:             const PetscInt bcidx = b * Nc + c;
2412:             const PetscInt qcidx = (q * 2 + s) * Nc + c;

2414:             elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2415:             for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2416:           }
2417:         }
2418:       }
2419:     } else {
2420:       // Integrating over endcaps of cohesive cell, so insert for correct side
2421:       for (PetscInt b = 0; b < Nb; ++b) {
2422:         for (PetscInt c = 0; c < Nc; ++c) {
2423:           const PetscInt bcidx = b * Nc + c;
2424:           const PetscInt qcidx = q * Nc + c;

2426:           elemVec[Nb * side + b] += tmpBasis[bcidx] * f0[qcidx];
2427:           for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * side + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2428:         }
2429:       }
2430:     }
2431:   }
2432:   return PETSC_SUCCESS;
2433: }

2435: 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[])
2436: {
2437:   const PetscInt   cdim      = TI->cdim;
2438:   const PetscInt   dE        = fegeom->dimEmbed;
2439:   const PetscInt   NqI       = TI->Np;
2440:   const PetscInt   NbI       = TI->Nb;
2441:   const PetscInt   NcI       = TI->Nc;
2442:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2443:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * cdim];
2444:   const PetscInt   NqJ       = TJ->Np;
2445:   const PetscInt   NbJ       = TJ->Nb;
2446:   const PetscInt   NcJ       = TJ->Nc;
2447:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2448:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * cdim];
2449:   PetscInt         f, fc, g, gc, df, dg;

2451:   for (f = 0; f < NbI; ++f) {
2452:     for (fc = 0; fc < NcI; ++fc) {
2453:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2455:       tmpBasisI[fidx] = basisI[fidx];
2456:       for (df = 0; df < cdim; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * cdim + df];
2457:     }
2458:   }
2459:   PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2460:   PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2461:   for (g = 0; g < NbJ; ++g) {
2462:     for (gc = 0; gc < NcJ; ++gc) {
2463:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2465:       tmpBasisJ[gidx] = basisJ[gidx];
2466:       for (dg = 0; dg < cdim; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * cdim + dg];
2467:     }
2468:   }
2469:   PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2470:   PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2471:   for (f = 0; f < NbI; ++f) {
2472:     for (fc = 0; fc < NcI; ++fc) {
2473:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2474:       const PetscInt i    = offsetI + f;  /* Element matrix row */
2475:       for (g = 0; g < NbJ; ++g) {
2476:         for (gc = 0; gc < NcJ; ++gc) {
2477:           const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2478:           const PetscInt j    = offsetJ + g;  /* Element matrix column */
2479:           const PetscInt fOff = eOffset + i * totDim + j;

2481:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2482:           for (df = 0; df < dE; ++df) {
2483:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2484:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2485:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2486:           }
2487:         }
2488:       }
2489:     }
2490:   }
2491:   return PETSC_SUCCESS;
2492: }

2494: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt s, PetscInt t, 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[])
2495: {
2496:   const PetscInt   dE        = TI->cdim;
2497:   const PetscInt   NqI       = TI->Np;
2498:   const PetscInt   NbI       = TI->Nb;
2499:   const PetscInt   NcI       = TI->Nc;
2500:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2501:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2502:   const PetscInt   NqJ       = TJ->Np;
2503:   const PetscInt   NbJ       = TJ->Nb;
2504:   const PetscInt   NcJ       = TJ->Nc;
2505:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2506:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2507:   const PetscInt   so        = isHybridI ? 0 : s;
2508:   const PetscInt   to        = isHybridJ ? 0 : t;
2509:   PetscInt         f, fc, g, gc, df, dg;

2511:   for (f = 0; f < NbI; ++f) {
2512:     for (fc = 0; fc < NcI; ++fc) {
2513:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2515:       tmpBasisI[fidx] = basisI[fidx];
2516:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2517:     }
2518:   }
2519:   PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2520:   PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2521:   for (g = 0; g < NbJ; ++g) {
2522:     for (gc = 0; gc < NcJ; ++gc) {
2523:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2525:       tmpBasisJ[gidx] = basisJ[gidx];
2526:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2527:     }
2528:   }
2529:   PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2530:   // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2531:   // PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2532:   for (f = 0; f < NbI; ++f) {
2533:     for (fc = 0; fc < NcI; ++fc) {
2534:       const PetscInt fidx = f * NcI + fc;           /* Test function basis index */
2535:       const PetscInt i    = offsetI + NbI * so + f; /* Element matrix row */
2536:       for (g = 0; g < NbJ; ++g) {
2537:         for (gc = 0; gc < NcJ; ++gc) {
2538:           const PetscInt gidx = g * NcJ + gc;           /* Trial function basis index */
2539:           const PetscInt j    = offsetJ + NbJ * to + g; /* Element matrix column */
2540:           const PetscInt fOff = eOffset + i * totDim + j;

2542:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2543:           for (df = 0; df < dE; ++df) {
2544:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2545:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2546:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2547:           }
2548:         }
2549:       }
2550:     }
2551:   }
2552:   return PETSC_SUCCESS;
2553: }

2555: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2556: {
2557:   PetscDualSpace  dsp;
2558:   DM              dm;
2559:   PetscQuadrature quadDef;
2560:   PetscInt        dim, cdim, Nq;

2562:   PetscFunctionBegin;
2563:   PetscCall(PetscFEGetDualSpace(fe, &dsp));
2564:   PetscCall(PetscDualSpaceGetDM(dsp, &dm));
2565:   PetscCall(DMGetDimension(dm, &dim));
2566:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2567:   PetscCall(PetscFEGetQuadrature(fe, &quadDef));
2568:   quad = quad ? quad : quadDef;
2569:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
2570:   PetscCall(PetscMalloc1(Nq * cdim, &cgeom->v));
2571:   PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->J));
2572:   PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ));
2573:   PetscCall(PetscMalloc1(Nq, &cgeom->detJ));
2574:   cgeom->dim       = dim;
2575:   cgeom->dimEmbed  = cdim;
2576:   cgeom->numCells  = 1;
2577:   cgeom->numPoints = Nq;
2578:   PetscCall(DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ));
2579:   PetscFunctionReturn(PETSC_SUCCESS);
2580: }

2582: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2583: {
2584:   PetscFunctionBegin;
2585:   PetscCall(PetscFree(cgeom->v));
2586:   PetscCall(PetscFree(cgeom->J));
2587:   PetscCall(PetscFree(cgeom->invJ));
2588:   PetscCall(PetscFree(cgeom->detJ));
2589:   PetscFunctionReturn(PETSC_SUCCESS);
2590: }