Actual source code: fe.c

petsc-master 2020-04-04
Report Typos and Errors
  1: /* Basis Jet Tabulation

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

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

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

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

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

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

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

 22:    \alpha = V^{-1}

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

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

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

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

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

 51: PetscClassId PETSCFE_CLASSID = 0;

 53: PetscFunctionList PetscFEList              = NULL;
 54: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

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

 59:   Not Collective

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

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

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

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

 83:   Level: advanced

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

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

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

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

100:   Collective on fem

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

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

109:   Level: intermediate

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

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

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

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

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

140:   Not Collective

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

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

148:   Level: intermediate

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

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

166: /*@C
167:    PetscFEViewFromOptions - View from Options

169:    Collective on PetscFE

171:    Input Parameters:
172: +  A - the PetscFE object
173: .  obj - Optional object
174: -  name - command line option

176:    Level: intermediate
177: .seealso:  PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
178: @*/
179: PetscErrorCode  PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
180: {

185:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
186:   return(0);
187: }

189: /*@C
190:   PetscFEView - Views a PetscFE

192:   Collective on fem

194:   Input Parameter:
195: + fem - the PetscFE object to view
196: - viewer   - the viewer

198:   Level: beginner

200: .seealso PetscFEDestroy()
201: @*/
202: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
203: {
204:   PetscBool      iascii;

210:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
211:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
212:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
213:   if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
214:   return(0);
215: }

