Actual source code: fe.c

petsc-master 2019-11-16
Report Typos and Errors
  1: /* Basis Jet Tabulation

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

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

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

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

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

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

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

 22:    \alpha = V^{-1}

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

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

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

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

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

 51: PetscClassId PETSCFE_CLASSID = 0;

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

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

 59:   Not Collective

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

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

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

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

 83:   Level: advanced

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

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

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

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

100:   Collective on fem

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

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

109:   Level: intermediate

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

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

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

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

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

140:   Not Collective

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

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

148:   Level: intermediate

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

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

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

169:    Collective on PetscFE

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

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

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

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

192:   Collective on fem

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

198:   Level: beginner

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

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

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

220:   Collective on fem

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

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

229:   Level: intermediate

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

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

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

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

271:   Collective on fem

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

276:   Level: intermediate

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

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

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

295:   Collective on fem

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

300:   Level: beginner

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

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

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

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

318:     PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
319:     for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
320:   }
321:   PetscFree((*fem)->subspaces);
322:   PetscFree((*fem)->invV);
323:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
324:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
325:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
326:   PetscSpaceDestroy(&(*fem)->basisSpace);
327:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
328:   PetscQuadratureDestroy(&(*fem)->quadrature);
329:   PetscQuadratureDestroy(&(*fem)->faceQuadrature);

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

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

339:   Collective

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

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

347:   Level: beginner

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

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

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

364:   f->basisSpace    = NULL;
365:   f->dualSpace     = NULL;
366:   f->numComponents = 1;
367:   f->subspaces     = NULL;
368:   f->invV          = NULL;
369:   f->B             = NULL;
370:   f->D             = NULL;
371:   f->H             = NULL;
372:   f->Bf            = NULL;
373:   f->Df            = NULL;
374:   f->Hf            = NULL;
375:   PetscArrayzero(&f->quadrature, 1);
376:   PetscArrayzero(&f->faceQuadrature, 1);
377:   f->blockSize     = 0;
378:   f->numBlocks     = 1;
379:   f->batchSize     = 0;
380:   f->numBatches    = 1;

382:   *fem = f;
383:   return(0);
384: }

386: /*@
387:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

389:   Not collective

391:   Input Parameter:
392: . fem - The PetscFE object

394:   Output Parameter:
395: . dim - The spatial dimension

397:   Level: intermediate

399: .seealso: PetscFECreate()
400: @*/
401: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
402: {
403:   DM             dm;

409:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
410:   DMGetDimension(dm, dim);
411:   return(0);
412: }

414: /*@
415:   PetscFESetNumComponents - Sets the number of components in the element

417:   Not collective

419:   Input Parameters:
420: + fem - The PetscFE object
421: - comp - The number of field components

423:   Level: intermediate

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

435: /*@
436:   PetscFEGetNumComponents - Returns the number of components in the element

438:   Not collective

440:   Input Parameter:
441: . fem - The PetscFE object

443:   Output Parameter:
444: . comp - The number of field components

446:   Level: intermediate

448: .seealso: PetscFECreate()
449: @*/
450: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
451: {
455:   *comp = fem->numComponents;
456:   return(0);
457: }

459: /*@
460:   PetscFESetTileSizes - Sets the tile sizes for evaluation

462:   Not collective

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

471:   Level: intermediate

473: .seealso: PetscFECreate()
474: @*/
475: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
476: {
479:   fem->blockSize  = blockSize;
480:   fem->numBlocks  = numBlocks;
481:   fem->batchSize  = batchSize;
482:   fem->numBatches = numBatches;
483:   return(0);
484: }

486: /*@
487:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

489:   Not collective

491:   Input Parameter:
492: . fem - The PetscFE object

494:   Output Parameters:
495: + blockSize - The number of elements in a block
496: . numBlocks - The number of blocks in a batch
497: . batchSize - The number of elements in a batch
498: - numBatches - The number of batches in a chunk

500:   Level: intermediate

502: .seealso: PetscFECreate()
503: @*/
504: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
505: {
512:   if (blockSize)  *blockSize  = fem->blockSize;
513:   if (numBlocks)  *numBlocks  = fem->numBlocks;
514:   if (batchSize)  *batchSize  = fem->batchSize;
515:   if (numBatches) *numBatches = fem->numBatches;
516:   return(0);
517: }

519: /*@
520:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

522:   Not collective

524:   Input Parameter:
525: . fem - The PetscFE object

527:   Output Parameter:
528: . sp - The PetscSpace object

530:   Level: intermediate

532: .seealso: PetscFECreate()
533: @*/
534: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
535: {
539:   *sp = fem->basisSpace;
540:   return(0);
541: }

