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:   if (fem->ops->destroy) {
126:     (*fem->ops->destroy)(fem);
127:     fem->ops->destroy = NULL;
128:   }
129:   (*r)(fem);
130:   PetscObjectChangeTypeName((PetscObject) fem, name);
131:   return 0;
132: }

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

137:   Not Collective

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

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

145:   Level: intermediate

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

160: /*@C
161:    PetscFEViewFromOptions - View from Options

163:    Collective on PetscFE

165:    Input Parameters:
166: +  A - the PetscFE object
167: .  obj - Optional object
168: -  name - command line option

170:    Level: intermediate
171: .seealso:  PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
172: @*/
173: PetscErrorCode  PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
174: {
176:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
177:   return 0;
178: }

180: /*@C
181:   PetscFEView - Views a PetscFE

183:   Collective on fem

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

189:   Level: beginner

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

199:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);
200:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
201:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
202:   if (fem->ops->view) (*fem->ops->view)(fem, viewer);
203:   return 0;
204: }

206: /*@
207:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

209:   Collective on fem

211:   Input Parameter:
212: . fem - the PetscFE object to set options for

214:   Options Database:
215: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
216: - -petscfe_num_batches - the number of cell batches to integrate serially

218:   Level: intermediate

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

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

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

256: /*@C
257:   PetscFESetUp - Construct data structures for the PetscFE

259:   Collective on fem

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

264:   Level: intermediate

266: .seealso PetscFEView(), PetscFEDestroy()
267: @*/
268: PetscErrorCode PetscFESetUp(PetscFE fem)
269: {
271:   if (fem->setupcalled) return 0;
272:   PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
273:   fem->setupcalled = PETSC_TRUE;
274:   if (fem->ops->setup) (*fem->ops->setup)(fem);
275:   PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
276:   return 0;
277: }

