Actual source code: fe.c

petsc-3.10.3 2018-12-18
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: - v   - the viewer

178:   Level: developer

180: .seealso PetscFEDestroy()
181: @*/
182: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
183: {

188:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);}
189:   if (fem->ops->view) {(*fem->ops->view)(fem, v);}
190:   return(0);
191: }

193: /*@
194:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

196:   Collective on PetscFE

198:   Input Parameter:
199: . fem - the PetscFE object to set options for

201:   Options Database:
202: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
203: . -petscfe_num_batches the number of cell batches to integrate serially

205:   Level: developer

207: .seealso PetscFEView()
208: @*/
209: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
210: {
211:   const char    *defaultType;
212:   char           name[256];
213:   PetscBool      flg;

218:   if (!((PetscObject) fem)->type_name) {
219:     defaultType = PETSCFEBASIC;
220:   } else {
221:     defaultType = ((PetscObject) fem)->type_name;
222:   }
223:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

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

244: /*@C
245:   PetscFESetUp - Construct data structures for the PetscFE

247:   Collective on PetscFE

249:   Input Parameter:
250: . fem - the PetscFE object to setup

252:   Level: developer

254: .seealso PetscFEView(), PetscFEDestroy()
255: @*/
256: PetscErrorCode PetscFESetUp(PetscFE fem)
257: {

262:   if (fem->setupcalled) return(0);
263:   fem->setupcalled = PETSC_TRUE;
264:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
265:   return(0);
266: }

268: /*@
269:   PetscFEDestroy - Destroys a PetscFE object

271:   Collective on PetscFE

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

276:   Level: developer

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

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

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

291:   if ((*fem)->subspaces) {
292:     PetscInt dim, d;

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

307:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
308:   PetscHeaderDestroy(fem);
309:   return(0);
310: }

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

315:   Collective on MPI_Comm

317:   Input Parameter:
318: . comm - The communicator for the PetscFE object

320:   Output Parameter:
321: . fem - The PetscFE object

323:   Level: beginner

325: .seealso: PetscFESetType(), PETSCFEGALERKIN
326: @*/
327: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
328: {
329:   PetscFE        f;

334:   PetscCitationsRegister(FECitation,&FEcite);
335:   *fem = NULL;
336:   PetscFEInitializePackage();

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

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

358:   *fem = f;
359:   return(0);
360: }

362: /*@
363:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

365:   Not collective

367:   Input Parameter:
368: . fem - The PetscFE object

370:   Output Parameter:
371: . dim - The spatial dimension

373:   Level: intermediate

375: .seealso: PetscFECreate()
376: @*/
377: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
378: {
379:   DM             dm;

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

390: /*@
391:   PetscFESetNumComponents - Sets the number of components in the element

393:   Not collective

395:   Input Parameters:
396: + fem - The PetscFE object
397: - comp - The number of field components

399:   Level: intermediate

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

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

414:   Not collective

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

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

422:   Level: intermediate

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

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

438:   Not collective

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

447:   Level: intermediate

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

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

465:   Not collective

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

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

476:   Level: intermediate

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

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

498:   Not collective

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

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

506:   Level: intermediate

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

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

522:   Not collective

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

528:   Level: intermediate

530: .seealso: PetscFECreate()
531: @*/
532: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
533: {

539:   PetscSpaceDestroy(&fem->basisSpace);
540:   fem->basisSpace = sp;
541:   PetscObjectReference((PetscObject) fem->basisSpace);
542:   return(0);
543: }

545: /*@
546:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

548:   Not collective

550:   Input Parameter:
551: . fem - The PetscFE object

553:   Output Parameter:
554: . sp - The PetscDualSpace object

556:   Level: intermediate

558: .seealso: PetscFECreate()
559: @*/
560: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
561: {
565:   *sp = fem->dualSpace;
566:   return(0);
567: }

569: /*@
570:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

572:   Not collective

574:   Input Parameters:
575: + fem - The PetscFE object
576: - sp - The PetscDualSpace object

578:   Level: intermediate

580: .seealso: PetscFECreate()
581: @*/
582: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
583: {

589:   PetscDualSpaceDestroy(&fem->dualSpace);
590:   fem->dualSpace = sp;
591:   PetscObjectReference((PetscObject) fem->dualSpace);
592:   return(0);
593: }

595: /*@
596:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

598:   Not collective

600:   Input Parameter:
601: . fem - The PetscFE object

603:   Output Parameter:
604: . q - The PetscQuadrature object

606:   Level: intermediate

608: .seealso: PetscFECreate()
609: @*/
610: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
611: {
615:   *q = fem->quadrature;
616:   return(0);
617: }

619: /*@
620:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

622:   Not collective

624:   Input Parameters:
625: + fem - The PetscFE object
626: - q - The PetscQuadrature object

628:   Level: intermediate

630: .seealso: PetscFECreate()
631: @*/
632: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
633: {
634:   PetscInt       Nc, qNc;

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

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

652:   Not collective

654:   Input Parameter:
655: . fem - The PetscFE object

657:   Output Parameter:
658: . q - The PetscQuadrature object

660:   Level: intermediate

662: .seealso: PetscFECreate()
663: @*/
664: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
665: {
669:   *q = fem->faceQuadrature;
670:   return(0);
671: }

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

676:   Not collective

678:   Input Parameters:
679: + fem - The PetscFE object
680: - q - The PetscQuadrature object

682:   Level: intermediate

684: .seealso: PetscFECreate()
685: @*/
686: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
687: {

692:   PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
693:   PetscQuadratureDestroy(&fem->faceQuadrature);
694:   fem->faceQuadrature = q;
695:   PetscObjectReference((PetscObject) q);
696:   return(0);
697: }

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

702:   Not collective

704:   Input Parameter:
705: . fem - The PetscFE object

707:   Output Parameter:
708: . numDof - Array with the number of dofs per dimension

710:   Level: intermediate

712: .seealso: PetscFECreate()
713: @*/
714: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
715: {

721:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
722:   return(0);
723: }

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

728:   Not collective

730:   Input Parameter:
731: . fem - The PetscFE object

733:   Output Parameters:
734: + B - The basis function values at quadrature points
735: . D - The basis function derivatives at quadrature points
736: - H - The basis function second derivatives at quadrature points

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

743:   Level: intermediate

745: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
746: @*/
747: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
748: {
749:   PetscInt         npoints;
750:   const PetscReal *points;
751:   PetscErrorCode   ierr;

758:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
759:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
760:   if (B) *B = fem->B;
761:   if (D) *D = fem->D;
762:   if (H) *H = fem->H;
763:   return(0);
764: }

766: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
767: {
768:   PetscErrorCode   ierr;

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

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

809: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
810: {
811:   PetscErrorCode   ierr;

816:   if (!fem->F) {
817:     PetscDualSpace  sp;
818:     DM              dm;
819:     const PetscInt *cone;
820:     PetscReal      *centroids;
821:     PetscInt        dim, numFaces, f;

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

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

840:   Not collective

842:   Input Parameters:
843: + fem     - The PetscFE object
844: . npoints - The number of tabulation points
845: - points  - The tabulation point coordinates

847:   Output Parameters:
848: + B - The basis function values at tabulation points
849: . D - The basis function derivatives at tabulation points
850: - H - The basis function second derivatives at tabulation points

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

857:   Level: intermediate

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

870:   if (!npoints) {
871:     if (B) *B = NULL;
872:     if (D) *D = NULL;
873:     if (H) *H = NULL;
874:     return(0);
875:   }
881:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
882:   DMGetDimension(dm, &dim);
883:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
884:   PetscFEGetNumComponents(fem, &comp);
885:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
886:   if (!dim) {
887:     if (D) *D = NULL;
888:     if (H) *H = NULL;
889:   } else {
890:     if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
891:     if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
892:   }
893:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
894:   return(0);
895: }

897: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
898: {
899:   DM             dm;

904:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
905:   if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
906:   if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
907:   if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
908:   return(0);
909: }

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

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

971: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
972: {
973:   PetscInt       hStart, hEnd;
974:   PetscDualSpace dsp;
975:   DM             dm;

981:   *trFE = NULL;
982:   PetscFEGetDualSpace(fe,&dsp);
983:   PetscDualSpaceGetDM(dsp,&dm);
984:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
985:   if (hEnd <= hStart) return(0);
986:   PetscFECreatePointTrace(fe,hStart,trFE);
987:   return(0);
988: }


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

994:   Not collective

996:   Input Parameter:
997: . fe - The PetscFE

999:   Output Parameter:
1000: . dim - The dimension

1002:   Level: intermediate

1004: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1005: @*/
1006: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1007: {

1013:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1014:   return(0);
1015: }

1017: /*
1018: Purpose: Compute element vector for chunk of elements

1020: Input:
1021:   Sizes:
1022:      Ne:  number of elements
1023:      Nf:  number of fields
1024:      PetscFE
1025:        dim: spatial dimension
1026:        Nb:  number of basis functions
1027:        Nc:  number of field components
1028:        PetscQuadrature
1029:          Nq:  number of quadrature points

1031:   Geometry:
1032:      PetscFEGeom[Ne] possibly *Nq
1033:        PetscReal v0s[dim]
1034:        PetscReal n[dim]
1035:        PetscReal jacobians[dim*dim]
1036:        PetscReal jacobianInverses[dim*dim]
1037:        PetscReal jacobianDeterminants
1038:   FEM:
1039:      PetscFE
1040:        PetscQuadrature
1041:          PetscReal   quadPoints[Nq*dim]
1042:          PetscReal   quadWeights[Nq]
1043:        PetscReal   basis[Nq*Nb*Nc]
1044:        PetscReal   basisDer[Nq*Nb*Nc*dim]
1045:      PetscScalar coefficients[Ne*Nb*Nc]
1046:      PetscScalar elemVec[Ne*Nb*Nc]

1048:   Problem:
1049:      PetscInt f: the active field
1050:      f0, f1

1052:   Work Space:
1053:      PetscFE
1054:        PetscScalar f0[Nq*dim];
1055:        PetscScalar f1[Nq*dim*dim];
1056:        PetscScalar u[Nc];
1057:        PetscScalar gradU[Nc*dim];
1058:        PetscReal   x[dim];
1059:        PetscScalar realSpaceDer[dim];

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

1063: Input:
1064:   Sizes:
1065:      N_cb: Number of serial cell batches

1067:   Geometry:
1068:      PetscReal v0s[Ne*dim]
1069:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1070:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1071:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1072:   FEM:
1073:      static PetscReal   quadPoints[Nq*dim]
1074:      static PetscReal   quadWeights[Nq]
1075:      static PetscReal   basis[Nq*Nb*Nc]
1076:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1077:      PetscScalar coefficients[Ne*Nb*Nc]
1078:      PetscScalar elemVec[Ne*Nb*Nc]

1080: ex62.c:
1081:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1082:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1083:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1084:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1086: ex52.c:
1087:   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)
1088:   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)

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

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

1097: ex52_integrateElementOpenCL.c:
1098: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1099:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1100:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1108:   Not collective

1110:   Input Parameters:
1111: + fem          - The PetscFE object for the field being integrated
1112: . prob         - The PetscDS specifying the discretizations and continuum functions
1113: . field        - The field being integrated
1114: . Ne           - The number of elements in the chunk
1115: . cgeom        - The cell geometry for each cell in the chunk
1116: . coefficients - The array of FEM basis coefficients for the elements
1117: . probAux      - The PetscDS specifying the auxiliary discretizations
1118: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1120:   Output Parameter
1121: . integral     - the integral for this field

1123:   Level: developer

1125: .seealso: PetscFEIntegrateResidual()
1126: @*/
1127: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1128:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1129: {

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

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

1142:   Not collective

1144:   Input Parameters:
1145: + fem          - The PetscFE object for the field being integrated
1146: . prob         - The PetscDS specifying the discretizations and continuum functions
1147: . field        - The field being integrated
1148: . obj_func     - The function to be integrated
1149: . Ne           - The number of elements in the chunk
1150: . fgeom        - The face geometry for each face in the chunk
1151: . coefficients - The array of FEM basis coefficients for the elements
1152: . probAux      - The PetscDS specifying the auxiliary discretizations
1153: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1155:   Output Parameter
1156: . integral     - the integral for this field

1158:   Level: developer

1160: .seealso: PetscFEIntegrateResidual()
1161: @*/
1162: PetscErrorCode PetscFEIntegrateBd(PetscFE fem, PetscDS prob, PetscInt field,
1163:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1164:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1165:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1166:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1167:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1168: {

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

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

1181:   Not collective

1183:   Input Parameters:
1184: + fem          - The PetscFE object for the field being integrated
1185: . prob         - The PetscDS specifying the discretizations and continuum functions
1186: . field        - The field being integrated
1187: . Ne           - The number of elements in the chunk
1188: . cgeom        - The cell geometry for each cell in the chunk
1189: . coefficients - The array of FEM basis coefficients for the elements
1190: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1191: . probAux      - The PetscDS specifying the auxiliary discretizations
1192: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1193: - t            - The time

1195:   Output Parameter
1196: . elemVec      - the element residual vectors from each element

1198:   Note:
1199: $ Loop over batch of elements (e):
1200: $   Loop over quadrature points (q):
1201: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1202: $     Call f_0 and f_1
1203: $   Loop over element vector entries (f,fc --> i):
1204: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1206:   Level: developer

1208: .seealso: PetscFEIntegrateResidual()
1209: @*/
1210: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1211:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1212: {

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

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

1225:   Not collective

1227:   Input Parameters:
1228: + fem          - The PetscFE object for the field being integrated
1229: . prob         - The PetscDS specifying the discretizations and continuum functions
1230: . field        - The field being integrated
1231: . Ne           - The number of elements in the chunk
1232: . fgeom        - The face geometry for each cell in the chunk
1233: . coefficients - The array of FEM basis coefficients for the elements
1234: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1235: . probAux      - The PetscDS specifying the auxiliary discretizations
1236: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1237: - t            - The time

1239:   Output Parameter
1240: . elemVec      - the element residual vectors from each element

1242:   Level: developer

1244: .seealso: PetscFEIntegrateResidual()
1245: @*/
1246: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1247:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1248: {

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

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

1260:   Not collective

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

1277:   Output Parameter
1278: . elemMat      - the element matrices for the Jacobian from each element

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

1291: .seealso: PetscFEIntegrateResidual()
1292: @*/
1293: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1294:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1295: {

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

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

1307:   Not collective

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

1323:   Output Parameter
1324: . elemMat              - the element matrices for the Jacobian from each element

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

1337: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1338: @*/
1339: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1340:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1341: {

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

1350: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1351: {
1352:   PetscSpace      P, subP;
1353:   PetscDualSpace  Q, subQ;
1354:   PetscQuadrature subq;
1355:   PetscFEType     fetype;
1356:   PetscInt        dim, Nc;
1357:   PetscErrorCode  ierr;

1362:   if (height == 0) {
1363:     *subfe = fe;
1364:     return(0);
1365:   }
1366:   PetscFEGetBasisSpace(fe, &P);
1367:   PetscFEGetDualSpace(fe, &Q);
1368:   PetscFEGetNumComponents(fe, &Nc);
1369:   PetscFEGetFaceQuadrature(fe, &subq);
1370:   PetscDualSpaceGetDimension(Q, &dim);
1371:   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);}
1372:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1373:   if (height <= dim) {
1374:     if (!fe->subspaces[height-1]) {
1375:       PetscFE sub;

1377:       PetscSpaceGetHeightSubspace(P, height, &subP);
1378:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1379:       PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1380:       PetscFEGetType(fe, &fetype);
1381:       PetscFESetType(sub, fetype);
1382:       PetscFESetBasisSpace(sub, subP);
1383:       PetscFESetDualSpace(sub, subQ);
1384:       PetscFESetNumComponents(sub, Nc);
1385:       PetscFESetUp(sub);
1386:       PetscFESetQuadrature(sub, subq);
1387:       fe->subspaces[height-1] = sub;
1388:     }
1389:     *subfe = fe->subspaces[height-1];
1390:   } else {
1391:     *subfe = NULL;
1392:   }
1393:   return(0);
1394: }

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

1401:   Collective on PetscFE

1403:   Input Parameter:
1404: . fe - The initial PetscFE

1406:   Output Parameter:
1407: . feRef - The refined PetscFE

1409:   Level: developer

1411: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1412: @*/
1413: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1414: {
1415:   PetscSpace       P, Pref;
1416:   PetscDualSpace   Q, Qref;
1417:   DM               K, Kref;
1418:   PetscQuadrature  q, qref;
1419:   const PetscReal *v0, *jac;
1420:   PetscInt         numComp, numSubelements;
1421:   PetscErrorCode   ierr;

1424:   PetscFEGetBasisSpace(fe, &P);
1425:   PetscFEGetDualSpace(fe, &Q);
1426:   PetscFEGetQuadrature(fe, &q);
1427:   PetscDualSpaceGetDM(Q, &K);
1428:   /* Create space */
1429:   PetscObjectReference((PetscObject) P);
1430:   Pref = P;
1431:   /* Create dual space */
1432:   PetscDualSpaceDuplicate(Q, &Qref);
1433:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1434:   PetscDualSpaceSetDM(Qref, Kref);
1435:   DMDestroy(&Kref);
1436:   PetscDualSpaceSetUp(Qref);
1437:   /* Create element */
1438:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1439:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1440:   PetscFESetBasisSpace(*feRef, Pref);
1441:   PetscFESetDualSpace(*feRef, Qref);
1442:   PetscFEGetNumComponents(fe,    &numComp);
1443:   PetscFESetNumComponents(*feRef, numComp);
1444:   PetscFESetUp(*feRef);
1445:   PetscSpaceDestroy(&Pref);
1446:   PetscDualSpaceDestroy(&Qref);
1447:   /* Create quadrature */
1448:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1449:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1450:   PetscFESetQuadrature(*feRef, qref);
1451:   PetscQuadratureDestroy(&qref);
1452:   return(0);
1453: }

1455: /*@C
1456:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1458:   Collective on DM

1460:   Input Parameters:
1461: + comm      - The MPI comm
1462: . dim       - The spatial dimension
1463: . Nc        - The number of components
1464: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1465: . prefix    - The options prefix, or NULL
1466: - qorder    - The quadrature order

1468:   Output Parameter:
1469: . fem - The PetscFE object

1471:   Level: beginner

1473: .keywords: PetscFE, finite element
1474: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1475: @*/
1476: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1477: {
1478:   PetscQuadrature q, fq;
1479:   DM              K;
1480:   PetscSpace      P;
1481:   PetscDualSpace  Q;
1482:   PetscInt        order, quadPointsPerEdge;
1483:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1484:   PetscErrorCode  ierr;

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