543: /*@
544:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

546:   Not collective

548:   Input Parameters:
549: + fem - The PetscFE object
550: - sp - The PetscSpace object

552:   Level: intermediate

554: .seealso: PetscFECreate()
555: @*/
556: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
557: {

563:   PetscSpaceDestroy(&fem->basisSpace);
564:   fem->basisSpace = sp;
565:   PetscObjectReference((PetscObject) fem->basisSpace);
566:   return(0);
567: }

569: /*@
570:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

572:   Not collective

574:   Input Parameter:
575: . fem - The PetscFE object

577:   Output Parameter:
578: . sp - The PetscDualSpace object

580:   Level: intermediate

582: .seealso: PetscFECreate()
583: @*/
584: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
585: {
589:   *sp = fem->dualSpace;
590:   return(0);
591: }

593: /*@
594:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

596:   Not collective

598:   Input Parameters:
599: + fem - The PetscFE object
600: - sp - The PetscDualSpace object

602:   Level: intermediate

604: .seealso: PetscFECreate()
605: @*/
606: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
607: {

613:   PetscDualSpaceDestroy(&fem->dualSpace);
614:   fem->dualSpace = sp;
615:   PetscObjectReference((PetscObject) fem->dualSpace);
616:   return(0);
617: }

619: /*@
620:   PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products

622:   Not collective

624:   Input Parameter:
625: . fem - The PetscFE object

627:   Output Parameter:
628: . q - The PetscQuadrature object

630:   Level: intermediate

632: .seealso: PetscFECreate()
633: @*/
634: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
635: {
639:   *q = fem->quadrature;
640:   return(0);
641: }

643: /*@
644:   PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products

646:   Not collective

648:   Input Parameters:
649: + fem - The PetscFE object
650: - q - The PetscQuadrature object

652:   Level: intermediate

654: .seealso: PetscFECreate()
655: @*/
656: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
657: {
658:   PetscInt       Nc, qNc;

663:   PetscFEGetNumComponents(fem, &Nc);
664:   PetscQuadratureGetNumComponents(q, &qNc);
665:   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);
666:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
667:   PetscQuadratureDestroy(&fem->quadrature);
668:   fem->quadrature = q;
669:   PetscObjectReference((PetscObject) q);
670:   return(0);
671: }

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

676:   Not collective

678:   Input Parameter:
679: . fem - The PetscFE object

681:   Output Parameter:
682: . q - The PetscQuadrature object

684:   Level: intermediate

686: .seealso: PetscFECreate()
687: @*/
688: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
689: {
693:   *q = fem->faceQuadrature;
694:   return(0);
695: }

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

700:   Not collective

702:   Input Parameters:
703: + fem - The PetscFE object
704: - q - The PetscQuadrature object

706:   Level: intermediate

708: .seealso: PetscFECreate()
709: @*/
710: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
711: {

716:   PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
717:   PetscQuadratureDestroy(&fem->faceQuadrature);
718:   fem->faceQuadrature = q;
719:   PetscObjectReference((PetscObject) q);
720:   return(0);
721: }

723: /*@
724:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

726:   Not collective

728:   Input Parameters:
729: + sfe - The PetscFE source for the quadratures
730: - tfe - The PetscFE target for the quadratures

732:   Level: intermediate

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

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

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

754:   Not collective

756:   Input Parameter:
757: . fem - The PetscFE object

759:   Output Parameter:
760: . numDof - Array with the number of dofs per dimension

762:   Level: intermediate

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

773:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
774:   return(0);
775: }

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

780:   Not collective

782:   Input Parameter:
783: . fem - The PetscFE object

785:   Output Parameters:
786: + B - The basis function values at quadrature points
787: . D - The basis function derivatives at quadrature points
788: - H - The basis function second derivatives at quadrature points

790:   Note:
791: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
792: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
793: $ H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

795:   Level: intermediate

797: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
798: @*/
799: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
800: {
801:   PetscInt         npoints;
802:   const PetscReal *points;
803:   PetscErrorCode   ierr;

810:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
811:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
812:   if (B) *B = fem->B;
813:   if (D) *D = fem->D;
814:   if (H) *H = fem->H;
815:   return(0);
816: }

