Actual source code: fe.c

  1: /* Basis Jet Tabulation

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

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

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

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

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

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

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

 22:    \alpha = V^{-1}

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

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

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

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

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

 51: PetscClassId PETSCFE_CLASSID = 0;

 53: PetscLogEvent PETSCFE_SetUp;

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

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

 61:   Not Collective

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

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

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

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

 85:   Level: advanced

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

 89: @*/
 90: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
 91: {
 92:   PetscFunctionListAdd(&PetscFEList, sname, function);
 93:   return 0;
 94: }

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

 99:   Collective on fem

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: `PetscFEGetType()`, `PetscFECreate()`
111: @*/
112: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
113: {
114:   PetscErrorCode (*r)(PetscFE);
115:   PetscBool match;

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

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

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

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

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

136:   Not Collective

138:   Input Parameter:
139: . fem  - The PetscFE

141:   Output Parameter:
142: . name - The PetscFE type name

144:   Level: intermediate

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

157: /*@C
158:    PetscFEViewFromOptions - View from Options

160:    Collective on PetscFE

162:    Input Parameters:
163: +  A - the PetscFE object
164: .  obj - Optional object
165: -  name - command line option

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

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

180:   Collective on fem

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

186:   Level: beginner

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

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

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

206:   Collective on fem

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

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

215:   Level: intermediate

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

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

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

250: /*@C
251:   PetscFESetUp - Construct data structures for the PetscFE

253:   Collective on fem

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

258:   Level: intermediate

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

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

276:   Collective on fem

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

281:   Level: beginner

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

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

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

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

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

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

324:   Collective

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

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

332:   Level: beginner

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

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

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

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

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

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

369:   Not collective

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

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

377:   Level: intermediate

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

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

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

395:   Not collective

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

401:   Level: intermediate

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

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

415:   Not collective

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

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

423:   Level: intermediate

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

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

438:   Not collective

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

447:   Level: intermediate

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

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

464:   Not collective

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

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

475:   Level: intermediate

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

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

496:   Not collective

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

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

504:   Level: intermediate

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

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

519:   Not collective

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

525:   Level: intermediate

527: .seealso: `PetscFECreate()`
528: @*/
529: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
530: {
533:   PetscSpaceDestroy(&fem->basisSpace);
534:   fem->basisSpace = sp;
535:   PetscObjectReference((PetscObject)fem->basisSpace);
536:   return 0;
537: }

539: /*@
540:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

542:   Not collective

544:   Input Parameter:
545: . fem - The PetscFE object

547:   Output Parameter:
548: . sp - The PetscDualSpace object

550:   Level: intermediate

552: .seealso: `PetscFECreate()`
553: @*/
554: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
555: {
558:   *sp = fem->dualSpace;
559:   return 0;
560: }

562: /*@
563:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

565:   Not collective

567:   Input Parameters:
568: + fem - The PetscFE object
569: - sp - The PetscDualSpace object

571:   Level: intermediate

573: .seealso: `PetscFECreate()`
574: @*/
575: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
576: {
579:   PetscDualSpaceDestroy(&fem->dualSpace);
580:   fem->dualSpace = sp;
581:   PetscObjectReference((PetscObject)fem->dualSpace);
582:   return 0;
583: }

585: /*@
586:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

588:   Not collective

590:   Input Parameter:
591: . fem - The PetscFE object

593:   Output Parameter:
594: . q - The PetscQuadrature object

596:   Level: intermediate

598: .seealso: `PetscFECreate()`
599: @*/
600: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
601: {
604:   *q = fem->quadrature;
605:   return 0;
606: }

608: /*@
609:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

611:   Not collective

613:   Input Parameters:
614: + fem - The PetscFE object
615: - q - The PetscQuadrature object

617:   Level: intermediate

619: .seealso: `PetscFECreate()`
620: @*/
621: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
622: {
623:   PetscInt Nc, qNc;

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

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

641:   Not collective

643:   Input Parameter:
644: . fem - The PetscFE object

646:   Output Parameter:
647: . q - The PetscQuadrature object

649:   Level: intermediate

651: .seealso: `PetscFECreate()`
652: @*/
653: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
654: {
657:   *q = fem->faceQuadrature;
658:   return 0;
659: }

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

664:   Not collective

666:   Input Parameters:
667: + fem - The PetscFE object
668: - q - The PetscQuadrature object

670:   Level: intermediate

672: .seealso: `PetscFECreate()`
673: @*/
674: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
675: {
676:   PetscInt Nc, qNc;

679:   PetscFEGetNumComponents(fem, &Nc);
680:   PetscQuadratureGetNumComponents(q, &qNc);
682:   PetscTabulationDestroy(&fem->Tf);
683:   PetscQuadratureDestroy(&fem->faceQuadrature);
684:   fem->faceQuadrature = q;
685:   PetscObjectReference((PetscObject)q);
686:   return 0;
687: }

689: /*@
690:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

692:   Not collective

694:   Input Parameters:
695: + sfe - The PetscFE source for the quadratures
696: - tfe - The PetscFE target for the quadratures

698:   Level: intermediate

700: .seealso: `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
701: @*/
702: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
703: {
704:   PetscQuadrature q;

708:   PetscFEGetQuadrature(sfe, &q);
709:   PetscFESetQuadrature(tfe, q);
710:   PetscFEGetFaceQuadrature(sfe, &q);
711:   PetscFESetFaceQuadrature(tfe, q);
712:   return 0;
713: }

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

718:   Not collective

720:   Input Parameter:
721: . fem - The PetscFE object

723:   Output Parameter:
724: . numDof - Array with the number of dofs per dimension

726:   Level: intermediate

728: .seealso: `PetscFECreate()`
729: @*/
730: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
731: {
734:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
735:   return 0;
736: }

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

741:   Not collective

743:   Input Parameters:
744: + fem - The PetscFE object
745: - k   - The highest derivative we need to tabulate, very often 1

747:   Output Parameter:
748: . T - The basis function values and derivatives at quadrature points

750:   Note:
751: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
752: $ 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
753: $ 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

755:   Level: intermediate

757: .seealso: `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
758: @*/
759: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
760: {
761:   PetscInt         npoints;
762:   const PetscReal *points;

766:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
767:   if (!fem->T) PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);
769:   *T = fem->T;
770:   return 0;
771: }

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

776:   Not collective

778:   Input Parameters:
779: + fem - The PetscFE object
780: - k   - The highest derivative we need to tabulate, very often 1

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

785:   Note:
786: $ 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
787: $ 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
788: $ 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

790:   Level: intermediate

792: .seealso: `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
793: @*/
794: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
795: {
798:   if (!fem->Tf) {
799:     const PetscReal  xi0[3] = {-1., -1., -1.};
800:     PetscReal        v0[3], J[9], detJ;
801:     PetscQuadrature  fq;
802:     PetscDualSpace   sp;
803:     DM               dm;
804:     const PetscInt  *faces;
805:     PetscInt         dim, numFaces, f, npoints, q;
806:     const PetscReal *points;
807:     PetscReal       *facePoints;

809:     PetscFEGetDualSpace(fem, &sp);
810:     PetscDualSpaceGetDM(sp, &dm);
811:     DMGetDimension(dm, &dim);
812:     DMPlexGetConeSize(dm, 0, &numFaces);
813:     DMPlexGetCone(dm, 0, &faces);
814:     PetscFEGetFaceQuadrature(fem, &fq);
815:     if (fq) {
816:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
817:       PetscMalloc1(numFaces * npoints * dim, &facePoints);
818:       for (f = 0; f < numFaces; ++f) {
819:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
820:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * npoints + q) * dim]);
821:       }
822:       PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
823:       PetscFree(facePoints);
824:     }
825:   }
827:   *Tf = fem->Tf;
828:   return 0;
829: }

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