279: /*@
280:   PetscFEDestroy - Destroys a PetscFE object

282:   Collective on fem

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

287:   Level: beginner

289: .seealso PetscFEView()
290: @*/
291: PetscErrorCode PetscFEDestroy(PetscFE *fem)
292: {
293:   if (!*fem) return 0;

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

299:   if ((*fem)->subspaces) {
300:     PetscInt dim, d;

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

319:   if ((*fem)->ops->destroy) (*(*fem)->ops->destroy)(*fem);
320:   PetscHeaderDestroy(fem);
321:   return 0;
322: }

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

327:   Collective

329:   Input Parameter:
330: . comm - The communicator for the PetscFE object

332:   Output Parameter:
333: . fem - The PetscFE object

335:   Level: beginner

337: .seealso: PetscFESetType(), PETSCFEGALERKIN
338: @*/
339: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
340: {
341:   PetscFE        f;

344:   PetscCitationsRegister(FECitation,&FEcite);
345:   *fem = NULL;
346:   PetscFEInitializePackage();

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

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

365:   *fem = f;
366:   return 0;
367: }

369: /*@
370:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

372:   Not collective

374:   Input Parameter:
375: . fem - The PetscFE object

377:   Output Parameter:
378: . dim - The spatial dimension

380:   Level: intermediate

382: .seealso: PetscFECreate()
383: @*/
384: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
385: {
386:   DM             dm;

390:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
391:   DMGetDimension(dm, dim);
392:   return 0;
393: }

395: /*@
396:   PetscFESetNumComponents - Sets the number of components in the element

398:   Not collective

400:   Input Parameters:
401: + fem - The PetscFE object
402: - comp - The number of field components

404:   Level: intermediate

406: .seealso: PetscFECreate()
407: @*/
408: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
409: {
411:   fem->numComponents = comp;
412:   return 0;
413: }

415: /*@
416:   PetscFEGetNumComponents - Returns the number of components in the element

418:   Not collective

420:   Input Parameter:
421: . fem - The PetscFE object

423:   Output Parameter:
424: . comp - The number of field components

426:   Level: intermediate

428: .seealso: PetscFECreate()
429: @*/
430: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
431: {
434:   *comp = fem->numComponents;
435:   return 0;
436: }

438: /*@
439:   PetscFESetTileSizes - Sets the tile sizes for evaluation

441:   Not collective

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

450:   Level: intermediate

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

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

467:   Not collective

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

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

478:   Level: intermediate

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

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

499:   Not collective

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

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

507:   Level: intermediate

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

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

522:   Not collective

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

528:   Level: intermediate

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

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

545:   Not collective

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

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

553:   Level: intermediate

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

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

568:   Not collective

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

574:   Level: intermediate

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

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

591:   Not collective

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

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

599:   Level: intermediate

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

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

614:   Not collective

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

620:   Level: intermediate

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

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

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

644:   Not collective

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

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

652:   Level: intermediate

654: .seealso: PetscFECreate()
655: @*/
656: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
657: {
660:   *q = fem->faceQuadrature;
661:   return 0;
662: }

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

667:   Not collective

669:   Input Parameters:
670: + fem - The PetscFE object
671: - q - The PetscQuadrature object

673:   Level: intermediate

675: .seealso: PetscFECreate()
676: @*/
677: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
678: {
679:   PetscInt       Nc, qNc;

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

692: /*@
693:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

695:   Not collective

697:   Input Parameters:
698: + sfe - The PetscFE source for the quadratures
699: - tfe - The PetscFE target for the quadratures

701:   Level: intermediate

703: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
704: @*/
705: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
706: {
707:   PetscQuadrature q;

711:   PetscFEGetQuadrature(sfe, &q);
712:   PetscFESetQuadrature(tfe,  q);
713:   PetscFEGetFaceQuadrature(sfe, &q);
714:   PetscFESetFaceQuadrature(tfe,  q);
715:   return 0;
716: }

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

721:   Not collective

723:   Input Parameter:
724: . fem - The PetscFE object

726:   Output Parameter:
727: . numDof - Array with the number of dofs per dimension

729:   Level: intermediate

731: .seealso: PetscFECreate()
732: @*/
733: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
734: {
737:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
738:   return 0;
739: }

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

744:   Not collective

746:   Input Parameters:
747: + fem - The PetscFE object
748: - k   - The highest derivative we need to tabulate, very often 1

750:   Output Parameter:
751: . T - The basis function values and derivatives at quadrature points

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

758:   Level: intermediate

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

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

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

779:   Not collective

781:   Input Parameters:
782: + fem - The PetscFE object
783: - k   - The highest derivative we need to tabulate, very often 1

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

788:   Note:
789: $ 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
790: $ 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
791: $ 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

793:   Level: intermediate

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

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

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

837:   Not collective

839:   Input Parameter:
840: . fem - The PetscFE object

842:   Output Parameters:
843: . Tc - The basis function values at face centroid points

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

848:   Level: intermediate

850: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
851: @*/
852: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
853: {
856:   if (!fem->Tc) {
857:     PetscDualSpace  sp;
858:     DM              dm;
859:     const PetscInt *cone;
860:     PetscReal      *centroids;
861:     PetscInt        dim, numFaces, f;

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

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

880:   Not collective

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

889:   Output Parameter:
890: . T - The basis function values and derivatives at tabulation points

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

897:   Level: intermediate

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

910:   if (!npoints || !fem->dualSpace || K < 0) {
911:     *T = NULL;
912:     return 0;
913:   }
917:   PetscFEGetDualSpace(fem, &Q);
918:   PetscDualSpaceGetDM(Q, &dm);
919:   DMGetDimension(dm, &cdim);
920:   PetscDualSpaceGetDimension(Q, &Nb);
921:   PetscFEGetNumComponents(fem, &Nc);
922:   PetscMalloc1(1, T);
923:   (*T)->K    = !cdim ? 0 : K;
924:   (*T)->Nr   = nrepl;
925:   (*T)->Np   = npoints;
926:   (*T)->Nb   = Nb;
927:   (*T)->Nc   = Nc;
928:   (*T)->cdim = cdim;
929:   PetscMalloc1((*T)->K+1, &(*T)->T);
930:   for (k = 0; k <= (*T)->K; ++k) {
931:     PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
932:   }
933:   (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
934:   return 0;
935: }

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

940:   Not collective

942:   Input Parameters:
943: + fem     - The PetscFE object
944: . npoints - The number of tabulation points
945: . points  - The tabulation point coordinates
946: . K       - The number of derivatives calculated
947: - T       - An existing tabulation object with enough allocated space

949:   Output Parameter:
950: . T - The basis function values and derivatives at tabulation points

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

957:   Level: intermediate

959: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
960: @*/
961: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
962: {
964:   if (!npoints || !fem->dualSpace || K < 0) return 0;
968:   if (PetscDefined(USE_DEBUG)) {
969:     DM               dm;
970:     PetscDualSpace   Q;
971:     PetscInt         Nb;   /* Dimension of FE space P */
972:     PetscInt         Nc;   /* Field components */
973:     PetscInt         cdim; /* Reference coordinate dimension */

975:     PetscFEGetDualSpace(fem, &Q);
976:     PetscDualSpaceGetDM(Q, &dm);
977:     DMGetDimension(dm, &cdim);
978:     PetscDualSpaceGetDimension(Q, &Nb);
979:     PetscFEGetNumComponents(fem, &Nc);
984:   }
985:   T->Nr = 1;
986:   T->Np = npoints;
987:   (*fem->ops->createtabulation)(fem, npoints, points, K, T);
988:   return 0;
989: }

991: /*@C
992:   PetscTabulationDestroy - Frees memory from the associated tabulation.

994:   Not collective

996:   Input Parameter:
997: . T - The tabulation

999:   Level: intermediate

1001: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1002: @*/
1003: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1004: {
1005:   PetscInt       k;

1008:   if (!T || !(*T)) return 0;
1009:   for (k = 0; k <= (*T)->K; ++k) PetscFree((*T)->T[k]);
1010:   PetscFree((*T)->T);
1011:   PetscFree(*T);
1012:   *T = NULL;
1013:   return 0;
1014: }

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

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

1077: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1078: {
1079:   PetscInt       hStart, hEnd;
1080:   PetscDualSpace dsp;
1081:   DM             dm;

1085:   *trFE = NULL;
1086:   PetscFEGetDualSpace(fe,&dsp);
1087:   PetscDualSpaceGetDM(dsp,&dm);
1088:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1089:   if (hEnd <= hStart) return 0;
1090:   PetscFECreatePointTrace(fe,hStart,trFE);
1091:   return 0;
1092: }

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

1097:   Not collective

1099:   Input Parameter:
1100: . fe - The PetscFE

1102:   Output Parameter:
1103: . dim - The dimension

1105:   Level: intermediate

1107: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1108: @*/
1109: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1110: {
1113:   if (fem->ops->getdimension) (*fem->ops->getdimension)(fem, dim);
1114:   return 0;
1115: }

1117: /*@C
1118:   PetscFEPushforward - Map the reference element function to real space

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

1126:   Output Parameter:
1127: . vals   - The transformed function values

1129:   Level: advanced

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

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

1135: .seealso: PetscDualSpacePushforward()
1136: @*/
1137: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1138: {
1140:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1141:   return 0;
1142: }

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

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

1153:   Output Parameter:
1154: . vals   - The transformed function gradient values

1156:   Level: advanced

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

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

1162: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1163: @*/
1164: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1165: {
1167:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1168:   return 0;
1169: }

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

1174:   Input Parameters:
1175: + fe     - The PetscFE
1176: . fegeom - The cell geometry
1177: . Nv     - The number of function Hessian values
1178: - vals   - The function Hessian values

1180:   Output Parameter:
1181: . vals   - The transformed function Hessian values

1183:   Level: advanced

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

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

1189: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward()
1190: @*/
1191: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1192: {
1194:   PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1195:   return 0;
1196: }

1198: /*
1199: Purpose: Compute element vector for chunk of elements

1201: Input:
1202:   Sizes:
1203:      Ne:  number of elements
1204:      Nf:  number of fields
1205:      PetscFE
1206:        dim: spatial dimension
1207:        Nb:  number of basis functions
1208:        Nc:  number of field components
1209:        PetscQuadrature
1210:          Nq:  number of quadrature points

1212:   Geometry:
1213:      PetscFEGeom[Ne] possibly *Nq
1214:        PetscReal v0s[dim]
1215:        PetscReal n[dim]
1216:        PetscReal jacobians[dim*dim]
1217:        PetscReal jacobianInverses[dim*dim]
1218:        PetscReal jacobianDeterminants
1219:   FEM:
1220:      PetscFE
1221:        PetscQuadrature
1222:          PetscReal   quadPoints[Nq*dim]
1223:          PetscReal   quadWeights[Nq]
1224:        PetscReal   basis[Nq*Nb*Nc]
1225:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1226:      PetscScalar coefficients[Ne*Nb*Nc]
1227:      PetscScalar elemVec[Ne*Nb*Nc]

1229:   Problem:
1230:      PetscInt f: the active field
1231:      f0, f1

1233:   Work Space:
1234:      PetscFE
1235:        PetscScalar f0[Nq*dim];
1236:        PetscScalar f1[Nq*dim*dim];
1237:        PetscScalar u[Nc];
1238:        PetscScalar gradU[Nc*dim];
1239:        PetscReal   x[dim];
1240:        PetscScalar realSpaceDer[dim];

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

1244: Input:
1245:   Sizes:
1246:      N_cb: Number of serial cell batches

1248:   Geometry:
1249:      PetscReal v0s[Ne*dim]
1250:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1251:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1252:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1253:   FEM:
1254:      static PetscReal   quadPoints[Nq*dim]
1255:      static PetscReal   quadWeights[Nq]
1256:      static PetscReal   basis[Nq*Nb*Nc]
1257:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1258:      PetscScalar coefficients[Ne*Nb*Nc]
1259:      PetscScalar elemVec[Ne*Nb*Nc]

1261: ex62.c:
1262:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1263:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1264:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1265:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1267: ex52.c:
1268:   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)
1269:   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)

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

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

1278: ex52_integrateElementOpenCL.c:
1279: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1280:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1281:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1289:   Not collective

1291:   Input Parameters:
1292: + prob         - The PetscDS specifying the discretizations and continuum functions
1293: . field        - The field being integrated
1294: . Ne           - The number of elements in the chunk
1295: . cgeom        - The cell geometry for each cell in the chunk
1296: . coefficients - The array of FEM basis coefficients for the elements
1297: . probAux      - The PetscDS specifying the auxiliary discretizations
1298: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1300:   Output Parameter:
1301: . integral     - the integral for this field

1303:   Level: intermediate

1305: .seealso: PetscFEIntegrateResidual()
1306: @*/
1307: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1308:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1309: {
1310:   PetscFE        fe;

1313:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1314:   if (fe->ops->integrate) (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);
1315:   return 0;
1316: }

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

1321:   Not collective

1323:   Input Parameters:
1324: + prob         - The PetscDS specifying the discretizations and continuum functions
1325: . field        - The field being integrated
1326: . obj_func     - The function to be integrated
1327: . Ne           - The number of elements in the chunk
1328: . fgeom        - The face geometry for each face in the chunk
1329: . coefficients - The array of FEM basis coefficients for the elements
1330: . probAux      - The PetscDS specifying the auxiliary discretizations
1331: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1333:   Output Parameter:
1334: . integral     - the integral for this field

1336:   Level: intermediate

1338: .seealso: PetscFEIntegrateResidual()
1339: @*/
1340: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1341:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1342:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1343:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1344:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1345:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1346: {
1347:   PetscFE        fe;

1350:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1351:   if (fe->ops->integratebd) (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);
1352:   return 0;
1353: }

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

1358:   Not collective

1360:   Input Parameters:
1361: + ds           - The PetscDS specifying the discretizations and continuum functions
1362: . key          - The (label+value, field) being integrated
1363: . Ne           - The number of elements in the chunk
1364: . cgeom        - The cell geometry for each cell in the chunk
1365: . coefficients - The array of FEM basis coefficients for the elements
1366: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1367: . probAux      - The PetscDS specifying the auxiliary discretizations
1368: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1369: - t            - The time

1371:   Output Parameter:
1372: . elemVec      - the element residual vectors from each element

1374:   Note:
1375: $ Loop over batch of elements (e):
1376: $   Loop over quadrature points (q):
1377: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1378: $     Call f_0 and f_1
1379: $   Loop over element vector entries (f,fc --> i):
1380: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1382:   Level: intermediate

1384: .seealso: PetscFEIntegrateResidual()
1385: @*/
1386: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1387:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1388: {
1389:   PetscFE        fe;

1393:   PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);
1394:   if (fe->ops->integrateresidual) (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1395:   return 0;
1396: }

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

1401:   Not collective

1403:   Input Parameters:
1404: + ds           - The PetscDS specifying the discretizations and continuum functions
1405: . wf           - The PetscWeakForm object holding the pointwise functions
1406: . key          - The (label+value, field) being integrated
1407: . Ne           - The number of elements in the chunk
1408: . fgeom        - The face geometry for each cell in the chunk
1409: . coefficients - The array of FEM basis coefficients for the elements
1410: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1411: . probAux      - The PetscDS specifying the auxiliary discretizations
1412: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1413: - t            - The time

1415:   Output Parameter:
1416: . elemVec      - the element residual vectors from each element

1418:   Level: intermediate

1420: .seealso: PetscFEIntegrateResidual()
1421: @*/
1422: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1423:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1424: {
1425:   PetscFE        fe;

1428:   PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);
1429:   if (fe->ops->integratebdresidual) (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1430:   return 0;
1431: }

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

1436:   Not collective

1438:   Input Parameters:
1439: + prob         - The PetscDS specifying the discretizations and continuum functions
1440: . key          - The (label+value, field) being integrated
1441: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1442: . Ne           - The number of elements in the chunk
1443: . fgeom        - The face geometry for each cell in the chunk
1444: . coefficients - The array of FEM basis coefficients for the elements
1445: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1446: . probAux      - The PetscDS specifying the auxiliary discretizations
1447: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1448: - t            - The time

1450:   Output Parameter
1451: . elemVec      - the element residual vectors from each element

1453:   Level: developer

1455: .seealso: PetscFEIntegrateResidual()
1456: @*/
1457: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom,
1458:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1459: {
1460:   PetscFE        fe;

1463:   PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe);
1464:   if (fe->ops->integratehybridresidual) (*fe->ops->integratehybridresidual)(prob, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1465:   return 0;
1466: }

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

1471:   Not collective

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

1487:   Output Parameter:
1488: . elemMat      - the element matrices for the Jacobian from each element

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

1501: .seealso: PetscFEIntegrateResidual()
1502: @*/
1503: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1504:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1505: {
1506:   PetscFE        fe;
1507:   PetscInt       Nf;

1510:   PetscDSGetNumFields(ds, &Nf);
1511:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1512:   if (fe->ops->integratejacobian) (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1513:   return 0;
1514: }

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

1519:   Not collective

1521:   Input Parameters:
1522: + ds           - The PetscDS specifying the discretizations and continuum functions
1523: . wf           - The PetscWeakForm holding the pointwise functions
1524: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1525: . Ne           - The number of elements in the chunk
1526: . fgeom        - The face geometry for each cell in the chunk
1527: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1528: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1529: . probAux      - The PetscDS specifying the auxiliary discretizations
1530: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1531: . t            - The time
1532: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1534:   Output Parameter:
1535: . elemMat              - the element matrices for the Jacobian from each element

1537:   Note:
1538: $ Loop over batch of elements (e):
1539: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1540: $     Loop over quadrature points (q):
1541: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1542: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1543: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1544: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1545: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1546:   Level: intermediate

1548: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1549: @*/
1550: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1551:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1552: {
1553:   PetscFE        fe;
1554:   PetscInt       Nf;

1557:   PetscDSGetNumFields(ds, &Nf);
1558:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1559:   if (fe->ops->integratebdjacobian) (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1560:   return 0;
1561: }

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

1566:   Not collective

1568:   Input Parameters:
1569: + ds           - The PetscDS specifying the discretizations and continuum functions
1570: . jtype        - The type of matrix pointwise functions that should be used
1571: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1572: . s            - The side of the cell being integrated, 0 for negative and 1 for positive
1573: . Ne           - The number of elements in the chunk
1574: . fgeom        - The face geometry for each cell in the chunk
1575: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1576: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1577: . probAux      - The PetscDS specifying the auxiliary discretizations
1578: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1579: . t            - The time
1580: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1582:   Output Parameter
1583: . elemMat              - the element matrices for the Jacobian from each element

1585:   Note:
1586: $ Loop over batch of elements (e):
1587: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1588: $     Loop over quadrature points (q):
1589: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1590: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1591: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1592: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1593: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1594:   Level: developer

1596: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1597: @*/
1598: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom,
1599:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1600: {
1601:   PetscFE        fe;
1602:   PetscInt       Nf;

1605:   PetscDSGetNumFields(ds, &Nf);
1606:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1607:   if (fe->ops->integratehybridjacobian) (*fe->ops->integratehybridjacobian)(ds, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1608:   return 0;
1609: }

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

1614:   Input Parameters:
1615: + fe     - The finite element space
1616: - height - The height of the Plex point

1618:   Output Parameter:
1619: . subfe  - The subspace of this FE space

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

1623:   Level: advanced

1625: .seealso: PetscFECreateDefault()
1626: @*/
1627: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1628: {
1629:   PetscSpace      P, subP;
1630:   PetscDualSpace  Q, subQ;
1631:   PetscQuadrature subq;
1632:   PetscFEType     fetype;
1633:   PetscInt        dim, Nc;

1637:   if (height == 0) {
1638:     *subfe = fe;
1639:     return 0;
1640:   }
1641:   PetscFEGetBasisSpace(fe, &P);
1642:   PetscFEGetDualSpace(fe, &Q);
1643:   PetscFEGetNumComponents(fe, &Nc);
1644:   PetscFEGetFaceQuadrature(fe, &subq);
1645:   PetscDualSpaceGetDimension(Q, &dim);
1647:   if (!fe->subspaces) PetscCalloc1(dim, &fe->subspaces);
1648:   if (height <= dim) {
1649:     if (!fe->subspaces[height-1]) {
1650:       PetscFE     sub = NULL;
1651:       const char *name;

1653:       PetscSpaceGetHeightSubspace(P, height, &subP);
1654:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1655:       if (subQ) {
1656:         PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1657:         PetscObjectGetName((PetscObject) fe,  &name);
1658:         PetscObjectSetName((PetscObject) sub,  name);
1659:         PetscFEGetType(fe, &fetype);
1660:         PetscFESetType(sub, fetype);
1661:         PetscFESetBasisSpace(sub, subP);
1662:         PetscFESetDualSpace(sub, subQ);
1663:         PetscFESetNumComponents(sub, Nc);
1664:         PetscFESetUp(sub);
1665:         PetscFESetQuadrature(sub, subq);
1666:       }
1667:       fe->subspaces[height-1] = sub;
1668:     }
1669:     *subfe = fe->subspaces[height-1];
1670:   } else {
1671:     *subfe = NULL;
1672:   }
1673:   return 0;
1674: }

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

1681:   Collective on fem

1683:   Input Parameter:
1684: . fe - The initial PetscFE

1686:   Output Parameter:
1687: . feRef - The refined PetscFE

1689:   Level: advanced

1691: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1692: @*/
1693: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1694: {
1695:   PetscSpace       P, Pref;
1696:   PetscDualSpace   Q, Qref;
1697:   DM               K, Kref;
1698:   PetscQuadrature  q, qref;
1699:   const PetscReal *v0, *jac;
1700:   PetscInt         numComp, numSubelements;
1701:   PetscInt         cStart, cEnd, c;
1702:   PetscDualSpace  *cellSpaces;

1704:   PetscFEGetBasisSpace(fe, &P);
1705:   PetscFEGetDualSpace(fe, &Q);
1706:   PetscFEGetQuadrature(fe, &q);
1707:   PetscDualSpaceGetDM(Q, &K);
1708:   /* Create space */
1709:   PetscObjectReference((PetscObject) P);
1710:   Pref = P;
1711:   /* Create dual space */
1712:   PetscDualSpaceDuplicate(Q, &Qref);
1713:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1714:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1715:   PetscDualSpaceSetDM(Qref, Kref);
1716:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1717:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1718:   /* TODO: fix for non-uniform refinement */
1719:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1720:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1721:   PetscFree(cellSpaces);
1722:   DMDestroy(&Kref);
1723:   PetscDualSpaceSetUp(Qref);
1724:   /* Create element */
1725:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1726:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1727:   PetscFESetBasisSpace(*feRef, Pref);
1728:   PetscFESetDualSpace(*feRef, Qref);
1729:   PetscFEGetNumComponents(fe,    &numComp);
1730:   PetscFESetNumComponents(*feRef, numComp);
1731:   PetscFESetUp(*feRef);
1732:   PetscSpaceDestroy(&Pref);
1733:   PetscDualSpaceDestroy(&Qref);
1734:   /* Create quadrature */
1735:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1736:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1737:   PetscFESetQuadrature(*feRef, qref);
1738:   PetscQuadratureDestroy(&qref);
1739:   return 0;
1740: }

1742: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1743: {
1744:   PetscQuadrature q, fq;
1745:   DM              K;
1746:   PetscSpace      P;
1747:   PetscDualSpace  Q;
1748:   PetscInt        quadPointsPerEdge;
1749:   PetscBool       tensor;
1750:   char            name[64];
1751:   PetscErrorCode  ierr;

1755:   switch (ct) {
1756:     case DM_POLYTOPE_SEGMENT:
1757:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1758:     case DM_POLYTOPE_QUADRILATERAL:
1759:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1760:     case DM_POLYTOPE_HEXAHEDRON:
1761:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1762:       tensor = PETSC_TRUE;
1763:       break;
1764:     default: tensor = PETSC_FALSE;
1765:   }
1766:   /* Create space */
1767:   PetscSpaceCreate(comm, &P);
1768:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1769:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1770:   PetscSpacePolynomialSetTensor(P, tensor);
1771:   PetscSpaceSetNumComponents(P, Nc);
1772:   PetscSpaceSetNumVariables(P, dim);
1773:   if (degree >= 0) {
1774:     PetscSpaceSetDegree(P, degree, PETSC_DETERMINE);
1775:     if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1776:       PetscSpace Pend, Pside;

1778:       PetscSpaceCreate(comm, &Pend);
1779:       PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL);
1780:       PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE);
1781:       PetscSpaceSetNumComponents(Pend, Nc);
1782:       PetscSpaceSetNumVariables(Pend, dim-1);
1783:       PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE);
1784:       PetscSpaceCreate(comm, &Pside);
1785:       PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL);
1786:       PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE);
1787:       PetscSpaceSetNumComponents(Pside, 1);
1788:       PetscSpaceSetNumVariables(Pside, 1);
1789:       PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE);
1790:       PetscSpaceSetType(P, PETSCSPACETENSOR);
1791:       PetscSpaceTensorSetNumSubspaces(P, 2);
1792:       PetscSpaceTensorSetSubspace(P, 0, Pend);
1793:       PetscSpaceTensorSetSubspace(P, 1, Pside);
1794:       PetscSpaceDestroy(&Pend);
1795:       PetscSpaceDestroy(&Pside);
1796:     }
1797:   }
1798:   if (setFromOptions) PetscSpaceSetFromOptions(P);
1799:   PetscSpaceSetUp(P);
1800:   PetscSpaceGetDegree(P, &degree, NULL);
1801:   PetscSpacePolynomialGetTensor(P, &tensor);
1802:   PetscSpaceGetNumComponents(P, &Nc);
1803:   /* Create dual space */
1804:   PetscDualSpaceCreate(comm, &Q);
1805:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1806:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1807:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1808:   PetscDualSpaceSetDM(Q, K);
1809:   DMDestroy(&K);
1810:   PetscDualSpaceSetNumComponents(Q, Nc);
1811:   PetscDualSpaceSetOrder(Q, degree);
1812:   /* TODO For some reason, we need a tensor dualspace with wedges */
1813:   PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE);
1814:   if (setFromOptions) PetscDualSpaceSetFromOptions(Q);
1815:   PetscDualSpaceSetUp(Q);
1816:   /* Create finite element */
1817:   PetscFECreate(comm, fem);
1818:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1819:   PetscFESetType(*fem, PETSCFEBASIC);
1820:   PetscFESetBasisSpace(*fem, P);
1821:   PetscFESetDualSpace(*fem, Q);
1822:   PetscFESetNumComponents(*fem, Nc);
1823:   if (setFromOptions) PetscFESetFromOptions(*fem);
1824:   PetscFESetUp(*fem);
1825:   PetscSpaceDestroy(&P);
1826:   PetscDualSpaceDestroy(&Q);
1827:   /* Create quadrature (with specified order if given) */
1828:   qorder = qorder >= 0 ? qorder : degree;
1829:   if (setFromOptions) {
1830:     PetscObjectOptionsBegin((PetscObject)*fem);
1831:     PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1832:     PetscOptionsEnd();
1833:   }
1834:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1835:   switch (ct) {
1836:     case DM_POLYTOPE_SEGMENT:
1837:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1838:     case DM_POLYTOPE_QUADRILATERAL:
1839:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1840:     case DM_POLYTOPE_HEXAHEDRON:
1841:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1842:       PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1843:       PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1844:       break;
1845:     case DM_POLYTOPE_TRIANGLE:
1846:     case DM_POLYTOPE_TETRAHEDRON:
1847:       PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1848:       PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1849:       break;
1850:     case DM_POLYTOPE_TRI_PRISM:
1851:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1852:       {
1853:         PetscQuadrature q1, q2;

1855:         PetscDTStroudConicalQuadrature(2, 1, quadPointsPerEdge, -1.0, 1.0, &q1);
1856:         PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2);
1857:         PetscDTTensorQuadratureCreate(q1, q2, &q);
1858:         PetscQuadratureDestroy(&q1);
1859:         PetscQuadratureDestroy(&q2);
1860:       }
1861:       PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1862:       /* TODO Need separate quadratures for each face */
1863:       break;
1864:     default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
1865:   }
1866:   PetscFESetQuadrature(*fem, q);
1867:   PetscFESetFaceQuadrature(*fem, fq);
1868:   PetscQuadratureDestroy(&q);
1869:   PetscQuadratureDestroy(&fq);
1870:   /* Set finite element name */
1871:   switch (ct) {
1872:     case DM_POLYTOPE_SEGMENT:
1873:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1874:     case DM_POLYTOPE_QUADRILATERAL:
1875:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1876:     case DM_POLYTOPE_HEXAHEDRON:
1877:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1878:       PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree);
1879:       break;
1880:     case DM_POLYTOPE_TRIANGLE:
1881:     case DM_POLYTOPE_TETRAHEDRON:
1882:       PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree);
1883:       break;
1884:     case DM_POLYTOPE_TRI_PRISM:
1885:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1886:       PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree);
1887:       break;
1888:     default:
1889:       PetscSNPrintf(name, sizeof(name), "FE");
1890:   }
1891:   PetscFESetName(*fem, name);
1892:   return 0;
1893: }