818: /*@C
819:   PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points

821:   Not collective

823:   Input Parameter:
824: . fem - The PetscFE object

826:   Output Parameters:
827: + B - The basis function values at face quadrature points
828: . D - The basis function derivatives at face quadrature points
829: - H - The basis function second derivatives at face quadrature points

831:   Note:
832: $ Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
833: $ Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d
834: $ Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e

836:   Level: intermediate

838: .seealso: PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
839: @*/
840: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
841: {
842:   PetscErrorCode   ierr;

849:   if (!fem->Bf) {
850:     const PetscReal  xi0[3] = {-1., -1., -1.};
851:     PetscReal        v0[3], J[9], detJ;
852:     PetscQuadrature  fq;
853:     PetscDualSpace   sp;
854:     DM               dm;
855:     const PetscInt  *faces;
856:     PetscInt         dim, numFaces, f, npoints, q;
857:     const PetscReal *points;
858:     PetscReal       *facePoints;

860:     PetscFEGetDualSpace(fem, &sp);
861:     PetscDualSpaceGetDM(sp, &dm);
862:     DMGetDimension(dm, &dim);
863:     DMPlexGetConeSize(dm, 0, &numFaces);
864:     DMPlexGetCone(dm, 0, &faces);
865:     PetscFEGetFaceQuadrature(fem, &fq);
866:     if (fq) {
867:       PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
868:       PetscMalloc1(numFaces*npoints*dim, &facePoints);
869:       for (f = 0; f < numFaces; ++f) {
870:         DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
871:         for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
872:       }
873:       PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
874:       PetscFree(facePoints);
875:     }
876:   }
877:   if (Bf) *Bf = fem->Bf;
878:   if (Df) *Df = fem->Df;
879:   if (Hf) *Hf = fem->Hf;
880:   return(0);
881: }

883: /*@C
884:   PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face centroid points

886:   Not collective

888:   Input Parameter:
889: . fem - The PetscFE object

891:   Output Parameters:
892: + B - The basis function values at face centroid points
893: . D - The basis function derivatives at face centroid points
894: - H - The basis function second derivatives at face centroid points

896:   Note:
897: $ Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
898: $ Df[((f*pdim + i)*Nc + c)*dim + d] is the derivative value at point f for basis function i, component c, in direction d
899: $ Hf[(((f*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f for basis function i, component c, in directions d and e

901:   Level: intermediate

903: .seealso: PetscFEGetFaceTabulation(), PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
904: @*/
905: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
906: {
907:   PetscErrorCode   ierr;

912:   if (!fem->F) {
913:     PetscDualSpace  sp;
914:     DM              dm;
915:     const PetscInt *cone;
916:     PetscReal      *centroids;
917:     PetscInt        dim, numFaces, f;

919:     PetscFEGetDualSpace(fem, &sp);
920:     PetscDualSpaceGetDM(sp, &dm);
921:     DMGetDimension(dm, &dim);
922:     DMPlexGetConeSize(dm, 0, &numFaces);
923:     DMPlexGetCone(dm, 0, &cone);
924:     PetscMalloc1(numFaces*dim, &centroids);
925:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
926:     PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
927:     PetscFree(centroids);
928:   }
929:   *F = fem->F;
930:   return(0);
931: }

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

936:   Not collective

938:   Input Parameters:
939: + fem     - The PetscFE object
940: . npoints - The number of tabulation points
941: - points  - The tabulation point coordinates

943:   Output Parameters:
944: + B - The basis function values at tabulation points
945: . D - The basis function derivatives at tabulation points
946: - H - The basis function second derivatives at tabulation points

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

953:   Level: intermediate

955: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
956: @*/
957: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
958: {
959:   DM               dm;
960:   PetscInt         pdim; /* Dimension of FE space P */
961:   PetscInt         dim;  /* Spatial dimension */
962:   PetscInt         comp; /* Field components */
963:   PetscErrorCode   ierr;

966:   if (!npoints) {
967:     if (B) *B = NULL;
968:     if (D) *D = NULL;
969:     if (H) *H = NULL;
970:     return(0);
971:   }
977:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
978:   DMGetDimension(dm, &dim);
979:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
980:   PetscFEGetNumComponents(fem, &comp);
981:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
982:   if (!dim) {
983:     if (D) *D = NULL;
984:     if (H) *H = NULL;
985:   } else {
986:     if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
987:     if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
988:   }
989:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
990:   return(0);
991: }

993: /*@C
994:   PetscFERestoreTabulation - Frees memory from the associated tabulation.

996:   Not collective

998:   Input Parameters:
999: + fem     - The PetscFE object
1000: . npoints - The number of tabulation points
1001: . points  - The tabulation point coordinates
1002: . B - The basis function values at tabulation points
1003: . D - The basis function derivatives at tabulation points
1004: - H - The basis function second derivatives at tabulation points

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

1011:   Level: intermediate

1013: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation()
1014: @*/
1015: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1016: {
1017:   DM             dm;

1022:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
1023:   if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
1024:   if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
1025:   if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
1026:   return(0);
1027: }