217: /*@
218:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

220:   Collective on fem

222:   Input Parameter:
223: . fem - the PetscFE object to set options for

225:   Options Database:
226: + -petscfe_num_blocks  - the number of cell blocks to integrate concurrently
227: - -petscfe_num_batches - the number of cell batches to integrate serially

229:   Level: intermediate

231: .seealso PetscFEView()
232: @*/
233: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
234: {
235:   const char    *defaultType;
236:   char           name[256];
237:   PetscBool      flg;

242:   if (!((PetscObject) fem)->type_name) {
243:     defaultType = PETSCFEBASIC;
244:   } else {
245:     defaultType = ((PetscObject) fem)->type_name;
246:   }
247:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

249:   PetscObjectOptionsBegin((PetscObject) fem);
250:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
251:   if (flg) {
252:     PetscFESetType(fem, name);
253:   } else if (!((PetscObject) fem)->type_name) {
254:     PetscFESetType(fem, defaultType);
255:   }
256:   PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
257:   PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
258:   if (fem->ops->setfromoptions) {
259:     (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
260:   }
261:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
262:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
263:   PetscOptionsEnd();
264:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
265:   return(0);
266: }

268: /*@C
269:   PetscFESetUp - Construct data structures for the PetscFE

271:   Collective on fem

273:   Input Parameter:
274: . fem - the PetscFE object to setup

276:   Level: intermediate

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

286:   if (fem->setupcalled) return(0);
287:   fem->setupcalled = PETSC_TRUE;
288:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
289:   return(0);
290: }

292: /*@
293:   PetscFEDestroy - Destroys a PetscFE object

295:   Collective on fem

297:   Input Parameter:
298: . fem - the PetscFE object to destroy

300:   Level: beginner

302: .seealso PetscFEView()
303: @*/
304: PetscErrorCode PetscFEDestroy(PetscFE *fem)
305: {

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

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

315:   if ((*fem)->subspaces) {
316:     PetscInt dim, d;

318:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
319:     for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
320:   }
321:   PetscFree((*fem)->subspaces);
322:   PetscFree((*fem)->invV);
323:   PetscTabulationDestroy(&(*fem)->T);
324:   PetscTabulationDestroy(&(*fem)->Tf);
325:   PetscTabulationDestroy(&(*fem)->Tc);
326:   PetscSpaceDestroy(&(*fem)->basisSpace);
327:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
328:   PetscQuadratureDestroy(&(*fem)->quadrature);
329:   PetscQuadratureDestroy(&(*fem)->faceQuadrature);

331:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
332:   PetscHeaderDestroy(fem);
333:   return(0);
334: }

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

339:   Collective

341:   Input Parameter:
342: . comm - The communicator for the PetscFE object

344:   Output Parameter:
345: . fem - The PetscFE object

347:   Level: beginner

349: .seealso: PetscFESetType(), PETSCFEGALERKIN
350: @*/
351: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
352: {
353:   PetscFE        f;

358:   PetscCitationsRegister(FECitation,&FEcite);
359:   *fem = NULL;
360:   PetscFEInitializePackage();

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

364:   f->basisSpace    = NULL;
365:   f->dualSpace     = NULL;
366:   f->numComponents = 1;
367:   f->subspaces     = NULL;
368:   f->invV          = NULL;
369:   f->T             = NULL;
370:   f->Tf            = NULL;
371:   f->Tc            = NULL;
372:   PetscArrayzero(&f->quadrature, 1);
373:   PetscArrayzero(&f->faceQuadrature, 1);
374:   f->blockSize     = 0;
375:   f->numBlocks     = 1;
376:   f->batchSize     = 0;
377:   f->numBatches    = 1;

379:   *fem = f;
380:   return(0);
381: }

383: /*@
384:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

386:   Not collective

388:   Input Parameter:
389: . fem - The PetscFE object

391:   Output Parameter:
392: . dim - The spatial dimension

394:   Level: intermediate

396: .seealso: PetscFECreate()
397: @*/
398: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
399: {
400:   DM             dm;

406:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
407:   DMGetDimension(dm, dim);
408:   return(0);
409: }

411: /*@
412:   PetscFESetNumComponents - Sets the number of components in the element

414:   Not collective

416:   Input Parameters:
417: + fem - The PetscFE object
418: - comp - The number of field components

420:   Level: intermediate

422: .seealso: PetscFECreate()
423: @*/
424: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
425: {
428:   fem->numComponents = comp;
429:   return(0);
430: }

432: /*@
433:   PetscFEGetNumComponents - Returns the number of components in the element

435:   Not collective

437:   Input Parameter:
438: . fem - The PetscFE object

440:   Output Parameter:
441: . comp - The number of field components

443:   Level: intermediate

445: .seealso: PetscFECreate()
446: @*/
447: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
448: {
452:   *comp = fem->numComponents;
453:   return(0);
454: }

456: /*@
457:   PetscFESetTileSizes - Sets the tile sizes for evaluation

459:   Not collective

461:   Input Parameters:
462: + fem - The PetscFE object
463: . blockSize - The number of elements in a block
464: . numBlocks - The number of blocks in a batch
465: . batchSize - The number of elements in a batch
466: - numBatches - The number of batches in a chunk

468:   Level: intermediate

470: .seealso: PetscFECreate()
471: @*/
472: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
473: {
476:   fem->blockSize  = blockSize;
477:   fem->numBlocks  = numBlocks;
478:   fem->batchSize  = batchSize;
479:   fem->numBatches = numBatches;
480:   return(0);
481: }

483: /*@
484:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

486:   Not collective

488:   Input Parameter:
489: . fem - The PetscFE object

491:   Output Parameters:
492: + blockSize - The number of elements in a block
493: . numBlocks - The number of blocks in a batch
494: . batchSize - The number of elements in a batch
495: - numBatches - The number of batches in a chunk

497:   Level: intermediate

499: .seealso: PetscFECreate()
500: @*/
501: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
502: {
509:   if (blockSize)  *blockSize  = fem->blockSize;
510:   if (numBlocks)  *numBlocks  = fem->numBlocks;
511:   if (batchSize)  *batchSize  = fem->batchSize;
512:   if (numBatches) *numBatches = fem->numBatches;
513:   return(0);
514: }

516: /*@
517:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

519:   Not collective

521:   Input Parameter:
522: . fem - The PetscFE object

524:   Output Parameter:
525: . sp - The PetscSpace object

527:   Level: intermediate

529: .seealso: PetscFECreate()
530: @*/
531: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
532: {
536:   *sp = fem->basisSpace;
537:   return(0);
538: }

540: /*@
541:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

543:   Not collective

545:   Input Parameters:
546: + fem - The PetscFE object
547: - sp - The PetscSpace object

549:   Level: intermediate

551: .seealso: PetscFECreate()
552: @*/
553: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
554: {

560:   PetscSpaceDestroy(&fem->basisSpace);
561:   fem->basisSpace = sp;
562:   PetscObjectReference((PetscObject) fem->basisSpace);
563:   return(0);
564: }

566: /*@
567:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

569:   Not collective

571:   Input Parameter:
572: . fem - The PetscFE object

574:   Output Parameter:
575: . sp - The PetscDualSpace object

577:   Level: intermediate

579: .seealso: PetscFECreate()
580: @*/
581: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
582: {
586:   *sp = fem->dualSpace;
587:   return(0);
588: }

590: /*@
591:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

593:   Not collective

595:   Input Parameters:
596: + fem - The PetscFE object
597: - sp - The PetscDualSpace object

599:   Level: intermediate

601: .seealso: PetscFECreate()
602: @*/
603: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
604: {

610:   PetscDualSpaceDestroy(&fem->dualSpace);
611:   fem->dualSpace = sp;
612:   PetscObjectReference((PetscObject) fem->dualSpace);
613:   return(0);
614: }

616: /*@
617:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

619:   Not collective

621:   Input Parameter:
622: . fem - The PetscFE object

624:   Output Parameter:
625: . q - The PetscQuadrature object

627:   Level: intermediate

629: .seealso: PetscFECreate()
630: @*/
631: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
632: {
636:   *q = fem->quadrature;
637:   return(0);
638: }

640: /*@
641:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

643:   Not collective

645:   Input Parameters:
646: + fem - The PetscFE object
647: - q - The PetscQuadrature object

649:   Level: intermediate

651: .seealso: PetscFECreate()
652: @*/
653: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
654: {
655:   PetscInt       Nc, qNc;

660:   PetscFEGetNumComponents(fem, &Nc);
661:   PetscQuadratureGetNumComponents(q, &qNc);
662:   if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
663:   PetscTabulationDestroy(&fem->T);
664:   PetscTabulationDestroy(&fem->Tc);
665:   PetscQuadratureDestroy(&fem->quadrature);
666:   fem->quadrature = q;
667:   PetscObjectReference((PetscObject) q);
668:   return(0);
669: }

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

674:   Not collective

676:   Input Parameter:
677: . fem - The PetscFE object

679:   Output Parameter:
680: . q - The PetscQuadrature object

682:   Level: intermediate

684: .seealso: PetscFECreate()
685: @*/
686: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
687: {
691:   *q = fem->faceQuadrature;
692:   return(0);
693: }

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

698:   Not collective

700:   Input Parameters:
701: + fem - The PetscFE object
702: - q - The PetscQuadrature object

704:   Level: intermediate

706: .seealso: PetscFECreate()
707: @*/
708: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
709: {
710:   PetscInt       Nc, qNc;

715:   PetscFEGetNumComponents(fem, &Nc);
716:   PetscQuadratureGetNumComponents(q, &qNc);
717:   if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
718:   PetscTabulationDestroy(&fem->Tf);
719:   PetscQuadratureDestroy(&fem->faceQuadrature);
720:   fem->faceQuadrature = q;
721:   PetscObjectReference((PetscObject) q);
722:   return(0);
723: }

725: /*@
726:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

728:   Not collective

730:   Input Parameters:
731: + sfe - The PetscFE source for the quadratures
732: - tfe - The PetscFE target for the quadratures

734:   Level: intermediate

736: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
737: @*/
738: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
739: {
740:   PetscQuadrature q;
741:   PetscErrorCode  ierr;

746:   PetscFEGetQuadrature(sfe, &q);
747:   PetscFESetQuadrature(tfe,  q);
748:   PetscFEGetFaceQuadrature(sfe, &q);
749:   PetscFESetFaceQuadrature(tfe,  q);
750:   return(0);
751: }

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

756:   Not collective

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

761:   Output Parameter:
762: . numDof - Array with the number of dofs per dimension

764:   Level: intermediate

766: .seealso: PetscFECreate()
767: @*/
768: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
769: {

775:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
776:   return(0);
777: }

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

782:   Not collective

784:   Input Parameter:
785: . fem - The PetscFE object

787:   Output Parameter:
788: . T - The basis function values and derivatives at quadrature points

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

795:   Level: intermediate

797: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
798: @*/
799: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T)
800: {
801:   PetscInt         npoints;
802:   const PetscReal *points;
803:   PetscErrorCode   ierr;

808:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
809:   if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);}
810:   *T = fem->T;
811:   return(0);
812: }

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