1895: /*@C
1896:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1898:   Collective

1900:   Input Parameters:
1901: + comm      - The MPI comm
1902: . dim       - The spatial dimension
1903: . Nc        - The number of components
1904: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1905: . prefix    - The options prefix, or NULL
1906: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1908:   Output Parameter:
1909: . fem - The PetscFE object

1911:   Note:
1912:   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.

1914:   Level: beginner

1916: .seealso: PetscFECreateLagrange(), PetscFECreateByCell(), PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1917: @*/
1918: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1919: {
1920:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1921:   return 0;
1922: }

1924: /*@C
1925:   PetscFECreateByCell - Create a PetscFE for basic FEM computation

1927:   Collective

1929:   Input Parameters:
1930: + comm   - The MPI comm
1931: . dim    - The spatial dimension
1932: . Nc     - The number of components
1933: . ct     - The celltype of the reference cell
1934: . prefix - The options prefix, or NULL
1935: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1937:   Output Parameter:
1938: . fem - The PetscFE object

1940:   Note:
1941:   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.

1943:   Level: beginner

1945: .seealso: PetscFECreateDefault(), PetscFECreateLagrange(), PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1946: @*/
1947: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
1948: {
1949:   PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1950:   return 0;
1951: }

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

1956:   Collective

1958:   Input Parameters:
1959: + comm      - The MPI comm
1960: . dim       - The spatial dimension
1961: . Nc        - The number of components
1962: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1963: . k         - The degree k of the space
1964: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1966:   Output Parameter:
1967: . fem       - The PetscFE object

