Actual source code: fe.c

petsc-3.11.3 2019-06-26
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: .keywords: PetscFE, register
 86: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

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

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

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

101:   Collective on PetscFE

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

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

110:   Level: intermediate

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

123:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
124:   if (match) return(0);

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

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

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

142:   Not Collective

144:   Input Parameter:
145: . fem  - The PetscFE

147:   Output Parameter:
148: . name - The PetscFE type name

150:   Level: intermediate

152: .keywords: PetscFE, get, type, name
153: .seealso: PetscFESetType(), PetscFECreate()
154: @*/
155: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
156: {

162:   if (!PetscFERegisterAllCalled) {
163:     PetscFERegisterAll();
164:   }
165:   *name = ((PetscObject) fem)->type_name;
166:   return(0);
167: }

169: /*@C
170:   PetscFEView - Views a PetscFE

172:   Collective on PetscFE

174:   Input Parameter:
175: + fem - the PetscFE object to view
176: - viewer   - the viewer

178:   Level: developer

180: .seealso PetscFEDestroy()
181: @*/
182: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
183: {
184:   PetscBool      iascii;

190:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
191:   PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
192:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
193:   if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
194:   return(0);
195: }

197: /*@
198:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

200:   Collective on PetscFE

202:   Input Parameter:
203: . fem - the PetscFE object to set options for

205:   Options Database:
206: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
207: . -petscfe_num_batches the number of cell batches to integrate serially

209:   Level: developer

211: .seealso PetscFEView()
212: @*/
213: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
214: {
215:   const char    *defaultType;
216:   char           name[256];
217:   PetscBool      flg;

222:   if (!((PetscObject) fem)->type_name) {
223:     defaultType = PETSCFEBASIC;
224:   } else {
225:     defaultType = ((PetscObject) fem)->type_name;
226:   }
227:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

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

248: /*@C
249:   PetscFESetUp - Construct data structures for the PetscFE

251:   Collective on PetscFE

253:   Input Parameter:
254: . fem - the PetscFE object to setup

256:   Level: developer

258: .seealso PetscFEView(), PetscFEDestroy()
259: @*/
260: PetscErrorCode PetscFESetUp(PetscFE fem)
261: {

266:   if (fem->setupcalled) return(0);
267:   fem->setupcalled = PETSC_TRUE;
268:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
269:   return(0);
270: }

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

275:   Collective on PetscFE

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

280:   Level: developer

282: .seealso PetscFEView()
283: @*/
284: PetscErrorCode PetscFEDestroy(PetscFE *fem)
285: {

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

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

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

298:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
299:     for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
300:   }
301:   PetscFree((*fem)->subspaces);
302:   PetscFree((*fem)->invV);
303:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
304:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
305:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
306:   PetscSpaceDestroy(&(*fem)->basisSpace);
307:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
308:   PetscQuadratureDestroy(&(*fem)->quadrature);
309:   PetscQuadratureDestroy(&(*fem)->faceQuadrature);

311:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
312:   PetscHeaderDestroy(fem);
313:   return(0);
314: }

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

319:   Collective on MPI_Comm

321:   Input Parameter:
322: . comm - The communicator for the PetscFE object

324:   Output Parameter:
325: . fem - The PetscFE object

327:   Level: beginner