834:   Not collective

836:   Input Parameter:
837: . fem - The PetscFE object

839:   Output Parameters:
840: . Tc - The basis function values at face centroid points

842:   Note:
843: $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c

845:   Level: intermediate

847: .seealso: `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
848: @*/
849: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
850: {
853:   if (!fem->Tc) {
854:     PetscDualSpace  sp;
855:     DM              dm;
856:     const PetscInt *cone;
857:     PetscReal      *centroids;
858:     PetscInt        dim, numFaces, f;

860:     PetscFEGetDualSpace(fem, &sp);
861:     PetscDualSpaceGetDM(sp, &dm);
862:     DMGetDimension(dm, &dim);
863:     DMPlexGetConeSize(dm, 0, &numFaces);
864:     DMPlexGetCone(dm, 0, &cone);
865:     PetscMalloc1(numFaces * dim, &centroids);
866:     for (f = 0; f < numFaces; ++f) DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f * dim], NULL);
867:     PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
868:     PetscFree(centroids);
869:   }
870:   *Tc = fem->Tc;
871:   return 0;
872: }

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

877:   Not collective

879:   Input Parameters:
880: + fem     - The PetscFE object
881: . nrepl   - The number of replicas
882: . npoints - The number of tabulation points in a replica
883: . points  - The tabulation point coordinates
884: - K       - The number of derivatives calculated

886:   Output Parameter:
887: . T - The basis function values and derivatives at tabulation points

889:   Note:
890: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
891: $ 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
892: $ 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

894:   Level: intermediate

896: .seealso: `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
897: @*/
898: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
899: {
900:   DM             dm;
901:   PetscDualSpace Q;
902:   PetscInt       Nb;   /* Dimension of FE space P */
903:   PetscInt       Nc;   /* Field components */
904:   PetscInt       cdim; /* Reference coordinate dimension */
905:   PetscInt       k;

907:   if (!npoints || !fem->dualSpace || K < 0) {
908:     *T = NULL;
909:     return 0;
910:   }
914:   PetscFEGetDualSpace(fem, &Q);
915:   PetscDualSpaceGetDM(Q, &dm);
916:   DMGetDimension(dm, &cdim);
917:   PetscDualSpaceGetDimension(Q, &Nb);
918:   PetscFEGetNumComponents(fem, &Nc);
919:   PetscMalloc1(1, T);
920:   (*T)->K    = !cdim ? 0 : K;
921:   (*T)->Nr   = nrepl;
922:   (*T)->Np   = npoints;
923:   (*T)->Nb   = Nb;
924:   (*T)->Nc   = Nc;
925:   (*T)->cdim = cdim;
926:   PetscMalloc1((*T)->K + 1, &(*T)->T);
927:   for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
928:   PetscUseTypeMethod(fem, createtabulation, nrepl * npoints, points, K, *T);
929:   return 0;
930: }

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

935:   Not collective

937:   Input Parameters:
938: + fem     - The PetscFE object
939: . npoints - The number of tabulation points
940: . points  - The tabulation point coordinates
941: . K       - The number of derivatives calculated
942: - T       - An existing tabulation object with enough allocated space

944:   Output Parameter:
945: . T - The basis function values and derivatives at tabulation points

947:   Note:
948: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
949: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
950: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

952:   Level: intermediate

954: .seealso: `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
955: @*/
956: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
957: {
959:   if (!npoints || !fem->dualSpace || K < 0) return 0;
963:   if (PetscDefined(USE_DEBUG)) {
964:     DM             dm;
965:     PetscDualSpace Q;
966:     PetscInt       Nb;   /* Dimension of FE space P */
967:     PetscInt       Nc;   /* Field components */
968:     PetscInt       cdim; /* Reference coordinate dimension */

970:     PetscFEGetDualSpace(fem, &Q);
971:     PetscDualSpaceGetDM(Q, &dm);
972:     DMGetDimension(dm, &cdim);
973:     PetscDualSpaceGetDimension(Q, &Nb);
974:     PetscFEGetNumComponents(fem, &Nc);
979:   }
980:   T->Nr = 1;
981:   T->Np = npoints;
982:   PetscUseTypeMethod(fem, createtabulation, npoints, points, K, T);
983:   return 0;
984: }

986: /*@C
987:   PetscTabulationDestroy - Frees memory from the associated tabulation.

989:   Not collective

991:   Input Parameter:
992: . T - The tabulation

994:   Level: intermediate

996: .seealso: `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
997: @*/
998: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
999: {
1000:   PetscInt k;

1003:   if (!T || !(*T)) return 0;
1004:   for (k = 0; k <= (*T)->K; ++k) PetscFree((*T)->T[k]);
1005:   PetscFree((*T)->T);
1006:   PetscFree(*T);
1007:   *T = NULL;
1008:   return 0;
1009: }

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

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

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