1029: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1030: {
1031:   PetscSpace     bsp, bsubsp;
1032:   PetscDualSpace dsp, dsubsp;
1033:   PetscInt       dim, depth, numComp, i, j, coneSize, order;
1034:   PetscFEType    type;
1035:   DM             dm;
1036:   DMLabel        label;
1037:   PetscReal      *xi, *v, *J, detJ;
1038:   const char     *name;
1039:   PetscQuadrature origin, fullQuad, subQuad;

1045:   PetscFEGetBasisSpace(fe,&bsp);
1046:   PetscFEGetDualSpace(fe,&dsp);
1047:   PetscDualSpaceGetDM(dsp,&dm);
1048:   DMGetDimension(dm,&dim);
1049:   DMPlexGetDepthLabel(dm,&label);
1050:   DMLabelGetValue(label,refPoint,&depth);
1051:   PetscCalloc1(depth,&xi);
1052:   PetscMalloc1(dim,&v);
1053:   PetscMalloc1(dim*dim,&J);
1054:   for (i = 0; i < depth; i++) xi[i] = 0.;
1055:   PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1056:   PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1057:   DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1058:   /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1059:   for (i = 1; i < dim; i++) {
1060:     for (j = 0; j < depth; j++) {
1061:       J[i * depth + j] = J[i * dim + j];
1062:     }
1063:   }
1064:   PetscQuadratureDestroy(&origin);
1065:   PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1066:   PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1067:   PetscSpaceSetUp(bsubsp);
1068:   PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1069:   PetscFEGetType(fe,&type);
1070:   PetscFESetType(*trFE,type);
1071:   PetscFEGetNumComponents(fe,&numComp);
1072:   PetscFESetNumComponents(*trFE,numComp);
1073:   PetscFESetBasisSpace(*trFE,bsubsp);
1074:   PetscFESetDualSpace(*trFE,dsubsp);
1075:   PetscObjectGetName((PetscObject) fe, &name);
1076:   if (name) {PetscFESetName(*trFE, name);}
1077:   PetscFEGetQuadrature(fe,&fullQuad);
1078:   PetscQuadratureGetOrder(fullQuad,&order);
1079:   DMPlexGetConeSize(dm,refPoint,&coneSize);
1080:   if (coneSize == 2 * depth) {
1081:     PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1082:   } else {
1083:     PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1084:   }
1085:   PetscFESetQuadrature(*trFE,subQuad);
1086:   PetscFESetUp(*trFE);
1087:   PetscQuadratureDestroy(&subQuad);
1088:   PetscSpaceDestroy(&bsubsp);
1089:   return(0);
1090: }

1092: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1093: {
1094:   PetscInt       hStart, hEnd;
1095:   PetscDualSpace dsp;
1096:   DM             dm;

1102:   *trFE = NULL;
1103:   PetscFEGetDualSpace(fe,&dsp);
1104:   PetscDualSpaceGetDM(dsp,&dm);
1105:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1106:   if (hEnd <= hStart) return(0);
1107:   PetscFECreatePointTrace(fe,hStart,trFE);
1108:   return(0);
1109: }


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

1115:   Not collective

1117:   Input Parameter:
1118: . fe - The PetscFE

1120:   Output Parameter:
1121: . dim - The dimension

1123:   Level: intermediate

1125: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1126: @*/
1127: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1128: {

1134:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1135:   return(0);
1136: }

1138: /*@C
1139:   PetscFEPushforward - Map the reference element function to real space

1141:   Input Parameters:
1142: + fe     - The PetscFE
1143: . fegeom - The cell geometry
1144: . Nv     - The number of function values
1145: - vals   - The function values

1147:   Output Parameter:
1148: . vals   - The transformed function values

1150:   Level: advanced

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

1154: .seealso: PetscDualSpacePushforward()
1155: @*/
1156: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1157: {

1161:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1162:   return(0);
1163: }

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

1168:   Input Parameters:
1169: + fe     - The PetscFE
1170: . fegeom - The cell geometry
1171: . Nv     - The number of function gradient values
1172: - vals   - The function gradient values

1174:   Output Parameter:
1175: . vals   - The transformed function gradient values

1177:   Level: advanced

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

1181: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1182: @*/
1183: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1184: {

1188:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1189:   return(0);
1190: }