329: .seealso: PetscFESetType(), PETSCFEGALERKIN
330: @*/
331: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
332: {
333:   PetscFE        f;

338:   PetscCitationsRegister(FECitation,&FEcite);
339:   *fem = NULL;
340:   PetscFEInitializePackage();

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

344:   f->basisSpace    = NULL;
345:   f->dualSpace     = NULL;
346:   f->numComponents = 1;
347:   f->subspaces     = NULL;
348:   f->invV          = NULL;
349:   f->B             = NULL;
350:   f->D             = NULL;
351:   f->H             = NULL;
352:   f->Bf            = NULL;
353:   f->Df            = NULL;
354:   f->Hf            = NULL;
355:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
356:   PetscMemzero(&f->faceQuadrature, sizeof(PetscQuadrature));
357:   f->blockSize     = 0;
358:   f->numBlocks     = 1;
359:   f->batchSize     = 0;
360:   f->numBatches    = 1;

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

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

369:   Not collective

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

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

377:   Level: intermediate

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

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

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

397:   Not collective

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

403:   Level: intermediate

405: .seealso: PetscFECreate()
406: @*/
407: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
408: {
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: {
435:   *comp = fem->numComponents;
436:   return(0);
437: }

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

442:   Not collective

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

451:   Level: intermediate

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

466: /*@
467:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

469:   Not collective

471:   Input Parameter:
472: . fem - The PetscFE object

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

480:   Level: intermediate

482: .seealso: PetscFECreate()
483: @*/
484: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
485: {
492:   if (blockSize)  *blockSize  = fem->blockSize;
493:   if (numBlocks)  *numBlocks  = fem->numBlocks;
494:   if (batchSize)  *batchSize  = fem->batchSize;
495:   if (numBatches) *numBatches = fem->numBatches;
496:   return(0);
497: }

499: /*@
500:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

502:   Not collective

504:   Input Parameter:
505: . fem - The PetscFE object

507:   Output Parameter:
508: . sp - The PetscSpace object

510:   Level: intermediate

512: .seealso: PetscFECreate()
513: @*/
514: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
515: {
519:   *sp = fem->basisSpace;
520:   return(0);
521: }

523: /*@
524:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

526:   Not collective

528:   Input Parameters:
529: + fem - The PetscFE object
530: - sp - The PetscSpace object

532:   Level: intermediate

534: .seealso: PetscFECreate()
535: @*/
536: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
537: {

543:   PetscSpaceDestroy(&fem->basisSpace);
544:   fem->basisSpace = sp;
545:   PetscObjectReference((PetscObject) fem->basisSpace);
546:   return(0);
547: }

549: /*@
550:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

552:   Not collective

554:   Input Parameter:
555: . fem - The PetscFE object

557:   Output Parameter:
558: . sp - The PetscDualSpace object

560:   Level: intermediate

562: .seealso: PetscFECreate()
563: @*/
564: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
565: {
569:   *sp = fem->dualSpace;
570:   return(0);
571: }

573: /*@
574:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

576:   Not collective

578:   Input Parameters:
579: + fem - The PetscFE object
580: - sp - The PetscDualSpace object

582:   Level: intermediate

584: .seealso: PetscFECreate()
585: @*/
586: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
587: {

593:   PetscDualSpaceDestroy(&fem->dualSpace);
594:   fem->dualSpace = sp;
595:   PetscObjectReference((PetscObject) fem->dualSpace);
596:   return(0);
597: }

599: /*@
600:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

602:   Not collective

604:   Input Parameter:
605: . fem - The PetscFE object

607:   Output Parameter:
608: . q - The PetscQuadrature object

610:   Level: intermediate

612: .seealso: PetscFECreate()
613: @*/
614: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
615: {
619:   *q = fem->quadrature;
620:   return(0);
621: }

623: /*@
624:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

626:   Not collective

628:   Input Parameters:
629: + fem - The PetscFE object
630: - q - The PetscQuadrature object

632:   Level: intermediate

634: .seealso: PetscFECreate()
635: @*/
636: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
637: {
638:   PetscInt       Nc, qNc;

643:   PetscFEGetNumComponents(fem, &Nc);
644:   PetscQuadratureGetNumComponents(q, &qNc);
645:   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);
646:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
647:   PetscQuadratureDestroy(&fem->quadrature);
648:   fem->quadrature = q;
649:   PetscObjectReference((PetscObject) q);
650:   return(0);
651: }

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

656:   Not collective

658:   Input Parameter:
659: . fem - The PetscFE object

661:   Output Parameter:
662: . q - The PetscQuadrature object

664:   Level: intermediate

666: .seealso: PetscFECreate()
667: @*/
668: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
669: {
673:   *q = fem->faceQuadrature;
674:   return(0);
675: }

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

680:   Not collective

682:   Input Parameters:
683: + fem - The PetscFE object
684: - q - The PetscQuadrature object

686:   Level: intermediate

688: .seealso: PetscFECreate()
689: @*/
690: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
691: {

696:   PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
697:   PetscQuadratureDestroy(&fem->faceQuadrature);
698:   fem->faceQuadrature = q;
699:   PetscObjectReference((PetscObject) q);
700:   return(0);
701: }

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

706:   Not collective

708:   Input Parameter:
709: . fem - The PetscFE object

711:   Output Parameter:
712: . numDof - Array with the number of dofs per dimension

714:   Level: intermediate

716: .seealso: PetscFECreate()
717: @*/
718: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
719: {

725:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
726:   return(0);
727: }

729: /*@C
730:   PetscFEGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points

732:   Not collective

734:   Input Parameter:
735: . fem - The PetscFE object

737:   Output Parameters:
738: + B - The basis function values at quadrature points
739: . D - The basis function derivatives at quadrature points
740: - H - The basis function second derivatives at quadrature points

742:   Note:
743: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
744: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
745: $ 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

747:   Level: intermediate

749: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
750: @*/
751: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
752: {
753:   PetscInt         npoints;
754:   const PetscReal *points;
755:   PetscErrorCode   ierr;

762:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
763:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
764:   if (B) *B = fem->B;
765:   if (D) *D = fem->D;
766:   if (H) *H = fem->H;
767:   return(0);
768: }

770: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
771: {
772:   PetscErrorCode   ierr;

779:   if (!fem->Bf) {
780:     const PetscReal  xi0[3] = {-1., -1., -1.};
781:     PetscReal        v0[3], J[9], detJ;
782:     PetscQuadrature  fq;
783:     PetscDualSpace   sp;
784:     DM               dm;
785:     const PetscInt  *faces;
786:     PetscInt         dim, numFaces, f, npoints, q;
787:     const PetscReal *points;
788:     PetscReal       *facePoints;

790:     PetscFEGetDualSpace(fem, &sp);
791:     PetscDualSpaceGetDM(sp, &dm);
792:     DMGetDimension(dm, &dim);
793:     DMPlexGetConeSize(dm, 0, &numFaces);
794:     DMPlexGetCone(dm, 0, &faces);
795:     PetscFEGetFaceQuadrature(fem, &fq);
796:     if (fq) {
797:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
798:       PetscMalloc1(numFaces*npoints*dim, &facePoints);
799:       for (f = 0; f < numFaces; ++f) {
800:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
801:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
802:       }
803:       PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
804:       PetscFree(facePoints);
805:     }
806:   }
807:   if (Bf) *Bf = fem->Bf;
808:   if (Df) *Df = fem->Df;
809:   if (Hf) *Hf = fem->Hf;
810:   return(0);
811: }

813: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
814: {
815:   PetscErrorCode   ierr;

820:   if (!fem->F) {
821:     PetscDualSpace  sp;
822:     DM              dm;
823:     const PetscInt *cone;
824:     PetscReal      *centroids;
825:     PetscInt        dim, numFaces, f;

827:     PetscFEGetDualSpace(fem, &sp);
828:     PetscDualSpaceGetDM(sp, &dm);
829:     DMGetDimension(dm, &dim);
830:     DMPlexGetConeSize(dm, 0, &numFaces);
831:     DMPlexGetCone(dm, 0, &cone);
832:     PetscMalloc1(numFaces*dim, &centroids);
833:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
834:     PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
835:     PetscFree(centroids);
836:   }
837:   *F = fem->F;
838:   return(0);
839: }

841: /*@C
842:   PetscFEGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

844:   Not collective

846:   Input Parameters:
847: + fem     - The PetscFE object
848: . npoints - The number of tabulation points
849: - points  - The tabulation point coordinates

851:   Output Parameters:
852: + B - The basis function values at tabulation points
853: . D - The basis function derivatives at tabulation points
854: - H - The basis function second derivatives at tabulation points

856:   Note:
857: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
858: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
859: $ 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

861:   Level: intermediate

863: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
864: @*/
865: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
866: {
867:   DM               dm;
868:   PetscInt         pdim; /* Dimension of FE space P */
869:   PetscInt         dim;  /* Spatial dimension */
870:   PetscInt         comp; /* Field components */
871:   PetscErrorCode   ierr;

874:   if (!npoints) {
875:     if (B) *B = NULL;
876:     if (D) *D = NULL;
877:     if (H) *H = NULL;
878:     return(0);
879:   }
885:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
886:   DMGetDimension(dm, &dim);
887:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
888:   PetscFEGetNumComponents(fem, &comp);
889:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
890:   if (!dim) {
891:     if (D) *D = NULL;
892:     if (H) *H = NULL;
893:   } else {
894:     if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
895:     if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
896:   }
897:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
898:   return(0);
899: }

901: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
902: {
903:   DM             dm;

908:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
909:   if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
910:   if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
911:   if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
912:   return(0);
913: }

915: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
916: {
917:   PetscSpace     bsp, bsubsp;
918:   PetscDualSpace dsp, dsubsp;
919:   PetscInt       dim, depth, numComp, i, j, coneSize, order;
920:   PetscFEType    type;
921:   DM             dm;
922:   DMLabel        label;
923:   PetscReal      *xi, *v, *J, detJ;
924:   const char     *name;
925:   PetscQuadrature origin, fullQuad, subQuad;

931:   PetscFEGetBasisSpace(fe,&bsp);
932:   PetscFEGetDualSpace(fe,&dsp);
933:   PetscDualSpaceGetDM(dsp,&dm);
934:   DMGetDimension(dm,&dim);
935:   DMPlexGetDepthLabel(dm,&label);
936:   DMLabelGetValue(label,refPoint,&depth);
937:   PetscCalloc1(depth,&xi);
938:   PetscMalloc1(dim,&v);
939:   PetscMalloc1(dim*dim,&J);
940:   for (i = 0; i < depth; i++) xi[i] = 0.;
941:   PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
942:   PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
943:   DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
944:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
945:   for (i = 1; i < dim; i++) {
946:     for (j = 0; j < depth; j++) {
947:       J[i * depth + j] = J[i * dim + j];
948:     }
949:   }
950:   PetscQuadratureDestroy(&origin);
951:   PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
952:   PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
953:   PetscSpaceSetUp(bsubsp);
954:   PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
955:   PetscFEGetType(fe,&type);
956:   PetscFESetType(*trFE,type);
957:   PetscFEGetNumComponents(fe,&numComp);
958:   PetscFESetNumComponents(*trFE,numComp);
959:   PetscFESetBasisSpace(*trFE,bsubsp);
960:   PetscFESetDualSpace(*trFE,dsubsp);
961:   PetscObjectGetName((PetscObject) fe, &name);
962:   if (name) {PetscFESetName(*trFE, name);}
963:   PetscFEGetQuadrature(fe,&fullQuad);
964:   PetscQuadratureGetOrder(fullQuad,&order);
965:   DMPlexGetConeSize(dm,refPoint,&coneSize);
966:   if (coneSize == 2 * depth) {
967:     PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
968:   } else {
969:     PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
970:   }
971:   PetscFESetQuadrature(*trFE,subQuad);
972:   PetscFESetUp(*trFE);
973:   PetscQuadratureDestroy(&subQuad);
974:   PetscSpaceDestroy(&bsubsp);
975:   return(0);
976: }

978: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
979: {
980:   PetscInt       hStart, hEnd;
981:   PetscDualSpace dsp;
982:   DM             dm;

988:   *trFE = NULL;
989:   PetscFEGetDualSpace(fe,&dsp);
990:   PetscDualSpaceGetDM(dsp,&dm);
991:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
992:   if (hEnd <= hStart) return(0);
993:   PetscFECreatePointTrace(fe,hStart,trFE);
994:   return(0);
995: }


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

1001:   Not collective

1003:   Input Parameter:
1004: . fe - The PetscFE

1006:   Output Parameter:
1007: . dim - The dimension

1009:   Level: intermediate

1011: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1012: @*/
1013: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1014: {

1020:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1021:   return(0);
1022: }

1024: /*
1025: Purpose: Compute element vector for chunk of elements

1027: Input:
1028:   Sizes:
1029:      Ne:  number of elements
1030:      Nf:  number of fields
1031:      PetscFE
1032:        dim: spatial dimension
1033:        Nb:  number of basis functions
1034:        Nc:  number of field components
1035:        PetscQuadrature
1036:          Nq:  number of quadrature points

1038:   Geometry:
1039:      PetscFEGeom[Ne] possibly *Nq
1040:        PetscReal v0s[dim]
1041:        PetscReal n[dim]
1042:        PetscReal jacobians[dim*dim]
1043:        PetscReal jacobianInverses[dim*dim]
1044:        PetscReal jacobianDeterminants
1045:   FEM:
1046:      PetscFE
1047:        PetscQuadrature
1048:          PetscReal   quadPoints[Nq*dim]
1049:          PetscReal   quadWeights[Nq]
1050:        PetscReal   basis[Nq*Nb*Nc]
1051:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1052:      PetscScalar coefficients[Ne*Nb*Nc]
1053:      PetscScalar elemVec[Ne*Nb*Nc]

1055:   Problem:
1056:      PetscInt f: the active field
1057:      f0, f1

1059:   Work Space:
1060:      PetscFE
1061:        PetscScalar f0[Nq*dim];
1062:        PetscScalar f1[Nq*dim*dim];
1063:        PetscScalar u[Nc];
1064:        PetscScalar gradU[Nc*dim];
1065:        PetscReal   x[dim];
1066:        PetscScalar realSpaceDer[dim];

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

1070: Input:
1071:   Sizes:
1072:      N_cb: Number of serial cell batches

1074:   Geometry:
1075:      PetscReal v0s[Ne*dim]
1076:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1077:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1078:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1079:   FEM:
1080:      static PetscReal   quadPoints[Nq*dim]
1081:      static PetscReal   quadWeights[Nq]
1082:      static PetscReal   basis[Nq*Nb*Nc]
1083:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1084:      PetscScalar coefficients[Ne*Nb*Nc]
1085:      PetscScalar elemVec[Ne*Nb*Nc]

1087: ex62.c:
1088:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1089:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1090:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1091:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1093: ex52.c:
1094:   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)
1095:   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)

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

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

1104: ex52_integrateElementOpenCL.c:
1105: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1106:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1107:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1115:   Not collective

1117:   Input Parameters:
1118: + fem          - The PetscFE object for the field being integrated
1119: . prob         - The PetscDS specifying the discretizations and continuum functions
1120: . field        - The field being integrated
1121: . Ne           - The number of elements in the chunk
1122: . cgeom        - The cell geometry for each cell in the chunk
1123: . coefficients - The array of FEM basis coefficients for the elements
1124: . probAux      - The PetscDS specifying the auxiliary discretizations
1125: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1127:   Output Parameter
1128: . integral     - the integral for this field

1130:   Level: developer

1132: .seealso: PetscFEIntegrateResidual()
1133: @*/
1134: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1135:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1136: {

1142:   if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1143:   return(0);
1144: }

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

1149:   Not collective

1151:   Input Parameters:
1152: + fem          - The PetscFE object for the field being integrated
1153: . prob         - The PetscDS specifying the discretizations and continuum functions
1154: . field        - The field being integrated
1155: . obj_func     - The function to be integrated
1156: . Ne           - The number of elements in the chunk
1157: . fgeom        - The face geometry for each face in the chunk
1158: . coefficients - The array of FEM basis coefficients for the elements
1159: . probAux      - The PetscDS specifying the auxiliary discretizations
1160: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1162:   Output Parameter
1163: . integral     - the integral for this field

1165:   Level: developer

1167: .seealso: PetscFEIntegrateResidual()
1168: @*/
1169: PetscErrorCode PetscFEIntegrateBd(PetscFE fem, PetscDS prob, PetscInt field,
1170:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1171:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1172:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1173:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1174:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1175: {

1181:   if (fem->ops->integratebd) {(*fem->ops->integratebd)(fem, prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1182:   return(0);
1183: }

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

1188:   Not collective

1190:   Input Parameters:
1191: + fem          - The PetscFE object for the field being integrated
1192: . prob         - The PetscDS specifying the discretizations and continuum functions
1193: . field        - The field being integrated
1194: . Ne           - The number of elements in the chunk
1195: . cgeom        - The cell geometry for each cell in the chunk
1196: . coefficients - The array of FEM basis coefficients for the elements
1197: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1198: . probAux      - The PetscDS specifying the auxiliary discretizations
1199: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1200: - t            - The time

1202:   Output Parameter
1203: . elemVec      - the element residual vectors from each element

1205:   Note:
1206: $ Loop over batch of elements (e):
1207: $   Loop over quadrature points (q):
1208: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1209: $     Call f_0 and f_1
1210: $   Loop over element vector entries (f,fc --> i):
1211: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1213:   Level: developer

1215: .seealso: PetscFEIntegrateResidual()
1216: @*/
1217: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1218:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1219: {

1225:   if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1226:   return(0);
1227: }

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

1232:   Not collective

1234:   Input Parameters:
1235: + fem          - The PetscFE object for the field being integrated
1236: . prob         - The PetscDS specifying the discretizations and continuum functions
1237: . field        - The field being integrated
1238: . Ne           - The number of elements in the chunk
1239: . fgeom        - The face geometry for each cell in the chunk
1240: . coefficients - The array of FEM basis coefficients for the elements
1241: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1242: . probAux      - The PetscDS specifying the auxiliary discretizations
1243: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1244: - t            - The time

1246:   Output Parameter
1247: . elemVec      - the element residual vectors from each element

1249:   Level: developer

1251: .seealso: PetscFEIntegrateResidual()
1252: @*/
1253: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1254:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1255: {

1260:   if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1261:   return(0);
1262: }

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

1267:   Not collective

1269:   Input Parameters:
1270: + fem          - The PetscFE object for the field being integrated
1271: . prob         - The PetscDS specifying the discretizations and continuum functions
1272: . jtype        - The type of matrix pointwise functions that should be used
1273: . fieldI       - The test field being integrated
1274: . fieldJ       - The basis field being integrated
1275: . Ne           - The number of elements in the chunk
1276: . cgeom        - The cell geometry for each cell in the chunk
1277: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1278: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1279: . probAux      - The PetscDS specifying the auxiliary discretizations
1280: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1281: . t            - The time
1282: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1284:   Output Parameter
1285: . elemMat      - the element matrices for the Jacobian from each element

1287:   Note:
1288: $ Loop over batch of elements (e):
1289: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1290: $     Loop over quadrature points (q):
1291: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1292: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1293: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1294: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1295: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1296:   Level: developer

1298: .seealso: PetscFEIntegrateResidual()
1299: @*/
1300: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1301:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1302: {

1307:   if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1308:   return(0);
1309: }

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

1314:   Not collective

1316:   Input Parameters:
1317: + fem          = The PetscFE object for the field being integrated
1318: . prob         - The PetscDS specifying the discretizations and continuum functions
1319: . fieldI       - The test field being integrated
1320: . fieldJ       - The basis field being integrated
1321: . Ne           - The number of elements in the chunk
1322: . fgeom        - The face geometry for each cell in the chunk
1323: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1324: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1325: . probAux      - The PetscDS specifying the auxiliary discretizations
1326: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1327: . t            - The time
1328: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1330:   Output Parameter
1331: . elemMat              - the element matrices for the Jacobian from each element

1333:   Note:
1334: $ Loop over batch of elements (e):
1335: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1336: $     Loop over quadrature points (q):
1337: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1338: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1339: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1340: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1341: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1342:   Level: developer

1344: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1345: @*/
1346: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1347:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1348: {

1353:   if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1354:   return(0);
1355: }

1357: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1358: {
1359:   PetscSpace      P, subP;
1360:   PetscDualSpace  Q, subQ;
1361:   PetscQuadrature subq;
1362:   PetscFEType     fetype;
1363:   PetscInt        dim, Nc;
1364:   PetscErrorCode  ierr;

1369:   if (height == 0) {
1370:     *subfe = fe;
1371:     return(0);
1372:   }
1373:   PetscFEGetBasisSpace(fe, &P);
1374:   PetscFEGetDualSpace(fe, &Q);
1375:   PetscFEGetNumComponents(fe, &Nc);
1376:   PetscFEGetFaceQuadrature(fe, &subq);
1377:   PetscDualSpaceGetDimension(Q, &dim);
1378:   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);}
1379:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1380:   if (height <= dim) {
1381:     if (!fe->subspaces[height-1]) {
1382:       PetscFE     sub;
1383:       const char *name;

1385:       PetscSpaceGetHeightSubspace(P, height, &subP);
1386:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1387:       PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1388:       PetscObjectGetName((PetscObject) fe,  &name);
1389:       PetscObjectSetName((PetscObject) sub,  name);
1390:       PetscFEGetType(fe, &fetype);
1391:       PetscFESetType(sub, fetype);
1392:       PetscFESetBasisSpace(sub, subP);
1393:       PetscFESetDualSpace(sub, subQ);
1394:       PetscFESetNumComponents(sub, Nc);
1395:       PetscFESetUp(sub);
1396:       PetscFESetQuadrature(sub, subq);
1397:       fe->subspaces[height-1] = sub;
1398:     }
1399:     *subfe = fe->subspaces[height-1];
1400:   } else {
1401:     *subfe = NULL;
1402:   }
1403:   return(0);
1404: }

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

1411:   Collective on PetscFE

1413:   Input Parameter:
1414: . fe - The initial PetscFE

1416:   Output Parameter:
1417: . feRef - The refined PetscFE

1419:   Level: developer

1421: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1422: @*/
1423: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1424: {
1425:   PetscSpace       P, Pref;
1426:   PetscDualSpace   Q, Qref;
1427:   DM               K, Kref;
1428:   PetscQuadrature  q, qref;
1429:   const PetscReal *v0, *jac;
1430:   PetscInt         numComp, numSubelements;
1431:   PetscErrorCode   ierr;

1434:   PetscFEGetBasisSpace(fe, &P);
1435:   PetscFEGetDualSpace(fe, &Q);
1436:   PetscFEGetQuadrature(fe, &q);
1437:   PetscDualSpaceGetDM(Q, &K);
1438:   /* Create space */
1439:   PetscObjectReference((PetscObject) P);
1440:   Pref = P;
1441:   /* Create dual space */
1442:   PetscDualSpaceDuplicate(Q, &Qref);
1443:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1444:   PetscDualSpaceSetDM(Qref, Kref);
1445:   DMDestroy(&Kref);
1446:   PetscDualSpaceSetUp(Qref);
1447:   /* Create element */
1448:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1449:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1450:   PetscFESetBasisSpace(*feRef, Pref);
1451:   PetscFESetDualSpace(*feRef, Qref);
1452:   PetscFEGetNumComponents(fe,    &numComp);
1453:   PetscFESetNumComponents(*feRef, numComp);
1454:   PetscFESetUp(*feRef);
1455:   PetscSpaceDestroy(&Pref);
1456:   PetscDualSpaceDestroy(&Qref);
1457:   /* Create quadrature */
1458:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1459:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1460:   PetscFESetQuadrature(*feRef, qref);
1461:   PetscQuadratureDestroy(&qref);
1462:   return(0);
1463: }

1465: /*@C
1466:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1468:   Collective on DM

1470:   Input Parameters:
1471: + comm      - The MPI comm
1472: . dim       - The spatial dimension
1473: . Nc        - The number of components
1474: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1475: . prefix    - The options prefix, or NULL
1476: - qorder    - The quadrature order

1478:   Output Parameter:
1479: . fem - The PetscFE object

1481:   Level: beginner

1483: .keywords: PetscFE, finite element
1484: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1485: @*/
1486: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1487: {
1488:   PetscQuadrature q, fq;
1489:   DM              K;
1490:   PetscSpace      P;
1491:   PetscDualSpace  Q;
1492:   PetscInt        order, quadPointsPerEdge;
1493:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1494:   PetscErrorCode  ierr;

1497:   /* Create space */
1498:   PetscSpaceCreate(comm, &P);
1499:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1500:   PetscSpacePolynomialSetTensor(P, tensor);
1501:   PetscSpaceSetNumComponents(P, Nc);
1502:   PetscSpaceSetNumVariables(P, dim);
1503:   PetscSpaceSetFromOptions(P);
1504:   PetscSpaceSetUp(P);
1505:   PetscSpaceGetDegree(P, &order, NULL);
1506:   PetscSpacePolynomialGetTensor(P, &tensor);
1507:   /* Create dual space */
1508:   PetscDualSpaceCreate(comm, &Q);
1509:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1510:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1511:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1512:   PetscDualSpaceSetDM(Q, K);
1513:   DMDestroy(&K);
1514:   PetscDualSpaceSetNumComponents(Q, Nc);
1515:   PetscDualSpaceSetOrder(Q, order);
1516:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1517:   PetscDualSpaceSetFromOptions(Q);
1518:   PetscDualSpaceSetUp(Q);
1519:   /* Create element */
1520:   PetscFECreate(comm, fem);
1521:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1522:   PetscFESetFromOptions(*fem);
1523:   PetscFESetBasisSpace(*fem, P);
1524:   PetscFESetDualSpace(*fem, Q);
1525:   PetscFESetNumComponents(*fem, Nc);
1526:   PetscFESetUp(*fem);
1527:   PetscSpaceDestroy(&P);
1528:   PetscDualSpaceDestroy(&Q);
1529:   /* Create quadrature (with specified order if given) */
1530:   qorder = qorder >= 0 ? qorder : order;
1531:   PetscObjectOptionsBegin((PetscObject)*fem);
1532:   PetscOptionsInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadture points per edge","PetscFECreateDefault",qorder,&qorder,NULL);
1533:   PetscOptionsEnd();
1534:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1535:   if (isSimplex) {
1536:     PetscDTGaussJacobiQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1537:     PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1538:   } else {
1539:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1540:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1541:   }
1542:   PetscFESetQuadrature(*fem, q);
1543:   PetscFESetFaceQuadrature(*fem, fq);
1544:   PetscQuadratureDestroy(&q);
1545:   PetscQuadratureDestroy(&fq);
1546:   return(0);
1547: }

1549: /*@C
1550:   PetscFESetName - Names the FE and its subobjects

1552:   Not collective

1554:   Input Parameters:
1555: + fe   - The PetscFE
1556: - name - The name

1558:   Level: beginner

1560: .keywords: PetscFE, finite element
1561: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1562: @*/
1563: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1564: {
1565:   PetscSpace     P;
1566:   PetscDualSpace Q;

1570:   PetscFEGetBasisSpace(fe, &P);
1571:   PetscFEGetDualSpace(fe, &Q);
1572:   PetscObjectSetName((PetscObject) fe, name);
1573:   PetscObjectSetName((PetscObject) P,  name);
1574:   PetscObjectSetName((PetscObject) Q,  name);
1575:   return(0);
1576: }