1075:   *trFE = NULL;
1076:   PetscFEGetDualSpace(fe, &dsp);
1077:   PetscDualSpaceGetDM(dsp, &dm);
1078:   DMPlexGetHeightStratum(dm, height, &hStart, &hEnd);
1079:   if (hEnd <= hStart) return 0;
1080:   PetscFECreatePointTrace(fe, hStart, trFE);
1081:   return 0;
1082: }

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

1087:   Not collective

1089:   Input Parameter:
1090: . fe - The PetscFE

1092:   Output Parameter:
1093: . dim - The dimension

1095:   Level: intermediate

1097: .seealso: `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1098: @*/
1099: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1100: {
1103:   PetscTryTypeMethod(fem, getdimension, dim);
1104:   return 0;
1105: }

1107: /*@C
1108:   PetscFEPushforward - Map the reference element function to real space

1110:   Input Parameters:
1111: + fe     - The PetscFE
1112: . fegeom - The cell geometry
1113: . Nv     - The number of function values
1114: - vals   - The function values

1116:   Output Parameter:
1117: . vals   - The transformed function values

1119:   Level: advanced

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

1123:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1125: .seealso: `PetscDualSpacePushforward()`
1126: @*/
1127: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1128: {
1130:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1131:   return 0;
1132: }

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

1137:   Input Parameters:
1138: + fe     - The PetscFE
1139: . fegeom - The cell geometry
1140: . Nv     - The number of function gradient values
1141: - vals   - The function gradient values

1143:   Output Parameter:
1144: . vals   - The transformed function gradient values

1146:   Level: advanced

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

1150:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1152: .seealso: `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1153: @*/
1154: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1155: {
1157:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1158:   return 0;
1159: }

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

1164:   Input Parameters:
1165: + fe     - The PetscFE
1166: . fegeom - The cell geometry
1167: . Nv     - The number of function Hessian values
1168: - vals   - The function Hessian values

1170:   Output Parameter:
1171: . vals   - The transformed function Hessian values

1173:   Level: advanced

1175:   Note: This just forwards the call onto PetscDualSpacePushforwardHessian().