1192: /*
1193: Purpose: Compute element vector for chunk of elements

1195: Input:
1196:   Sizes:
1197:      Ne:  number of elements
1198:      Nf:  number of fields
1199:      PetscFE
1200:        dim: spatial dimension
1201:        Nb:  number of basis functions
1202:        Nc:  number of field components
1203:        PetscQuadrature
1204:          Nq:  number of quadrature points

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

1223:   Problem:
1224:      PetscInt f: the active field
1225:      f0, f1

1227:   Work Space:
1228:      PetscFE
1229:        PetscScalar f0[Nq*dim];
1230:        PetscScalar f1[Nq*dim*dim];
1231:        PetscScalar u[Nc];
1232:        PetscScalar gradU[Nc*dim];
1233:        PetscReal   x[dim];
1234:        PetscScalar realSpaceDer[dim];

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

1238: Input:
1239:   Sizes:
1240:      N_cb: Number of serial cell batches

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

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

1261: ex52.c:
1262:   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)
1263:   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)

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

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

1272: ex52_integrateElementOpenCL.c:
1273: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1274:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1275:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1283:   Not collective

1285:   Input Parameters:
1286: + fem          - The PetscFE object for the field being integrated
1287: . prob         - The PetscDS specifying the discretizations and continuum functions
1288: . field        - The field being integrated
1289: . Ne           - The number of elements in the chunk
1290: . cgeom        - The cell geometry for each cell in the chunk
1291: . coefficients - The array of FEM basis coefficients for the elements
1292: . probAux      - The PetscDS specifying the auxiliary discretizations
1293: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1295:   Output Parameter
1296: . integral     - the integral for this field

1298:   Level: intermediate

1300: .seealso: PetscFEIntegrateResidual()
1301: @*/
1302: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1303:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1304: {
1305:   PetscFE        fe;

1310:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1311:   if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1312:   return(0);
1313: }

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

1318:   Not collective

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

1331:   Output Parameter
1332: . integral     - the integral for this field

1334:   Level: intermediate

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

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

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

1358:   Not collective

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

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

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

1383:   Level: intermediate

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

1395:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1396:   if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1397:   return(0);
1398: }

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

1403:   Not collective

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

1417:   Output Parameter
1418: . elemVec      - the element residual vectors from each element

1420:   Level: intermediate

1422: .seealso: PetscFEIntegrateResidual()
1423: @*/
1424: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1425:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1426: {
1427:   PetscFE        fe;

1432:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1433:   if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1434:   return(0);
1435: }

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

1440:   Not collective

1442:   Input Parameters:
1443: + fem          - The PetscFE object for the field being integrated
1444: . prob         - The PetscDS specifying the discretizations and continuum functions
1445: . jtype        - The type of matrix pointwise functions that should be used
1446: . fieldI       - The test field being integrated
1447: . fieldJ       - The basis field being integrated
1448: . Ne           - The number of elements in the chunk
1449: . cgeom        - The cell geometry for each cell in the chunk
1450: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1451: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1452: . probAux      - The PetscDS specifying the auxiliary discretizations
1453: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1454: . t            - The time
1455: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1457:   Output Parameter
1458: . elemMat      - the element matrices for the Jacobian from each element

1460:   Note:
1461: $ Loop over batch of elements (e):
1462: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1463: $     Loop over quadrature points (q):
1464: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1465: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1466: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1467: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1468: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1469:   Level: intermediate