817:   Not collective

819:   Input Parameter:
820: . fem - The PetscFE object

822:   Output Parameters:
823: . Tf - The basis function values and derviatives at face quadrature points

825:   Note:
826: $ 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
827: $ 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
828: $ 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

830:   Level: intermediate

832: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
833: @*/
834: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf)
835: {
836:   PetscErrorCode   ierr;

841:   if (!fem->Tf) {
842:     const PetscReal  xi0[3] = {-1., -1., -1.};
843:     PetscReal        v0[3], J[9], detJ;
844:     PetscQuadrature  fq;
845:     PetscDualSpace   sp;
846:     DM               dm;
847:     const PetscInt  *faces;
848:     PetscInt         dim, numFaces, f, npoints, q;
849:     const PetscReal *points;
850:     PetscReal       *facePoints;

852:     PetscFEGetDualSpace(fem, &sp);
853:     PetscDualSpaceGetDM(sp, &dm);
854:     DMGetDimension(dm, &dim);
855:     DMPlexGetConeSize(dm, 0, &numFaces);
856:     DMPlexGetCone(dm, 0, &faces);
857:     PetscFEGetFaceQuadrature(fem, &fq);
858:     if (fq) {
859:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
860:       PetscMalloc1(numFaces*npoints*dim, &facePoints);
861:       for (f = 0; f < numFaces; ++f) {
862:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
863:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
864:       }
865:       PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);
866:       PetscFree(facePoints);
867:     }
868:   }
869:   *Tf = fem->Tf;
870:   return(0);
871: }

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

876:   Not collective

878:   Input Parameter:
879: . fem - The PetscFE object

881:   Output Parameters:
882: . Tc - The basis function values at face centroid points

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

887:   Level: intermediate

889: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
890: @*/
891: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
892: {
893:   PetscErrorCode   ierr;

898:   if (!fem->Tc) {
899:     PetscDualSpace  sp;
900:     DM              dm;
901:     const PetscInt *cone;
902:     PetscReal      *centroids;
903:     PetscInt        dim, numFaces, f;

905:     PetscFEGetDualSpace(fem, &sp);
906:     PetscDualSpaceGetDM(sp, &dm);
907:     DMGetDimension(dm, &dim);
908:     DMPlexGetConeSize(dm, 0, &numFaces);
909:     DMPlexGetCone(dm, 0, &cone);
910:     PetscMalloc1(numFaces*dim, &centroids);
911:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
912:     PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
913:     PetscFree(centroids);
914:   }
915:   *Tc = fem->Tc;
916:   return(0);
917: }

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

922:   Not collective

924:   Input Parameters:
925: + fem     - The PetscFE object
926: . nrepl   - The number of replicas
927: . npoints - The number of tabulation points in a replica
928: . points  - The tabulation point coordinates
929: - K       - The number of derivatives calculated

931:   Output Parameter:
932: . T - The basis function values and derivatives at tabulation points

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

939:   Level: intermediate

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

954:   if (!npoints || !fem->dualSpace || K < 0) {
955:     *T = NULL;
956:     return(0);
957:   }
961:   PetscFEGetDualSpace(fem, &Q);
962:   PetscDualSpaceGetDM(Q, &dm);
963:   DMGetDimension(dm, &cdim);
964:   PetscDualSpaceGetDimension(Q, &Nb);
965:   PetscFEGetNumComponents(fem, &Nc);
966:   PetscMalloc1(1, T);
967:   (*T)->K    = !cdim ? 0 : K;
968:   (*T)->Nr   = nrepl;
969:   (*T)->Np   = npoints;
970:   (*T)->Nb   = Nb;
971:   (*T)->Nc   = Nc;
972:   (*T)->cdim = cdim;
973:   PetscMalloc1((*T)->K+1, &(*T)->T);
974:   for (k = 0; k <= (*T)->K; ++k) {
975:     PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
976:   }
977:   (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
978:   return(0);
979: }

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

984:   Not collective

986:   Input Parameters:
987: + fem     - The PetscFE object
988: . npoints - The number of tabulation points
989: . points  - The tabulation point coordinates
990: . K       - The number of derivatives calculated
991: - T       - An existing tabulation object with enough allocated space

993:   Output Parameter:
994: . T - The basis function values and derivatives at tabulation points

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

1001:   Level: intermediate