1177:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1179: .seealso: `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1180: @*/
1181: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1182: {
1184:   PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1185:   return 0;
1186: }

1188: /*
1189: Purpose: Compute element vector for chunk of elements

1191: Input:
1192:   Sizes:
1193:      Ne:  number of elements
1194:      Nf:  number of fields
1195:      PetscFE
1196:        dim: spatial dimension
1197:        Nb:  number of basis functions
1198:        Nc:  number of field components
1199:        PetscQuadrature
1200:          Nq:  number of quadrature points

1202:   Geometry:
1203:      PetscFEGeom[Ne] possibly *Nq
1204:        PetscReal v0s[dim]
1205:        PetscReal n[dim]
1206:        PetscReal jacobians[dim*dim]
1207:        PetscReal jacobianInverses[dim*dim]
1208:        PetscReal jacobianDeterminants
1209:   FEM:
1210:      PetscFE
1211:        PetscQuadrature
1212:          PetscReal   quadPoints[Nq*dim]
1213:          PetscReal   quadWeights[Nq]
1214:        PetscReal   basis[Nq*Nb*Nc]
1215:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1216:      PetscScalar coefficients[Ne*Nb*Nc]
1217:      PetscScalar elemVec[Ne*Nb*Nc]

1219:   Problem:
1220:      PetscInt f: the active field
1221:      f0, f1

1223:   Work Space:
1224:      PetscFE
1225:        PetscScalar f0[Nq*dim];
1226:        PetscScalar f1[Nq*dim*dim];
1227:        PetscScalar u[Nc];
1228:        PetscScalar gradU[Nc*dim];
1229:        PetscReal   x[dim];
1230:        PetscScalar realSpaceDer[dim];

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

1234: Input:
1235:   Sizes:
1236:      N_cb: Number of serial cell batches

1238:   Geometry:
1239:      PetscReal v0s[Ne*dim]
1240:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1241:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1242:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1243:   FEM:
1244:      static PetscReal   quadPoints[Nq*dim]
1245:      static PetscReal   quadWeights[Nq]
1246:      static PetscReal   basis[Nq*Nb*Nc]
1247:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1248:      PetscScalar coefficients[Ne*Nb*Nc]
1249:      PetscScalar elemVec[Ne*Nb*Nc]

1251: ex62.c:
1252:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1253:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1254:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1255:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1257: ex52.c:
1258:   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)
1259:   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)

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

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

1268: ex52_integrateElementOpenCL.c:
1269: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1270:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1271:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1279:   Not collective

1281:   Input Parameters:
1282: + prob         - The PetscDS specifying the discretizations and continuum functions
1283: . field        - The field being integrated
1284: . Ne           - The number of elements in the chunk
1285: . cgeom        - The cell geometry for each cell in the chunk
1286: . coefficients - The array of FEM basis coefficients for the elements
1287: . probAux      - The PetscDS specifying the auxiliary discretizations
1288: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1290:   Output Parameter:
1291: . integral     - the integral for this field

1293:   Level: intermediate

1295: .seealso: `PetscFEIntegrateResidual()`
1296: @*/
1297: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1298: {
1299:   PetscFE fe;

1302:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1303:   if (fe->ops->integrate) (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);
1304:   return 0;
1305: }

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

1310:   Not collective

1312:   Input Parameters:
1313: + prob         - The PetscDS specifying the discretizations and continuum functions
1314: . field        - The field being integrated
1315: . obj_func     - The function to be integrated
1316: . Ne           - The number of elements in the chunk
1317: . fgeom        - The face geometry for each face in the chunk
1318: . coefficients - The array of FEM basis coefficients for the elements
1319: . probAux      - The PetscDS specifying the auxiliary discretizations
1320: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1322:   Output Parameter:
1323: . integral     - the integral for this field

1325:   Level: intermediate

1327: .seealso: `PetscFEIntegrateResidual()`
1328: @*/
1329: 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[])
1330: {
1331:   PetscFE fe;

1334:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1335:   if (fe->ops->integratebd) (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);
1336:   return 0;
1337: }

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

1342:   Not collective

1344:   Input Parameters:
1345: + ds           - The PetscDS specifying the discretizations and continuum functions
1346: . key          - The (label+value, field) being integrated
1347: . Ne           - The number of elements in the chunk
1348: . cgeom        - The cell geometry for each cell in the chunk
1349: . coefficients - The array of FEM basis coefficients for the elements
1350: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1351: . probAux      - The PetscDS specifying the auxiliary discretizations
1352: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1353: - t            - The time

1355:   Output Parameter:
1356: . elemVec      - the element residual vectors from each element

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

1366:   Level: intermediate

1368: .seealso: `PetscFEIntegrateResidual()`
1369: @*/
1370: 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[])
1371: {
1372:   PetscFE fe;

1376:   PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1377:   if (fe->ops->integrateresidual) (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1378:   return 0;
1379: }

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

1384:   Not collective

1386:   Input Parameters:
1387: + ds           - The PetscDS specifying the discretizations and continuum functions
1388: . wf           - The PetscWeakForm object holding the pointwise functions
1389: . key          - The (label+value, field) being integrated
1390: . Ne           - The number of elements in the chunk
1391: . fgeom        - The face geometry for each cell in the chunk
1392: . coefficients - The array of FEM basis coefficients for the elements
1393: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1394: . probAux      - The PetscDS specifying the auxiliary discretizations
1395: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1396: - t            - The time

1398:   Output Parameter:
1399: . elemVec      - the element residual vectors from each element

1401:   Level: intermediate

1403: .seealso: `PetscFEIntegrateResidual()`
1404: @*/
1405: 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[])
1406: {
1407:   PetscFE fe;

1410:   PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1411:   if (fe->ops->integratebdresidual) (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1412:   return 0;
1413: }

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

1418:   Not collective

1420:   Input Parameters:
1421: + prob         - The PetscDS specifying the discretizations and continuum functions
1422: . key          - The (label+value, field) being integrated
1423: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1424: . Ne           - The number of elements in the chunk
1425: . fgeom        - The face geometry for each cell in the chunk
1426: . coefficients - The array of FEM basis coefficients for the elements
1427: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1428: . probAux      - The PetscDS specifying the auxiliary discretizations
1429: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1430: - t            - The time

1432:   Output Parameter
1433: . elemVec      - the element residual vectors from each element

1435:   Level: developer

1437: .seealso: `PetscFEIntegrateResidual()`
1438: @*/
1439: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1440: {
1441:   PetscFE fe;

1444:   PetscDSGetDiscretization(prob, key.field, (PetscObject *)&fe);
1445:   if (fe->ops->integratehybridresidual) (*fe->ops->integratehybridresidual)(prob, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1446:   return 0;
1447: }

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

1452:   Not collective

1454:   Input Parameters:
1455: + ds           - The PetscDS specifying the discretizations and continuum functions
1456: . jtype        - The type of matrix pointwise functions that should be used
1457: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1458: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1459: . Ne           - The number of elements in the chunk
1460: . cgeom        - The cell geometry for each cell in the chunk
1461: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1462: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1463: . probAux      - The PetscDS specifying the auxiliary discretizations
1464: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1465: . t            - The time
1466: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1468:   Output Parameter:
1469: . elemMat      - the element matrices for the Jacobian from each element

1471:   Note:
1472: $ Loop over batch of elements (e):
1473: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1474: $     Loop over quadrature points (q):
1475: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1476: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1477: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1478: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1479: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1480:   Level: intermediate

1482: .seealso: `PetscFEIntegrateResidual()`
1483: @*/
1484: 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[])
1485: {
1486:   PetscFE  fe;
1487:   PetscInt Nf;

1490:   PetscDSGetNumFields(ds, &Nf);
1491:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1492:   if (fe->ops->integratejacobian) (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1493:   return 0;
1494: }

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

1499:   Not collective

1501:   Input Parameters:
1502: + ds           - The PetscDS specifying the discretizations and continuum functions
1503: . wf           - The PetscWeakForm holding the pointwise functions
1504: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1505: . Ne           - The number of elements in the chunk
1506: . fgeom        - The face geometry for each cell in the chunk
1507: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1508: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1509: . probAux      - The PetscDS specifying the auxiliary discretizations
1510: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1511: . t            - The time
1512: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1514:   Output Parameter:
1515: . elemMat              - the element matrices for the Jacobian from each element

1517:   Note:
1518: $ Loop over batch of elements (e):
1519: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1520: $     Loop over quadrature points (q):
1521: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1522: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1523: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1524: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1525: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1526:   Level: intermediate

1528: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1529: @*/
1530: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1531: {
1532:   PetscFE  fe;
1533:   PetscInt Nf;

1536:   PetscDSGetNumFields(ds, &Nf);
1537:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1538:   if (fe->ops->integratebdjacobian) (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1539:   return 0;
1540: }

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

1545:   Not collective

1547:   Input Parameters:
1548: + ds           - The PetscDS specifying the discretizations and continuum functions
1549: . jtype        - The type of matrix pointwise functions that should be used
1550: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1551: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1552: . Ne           - The number of elements in the chunk
1553: . fgeom        - The face geometry for each cell in the chunk
1554: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1555: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1556: . probAux      - The PetscDS specifying the auxiliary discretizations
1557: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1558: . t            - The time
1559: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1561:   Output Parameter
1562: . elemMat              - the element matrices for the Jacobian from each element

1564:   Note:
1565: $ Loop over batch of elements (e):
1566: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1567: $     Loop over quadrature points (q):
1568: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1569: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1570: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1571: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1572: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1573:   Level: developer

1575: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1576: @*/
1577: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1578: {
1579:   PetscFE  fe;
1580:   PetscInt Nf;

1583:   PetscDSGetNumFields(ds, &Nf);
1584:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1585:   if (fe->ops->integratehybridjacobian) (*fe->ops->integratehybridjacobian)(ds, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1586:   return 0;
1587: }

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

1592:   Input Parameters:
1593: + fe     - The finite element space
1594: - height - The height of the Plex point

1596:   Output Parameter:
1597: . subfe  - The subspace of this FE space

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

1601:   Level: advanced

1603: .seealso: `PetscFECreateDefault()`
1604: @*/
1605: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1606: {
1607:   PetscSpace      P, subP;
1608:   PetscDualSpace  Q, subQ;
1609:   PetscQuadrature subq;
1610:   PetscFEType     fetype;
1611:   PetscInt        dim, Nc;

1615:   if (height == 0) {
1616:     *subfe = fe;
1617:     return 0;
1618:   }
1619:   PetscFEGetBasisSpace(fe, &P);
1620:   PetscFEGetDualSpace(fe, &Q);
1621:   PetscFEGetNumComponents(fe, &Nc);
1622:   PetscFEGetFaceQuadrature(fe, &subq);
1623:   PetscDualSpaceGetDimension(Q, &dim);
1625:   if (!fe->subspaces) PetscCalloc1(dim, &fe->subspaces);
1626:   if (height <= dim) {
1627:     if (!fe->subspaces[height - 1]) {
1628:       PetscFE     sub = NULL;
1629:       const char *name;

1631:       PetscSpaceGetHeightSubspace(P, height, &subP);
1632:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1633:       if (subQ) {
1634:         PetscFECreate(PetscObjectComm((PetscObject)fe), &sub);
1635:         PetscObjectGetName((PetscObject)fe, &name);
1636:         PetscObjectSetName((PetscObject)sub, name);
1637:         PetscFEGetType(fe, &fetype);
1638:         PetscFESetType(sub, fetype);
1639:         PetscFESetBasisSpace(sub, subP);
1640:         PetscFESetDualSpace(sub, subQ);
1641:         PetscFESetNumComponents(sub, Nc);
1642:         PetscFESetUp(sub);
1643:         PetscFESetQuadrature(sub, subq);
1644:       }
1645:       fe->subspaces[height - 1] = sub;
1646:     }
1647:     *subfe = fe->subspaces[height - 1];
1648:   } else {
1649:     *subfe = NULL;
1650:   }
1651:   return 0;
1652: }

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

1659:   Collective on fem

1661:   Input Parameter:
1662: . fe - The initial PetscFE

1664:   Output Parameter:
1665: . feRef - The refined PetscFE

1667:   Level: advanced

1669: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1670: @*/
1671: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1672: {
1673:   PetscSpace       P, Pref;
1674:   PetscDualSpace   Q, Qref;
1675:   DM               K, Kref;
1676:   PetscQuadrature  q, qref;
1677:   const PetscReal *v0, *jac;
1678:   PetscInt         numComp, numSubelements;
1679:   PetscInt         cStart, cEnd, c;
1680:   PetscDualSpace  *cellSpaces;

1682:   PetscFEGetBasisSpace(fe, &P);
1683:   PetscFEGetDualSpace(fe, &Q);
1684:   PetscFEGetQuadrature(fe, &q);
1685:   PetscDualSpaceGetDM(Q, &K);
1686:   /* Create space */
1687:   PetscObjectReference((PetscObject)P);
1688:   Pref = P;
1689:   /* Create dual space */
1690:   PetscDualSpaceDuplicate(Q, &Qref);
1691:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1692:   DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref);
1693:   PetscDualSpaceSetDM(Qref, Kref);
1694:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1695:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1696:   /* TODO: fix for non-uniform refinement */
1697:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1698:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1699:   PetscFree(cellSpaces);
1700:   DMDestroy(&Kref);
1701:   PetscDualSpaceSetUp(Qref);
1702:   /* Create element */
1703:   PetscFECreate(PetscObjectComm((PetscObject)fe), feRef);
1704:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1705:   PetscFESetBasisSpace(*feRef, Pref);
1706:   PetscFESetDualSpace(*feRef, Qref);
1707:   PetscFEGetNumComponents(fe, &numComp);
1708:   PetscFESetNumComponents(*feRef, numComp);
1709:   PetscFESetUp(*feRef);
1710:   PetscSpaceDestroy(&Pref);
1711:   PetscDualSpaceDestroy(&Qref);
1712:   /* Create quadrature */
1713:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1714:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1715:   PetscFESetQuadrature(*feRef, qref);
1716:   PetscQuadratureDestroy(&qref);
1717:   return 0;
1718: }

1720: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1721: {
1722:   PetscSpace     P;
1723:   PetscDualSpace Q;
1724:   DM             K;
1725:   DMPolytopeType ct;
1726:   PetscInt       degree;
1727:   char           name[64];

1729:   PetscFEGetBasisSpace(fe, &P);
1730:   PetscSpaceGetDegree(P, &degree, NULL);
1731:   PetscFEGetDualSpace(fe, &Q);
1732:   PetscDualSpaceGetDM(Q, &K);
1733:   DMPlexGetCellType(K, 0, &ct);
1734:   switch (ct) {
1735:   case DM_POLYTOPE_SEGMENT:
1736:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1737:   case DM_POLYTOPE_QUADRILATERAL:
1738:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1739:   case DM_POLYTOPE_HEXAHEDRON:
1740:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1741:     PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree);
1742:     break;
1743:   case DM_POLYTOPE_TRIANGLE:
1744:   case DM_POLYTOPE_TETRAHEDRON:
1745:     PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree);
1746:     break;
1747:   case DM_POLYTOPE_TRI_PRISM:
1748:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1749:     PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree);
1750:     break;
1751:   default:
1752:     PetscSNPrintf(name, sizeof(name), "FE");
1753:   }
1754:   PetscFESetName(fe, name);
1755:   return 0;
1756: }

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

1762:   switch (ct) {
1763:   case DM_POLYTOPE_SEGMENT:
1764:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1765:   case DM_POLYTOPE_QUADRILATERAL:
1766:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1767:   case DM_POLYTOPE_HEXAHEDRON:
1768:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1769:     PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1770:     PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1771:     break;
1772:   case DM_POLYTOPE_TRIANGLE:
1773:   case DM_POLYTOPE_TETRAHEDRON:
1774:     PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1775:     PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1776:     break;
1777:   case DM_POLYTOPE_TRI_PRISM:
1778:   case DM_POLYTOPE_TRI_PRISM_TENSOR: {
1779:     PetscQuadrature q1, q2;

1781:     PetscDTStroudConicalQuadrature(2, 1, quadPointsPerEdge, -1.0, 1.0, &q1);
1782:     PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2);
1783:     PetscDTTensorQuadratureCreate(q1, q2, q);
1784:     PetscQuadratureDestroy(&q1);
1785:     PetscQuadratureDestroy(&q2);
1786:   }
1787:     PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1788:     /* TODO Need separate quadratures for each face */
1789:     break;
1790:   default:
1791:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
1792:   }
1793:   return 0;
1794: }

1796: /*@
1797:   PetscFECreateFromSpaces - Create a PetscFE from the basis and dual spaces

1799:   Collective

1801:   Input Parameters:
1802: + P  - The basis space
1803: . Q  - The dual space
1804: . q  - The cell quadrature
1805: - fq - The face quadrature

1807:   Output Parameter:
1808: . fem    - The PetscFE object

1810:   Note:
1811:   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.

1813:   Level: beginner

1815: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1816: @*/
1817: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1818: {
1819:   PetscInt    Nc;
1820:   const char *prefix;

1822:   PetscFECreate(PetscObjectComm((PetscObject)P), fem);
1823:   PetscObjectGetOptionsPrefix((PetscObject)P, &prefix);
1824:   PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix);
1825:   PetscFESetType(*fem, PETSCFEBASIC);
1826:   PetscFESetBasisSpace(*fem, P);
1827:   PetscFESetDualSpace(*fem, Q);
1828:   PetscSpaceGetNumComponents(P, &Nc);
1829:   PetscFESetNumComponents(*fem, Nc);
1830:   PetscFESetUp(*fem);
1831:   PetscSpaceDestroy(&P);
1832:   PetscDualSpaceDestroy(&Q);
1833:   PetscFESetQuadrature(*fem, q);
1834:   PetscFESetFaceQuadrature(*fem, fq);
1835:   PetscQuadratureDestroy(&q);
1836:   PetscQuadratureDestroy(&fq);
1837:   PetscFESetDefaultName_Private(*fem);
1838:   return 0;
1839: }

1841: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1842: {
1843:   DM              K;
1844:   PetscSpace      P;
1845:   PetscDualSpace  Q;
1846:   PetscQuadrature q, fq;
1847:   PetscBool       tensor;

1851:   switch (ct) {
1852:   case DM_POLYTOPE_SEGMENT:
1853:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1854:   case DM_POLYTOPE_QUADRILATERAL:
1855:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1856:   case DM_POLYTOPE_HEXAHEDRON:
1857:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1858:     tensor = PETSC_TRUE;
1859:     break;
1860:   default:
1861:     tensor = PETSC_FALSE;
1862:   }
1863:   /* Create space */
1864:   PetscSpaceCreate(comm, &P);
1865:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1866:   PetscObjectSetOptionsPrefix((PetscObject)P, prefix);
1867:   PetscSpacePolynomialSetTensor(P, tensor);
1868:   PetscSpaceSetNumComponents(P, Nc);
1869:   PetscSpaceSetNumVariables(P, dim);
1870:   if (degree >= 0) {
1871:     PetscSpaceSetDegree(P, degree, PETSC_DETERMINE);
1872:     if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1873:       PetscSpace Pend, Pside;

1875:       PetscSpaceCreate(comm, &Pend);
1876:       PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL);
1877:       PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE);
1878:       PetscSpaceSetNumComponents(Pend, Nc);
1879:       PetscSpaceSetNumVariables(Pend, dim - 1);
1880:       PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE);
1881:       PetscSpaceCreate(comm, &Pside);
1882:       PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL);
1883:       PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE);
1884:       PetscSpaceSetNumComponents(Pside, 1);
1885:       PetscSpaceSetNumVariables(Pside, 1);
1886:       PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE);
1887:       PetscSpaceSetType(P, PETSCSPACETENSOR);
1888:       PetscSpaceTensorSetNumSubspaces(P, 2);
1889:       PetscSpaceTensorSetSubspace(P, 0, Pend);
1890:       PetscSpaceTensorSetSubspace(P, 1, Pside);
1891:       PetscSpaceDestroy(&Pend);
1892:       PetscSpaceDestroy(&Pside);
1893:     }
1894:   }
1895:   if (setFromOptions) PetscSpaceSetFromOptions(P);
1896:   PetscSpaceSetUp(P);
1897:   PetscSpaceGetDegree(P, &degree, NULL);
1898:   PetscSpacePolynomialGetTensor(P, &tensor);
1899:   PetscSpaceGetNumComponents(P, &Nc);
1900:   /* Create dual space */
1901:   PetscDualSpaceCreate(comm, &Q);
1902:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1903:   PetscObjectSetOptionsPrefix((PetscObject)Q, prefix);
1904:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1905:   PetscDualSpaceSetDM(Q, K);
1906:   DMDestroy(&K);
1907:   PetscDualSpaceSetNumComponents(Q, Nc);
1908:   PetscDualSpaceSetOrder(Q, degree);
1909:   /* TODO For some reason, we need a tensor dualspace with wedges */
1910:   PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE);
1911:   if (setFromOptions) PetscDualSpaceSetFromOptions(Q);
1912:   PetscDualSpaceSetUp(Q);
1913:   /* Create quadrature */
1914:   qorder = qorder >= 0 ? qorder : degree;
1915:   if (setFromOptions) {
1916:     PetscObjectOptionsBegin((PetscObject)P);
1917:     PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0);
1918:     PetscOptionsEnd();
1919:   }
1920:   PetscFECreateDefaultQuadrature_Private(dim, ct, qorder, &q, &fq);
1921:   /* Create finite element */
1922:   PetscFECreateFromSpaces(P, Q, q, fq, fem);
1923:   if (setFromOptions) PetscFESetFromOptions(*fem);
1924:   return 0;
1925: }

1927: /*@C
1928:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1930:   Collective

1932:   Input Parameters:
1933: + comm      - The MPI comm
1934: . dim       - The spatial dimension
1935: . Nc        - The number of components
1936: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1937: . prefix    - The options prefix, or NULL
1938: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1940:   Output Parameter:
1941: . fem - The PetscFE object

1943:   Note:
1944:   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.

1946:   Level: beginner

1948: .seealso: `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1949: @*/
1950: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1951: {
1952:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1953:   return 0;
1954: }