1969:   Level: beginner

1971:   Notes:
1972:   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.

1974: .seealso: PetscFECreateLagrangeByCell(), PetscFECreateDefault(), PetscFECreateByCell(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1975: @*/
1976: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1977: {
1978:   PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem);
1979:   return 0;
1980: }

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

1985:   Collective

1987:   Input Parameters:
1988: + comm      - The MPI comm
1989: . dim       - The spatial dimension
1990: . Nc        - The number of components
1991: . ct        - The celltype of the reference cell
1992: . k         - The degree k of the space
1993: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1995:   Output Parameter:
1996: . fem       - The PetscFE object

1998:   Level: beginner

2000:   Notes:
2001:   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.

2003: .seealso: PetscFECreateLagrange(), PetscFECreateDefault(), PetscFECreateByCell(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
2004: @*/
2005: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2006: {
2007:   PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem);
2008:   return 0;
2009: }

2011: /*@C
2012:   PetscFESetName - Names the FE and its subobjects

2014:   Not collective

2016:   Input Parameters:
2017: + fe   - The PetscFE
2018: - name - The name

2020:   Level: intermediate

2022: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
2023: @*/
2024: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2025: {
2026:   PetscSpace     P;
2027:   PetscDualSpace Q;

2029:   PetscFEGetBasisSpace(fe, &P);
2030:   PetscFEGetDualSpace(fe, &Q);
2031:   PetscObjectSetName((PetscObject) fe, name);
2032:   PetscObjectSetName((PetscObject) P,  name);
2033:   PetscObjectSetName((PetscObject) Q,  name);
2034:   return 0;
2035: }