1003: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1004: @*/
1005: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1006: {

1010:   if (!npoints || !fem->dualSpace || K < 0) return(0);
1014: #ifdef PETSC_USE_DEBUG
1015:   {
1016:     DM               dm;
1017:     PetscDualSpace   Q;
1018:     PetscInt         Nb;   /* Dimension of FE space P */
1019:     PetscInt         Nc;   /* Field components */
1020:     PetscInt         cdim; /* Reference coordinate dimension */

1022:     PetscFEGetDualSpace(fem, &Q);
1023:     PetscDualSpaceGetDM(Q, &dm);
1024:     DMGetDimension(dm, &cdim);
1025:     PetscDualSpaceGetDimension(Q, &Nb);
1026:     PetscFEGetNumComponents(fem, &Nc);
1027:     if (T->K    != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K);
1028:     if (T->Nb   != Nb)              SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1029:     if (T->Nc   != Nc)              SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1030:     if (T->cdim != cdim)            SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1031:   }
1032: #endif
1033:   T->Nr = 1;
1034:   T->Np = npoints;
1035:   (*fem->ops->createtabulation)(fem, npoints, points, K, T);
1036:   return(0);
1037: }

1039: /*@C
1040:   PetscTabulationDestroy - Frees memory from the associated tabulation.

1042:   Not collective

1044:   Input Parameter:
1045: . T - The tabulation

1047:   Level: intermediate

1049: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1050: @*/
1051: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1052: {
1053:   PetscInt       k;

1058:   if (!T || !(*T)) return(0);
1059:   for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);}
1060:   PetscFree((*T)->T);
1061:   PetscFree(*T);
1062:   *T = NULL;
1063:   return(0);
1064: }

1066: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1067: {
1068:   PetscSpace     bsp, bsubsp;
1069:   PetscDualSpace dsp, dsubsp;
1070:   PetscInt       dim, depth, numComp, i, j, coneSize, order;
1071:   PetscFEType    type;
1072:   DM             dm;
1073:   DMLabel        label;
1074:   PetscReal      *xi, *v, *J, detJ;
1075:   const char     *name;
1076:   PetscQuadrature origin, fullQuad, subQuad;

1082:   PetscFEGetBasisSpace(fe,&bsp);
1083:   PetscFEGetDualSpace(fe,&dsp);
1084:   PetscDualSpaceGetDM(dsp,&dm);
1085:   DMGetDimension(dm,&dim);
1086:   DMPlexGetDepthLabel(dm,&label);
1087:   DMLabelGetValue(label,refPoint,&depth);
1088:   PetscCalloc1(depth,&xi);
1089:   PetscMalloc1(dim,&v);
1090:   PetscMalloc1(dim*dim,&J);
1091:   for (i = 0; i < depth; i++) xi[i] = 0.;
1092:   PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1093:   PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1094:   DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1095:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1096:   for (i = 1; i < dim; i++) {
1097:     for (j = 0; j < depth; j++) {
1098:       J[i * depth + j] = J[i * dim + j];
1099:     }
1100:   }
1101:   PetscQuadratureDestroy(&origin);
1102:   PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1103:   PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1104:   PetscSpaceSetUp(bsubsp);
1105:   PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1106:   PetscFEGetType(fe,&type);
1107:   PetscFESetType(*trFE,type);
1108:   PetscFEGetNumComponents(fe,&numComp);
1109:   PetscFESetNumComponents(*trFE,numComp);
1110:   PetscFESetBasisSpace(*trFE,bsubsp);
1111:   PetscFESetDualSpace(*trFE,dsubsp);
1112:   PetscObjectGetName((PetscObject) fe, &name);
1113:   if (name) {PetscFESetName(*trFE, name);}
1114:   PetscFEGetQuadrature(fe,&fullQuad);
1115:   PetscQuadratureGetOrder(fullQuad,&order);
1116:   DMPlexGetConeSize(dm,refPoint,&coneSize);
1117:   if (coneSize == 2 * depth) {
1118:     PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1119:   } else {
1120:     PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1121:   }
1122:   PetscFESetQuadrature(*trFE,subQuad);
1123:   PetscFESetUp(*trFE);
1124:   PetscQuadratureDestroy(&subQuad);
1125:   PetscSpaceDestroy(&bsubsp);
1126:   return(0);
1127: }

1129: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1130: {
1131:   PetscInt       hStart, hEnd;
1132:   PetscDualSpace dsp;
1133:   DM             dm;

1139:   *trFE = NULL;
1140:   PetscFEGetDualSpace(fe,&dsp);
1141:   PetscDualSpaceGetDM(dsp,&dm);
1142:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1143:   if (hEnd <= hStart) return(0);
1144:   PetscFECreatePointTrace(fe,hStart,trFE);
1145:   return(0);
1146: }


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

1152:   Not collective

1154:   Input Parameter:
1155: . fe - The PetscFE

1157:   Output Parameter:
1158: . dim - The dimension

1160:   Level: intermediate

1162: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1163: @*/
1164: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1165: {

1171:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1172:   return(0);
1173: }

1175: /*@C
1176:   PetscFEPushforward - Map the reference element function to real space

1178:   Input Parameters:
1179: + fe     - The PetscFE
1180: . fegeom - The cell geometry
1181: . Nv     - The number of function values
1182: - vals   - The function values

1184:   Output Parameter:
1185: . vals   - The transformed function values

1187:   Level: advanced

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

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

1193: .seealso: PetscDualSpacePushforward()
1194: @*/
1195: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1196: {

1200:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1201:   return(0);
1202: }

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

1207:   Input Parameters:
1208: + fe     - The PetscFE
1209: . fegeom - The cell geometry
1210: . Nv     - The number of function gradient values
1211: - vals   - The function gradient values

1213:   Output Parameter:
1214: . vals   - The transformed function gradient values

1216:   Level: advanced

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

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

1222: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1223: @*/
1224: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1225: {

1229:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1230:   return(0);
1231: }