1956: /*@C
1957:   PetscFECreateByCell - Create a PetscFE for basic FEM computation

1959:   Collective

1961:   Input Parameters:
1962: + comm   - The MPI comm
1963: . dim    - The spatial dimension
1964: . Nc     - The number of components
1965: . ct     - The celltype of the reference cell
1966: . prefix - The options prefix, or NULL
1967: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1969:   Output Parameter:
1970: . fem - The PetscFE object

1972:   Note:
1973:   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.

1975:   Level: beginner

1977: .seealso: `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1978: @*/
1979: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
1980: {
1981:   PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1982:   return 0;
1983: }

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

1988:   Collective

1990:   Input Parameters:
1991: + comm      - The MPI comm
1992: . dim       - The spatial dimension
1993: . Nc        - The number of components
1994: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1995: . k         - The degree k of the space
1996: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1998:   Output Parameter:
1999: . fem       - The PetscFE object

2001:   Level: beginner

2003:   Notes:
2004:   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.

2006: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2007: @*/
2008: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2009: {
2010:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem);
2011:   return 0;
2012: }

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

2017:   Collective

2019:   Input Parameters:
2020: + comm      - The MPI comm
2021: . dim       - The spatial dimension
2022: . Nc        - The number of components
2023: . ct        - The celltype of the reference cell
2024: . k         - The degree k of the space
2025: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