2037: 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[])
2038: {
2039:   PetscInt       dOffset = 0, fOffset = 0, f, g;

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

2053:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2054:     for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2055:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2056:     for (b = 0; b < Nbf; ++b) {
2057:       for (c = 0; c < Ncf; ++c) {
2058:         const PetscInt cidx = b*Ncf+c;

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

2071:           for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b];
2072:         }
2073:       }
2074:       PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);
2075:     }
2076:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2077:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2078:     if (u_t) {
2079:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2080:       for (b = 0; b < Nbf; ++b) {
2081:         for (c = 0; c < Ncf; ++c) {
2082:           const PetscInt cidx = b*Ncf+c;

2084:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2085:         }
2086:       }
2087:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2088:     }
2089:     fOffset += Ncf;
2090:     dOffset += Nbf;
2091:   }
2092:   return 0;
2093: }

2095: 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[])
2096: {
2097:   PetscInt       dOffset = 0, fOffset = 0, f, g;

2099:   /* f is the field number in the DS, g is the field number in u[] */
2100:   for (f = 0, g = 0; f < Nf; ++f) {
2101:     PetscFE          fe   = (PetscFE) ds->disc[f];
2102:     const PetscInt   dEt  = T[f]->cdim;
2103:     const PetscInt   dE   = fegeom->dimEmbed;
2104:     const PetscInt   Nq   = T[f]->Np;
2105:     const PetscInt   Nbf  = T[f]->Nb;
2106:     const PetscInt   Ncf  = T[f]->Nc;
2107:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2108:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*dEt];
2109:     PetscBool        isCohesive;
2110:     PetscInt         Ns, s;

2112:     if (!T[f]) continue;
2113:     PetscDSGetCohesive(ds, f, &isCohesive);
2114:     Ns   = isCohesive ? 1 : 2;
2115:     for (s = 0; s < Ns; ++s, ++g) {
2116:       PetscInt b, c, d;

2118:       for (c = 0; c < Ncf; ++c)    u[fOffset+c]      = 0.0;
2119:       for (d = 0; d < dE*Ncf; ++d) u_x[fOffset*dE+d] = 0.0;
2120:       for (b = 0; b < Nbf; ++b) {
2121:         for (c = 0; c < Ncf; ++c) {
2122:           const PetscInt cidx = b*Ncf+c;

2124:           u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2125:           for (d = 0; d < dEt; ++d) u_x[(fOffset+c)*dE+d] += Dq[cidx*dEt+d]*coefficients[dOffset+b];
2126:         }
2127:       }
2128:       PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2129:       PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*dE]);
2130:       if (u_t) {
2131:         for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2132:         for (b = 0; b < Nbf; ++b) {
2133:           for (c = 0; c < Ncf; ++c) {
2134:             const PetscInt cidx = b*Ncf+c;

2136:             u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2137:           }
2138:         }
2139:         PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2140:       }
2141:       fOffset += Ncf;
2142:       dOffset += Nbf;
2143:     }
2144:   }
2145:   return 0;
2146: }