1233: /*
1234: Purpose: Compute element vector for chunk of elements

1236: Input:
1237:   Sizes:
1238:      Ne:  number of elements
1239:      Nf:  number of fields
1240:      PetscFE
1241:        dim: spatial dimension
1242:        Nb:  number of basis functions
1243:        Nc:  number of field components
1244:        PetscQuadrature
1245:          Nq:  number of quadrature points

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

1264:   Problem:
1265:      PetscInt f: the active field
1266:      f0, f1

1268:   Work Space:
1269:      PetscFE
1270:        PetscScalar f0[Nq*dim];
1271:        PetscScalar f1[Nq*dim*dim];
1272:        PetscScalar u[Nc];
1273:        PetscScalar gradU[Nc*dim];
1274:        PetscReal   x[dim];
1275:        PetscScalar realSpaceDer[dim];

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

1279: Input:
1280:   Sizes:
1281:      N_cb: Number of serial cell batches

1283:   Geometry:
1284:      PetscReal v0s[Ne*dim]
1285:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1286:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1287:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1288:   FEM:
1289:      static PetscReal   quadPoints[Nq*dim]
1290:      static PetscReal   quadWeights[Nq]
1291:      static PetscReal   basis[Nq*Nb*Nc]
1292:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1293:      PetscScalar coefficients[Ne*Nb*Nc]
1294:      PetscScalar elemVec[Ne*Nb*Nc]

1296: ex62.c:
1297:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1298:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1299:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1300:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1302: ex52.c:
1303:   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)
1304:   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)

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

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

1313: ex52_integrateElementOpenCL.c:
1314: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1315:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1316:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1324:   Not collective

1326:   Input Parameters:
1327: + fem          - The PetscFE object for the field being integrated
1328: . prob         - The PetscDS specifying the discretizations and continuum functions
1329: . field        - The field being integrated
1330: . Ne           - The number of elements in the chunk
1331: . cgeom        - The cell geometry for each cell in the chunk
1332: . coefficients - The array of FEM basis coefficients for the elements
1333: . probAux      - The PetscDS specifying the auxiliary discretizations
1334: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1336:   Output Parameter:
1337: . integral     - the integral for this field

1339:   Level: intermediate

1341: .seealso: PetscFEIntegrateResidual()
1342: @*/
1343: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1344:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1345: {
1346:   PetscFE        fe;

1351:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1352:   if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1353:   return(0);
1354: }

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

1359:   Not collective

1361:   Input Parameters:
1362: + fem          - The PetscFE object for the field being integrated
1363: . prob         - The PetscDS specifying the discretizations and continuum functions
1364: . field        - The field being integrated
1365: . obj_func     - The function to be integrated
1366: . Ne           - The number of elements in the chunk
1367: . fgeom        - The face geometry for each face in the chunk
1368: . coefficients - The array of FEM basis coefficients for the elements
1369: . probAux      - The PetscDS specifying the auxiliary discretizations
1370: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1372:   Output Parameter:
1373: . integral     - the integral for this field

1375:   Level: intermediate

1377: .seealso: PetscFEIntegrateResidual()
1378: @*/
1379: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1380:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1381:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1382:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1383:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1384:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1385: {
1386:   PetscFE        fe;

1391:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1392:   if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1393:   return(0);
1394: }

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

1399:   Not collective

1401:   Input Parameters:
1402: + fem          - The PetscFE object for the field being integrated
1403: . prob         - The PetscDS specifying the discretizations and continuum functions
1404: . field        - The field being integrated
1405: . Ne           - The number of elements in the chunk
1406: . cgeom        - The cell geometry for each cell in the chunk
1407: . coefficients - The array of FEM basis coefficients for the elements
1408: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1409: . probAux      - The PetscDS specifying the auxiliary discretizations
1410: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1411: - t            - The time

1413:   Output Parameter:
1414: . elemVec      - the element residual vectors from each element

1416:   Note:
1417: $ Loop over batch of elements (e):
1418: $   Loop over quadrature points (q):
1419: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1420: $     Call f_0 and f_1
1421: $   Loop over element vector entries (f,fc --> i):
1422: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1424:   Level: intermediate

1426: .seealso: PetscFEIntegrateResidual()
1427: @*/
1428: PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1429:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1430: {
1431:   PetscFE        fe;

1436:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1437:   if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1438:   return(0);
1439: }

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

1444:   Not collective

1446:   Input Parameters:
1447: + fem          - The PetscFE object for the field being integrated
1448: . prob         - The PetscDS specifying the discretizations and continuum functions
1449: . field        - The field being integrated
1450: . Ne           - The number of elements in the chunk
1451: . fgeom        - The face geometry for each cell in the chunk
1452: . coefficients - The array of FEM basis coefficients for the elements
1453: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1454: . probAux      - The PetscDS specifying the auxiliary discretizations
1455: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1456: - t            - The time

1458:   Output Parameter:
1459: . elemVec      - the element residual vectors from each element

1461:   Level: intermediate