2027:   Output Parameter:
2028: . fem       - The PetscFE object

2030:   Level: beginner

2032:   Notes:
2033:   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.

2035: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2036: @*/
2037: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2038: {
2039:   PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem);
2040:   return 0;
2041: }

2043: /*@C
2044:   PetscFESetName - Names the FE and its subobjects

2046:   Not collective

2048:   Input Parameters:
2049: + fe   - The PetscFE
2050: - name - The name

2052:   Level: intermediate

2054: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2055: @*/
2056: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2057: {
2058:   PetscSpace     P;
2059:   PetscDualSpace Q;

2061:   PetscFEGetBasisSpace(fe, &P);
2062:   PetscFEGetDualSpace(fe, &Q);
2063:   PetscObjectSetName((PetscObject)fe, name);
2064:   PetscObjectSetName((PetscObject)P, name);
2065:   PetscObjectSetName((PetscObject)Q, name);
2066:   return 0;
2067: }

2069: 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[])
2070: {
2071:   PetscInt dOffset = 0, fOffset = 0, f, g;

2073:   for (f = 0; f < Nf; ++f) {
2074:     PetscFE          fe;
2075:     const PetscInt   k       = ds->jetDegree[f];
2076:     const PetscInt   cdim    = T[f]->cdim;
2077:     const PetscInt   Nq      = T[f]->Np;
2078:     const PetscInt   Nbf     = T[f]->Nb;
2079:     const PetscInt   Ncf     = T[f]->Nc;
2080:     const PetscReal *Bq      = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2081:     const PetscReal *Dq      = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2082:     const PetscReal *Hq      = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2083:     PetscInt         hOffset = 0, b, c, d;

2085:     PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
2086:     for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2087:     for (d = 0; d < cdim * Ncf; ++d) u_x[fOffset * cdim + d] = 0.0;
2088:     for (b = 0; b < Nbf; ++b) {
2089:       for (c = 0; c < Ncf; ++c) {
2090:         const PetscInt cidx = b * Ncf + c;

2092:         u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2093:         for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * cdim + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2094:       }
2095:     }
2096:     if (k > 1) {
2097:       for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * cdim;
2098:       for (d = 0; d < cdim * cdim * Ncf; ++d) u_x[hOffset + fOffset * cdim * cdim + d] = 0.0;
2099:       for (b = 0; b < Nbf; ++b) {
2100:         for (c = 0; c < Ncf; ++c) {
2101:           const PetscInt cidx = b * Ncf + c;

2103:           for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * cdim * cdim + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2104:         }
2105:       }
2106:       PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * cdim * cdim]);
2107:     }
2108:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2109:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * cdim]);
2110:     if (u_t) {
2111:       for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2112:       for (b = 0; b < Nbf; ++b) {
2113:         for (c = 0; c < Ncf; ++c) {
2114:           const PetscInt cidx = b * Ncf + c;

2116:           u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2117:         }
2118:       }
2119:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2120:     }
2121:     fOffset += Ncf;
2122:     dOffset += Nbf;
2123:   }
2124:   return 0;
2125: }

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