1471: .seealso: PetscFEIntegrateResidual()
1472: @*/
1473: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1474:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1475: {
1476:   PetscFE        fe;

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

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

1489:   Not collective

1491:   Input Parameters:
1492: . prob         - The PetscDS specifying the discretizations and continuum functions
1493: . fieldI       - The test field being integrated
1494: . fieldJ       - The basis field being integrated
1495: . Ne           - The number of elements in the chunk
1496: . fgeom        - The face geometry for each cell in the chunk
1497: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1498: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1499: . probAux      - The PetscDS specifying the auxiliary discretizations
1500: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1501: . t            - The time
1502: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1504:   Output Parameter
1505: . elemMat              - the element matrices for the Jacobian from each element

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

1518: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1519: @*/
1520: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1521:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1522: {
1523:   PetscFE        fe;

1528:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1529:   if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1530:   return(0);
1531: }

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

1536:   Input Parameters:
1537: + fe     - The finite element space
1538: - height - The height of the Plex point

1540:   Output Parameter:
1541: . subfe  - The subspace of this FE space

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

1545:   Level: advanced

1547: .seealso: PetscFECreateDefault()
1548: @*/
1549: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1550: {
1551:   PetscSpace      P, subP;
1552:   PetscDualSpace  Q, subQ;
1553:   PetscQuadrature subq;
1554:   PetscFEType     fetype;
1555:   PetscInt        dim, Nc;
1556:   PetscErrorCode  ierr;

1561:   if (height == 0) {
1562:     *subfe = fe;
1563:     return(0);
1564:   }
1565:   PetscFEGetBasisSpace(fe, &P);
1566:   PetscFEGetDualSpace(fe, &Q);
1567:   PetscFEGetNumComponents(fe, &Nc);
1568:   PetscFEGetFaceQuadrature(fe, &subq);
1569:   PetscDualSpaceGetDimension(Q, &dim);
1570:   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);}
1571:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1572:   if (height <= dim) {
1573:     if (!fe->subspaces[height-1]) {
1574:       PetscFE     sub;
1575:       const char *name;

1577:       PetscSpaceGetHeightSubspace(P, height, &subP);
1578:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1579:       PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1580:       PetscObjectGetName((PetscObject) fe,  &name);
1581:       PetscObjectSetName((PetscObject) sub,  name);
1582:       PetscFEGetType(fe, &fetype);
1583:       PetscFESetType(sub, fetype);
1584:       PetscFESetBasisSpace(sub, subP);
1585:       PetscFESetDualSpace(sub, subQ);
1586:       PetscFESetNumComponents(sub, Nc);
1587:       PetscFESetUp(sub);
1588:       PetscFESetQuadrature(sub, subq);
1589:       fe->subspaces[height-1] = sub;
1590:     }
1591:     *subfe = fe->subspaces[height-1];
1592:   } else {
1593:     *subfe = NULL;
1594:   }
1595:   return(0);
1596: }

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

1603:   Collective on fem

1605:   Input Parameter:
1606: . fe - The initial PetscFE

1608:   Output Parameter:
1609: . feRef - The refined PetscFE

1611:   Level: advanced

1613: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1614: @*/
1615: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1616: {
1617:   PetscSpace       P, Pref;
1618:   PetscDualSpace   Q, Qref;
1619:   DM               K, Kref;
1620:   PetscQuadrature  q, qref;
1621:   const PetscReal *v0, *jac;
1622:   PetscInt         numComp, numSubelements;
1623:   PetscErrorCode   ierr;

1626:   PetscFEGetBasisSpace(fe, &P);
1627:   PetscFEGetDualSpace(fe, &Q);
1628:   PetscFEGetQuadrature(fe, &q);
1629:   PetscDualSpaceGetDM(Q, &K);
1630:   /* Create space */
1631:   PetscObjectReference((PetscObject) P);
1632:   Pref = P;
1633:   /* Create dual space */
1634:   PetscDualSpaceDuplicate(Q, &Qref);
1635:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1636:   PetscDualSpaceSetDM(Qref, Kref);
1637:   DMDestroy(&Kref);
1638:   PetscDualSpaceSetUp(Qref);
1639:   /* Create element */
1640:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1641:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1642:   PetscFESetBasisSpace(*feRef, Pref);
1643:   PetscFESetDualSpace(*feRef, Qref);
1644:   PetscFEGetNumComponents(fe,    &numComp);
1645:   PetscFESetNumComponents(*feRef, numComp);
1646:   PetscFESetUp(*feRef);
1647:   PetscSpaceDestroy(&Pref);
1648:   PetscDualSpaceDestroy(&Qref);
1649:   /* Create quadrature */
1650:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1651:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1652:   PetscFESetQuadrature(*feRef, qref);
1653:   PetscQuadratureDestroy(&qref);
1654:   return(0);
1655: }