1463: .seealso: PetscFEIntegrateResidual()
1464: @*/
1465: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1466:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1467: {
1468:   PetscFE        fe;

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

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

1481:   Not collective

1483:   Input Parameters:
1484: + fem          - The PetscFE object for the field being integrated
1485: . prob         - The PetscDS specifying the discretizations and continuum functions
1486: . jtype        - The type of matrix pointwise functions that should be used
1487: . fieldI       - The test field being integrated
1488: . fieldJ       - The basis field being integrated
1489: . Ne           - The number of elements in the chunk
1490: . cgeom        - The cell geometry for each cell in the chunk
1491: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1492: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1493: . probAux      - The PetscDS specifying the auxiliary discretizations
1494: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1495: . t            - The time
1496: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

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

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

1512: .seealso: PetscFEIntegrateResidual()
1513: @*/
1514: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1515:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1516: {
1517:   PetscFE        fe;

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

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

1530:   Not collective

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

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

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

1559: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1560: @*/
1561: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1562:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1563: {
1564:   PetscFE        fe;

1569:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1570:   if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1571:   return(0);
1572: }

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

1577:   Input Parameters:
1578: + fe     - The finite element space
1579: - height - The height of the Plex point

1581:   Output Parameter:
1582: . subfe  - The subspace of this FE space

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

1586:   Level: advanced

1588: .seealso: PetscFECreateDefault()
1589: @*/
1590: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1591: {
1592:   PetscSpace      P, subP;
1593:   PetscDualSpace  Q, subQ;
1594:   PetscQuadrature subq;
1595:   PetscFEType     fetype;
1596:   PetscInt        dim, Nc;
1597:   PetscErrorCode  ierr;

1602:   if (height == 0) {
1603:     *subfe = fe;
1604:     return(0);
1605:   }
1606:   PetscFEGetBasisSpace(fe, &P);
1607:   PetscFEGetDualSpace(fe, &Q);
1608:   PetscFEGetNumComponents(fe, &Nc);
1609:   PetscFEGetFaceQuadrature(fe, &subq);
1610:   PetscDualSpaceGetDimension(Q, &dim);
1611:   if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);}
1612:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1613:   if (height <= dim) {
1614:     if (!fe->subspaces[height-1]) {
1615:       PetscFE     sub;
1616:       const char *name;

1618:       PetscSpaceGetHeightSubspace(P, height, &subP);
1619:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1620:       PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1621:       PetscObjectGetName((PetscObject) fe,  &name);
1622:       PetscObjectSetName((PetscObject) sub,  name);
1623:       PetscFEGetType(fe, &fetype);
1624:       PetscFESetType(sub, fetype);
1625:       PetscFESetBasisSpace(sub, subP);
1626:       PetscFESetDualSpace(sub, subQ);
1627:       PetscFESetNumComponents(sub, Nc);
1628:       PetscFESetUp(sub);
1629:       PetscFESetQuadrature(sub, subq);
1630:       fe->subspaces[height-1] = sub;
1631:     }
1632:     *subfe = fe->subspaces[height-1];
1633:   } else {
1634:     *subfe = NULL;
1635:   }
1636:   return(0);
1637: }

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

1644:   Collective on fem

1646:   Input Parameter:
1647: . fe - The initial PetscFE

1649:   Output Parameter:
1650: . feRef - The refined PetscFE

1652:   Level: advanced

1654: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1655: @*/
1656: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1657: {
1658:   PetscSpace       P, Pref;
1659:   PetscDualSpace   Q, Qref;
1660:   DM               K, Kref;
1661:   PetscQuadrature  q, qref;
1662:   const PetscReal *v0, *jac;
1663:   PetscInt         numComp, numSubelements;
1664:   PetscInt         cStart, cEnd, c;
1665:   PetscDualSpace  *cellSpaces;
1666:   PetscErrorCode   ierr;

1669:   PetscFEGetBasisSpace(fe, &P);
1670:   PetscFEGetDualSpace(fe, &Q);
1671:   PetscFEGetQuadrature(fe, &q);
1672:   PetscDualSpaceGetDM(Q, &K);
1673:   /* Create space */
1674:   PetscObjectReference((PetscObject) P);
1675:   Pref = P;
1676:   /* Create dual space */
1677:   PetscDualSpaceDuplicate(Q, &Qref);
1678:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1679:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1680:   PetscDualSpaceSetDM(Qref, Kref);
1681:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1682:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1683:   /* TODO: fix for non-uniform refinement */
1684:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1685:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1686:   PetscFree(cellSpaces);
1687:   DMDestroy(&Kref);
1688:   PetscDualSpaceSetUp(Qref);
1689:   /* Create element */
1690:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1691:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1692:   PetscFESetBasisSpace(*feRef, Pref);
1693:   PetscFESetDualSpace(*feRef, Qref);
1694:   PetscFEGetNumComponents(fe,    &numComp);
1695:   PetscFESetNumComponents(*feRef, numComp);
1696:   PetscFESetUp(*feRef);
1697:   PetscSpaceDestroy(&Pref);
1698:   PetscDualSpaceDestroy(&Qref);
1699:   /* Create quadrature */
1700:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1701:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1702:   PetscFESetQuadrature(*feRef, qref);
1703:   PetscQuadratureDestroy(&qref);
1704:   return(0);
1705: }

1707: /*@C
1708:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1710:   Collective

1712:   Input Parameters:
1713: + comm      - The MPI comm
1714: . dim       - The spatial dimension
1715: . Nc        - The number of components
1716: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1717: . prefix    - The options prefix, or NULL
1718: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1720:   Output Parameter:
1721: . fem - The PetscFE object

1723:   Note:
1724:   Each object is SetFromOption() during creation, so that the object may be customized from the command line.

1726:   Level: beginner

1728: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1729: @*/
1730: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1731: {
1732:   PetscQuadrature q, fq;
1733:   DM              K;
1734:   PetscSpace      P;
1735:   PetscDualSpace  Q;
1736:   PetscInt        order, quadPointsPerEdge;
1737:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1738:   PetscErrorCode  ierr;

1741:   /* Create space */
1742:   PetscSpaceCreate(comm, &P);
1743:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1744:   PetscSpacePolynomialSetTensor(P, tensor);
1745:   PetscSpaceSetNumComponents(P, Nc);
1746:   PetscSpaceSetNumVariables(P, dim);
1747:   PetscSpaceSetFromOptions(P);
1748:   PetscSpaceSetUp(P);
1749:   PetscSpaceGetDegree(P, &order, NULL);
1750:   PetscSpacePolynomialGetTensor(P, &tensor);
1751:   /* Create dual space */
1752:   PetscDualSpaceCreate(comm, &Q);
1753:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1754:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1755:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1756:   PetscDualSpaceSetDM(Q, K);
1757:   DMDestroy(&K);
1758:   PetscDualSpaceSetNumComponents(Q, Nc);
1759:   PetscDualSpaceSetOrder(Q, order);
1760:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1761:   PetscDualSpaceSetFromOptions(Q);
1762:   PetscDualSpaceSetUp(Q);
1763:   /* Create element */
1764:   PetscFECreate(comm, fem);
1765:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1766:   PetscFESetBasisSpace(*fem, P);
1767:   PetscFESetDualSpace(*fem, Q);
1768:   PetscFESetNumComponents(*fem, Nc);
1769:   PetscFESetFromOptions(*fem);
1770:   PetscFESetUp(*fem);
1771:   PetscSpaceDestroy(&P);
1772:   PetscDualSpaceDestroy(&Q);
1773:   /* Create quadrature (with specified order if given) */
1774:   qorder = qorder >= 0 ? qorder : order;
1775:   PetscObjectOptionsBegin((PetscObject)*fem);
1776:   PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1777:   PetscOptionsEnd();
1778:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1779:   if (isSimplex) {
1780:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1781:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1782:   } else {
1783:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1784:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1785:   }
1786:   PetscFESetQuadrature(*fem, q);
1787:   PetscFESetFaceQuadrature(*fem, fq);
1788:   PetscQuadratureDestroy(&q);
1789:   PetscQuadratureDestroy(&fq);
1790:   return(0);
1791: }

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