2148: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2149: {
2150:   PetscFE         fe;
2151:   PetscTabulation Tc;
2152:   PetscInt        b, c;

2154:   if (!prob) return 0;
2155:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2156:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2157:   {
2158:     const PetscReal *faceBasis = Tc->T[0];
2159:     const PetscInt   Nb        = Tc->Nb;
2160:     const PetscInt   Nc        = Tc->Nc;

2162:     for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2163:     for (b = 0; b < Nb; ++b) {
2164:       for (c = 0; c < Nc; ++c) {
2165:         u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2166:       }
2167:     }
2168:   }
2169:   return 0;
2170: }

2172: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2173: {
2174:   PetscFEGeom      pgeom;
2175:   const PetscInt   dEt      = T->cdim;
2176:   const PetscInt   dE       = fegeom->dimEmbed;
2177:   const PetscInt   Nq       = T->Np;
2178:   const PetscInt   Nb       = T->Nb;
2179:   const PetscInt   Nc       = T->Nc;
2180:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2181:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dEt];
2182:   PetscInt         q, b, c, d;

2184:   for (q = 0; q < Nq; ++q) {
2185:     for (b = 0; b < Nb; ++b) {
2186:       for (c = 0; c < Nc; ++c) {
2187:         const PetscInt bcidx = b*Nc+c;

2189:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2190:         for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dEt+bcidx*dEt+d];
2191:         for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = 0.0;
2192:       }
2193:     }
2194:     PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);
2195:     PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);
2196:     PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer);
2197:     for (b = 0; b < Nb; ++b) {
2198:       for (c = 0; c < Nc; ++c) {
2199:         const PetscInt bcidx = b*Nc+c;
2200:         const PetscInt qcidx = q*Nc+c;

2202:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2203:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2204:       }
2205:     }
2206:   }
2207:   return(0);
2208: }