1657: /*@C
1658:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1660:   Collective

1662:   Input Parameters:
1663: + comm      - The MPI comm
1664: . dim       - The spatial dimension
1665: . Nc        - The number of components
1666: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1667: . prefix    - The options prefix, or NULL
1668: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1670:   Output Parameter:
1671: . fem - The PetscFE object

1673:   Level: beginner

1675: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1676: @*/
1677: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1678: {
1679:   PetscQuadrature q, fq;
1680:   DM              K;
1681:   PetscSpace      P;
1682:   PetscDualSpace  Q;
1683:   PetscInt        order, quadPointsPerEdge;
1684:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1685:   PetscErrorCode  ierr;

1688:   /* Create space */
1689:   PetscSpaceCreate(comm, &P);
1690:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1691:   PetscSpacePolynomialSetTensor(P, tensor);
1692:   PetscSpaceSetNumComponents(P, Nc);
1693:   PetscSpaceSetNumVariables(P, dim);
1694:   PetscSpaceSetFromOptions(P);
1695:   PetscSpaceSetUp(P);
1696:   PetscSpaceGetDegree(P, &order, NULL);
1697:   PetscSpacePolynomialGetTensor(P, &tensor);
1698:   /* Create dual space */
1699:   PetscDualSpaceCreate(comm, &Q);
1700:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1701:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1702:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1703:   PetscDualSpaceSetDM(Q, K);
1704:   DMDestroy(&K);
1705:   PetscDualSpaceSetNumComponents(Q, Nc);
1706:   PetscDualSpaceSetOrder(Q, order);
1707:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1708:   PetscDualSpaceSetFromOptions(Q);
1709:   PetscDualSpaceSetUp(Q);
1710:   /* Create element */
1711:   PetscFECreate(comm, fem);
1712:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1713:   PetscFESetBasisSpace(*fem, P);
1714:   PetscFESetDualSpace(*fem, Q);
1715:   PetscFESetNumComponents(*fem, Nc);
1716:   PetscFESetFromOptions(*fem);
1717:   PetscFESetUp(*fem);
1718:   PetscSpaceDestroy(&P);
1719:   PetscDualSpaceDestroy(&Q);
1720:   /* Create quadrature (with specified order if given) */
1721:   qorder = qorder >= 0 ? qorder : order;
1722:   PetscObjectOptionsBegin((PetscObject)*fem);
1723:   PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1724:   PetscOptionsEnd();
1725:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1726:   if (isSimplex) {
1727:     PetscDTGaussJacobiQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1728:     PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1729:   } else {
1730:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1731:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1732:   }
1733:   PetscFESetQuadrature(*fem, q);
1734:   PetscFESetFaceQuadrature(*fem, fq);
1735:   PetscQuadratureDestroy(&q);
1736:   PetscQuadratureDestroy(&fq);
1737:   return(0);
1738: }

1740: /*@C
1741:   PetscFESetName - Names the FE and its subobjects

1743:   Not collective

1745:   Input Parameters:
1746: + fe   - The PetscFE
1747: - name - The name

1749:   Level: intermediate

1751: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1752: @*/
1753: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1754: {
1755:   PetscSpace     P;
1756:   PetscDualSpace Q;

1760:   PetscFEGetBasisSpace(fe, &P);
1761:   PetscFEGetDualSpace(fe, &Q);
1762:   PetscObjectSetName((PetscObject) fe, name);
1763:   PetscObjectSetName((PetscObject) P,  name);
1764:   PetscObjectSetName((PetscObject) Q,  name);
1765:   return(0);
1766: }

1768: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1769: {
1770:   PetscInt       dOffset = 0, fOffset = 0, f;

1773:   for (f = 0; f < Nf; ++f) {
1774:     PetscFE          fe;
1775:     const PetscInt   Nbf = Nb[f], Ncf = Nc[f];
1776:     const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
1777:     const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
1778:     PetscInt         b, c, d;

1780:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
1781:     for (c = 0; c < Ncf; ++c)     u[fOffset+c] = 0.0;
1782:     for (d = 0; d < dim*Ncf; ++d) u_x[fOffset*dim+d] = 0.0;
1783:     for (b = 0; b < Nbf; ++b) {
1784:       for (c = 0; c < Ncf; ++c) {
1785:         const PetscInt cidx = b*Ncf+c;

1787:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
1788:         for (d = 0; d < dim; ++d) u_x[(fOffset+c)*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
1789:       }
1790:     }
1791:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
1792:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*dim]);
1793:     if (u_t) {
1794:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
1795:       for (b = 0; b < Nbf; ++b) {
1796:         for (c = 0; c < Ncf; ++c) {
1797:           const PetscInt cidx = b*Ncf+c;

1799:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
1800:         }
1801:       }
1802:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
1803:     }
1804:     fOffset += Ncf;
1805:     dOffset += Nbf;
1806:   }
1807:   return 0;
1808: }

1810: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
1811: {
1812:   PetscFE        fe;
1813:   PetscReal     *faceBasis;
1814:   PetscInt       Nb, Nc, b, c;

1817:   if (!prob) return 0;
1818:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1819:   PetscFEGetDimension(fe, &Nb);
1820:   PetscFEGetNumComponents(fe, &Nc);
1821:   PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
1822:   for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
1823:   for (b = 0; b < Nb; ++b) {
1824:     for (c = 0; c < Nc; ++c) {
1825:       const PetscInt cidx = b*Nc+c;

1827:       u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
1828:     }
1829:   }
1830:   return 0;
1831: }