1796:   Collective

1798:   Input Parameters:
1799: + comm      - The MPI comm
1800: . dim       - The spatial dimension
1801: . Nc        - The number of components
1802: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1803: . k         - The degree k of the space
1804: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

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

1809:   Level: beginner

1811:   Notes:
1812:   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.

1814: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1815: @*/
1816: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1817: {
1818:   PetscQuadrature q, fq;
1819:   DM              K;
1820:   PetscSpace      P;
1821:   PetscDualSpace  Q;
1822:   PetscInt        quadPointsPerEdge;
1823:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1824:   char            name[64];
1825:   PetscErrorCode  ierr;

1828:   /* Create space */
1829:   PetscSpaceCreate(comm, &P);
1830:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1831:   PetscSpacePolynomialSetTensor(P, tensor);
1832:   PetscSpaceSetNumComponents(P, Nc);
1833:   PetscSpaceSetNumVariables(P, dim);
1834:   PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1835:   PetscSpaceSetUp(P);
1836:   /* Create dual space */
1837:   PetscDualSpaceCreate(comm, &Q);
1838:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1839:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1840:   PetscDualSpaceSetDM(Q, K);
1841:   DMDestroy(&K);
1842:   PetscDualSpaceSetNumComponents(Q, Nc);
1843:   PetscDualSpaceSetOrder(Q, k);
1844:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1845:   PetscDualSpaceSetUp(Q);
1846:   /* Create element */
1847:   PetscFECreate(comm, fem);
1848:   PetscSNPrintf(name, 64, "P%d", (int) k);
1849:   PetscObjectSetName((PetscObject) *fem, name);
1850:   PetscFESetType(*fem, PETSCFEBASIC);
1851:   PetscFESetBasisSpace(*fem, P);
1852:   PetscFESetDualSpace(*fem, Q);
1853:   PetscFESetNumComponents(*fem, Nc);
1854:   PetscFESetUp(*fem);
1855:   PetscSpaceDestroy(&P);
1856:   PetscDualSpaceDestroy(&Q);
1857:   /* Create quadrature (with specified order if given) */
1858:   qorder = qorder >= 0 ? qorder : k;
1859:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1860:   if (isSimplex) {
1861:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1862:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1863:   } else {
1864:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1865:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1866:   }
1867:   PetscFESetQuadrature(*fem, q);
1868:   PetscFESetFaceQuadrature(*fem, fq);
1869:   PetscQuadratureDestroy(&q);
1870:   PetscQuadratureDestroy(&fq);
1871:   return(0);
1872: }

1874: /*@C
1875:   PetscFESetName - Names the FE and its subobjects

1877:   Not collective

1879:   Input Parameters:
1880: + fe   - The PetscFE
1881: - name - The name

1883:   Level: intermediate

1885: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1886: @*/
1887: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1888: {
1889:   PetscSpace     P;
1890:   PetscDualSpace Q;

1894:   PetscFEGetBasisSpace(fe, &P);
1895:   PetscFEGetDualSpace(fe, &Q);
1896:   PetscObjectSetName((PetscObject) fe, name);
1897:   PetscObjectSetName((PetscObject) P,  name);
1898:   PetscObjectSetName((PetscObject) Q,  name);
1899:   return(0);
1900: }

1902: 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[])
1903: {
1904:   PetscInt       dOffset = 0, fOffset = 0, f;

1907:   for (f = 0; f < Nf; ++f) {
1908:     PetscFE          fe;
1909:     const PetscInt   cdim = T[f]->cdim;
1910:     const PetscInt   Nq   = T[f]->Np;
1911:     const PetscInt   Nbf  = T[f]->Nb;
1912:     const PetscInt   Ncf  = T[f]->Nc;
1913:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
1914:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
1915:     PetscInt         b, c, d;

1917:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
1918:     for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
1919:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
1920:     for (b = 0; b < Nbf; ++b) {
1921:       for (c = 0; c < Ncf; ++c) {
1922:         const PetscInt cidx = b*Ncf+c;

1924:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
1925:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
1926:       }
1927:     }
1928:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
1929:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
1930:     if (u_t) {
1931:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
1932:       for (b = 0; b < Nbf; ++b) {
1933:         for (c = 0; c < Ncf; ++c) {
1934:           const PetscInt cidx = b*Ncf+c;

1936:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
1937:         }
1938:       }
1939:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
1940:     }
1941:     fOffset += Ncf;
1942:     dOffset += Nbf;
1943:   }
1944:   return 0;
1945: }