2210: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2211: {
2212:   const PetscInt   dE       = T->cdim;
2213:   const PetscInt   Nq       = T->Np;
2214:   const PetscInt   Nb       = T->Nb;
2215:   const PetscInt   Nc       = T->Nc;
2216:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2217:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2218:   PetscInt         q, b, c, d;

2220:   for (q = 0; q < Nq; ++q) {
2221:     for (b = 0; b < Nb; ++b) {
2222:       for (c = 0; c < Nc; ++c) {
2223:         const PetscInt bcidx = b*Nc+c;

2225:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2226:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2227:       }
2228:     }
2229:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2230:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2231:     for (b = 0; b < Nb; ++b) {
2232:       for (c = 0; c < Nc; ++c) {
2233:         const PetscInt bcidx = b*Nc+c;
2234:         const PetscInt qcidx = q*Nc+c;

2236:         elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2237:         for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2238:       }
2239:     }
2240:   }
2241:   return(0);
2242: }

2244: 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[])
2245: {
2246:   const PetscInt   dE        = TI->cdim;
2247:   const PetscInt   NqI       = TI->Np;
2248:   const PetscInt   NbI       = TI->Nb;
2249:   const PetscInt   NcI       = TI->Nc;
2250:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2251:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2252:   const PetscInt   NqJ       = TJ->Np;
2253:   const PetscInt   NbJ       = TJ->Nb;
2254:   const PetscInt   NcJ       = TJ->Nc;
2255:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2256:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2257:   PetscInt         f, fc, g, gc, df, dg;

2259:   for (f = 0; f < NbI; ++f) {
2260:     for (fc = 0; fc < NcI; ++fc) {
2261:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

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

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

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

2304: 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[])
2305: {
2306:   const PetscInt   dE        = TI->cdim;
2307:   const PetscInt   NqI       = TI->Np;
2308:   const PetscInt   NbI       = TI->Nb;
2309:   const PetscInt   NcI       = TI->Nc;
2310:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2311:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2312:   const PetscInt   NqJ       = TJ->Np;
2313:   const PetscInt   NbJ       = TJ->Nb;
2314:   const PetscInt   NcJ       = TJ->Nc;
2315:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2316:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2317:   const PetscInt   so        = isHybridI ? 0 : s;
2318:   const PetscInt   to        = isHybridJ ? 0 : s;
2319:   PetscInt         f, fc, g, gc, df, dg;

2321:   for (f = 0; f < NbI; ++f) {
2322:     for (fc = 0; fc < NcI; ++fc) {
2323:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

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

2335:       tmpBasisJ[gidx] = basisJ[gidx];
2336:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2337:     }
2338:   }
2339:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2340:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2341:   for (f = 0; f < NbI; ++f) {
2342:     for (fc = 0; fc < NcI; ++fc) {
2343:       const PetscInt fidx = f*NcI+fc;         /* Test function basis index */
2344:       const PetscInt i    = offsetI+NbI*so+f; /* Element matrix row */
2345:       for (g = 0; g < NbJ; ++g) {
2346:         for (gc = 0; gc < NcJ; ++gc) {
2347:           const PetscInt gidx = g*NcJ+gc;         /* Trial function basis index */
2348:           const PetscInt j    = offsetJ+NbJ*to+g; /* Element matrix column */
2349:           const PetscInt fOff = eOffset+i*totDim+j;

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

2366: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2367: {
2368:   PetscDualSpace  dsp;
2369:   DM              dm;
2370:   PetscQuadrature quadDef;
2371:   PetscInt        dim, cdim, Nq;

2373:   PetscFEGetDualSpace(fe, &dsp);
2374:   PetscDualSpaceGetDM(dsp, &dm);
2375:   DMGetDimension(dm, &dim);
2376:   DMGetCoordinateDim(dm, &cdim);
2377:   PetscFEGetQuadrature(fe, &quadDef);
2378:   quad = quad ? quad : quadDef;
2379:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2380:   PetscMalloc1(Nq*cdim,      &cgeom->v);
2381:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2382:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2383:   PetscMalloc1(Nq,           &cgeom->detJ);
2384:   cgeom->dim       = dim;
2385:   cgeom->dimEmbed  = cdim;
2386:   cgeom->numCells  = 1;
2387:   cgeom->numPoints = Nq;
2388:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2389:   return 0;
2390: }

2392: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2393: {
2394:   PetscFree(cgeom->v);
2395:   PetscFree(cgeom->J);
2396:   PetscFree(cgeom->invJ);
2397:   PetscFree(cgeom->detJ);
2398:   return 0;
2399: }