1833: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
1834: {
1835:   PetscInt       q, b, c, d;

1838:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
1839:   for (q = 0; q < Nq; ++q) {
1840:     for (b = 0; b < Nb; ++b) {
1841:       for (c = 0; c < Nc; ++c) {
1842:         const PetscInt bcidx = b*Nc+c;

1844:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
1845:         for (d = 0; d < dim; ++d) tmpBasisDer[bcidx*dim+d] = basisDer[q*Nb*Nc*dim+bcidx*dim+d];
1846:       }
1847:     }
1848:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
1849:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
1850:     for (b = 0; b < Nb; ++b) {
1851:       for (c = 0; c < Nc; ++c) {
1852:         const PetscInt bcidx = b*Nc+c;
1853:         const PetscInt qcidx = q*Nc+c;

1855:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
1856:         for (d = 0; d < dim; ++d) elemVec[b] += tmpBasisDer[bcidx*dim+d]*f1[qcidx*dim+d];
1857:       }
1858:     }
1859:   }
1860:   return(0);
1861: }

1863: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt dim, PetscInt NbI, PetscInt NcI, const PetscReal basisI[], const PetscReal basisDerI[], PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscInt NbJ, PetscInt NcJ, const PetscReal basisJ[], const PetscReal basisDerJ[], PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
1864: {
1865:   PetscInt       f, fc, g, gc, df, dg;

1868:   for (f = 0; f < NbI; ++f) {
1869:     for (fc = 0; fc < NcI; ++fc) {
1870:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

1872:       tmpBasisI[fidx] = basisI[fidx];
1873:       for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dim+df] = basisDerI[fidx*dim+df];
1874:     }
1875:   }
1876:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
1877:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
1878:   for (g = 0; g < NbJ; ++g) {
1879:     for (gc = 0; gc < NcJ; ++gc) {
1880:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

1882:       tmpBasisJ[gidx] = basisJ[gidx];
1883:       for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dim+dg] = basisDerJ[gidx*dim+dg];
1884:     }
1885:   }
1886:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
1887:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
1888:   for (f = 0; f < NbI; ++f) {
1889:     for (fc = 0; fc < NcI; ++fc) {
1890:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
1891:       const PetscInt i    = offsetI+f; /* Element matrix row */
1892:       for (g = 0; g < NbJ; ++g) {
1893:         for (gc = 0; gc < NcJ; ++gc) {
1894:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
1895:           const PetscInt j    = offsetJ+g; /* Element matrix column */
1896:           const PetscInt fOff = eOffset+i*totDim+j;

1898:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
1899:           for (df = 0; df < dim; ++df) {
1900:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dim+df]*tmpBasisDerJ[gidx*dim+df];
1901:             elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g2[(fc*NcJ+gc)*dim+df]*tmpBasisJ[gidx];
1902:             for (dg = 0; dg < dim; ++dg) {
1903:               elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g3[((fc*NcJ+gc)*dim+df)*dim+dg]*tmpBasisDerJ[gidx*dim+dg];
1904:             }
1905:           }
1906:         }
1907:       }
1908:     }
1909:   }
1910:   return(0);
1911: }

1913: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
1914: {
1915:   PetscDualSpace  dsp;
1916:   DM              dm;
1917:   PetscQuadrature quadDef;
1918:   PetscInt        dim, cdim, Nq;
1919:   PetscErrorCode  ierr;

1922:   PetscFEGetDualSpace(fe, &dsp);
1923:   PetscDualSpaceGetDM(dsp, &dm);
1924:   DMGetDimension(dm, &dim);
1925:   DMGetCoordinateDim(dm, &cdim);
1926:   PetscFEGetQuadrature(fe, &quadDef);
1927:   quad = quad ? quad : quadDef;
1928:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
1929:   PetscMalloc1(Nq*cdim,      &cgeom->v);
1930:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
1931:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
1932:   PetscMalloc1(Nq,           &cgeom->detJ);
1933:   cgeom->dim       = dim;
1934:   cgeom->dimEmbed  = cdim;
1935:   cgeom->numCells  = 1;
1936:   cgeom->numPoints = Nq;
1937:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
1938:   return(0);
1939: }

1941: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
1942: {

1946:   PetscFree(cgeom->v);
1947:   PetscFree(cgeom->J);
1948:   PetscFree(cgeom->invJ);
1949:   PetscFree(cgeom->detJ);
1950:   return(0);
1951: }