1947: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
1948: {
1949:   PetscFE         fe;
1950:   PetscTabulation Tc;
1951:   PetscInt        b, c;
1952:   PetscErrorCode  ierr;

1954:   if (!prob) return 0;
1955:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1956:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
1957:   {
1958:     const PetscReal *faceBasis = Tc->T[0];
1959:     const PetscInt   Nb        = Tc->Nb;
1960:     const PetscInt   Nc        = Tc->Nc;

1962:     for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
1963:     for (b = 0; b < Nb; ++b) {
1964:       for (c = 0; c < Nc; ++c) {
1965:         const PetscInt cidx = b*Nc+c;

1967:         u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
1968:       }
1969:     }
1970:   }
1971:   return 0;
1972: }

1974: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
1975: {
1976:   const PetscInt   dim      = T->cdim;
1977:   const PetscInt   Nq       = T->Np;
1978:   const PetscInt   Nb       = T->Nb;
1979:   const PetscInt   Nc       = T->Nc;
1980:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
1981:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dim];
1982:   PetscInt         q, b, c, d;
1983:   PetscErrorCode   ierr;

1985:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
1986:   for (q = 0; q < Nq; ++q) {
1987:     for (b = 0; b < Nb; ++b) {
1988:       for (c = 0; c < Nc; ++c) {
1989:         const PetscInt bcidx = b*Nc+c;

1991:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
1992:         for (d = 0; d < dim; ++d) tmpBasisDer[bcidx*dim+d] = basisDer[q*Nb*Nc*dim+bcidx*dim+d];
1993:       }
1994:     }
1995:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
1996:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
1997:     for (b = 0; b < Nb; ++b) {
1998:       for (c = 0; c < Nc; ++c) {
1999:         const PetscInt bcidx = b*Nc+c;
2000:         const PetscInt qcidx = q*Nc+c;

2002:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2003:         for (d = 0; d < dim; ++d) elemVec[b] += tmpBasisDer[bcidx*dim+d]*f1[qcidx*dim+d];
2004:       }
2005:     }
2006:   }
2007:   return(0);
2008: }

2010: 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[])
2011: {
2012:   const PetscInt   dim       = TI->cdim;
2013:   const PetscInt   NqI       = TI->Np;
2014:   const PetscInt   NbI       = TI->Nb;
2015:   const PetscInt   NcI       = TI->Nc;
2016:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2017:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dim];
2018:   const PetscInt   NqJ       = TJ->Np;
2019:   const PetscInt   NbJ       = TJ->Nb;
2020:   const PetscInt   NcJ       = TJ->Nc;
2021:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2022:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dim];
2023:   PetscInt         f, fc, g, gc, df, dg;
2024:   PetscErrorCode   ierr;

2026:   for (f = 0; f < NbI; ++f) {
2027:     for (fc = 0; fc < NcI; ++fc) {
2028:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2030:       tmpBasisI[fidx] = basisI[fidx];
2031:       for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dim+df] = basisDerI[fidx*dim+df];
2032:     }
2033:   }
2034:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2035:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2036:   for (g = 0; g < NbJ; ++g) {
2037:     for (gc = 0; gc < NcJ; ++gc) {
2038:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2040:       tmpBasisJ[gidx] = basisJ[gidx];
2041:       for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dim+dg] = basisDerJ[gidx*dim+dg];
2042:     }
2043:   }
2044:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2045:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2046:   for (f = 0; f < NbI; ++f) {
2047:     for (fc = 0; fc < NcI; ++fc) {
2048:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2049:       const PetscInt i    = offsetI+f; /* Element matrix row */
2050:       for (g = 0; g < NbJ; ++g) {
2051:         for (gc = 0; gc < NcJ; ++gc) {
2052:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2053:           const PetscInt j    = offsetJ+g; /* Element matrix column */
2054:           const PetscInt fOff = eOffset+i*totDim+j;

2056:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2057:           for (df = 0; df < dim; ++df) {
2058:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dim+df]*tmpBasisDerJ[gidx*dim+df];
2059:             elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g2[(fc*NcJ+gc)*dim+df]*tmpBasisJ[gidx];
2060:             for (dg = 0; dg < dim; ++dg) {
2061:               elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g3[((fc*NcJ+gc)*dim+df)*dim+dg]*tmpBasisDerJ[gidx*dim+dg];
2062:             }
2063:           }
2064:         }
2065:       }
2066:     }
2067:   }
2068:   return(0);
2069: }

2071: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2072: {
2073:   PetscDualSpace  dsp;
2074:   DM              dm;
2075:   PetscQuadrature quadDef;
2076:   PetscInt        dim, cdim, Nq;
2077:   PetscErrorCode  ierr;

2080:   PetscFEGetDualSpace(fe, &dsp);
2081:   PetscDualSpaceGetDM(dsp, &dm);
2082:   DMGetDimension(dm, &dim);
2083:   DMGetCoordinateDim(dm, &cdim);
2084:   PetscFEGetQuadrature(fe, &quadDef);
2085:   quad = quad ? quad : quadDef;
2086:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2087:   PetscMalloc1(Nq*cdim,      &cgeom->v);
2088:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2089:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2090:   PetscMalloc1(Nq,           &cgeom->detJ);
2091:   cgeom->dim       = dim;
2092:   cgeom->dimEmbed  = cdim;
2093:   cgeom->numCells  = 1;
2094:   cgeom->numPoints = Nq;
2095:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2096:   return(0);
2097: }

2099: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2100: {

2104:   PetscFree(cgeom->v);
2105:   PetscFree(cgeom->J);
2106:   PetscFree(cgeom->invJ);
2107:   PetscFree(cgeom->detJ);
2108:   return(0);
2109: }