2131:   /* f is the field number in the DS, g is the field number in u[] */
2132:   for (f = 0, g = 0; f < Nf; ++f) {
2133:     PetscFE          fe  = (PetscFE)ds->disc[f];
2134:     const PetscInt   dEt = T[f]->cdim;
2135:     const PetscInt   dE  = fegeom->dimEmbed;
2136:     const PetscInt   Nq  = T[f]->Np;
2137:     const PetscInt   Nbf = T[f]->Nb;
2138:     const PetscInt   Ncf = T[f]->Nc;
2139:     const PetscReal *Bq  = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2140:     const PetscReal *Dq  = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2141:     PetscBool        isCohesive;
2142:     PetscInt         Ns, s;

2144:     if (!T[f]) continue;
2145:     PetscDSGetCohesive(ds, f, &isCohesive);
2146:     Ns = isCohesive ? 1 : 2;
2147:     for (s = 0; s < Ns; ++s, ++g) {
2148:       PetscInt b, c, d;

2150:       for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2151:       for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2152:       for (b = 0; b < Nbf; ++b) {
2153:         for (c = 0; c < Ncf; ++c) {
2154:           const PetscInt cidx = b * Ncf + c;

2156:           u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2157:           for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2158:         }
2159:       }
2160:       PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2161:       PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]);
2162:       if (u_t) {
2163:         for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2164:         for (b = 0; b < Nbf; ++b) {
2165:           for (c = 0; c < Ncf; ++c) {
2166:             const PetscInt cidx = b * Ncf + c;

2168:             u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2169:           }
2170:         }
2171:         PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2172:       }
2173:       fOffset += Ncf;
2174:       dOffset += Nbf;
2175:     }
2176:   }
2177:   return 0;
2178: }

2180: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2181: {
2182:   PetscFE         fe;
2183:   PetscTabulation Tc;
2184:   PetscInt        b, c;

2186:   if (!prob) return 0;
2187:   PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
2188:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2189:   {
2190:     const PetscReal *faceBasis = Tc->T[0];
2191:     const PetscInt   Nb        = Tc->Nb;
2192:     const PetscInt   Nc        = Tc->Nc;

2194:     for (c = 0; c < Nc; ++c) u[c] = 0.0;
2195:     for (b = 0; b < Nb; ++b) {
2196:       for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2197:     }
2198:   }
2199:   return 0;
2200: }

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

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

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

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

2240: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2241: {
2242:   const PetscInt   dE       = T->cdim;
2243:   const PetscInt   Nq       = T->Np;
2244:   const PetscInt   Nb       = T->Nb;
2245:   const PetscInt   Nc       = T->Nc;
2246:   const PetscReal *basis    = &T->T[0][r * Nq * Nb * Nc];
2247:   const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];
2248:   PetscInt         q, b, c, d;

2250:   for (q = 0; q < Nq; ++q) {
2251:     for (b = 0; b < Nb; ++b) {
2252:       for (c = 0; c < Nc; ++c) {
2253:         const PetscInt bcidx = b * Nc + c;

2255:         tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2256:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2257:       }
2258:     }
2259:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2260:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2261:     for (b = 0; b < Nb; ++b) {
2262:       for (c = 0; c < Nc; ++c) {
2263:         const PetscInt bcidx = b * Nc + c;
2264:         const PetscInt qcidx = q * Nc + c;

2266:         elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2267:         for (d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2268:       }
2269:     }
2270:   }
2271:   return (0);
2272: }

2274: 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[])
2275: {
2276:   const PetscInt   dE        = TI->cdim;
2277:   const PetscInt   NqI       = TI->Np;
2278:   const PetscInt   NbI       = TI->Nb;
2279:   const PetscInt   NcI       = TI->Nc;
2280:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2281:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2282:   const PetscInt   NqJ       = TJ->Np;
2283:   const PetscInt   NbJ       = TJ->Nb;
2284:   const PetscInt   NcJ       = TJ->Nc;
2285:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2286:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2287:   PetscInt         f, fc, g, gc, df, dg;

2289:   for (f = 0; f < NbI; ++f) {
2290:     for (fc = 0; fc < NcI; ++fc) {
2291:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2293:       tmpBasisI[fidx] = basisI[fidx];
2294:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2295:     }
2296:   }
2297:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2298:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2299:   for (g = 0; g < NbJ; ++g) {
2300:     for (gc = 0; gc < NcJ; ++gc) {
2301:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2303:       tmpBasisJ[gidx] = basisJ[gidx];
2304:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2305:     }
2306:   }
2307:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2308:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2309:   for (f = 0; f < NbI; ++f) {
2310:     for (fc = 0; fc < NcI; ++fc) {
2311:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2312:       const PetscInt i    = offsetI + f;  /* Element matrix row */
2313:       for (g = 0; g < NbJ; ++g) {
2314:         for (gc = 0; gc < NcJ; ++gc) {
2315:           const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2316:           const PetscInt j    = offsetJ + g;  /* Element matrix column */
2317:           const PetscInt fOff = eOffset + i * totDim + j;

2319:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2320:           for (df = 0; df < dE; ++df) {
2321:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2322:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2323:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2324:           }
2325:         }
2326:       }
2327:     }
2328:   }
2329:   return (0);
2330: }

2332: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt s, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2333: {
2334:   const PetscInt   dE        = TI->cdim;
2335:   const PetscInt   NqI       = TI->Np;
2336:   const PetscInt   NbI       = TI->Nb;
2337:   const PetscInt   NcI       = TI->Nc;
2338:   const PetscReal *basisI    = &TI->T[0][(r * NqI + q) * NbI * NcI];
2339:   const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2340:   const PetscInt   NqJ       = TJ->Np;
2341:   const PetscInt   NbJ       = TJ->Nb;
2342:   const PetscInt   NcJ       = TJ->Nc;
2343:   const PetscReal *basisJ    = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2344:   const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2345:   const PetscInt   so        = isHybridI ? 0 : s;
2346:   const PetscInt   to        = isHybridJ ? 0 : s;
2347:   PetscInt         f, fc, g, gc, df, dg;

2349:   for (f = 0; f < NbI; ++f) {
2350:     for (fc = 0; fc < NcI; ++fc) {
2351:       const PetscInt fidx = f * NcI + fc; /* Test function basis index */

2353:       tmpBasisI[fidx] = basisI[fidx];
2354:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2355:     }
2356:   }
2357:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2358:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2359:   for (g = 0; g < NbJ; ++g) {
2360:     for (gc = 0; gc < NcJ; ++gc) {
2361:       const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */

2363:       tmpBasisJ[gidx] = basisJ[gidx];
2364:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2365:     }
2366:   }
2367:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2368:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2369:   for (f = 0; f < NbI; ++f) {
2370:     for (fc = 0; fc < NcI; ++fc) {
2371:       const PetscInt fidx = f * NcI + fc;           /* Test function basis index */
2372:       const PetscInt i    = offsetI + NbI * so + f; /* Element matrix row */
2373:       for (g = 0; g < NbJ; ++g) {
2374:         for (gc = 0; gc < NcJ; ++gc) {
2375:           const PetscInt gidx = g * NcJ + gc;           /* Trial function basis index */
2376:           const PetscInt j    = offsetJ + NbJ * to + g; /* Element matrix column */
2377:           const PetscInt fOff = eOffset + i * totDim + j;

2379:           elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2380:           for (df = 0; df < dE; ++df) {
2381:             elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2382:             elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2383:             for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2384:           }
2385:         }
2386:       }
2387:     }
2388:   }
2389:   return (0);
2390: }

2392: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2393: {
2394:   PetscDualSpace  dsp;
2395:   DM              dm;
2396:   PetscQuadrature quadDef;
2397:   PetscInt        dim, cdim, Nq;

2399:   PetscFEGetDualSpace(fe, &dsp);
2400:   PetscDualSpaceGetDM(dsp, &dm);
2401:   DMGetDimension(dm, &dim);
2402:   DMGetCoordinateDim(dm, &cdim);
2403:   PetscFEGetQuadrature(fe, &quadDef);
2404:   quad = quad ? quad : quadDef;
2405:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2406:   PetscMalloc1(Nq * cdim, &cgeom->v);
2407:   PetscMalloc1(Nq * cdim * cdim, &cgeom->J);
2408:   PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ);
2409:   PetscMalloc1(Nq, &cgeom->detJ);
2410:   cgeom->dim       = dim;
2411:   cgeom->dimEmbed  = cdim;
2412:   cgeom->numCells  = 1;
2413:   cgeom->numPoints = Nq;
2414:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2415:   return 0;
2416: }

2418: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2419: {
2420:   PetscFree(cgeom->v);
2421:   PetscFree(cgeom->J);
2422:   PetscFree(cgeom->invJ);
2423:   PetscFree(cgeom->detJ);
2424:   return 0;
2425: }