Actual source code: fe.c

petsc-master 2020-10-23
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: PetscLogEvent PETSCFE_SetUp;

 55: PetscFunctionList PetscFEList              = NULL;
 56: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

 58: /*@C
 59:   PetscFERegister - Adds a new PetscFE implementation

 61:   Not Collective

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

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

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

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

 85:   Level: advanced

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

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

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

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

102:   Collective on fem

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

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

111:   Level: intermediate

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: .seealso: PetscFESetType(), PetscFECreate()
153: @*/
154: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
155: {

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

168: /*@C
169:    PetscFEViewFromOptions - View from Options

171:    Collective on PetscFE

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

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

187:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
188:   return(0);
189: }

191: /*@C
192:   PetscFEView - Views a PetscFE

194:   Collective on fem

196:   Input Parameter:
197: + fem - the PetscFE object to view
198: - viewer   - the viewer

200:   Level: beginner

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

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

219: /*@
220:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

222:   Collective on fem

224:   Input Parameter:
225: . fem - the PetscFE object to set options for

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

231:   Level: intermediate

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

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

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

270: /*@C
271:   PetscFESetUp - Construct data structures for the PetscFE

273:   Collective on fem

275:   Input Parameter:
276: . fem - the PetscFE object to setup

278:   Level: intermediate

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

288:   if (fem->setupcalled) return(0);
289:   PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
290:   fem->setupcalled = PETSC_TRUE;
291:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
292:   PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
293:   return(0);
294: }

296: /*@
297:   PetscFEDestroy - Destroys a PetscFE object

299:   Collective on fem

301:   Input Parameter:
302: . fem - the PetscFE object to destroy

304:   Level: beginner

306: .seealso PetscFEView()
307: @*/
308: PetscErrorCode PetscFEDestroy(PetscFE *fem)
309: {

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

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

319:   if ((*fem)->subspaces) {
320:     PetscInt dim, d;

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

335:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
336:   PetscHeaderDestroy(fem);
337:   return(0);
338: }

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

343:   Collective

345:   Input Parameter:
346: . comm - The communicator for the PetscFE object

348:   Output Parameter:
349: . fem - The PetscFE object

351:   Level: beginner

353: .seealso: PetscFESetType(), PETSCFEGALERKIN
354: @*/
355: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
356: {
357:   PetscFE        f;

362:   PetscCitationsRegister(FECitation,&FEcite);
363:   *fem = NULL;
364:   PetscFEInitializePackage();

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

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

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

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

390:   Not collective

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

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

398:   Level: intermediate

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

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

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

418:   Not collective

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

424:   Level: intermediate

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

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

439:   Not collective

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

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

447:   Level: intermediate

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

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

463:   Not collective

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

472:   Level: intermediate

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

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

490:   Not collective

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

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

501:   Level: intermediate

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

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

523:   Not collective

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

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

531:   Level: intermediate

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

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

547:   Not collective

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

553:   Level: intermediate

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

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

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

573:   Not collective

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

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

581:   Level: intermediate

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

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

597:   Not collective

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

603:   Level: intermediate

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

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

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

623:   Not collective

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

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

631:   Level: intermediate

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

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

647:   Not collective

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

653:   Level: intermediate

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

664:   if (q == fem->quadrature) return(0);
665:   PetscFEGetNumComponents(fem, &Nc);
666:   PetscQuadratureGetNumComponents(q, &qNc);
667:   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);
668:   PetscTabulationDestroy(&fem->T);
669:   PetscTabulationDestroy(&fem->Tc);
670:   PetscObjectReference((PetscObject) q);
671:   PetscQuadratureDestroy(&fem->quadrature);
672:   fem->quadrature = q;
673:   return(0);
674: }

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

679:   Not collective

681:   Input Parameter:
682: . fem - The PetscFE object

684:   Output Parameter:
685: . q - The PetscQuadrature object

687:   Level: intermediate

689: .seealso: PetscFECreate()
690: @*/
691: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
692: {
696:   *q = fem->faceQuadrature;
697:   return(0);
698: }

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

703:   Not collective

705:   Input Parameters:
706: + fem - The PetscFE object
707: - q - The PetscQuadrature object

709:   Level: intermediate

711: .seealso: PetscFECreate()
712: @*/
713: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
714: {
715:   PetscInt       Nc, qNc;

720:   PetscFEGetNumComponents(fem, &Nc);
721:   PetscQuadratureGetNumComponents(q, &qNc);
722:   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);
723:   PetscTabulationDestroy(&fem->Tf);
724:   PetscQuadratureDestroy(&fem->faceQuadrature);
725:   fem->faceQuadrature = q;
726:   PetscObjectReference((PetscObject) q);
727:   return(0);
728: }

730: /*@
731:   PetscFECopyQuadrature - Copy both volumetric and surface quadrature

733:   Not collective

735:   Input Parameters:
736: + sfe - The PetscFE source for the quadratures
737: - tfe - The PetscFE target for the quadratures

739:   Level: intermediate

741: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
742: @*/
743: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
744: {
745:   PetscQuadrature q;
746:   PetscErrorCode  ierr;

751:   PetscFEGetQuadrature(sfe, &q);
752:   PetscFESetQuadrature(tfe,  q);
753:   PetscFEGetFaceQuadrature(sfe, &q);
754:   PetscFESetFaceQuadrature(tfe,  q);
755:   return(0);
756: }

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

761:   Not collective

763:   Input Parameter:
764: . fem - The PetscFE object

766:   Output Parameter:
767: . numDof - Array with the number of dofs per dimension

769:   Level: intermediate

771: .seealso: PetscFECreate()
772: @*/
773: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
774: {

780:   PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
781:   return(0);
782: }

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

787:   Not collective

789:   Input Parameter:
790: . fem - The PetscFE object

792:   Output Parameter:
793: . T - The basis function values and derivatives at quadrature points

795:   Note:
796: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
797: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
798: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

800:   Level: intermediate

802: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
803: @*/
804: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T)
805: {
806:   PetscInt         npoints;
807:   const PetscReal *points;
808:   PetscErrorCode   ierr;

813:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
814:   if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);}
815:   *T = fem->T;
816:   return(0);
817: }

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

822:   Not collective

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

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

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

835:   Level: intermediate

837: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
838: @*/
839: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf)
840: {
841:   PetscErrorCode   ierr;

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

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

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

881:   Not collective

883:   Input Parameter:
884: . fem - The PetscFE object

886:   Output Parameters:
887: . Tc - The basis function values at face centroid points

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

892:   Level: intermediate

894: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
895: @*/
896: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
897: {
898:   PetscErrorCode   ierr;

903:   if (!fem->Tc) {
904:     PetscDualSpace  sp;
905:     DM              dm;
906:     const PetscInt *cone;
907:     PetscReal      *centroids;
908:     PetscInt        dim, numFaces, f;

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

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

927:   Not collective

929:   Input Parameters:
930: + fem     - The PetscFE object
931: . nrepl   - The number of replicas
932: . npoints - The number of tabulation points in a replica
933: . points  - The tabulation point coordinates
934: - K       - The number of derivatives calculated

936:   Output Parameter:
937: . T - The basis function values and derivatives at tabulation points

939:   Note:
940: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
941: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
942: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

944:   Level: intermediate

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

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

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

989:   Not collective

991:   Input Parameters:
992: + fem     - The PetscFE object
993: . npoints - The number of tabulation points
994: . points  - The tabulation point coordinates
995: . K       - The number of derivatives calculated
996: - T       - An existing tabulation object with enough allocated space

998:   Output Parameter:
999: . T - The basis function values and derivatives at tabulation points

1001:   Note:
1002: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1003: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1004: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

1006:   Level: intermediate

1008: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1009: @*/
1010: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1011: {

1015:   if (!npoints || !fem->dualSpace || K < 0) return(0);
1019:   if (PetscDefined(USE_DEBUG)) {
1020:     DM               dm;
1021:     PetscDualSpace   Q;
1022:     PetscInt         Nb;   /* Dimension of FE space P */
1023:     PetscInt         Nc;   /* Field components */
1024:     PetscInt         cdim; /* Reference coordinate dimension */

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

1042: /*@C
1043:   PetscTabulationDestroy - Frees memory from the associated tabulation.

1045:   Not collective

1047:   Input Parameter:
1048: . T - The tabulation

1050:   Level: intermediate

1052: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1053: @*/
1054: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1055: {
1056:   PetscInt       k;

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

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

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

1132: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1133: {
1134:   PetscInt       hStart, hEnd;
1135:   PetscDualSpace dsp;
1136:   DM             dm;

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


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

1155:   Not collective

1157:   Input Parameter:
1158: . fe - The PetscFE

1160:   Output Parameter:
1161: . dim - The dimension

1163:   Level: intermediate

1165: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1166: @*/
1167: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1168: {

1174:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1175:   return(0);
1176: }

1178: /*@C
1179:   PetscFEPushforward - Map the reference element function to real space

1181:   Input Parameters:
1182: + fe     - The PetscFE
1183: . fegeom - The cell geometry
1184: . Nv     - The number of function values
1185: - vals   - The function values

1187:   Output Parameter:
1188: . vals   - The transformed function values

1190:   Level: advanced

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

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

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

1203:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1204:   return(0);
1205: }

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

1210:   Input Parameters:
1211: + fe     - The PetscFE
1212: . fegeom - The cell geometry
1213: . Nv     - The number of function gradient values
1214: - vals   - The function gradient values

1216:   Output Parameter:
1217: . vals   - The transformed function gradient values

1219:   Level: advanced

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

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

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

1232:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1233:   return(0);
1234: }

1236: /*
1237: Purpose: Compute element vector for chunk of elements

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

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

1267:   Problem:
1268:      PetscInt f: the active field
1269:      f0, f1

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

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

1282: Input:
1283:   Sizes:
1284:      N_cb: Number of serial cell batches

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

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

1305: ex52.c:
1306:   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)
1307:   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)

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

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

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

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

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

1327:   Not collective

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

1338:   Output Parameter:
1339: . integral     - the integral for this field

1341:   Level: intermediate

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

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

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

1361:   Not collective

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

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

1376:   Level: intermediate

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

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

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

1400:   Not collective

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

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

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

1424:   Level: intermediate

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

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

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

1444:   Not collective

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

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

1460:   Level: intermediate

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

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

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

1480:   Not collective

1482:   Input Parameters:
1483: + prob         - The PetscDS specifying the discretizations and continuum functions
1484: . field        - The field being integrated
1485: . Ne           - The number of elements in the chunk
1486: . fgeom        - The face geometry for each cell in the chunk
1487: . coefficients - The array of FEM basis coefficients for the elements
1488: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1489: . probAux      - The PetscDS specifying the auxiliary discretizations
1490: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1491: - t            - The time

1493:   Output Parameter
1494: . elemVec      - the element residual vectors from each element

1496:   Level: developer

1498: .seealso: PetscFEIntegrateResidual()
1499: @*/
1500: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1501:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1502: {
1503:   PetscFE        fe;

1508:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1509:   if (fe->ops->integratehybridresidual) {(*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1510:   return(0);
1511: }

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

1516:   Not collective

1518:   Input Parameters:
1519: + prob         - The PetscDS specifying the discretizations and continuum functions
1520: . jtype        - The type of matrix pointwise functions that should be used
1521: . fieldI       - The test field being integrated
1522: . fieldJ       - The basis field being integrated
1523: . Ne           - The number of elements in the chunk
1524: . cgeom        - The cell geometry for each cell in the chunk
1525: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1526: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1527: . probAux      - The PetscDS specifying the auxiliary discretizations
1528: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1529: . t            - The time
1530: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1532:   Output Parameter:
1533: . elemMat      - the element matrices for the Jacobian from each element

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

1546: .seealso: PetscFEIntegrateResidual()
1547: @*/
1548: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1549:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1550: {
1551:   PetscFE        fe;

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

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

1564:   Not collective

1566:   Input Parameters:
1567: + prob         - The PetscDS specifying the discretizations and continuum functions
1568: . fieldI       - The test field being integrated
1569: . fieldJ       - The basis field being integrated
1570: . Ne           - The number of elements in the chunk
1571: . fgeom        - The face geometry for each cell in the chunk
1572: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1573: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1574: . probAux      - The PetscDS specifying the auxiliary discretizations
1575: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1576: . t            - The time
1577: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1579:   Output Parameter:
1580: . elemMat              - the element matrices for the Jacobian from each element

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

1593: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1594: @*/
1595: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1596:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1597: {
1598:   PetscFE        fe;

1603:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1604:   if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1605:   return(0);
1606: }

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

1611:   Not collective

1613:   Input Parameters:
1614: + prob         - The PetscDS specifying the discretizations and continuum functions
1615: . jtype        - The type of matrix pointwise functions that should be used
1616: . fieldI       - The test field being integrated
1617: . fieldJ       - The basis field being integrated
1618: . Ne           - The number of elements in the chunk
1619: . fgeom        - The face geometry for each cell in the chunk
1620: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1621: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1622: . probAux      - The PetscDS specifying the auxiliary discretizations
1623: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1624: . t            - The time
1625: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1627:   Output Parameter
1628: . elemMat              - the element matrices for the Jacobian from each element

1630:   Note:
1631: $ Loop over batch of elements (e):
1632: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1633: $     Loop over quadrature points (q):
1634: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1635: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1636: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1637: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1638: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1639:   Level: developer

1641: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1642: @*/
1643: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1644:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1645: {
1646:   PetscFE        fe;

1651:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1652:   if (fe->ops->integratehybridjacobian) {(*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1653:   return(0);
1654: }

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

1659:   Input Parameters:
1660: + fe     - The finite element space
1661: - height - The height of the Plex point

1663:   Output Parameter:
1664: . subfe  - The subspace of this FE space

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

1668:   Level: advanced

1670: .seealso: PetscFECreateDefault()
1671: @*/
1672: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1673: {
1674:   PetscSpace      P, subP;
1675:   PetscDualSpace  Q, subQ;
1676:   PetscQuadrature subq;
1677:   PetscFEType     fetype;
1678:   PetscInt        dim, Nc;
1679:   PetscErrorCode  ierr;

1684:   if (height == 0) {
1685:     *subfe = fe;
1686:     return(0);
1687:   }
1688:   PetscFEGetBasisSpace(fe, &P);
1689:   PetscFEGetDualSpace(fe, &Q);
1690:   PetscFEGetNumComponents(fe, &Nc);
1691:   PetscFEGetFaceQuadrature(fe, &subq);
1692:   PetscDualSpaceGetDimension(Q, &dim);
1693:   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);
1694:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1695:   if (height <= dim) {
1696:     if (!fe->subspaces[height-1]) {
1697:       PetscFE     sub = NULL;
1698:       const char *name;

1700:       PetscSpaceGetHeightSubspace(P, height, &subP);
1701:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1702:       if (subQ) {
1703:         PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1704:         PetscObjectGetName((PetscObject) fe,  &name);
1705:         PetscObjectSetName((PetscObject) sub,  name);
1706:         PetscFEGetType(fe, &fetype);
1707:         PetscFESetType(sub, fetype);
1708:         PetscFESetBasisSpace(sub, subP);
1709:         PetscFESetDualSpace(sub, subQ);
1710:         PetscFESetNumComponents(sub, Nc);
1711:         PetscFESetUp(sub);
1712:         PetscFESetQuadrature(sub, subq);
1713:       }
1714:       fe->subspaces[height-1] = sub;
1715:     }
1716:     *subfe = fe->subspaces[height-1];
1717:   } else {
1718:     *subfe = NULL;
1719:   }
1720:   return(0);
1721: }

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

1728:   Collective on fem

1730:   Input Parameter:
1731: . fe - The initial PetscFE

1733:   Output Parameter:
1734: . feRef - The refined PetscFE

1736:   Level: advanced

1738: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1739: @*/
1740: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1741: {
1742:   PetscSpace       P, Pref;
1743:   PetscDualSpace   Q, Qref;
1744:   DM               K, Kref;
1745:   PetscQuadrature  q, qref;
1746:   const PetscReal *v0, *jac;
1747:   PetscInt         numComp, numSubelements;
1748:   PetscInt         cStart, cEnd, c;
1749:   PetscDualSpace  *cellSpaces;
1750:   PetscErrorCode   ierr;

1753:   PetscFEGetBasisSpace(fe, &P);
1754:   PetscFEGetDualSpace(fe, &Q);
1755:   PetscFEGetQuadrature(fe, &q);
1756:   PetscDualSpaceGetDM(Q, &K);
1757:   /* Create space */
1758:   PetscObjectReference((PetscObject) P);
1759:   Pref = P;
1760:   /* Create dual space */
1761:   PetscDualSpaceDuplicate(Q, &Qref);
1762:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1763:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1764:   PetscDualSpaceSetDM(Qref, Kref);
1765:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1766:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1767:   /* TODO: fix for non-uniform refinement */
1768:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1769:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1770:   PetscFree(cellSpaces);
1771:   DMDestroy(&Kref);
1772:   PetscDualSpaceSetUp(Qref);
1773:   /* Create element */
1774:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1775:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1776:   PetscFESetBasisSpace(*feRef, Pref);
1777:   PetscFESetDualSpace(*feRef, Qref);
1778:   PetscFEGetNumComponents(fe,    &numComp);
1779:   PetscFESetNumComponents(*feRef, numComp);
1780:   PetscFESetUp(*feRef);
1781:   PetscSpaceDestroy(&Pref);
1782:   PetscDualSpaceDestroy(&Qref);
1783:   /* Create quadrature */
1784:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1785:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1786:   PetscFESetQuadrature(*feRef, qref);
1787:   PetscQuadratureDestroy(&qref);
1788:   return(0);
1789: }

1791: /*@C
1792:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1794:   Collective

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

1804:   Output Parameter:
1805: . fem - The PetscFE object

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

1810:   Level: beginner

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

1825:   /* Create space */
1826:   PetscSpaceCreate(comm, &P);
1827:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1828:   PetscSpacePolynomialSetTensor(P, tensor);
1829:   PetscSpaceSetNumComponents(P, Nc);
1830:   PetscSpaceSetNumVariables(P, dim);
1831:   PetscSpaceSetFromOptions(P);
1832:   PetscSpaceSetUp(P);
1833:   PetscSpaceGetDegree(P, &order, NULL);
1834:   PetscSpacePolynomialGetTensor(P, &tensor);
1835:   /* Create dual space */
1836:   PetscDualSpaceCreate(comm, &Q);
1837:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1838:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1839:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1840:   PetscDualSpaceSetDM(Q, K);
1841:   DMDestroy(&K);
1842:   PetscDualSpaceSetNumComponents(Q, Nc);
1843:   PetscDualSpaceSetOrder(Q, order);
1844:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1845:   PetscDualSpaceSetFromOptions(Q);
1846:   PetscDualSpaceSetUp(Q);
1847:   /* Create element */
1848:   PetscFECreate(comm, fem);
1849:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1850:   PetscFESetBasisSpace(*fem, P);
1851:   PetscFESetDualSpace(*fem, Q);
1852:   PetscFESetNumComponents(*fem, Nc);
1853:   PetscFESetFromOptions(*fem);
1854:   PetscFESetUp(*fem);
1855:   PetscSpaceDestroy(&P);
1856:   PetscDualSpaceDestroy(&Q);
1857:   /* Create quadrature (with specified order if given) */
1858:   qorder = qorder >= 0 ? qorder : order;
1859:   PetscObjectOptionsBegin((PetscObject)*fem);
1860:   PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1861:   PetscOptionsEnd();
1862:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1863:   if (isSimplex) {
1864:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1865:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1866:   } else {
1867:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1868:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1869:   }
1870:   PetscFESetQuadrature(*fem, q);
1871:   PetscFESetFaceQuadrature(*fem, fq);
1872:   PetscQuadratureDestroy(&q);
1873:   PetscQuadratureDestroy(&fq);
1874:   return(0);
1875: }

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

1880:   Collective

1882:   Input Parameters:
1883: + comm      - The MPI comm
1884: . dim       - The spatial dimension
1885: . Nc        - The number of components
1886: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1887: . k         - The degree k of the space
1888: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1890:   Output Parameter:
1891: . fem       - The PetscFE object

1893:   Level: beginner

1895:   Notes:
1896:   For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k.

1898: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1899: @*/
1900: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1901: {
1902:   PetscQuadrature q, fq;
1903:   DM              K;
1904:   PetscSpace      P;
1905:   PetscDualSpace  Q;
1906:   PetscInt        quadPointsPerEdge;
1907:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1908:   char            name[64];
1909:   PetscErrorCode  ierr;

1912:   /* Create space */
1913:   PetscSpaceCreate(comm, &P);
1914:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1915:   PetscSpacePolynomialSetTensor(P, tensor);
1916:   PetscSpaceSetNumComponents(P, Nc);
1917:   PetscSpaceSetNumVariables(P, dim);
1918:   PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1919:   PetscSpaceSetUp(P);
1920:   /* Create dual space */
1921:   PetscDualSpaceCreate(comm, &Q);
1922:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1923:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1924:   PetscDualSpaceSetDM(Q, K);
1925:   DMDestroy(&K);
1926:   PetscDualSpaceSetNumComponents(Q, Nc);
1927:   PetscDualSpaceSetOrder(Q, k);
1928:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1929:   PetscDualSpaceSetUp(Q);
1930:   /* Create finite element */
1931:   PetscFECreate(comm, fem);
1932:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1933:   PetscObjectSetName((PetscObject) *fem, name);
1934:   PetscFESetType(*fem, PETSCFEBASIC);
1935:   PetscFESetBasisSpace(*fem, P);
1936:   PetscFESetDualSpace(*fem, Q);
1937:   PetscFESetNumComponents(*fem, Nc);
1938:   PetscFESetUp(*fem);
1939:   PetscSpaceDestroy(&P);
1940:   PetscDualSpaceDestroy(&Q);
1941:   /* Create quadrature (with specified order if given) */
1942:   qorder = qorder >= 0 ? qorder : k;
1943:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1944:   if (isSimplex) {
1945:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1946:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1947:   } else {
1948:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1949:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1950:   }
1951:   PetscFESetQuadrature(*fem, q);
1952:   PetscFESetFaceQuadrature(*fem, fq);
1953:   PetscQuadratureDestroy(&q);
1954:   PetscQuadratureDestroy(&fq);
1955:   /* Set finite element name */
1956:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1957:   PetscFESetName(*fem, name);
1958:   return(0);
1959: }

1961: /*@C
1962:   PetscFESetName - Names the FE and its subobjects

1964:   Not collective

1966:   Input Parameters:
1967: + fe   - The PetscFE
1968: - name - The name

1970:   Level: intermediate

1972: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1973: @*/
1974: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1975: {
1976:   PetscSpace     P;
1977:   PetscDualSpace Q;

1981:   PetscFEGetBasisSpace(fe, &P);
1982:   PetscFEGetDualSpace(fe, &Q);
1983:   PetscObjectSetName((PetscObject) fe, name);
1984:   PetscObjectSetName((PetscObject) P,  name);
1985:   PetscObjectSetName((PetscObject) Q,  name);
1986:   return(0);
1987: }

1989: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1990: {
1991:   PetscInt       dOffset = 0, fOffset = 0, f;

1994:   for (f = 0; f < Nf; ++f) {
1995:     PetscFE          fe;
1996:     const PetscInt   cdim = T[f]->cdim;
1997:     const PetscInt   Nq   = T[f]->Np;
1998:     const PetscInt   Nbf  = T[f]->Nb;
1999:     const PetscInt   Ncf  = T[f]->Nc;
2000:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2001:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2002:     PetscInt         b, c, d;

2004:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2005:     for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2006:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2007:     for (b = 0; b < Nbf; ++b) {
2008:       for (c = 0; c < Ncf; ++c) {
2009:         const PetscInt cidx = b*Ncf+c;

2011:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2012:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2013:       }
2014:     }
2015:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2016:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2017:     if (u_t) {
2018:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2019:       for (b = 0; b < Nbf; ++b) {
2020:         for (c = 0; c < Ncf; ++c) {
2021:           const PetscInt cidx = b*Ncf+c;

2023:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2024:         }
2025:       }
2026:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2027:     }
2028:     fOffset += Ncf;
2029:     dOffset += Nbf;
2030:   }
2031:   return 0;
2032: }

2034: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2035: {
2036:   PetscInt       dOffset = 0, fOffset = 0, g;

2039:   for (g = 0; g < 2*Nf-1; ++g) {
2040:     if (!T[g/2]) continue;
2041:     {
2042:     PetscFE          fe;
2043:     const PetscInt   f   = g/2;
2044:     const PetscInt   cdim = T[f]->cdim;
2045:     const PetscInt   Nq   = T[f]->Np;
2046:     const PetscInt   Nbf  = T[f]->Nb;
2047:     const PetscInt   Ncf  = T[f]->Nc;
2048:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2049:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2050:     PetscInt         b, c, d;

2052:     fe = (PetscFE) ds->disc[f];
2053:     for (c = 0; c < Ncf; ++c)      u[fOffset+c] = 0.0;
2054:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2055:     for (b = 0; b < Nbf; ++b) {
2056:       for (c = 0; c < Ncf; ++c) {
2057:         const PetscInt cidx = b*Ncf+c;

2059:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2060:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2061:       }
2062:     }
2063:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2064:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2065:     if (u_t) {
2066:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2067:       for (b = 0; b < Nbf; ++b) {
2068:         for (c = 0; c < Ncf; ++c) {
2069:           const PetscInt cidx = b*Ncf+c;

2071:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2072:         }
2073:       }
2074:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2075:     }
2076:     fOffset += Ncf;
2077:     dOffset += Nbf;
2078:   }
2079:   }
2080:   return 0;
2081: }

2083: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2084: {
2085:   PetscFE         fe;
2086:   PetscTabulation Tc;
2087:   PetscInt        b, c;
2088:   PetscErrorCode  ierr;

2090:   if (!prob) return 0;
2091:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2092:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2093:   {
2094:     const PetscReal *faceBasis = Tc->T[0];
2095:     const PetscInt   Nb        = Tc->Nb;
2096:     const PetscInt   Nc        = Tc->Nc;

2098:     for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2099:     for (b = 0; b < Nb; ++b) {
2100:       for (c = 0; c < Nc; ++c) {
2101:         u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2102:       }
2103:     }
2104:   }
2105:   return 0;
2106: }

2108: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2109: {
2110:   const PetscInt   dE       = T->cdim; /* fegeom->dimEmbed */
2111:   const PetscInt   Nq       = T->Np;
2112:   const PetscInt   Nb       = T->Nb;
2113:   const PetscInt   Nc       = T->Nc;
2114:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2115:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2116:   PetscInt         q, b, c, d;
2117:   PetscErrorCode   ierr;

2119:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2120:   for (q = 0; q < Nq; ++q) {
2121:     for (b = 0; b < Nb; ++b) {
2122:       for (c = 0; c < Nc; ++c) {
2123:         const PetscInt bcidx = b*Nc+c;

2125:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2126:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2127:       }
2128:     }
2129:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2130:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2131:     for (b = 0; b < Nb; ++b) {
2132:       for (c = 0; c < Nc; ++c) {
2133:         const PetscInt bcidx = b*Nc+c;
2134:         const PetscInt qcidx = q*Nc+c;

2136:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2137:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2138:       }
2139:     }
2140:   }
2141:   return(0);
2142: }

2144: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2145: {
2146:   const PetscInt   dE       = T->cdim;
2147:   const PetscInt   Nq       = T->Np;
2148:   const PetscInt   Nb       = T->Nb;
2149:   const PetscInt   Nc       = T->Nc;
2150:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2151:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2152:   PetscInt         q, b, c, d, s;
2153:   PetscErrorCode   ierr;

2155:   for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0;
2156:   for (q = 0; q < Nq; ++q) {
2157:     for (b = 0; b < Nb; ++b) {
2158:       for (c = 0; c < Nc; ++c) {
2159:         const PetscInt bcidx = b*Nc+c;

2161:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2162:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2163:       }
2164:     }
2165:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2166:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2167:     for (s = 0; s < 2; ++s) {
2168:       for (b = 0; b < Nb; ++b) {
2169:         for (c = 0; c < Nc; ++c) {
2170:           const PetscInt bcidx = b*Nc+c;
2171:           const PetscInt qcidx = (q*2+s)*Nc+c;

2173:           elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2174:           for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2175:         }
2176:       }
2177:     }
2178:   }
2179:   return(0);
2180: }

2182: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2183: {
2184:   const PetscInt   dE        = TI->cdim;
2185:   const PetscInt   NqI       = TI->Np;
2186:   const PetscInt   NbI       = TI->Nb;
2187:   const PetscInt   NcI       = TI->Nc;
2188:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2189:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2190:   const PetscInt   NqJ       = TJ->Np;
2191:   const PetscInt   NbJ       = TJ->Nb;
2192:   const PetscInt   NcJ       = TJ->Nc;
2193:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2194:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2195:   PetscInt         f, fc, g, gc, df, dg;
2196:   PetscErrorCode   ierr;

2198:   for (f = 0; f < NbI; ++f) {
2199:     for (fc = 0; fc < NcI; ++fc) {
2200:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2202:       tmpBasisI[fidx] = basisI[fidx];
2203:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2204:     }
2205:   }
2206:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2207:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2208:   for (g = 0; g < NbJ; ++g) {
2209:     for (gc = 0; gc < NcJ; ++gc) {
2210:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2212:       tmpBasisJ[gidx] = basisJ[gidx];
2213:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2214:     }
2215:   }
2216:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2217:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2218:   for (f = 0; f < NbI; ++f) {
2219:     for (fc = 0; fc < NcI; ++fc) {
2220:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2221:       const PetscInt i    = offsetI+f; /* Element matrix row */
2222:       for (g = 0; g < NbJ; ++g) {
2223:         for (gc = 0; gc < NcJ; ++gc) {
2224:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2225:           const PetscInt j    = offsetJ+g; /* Element matrix column */
2226:           const PetscInt fOff = eOffset+i*totDim+j;

2228:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2229:           for (df = 0; df < dE; ++df) {
2230:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2231:             elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2232:             for (dg = 0; dg < dE; ++dg) {
2233:               elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2234:             }
2235:           }
2236:         }
2237:       }
2238:     }
2239:   }
2240:   return(0);
2241: }

2243: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2244: {
2245:   const PetscInt   dE        = TI->cdim;
2246:   const PetscInt   NqI       = TI->Np;
2247:   const PetscInt   NbI       = TI->Nb;
2248:   const PetscInt   NcI       = TI->Nc;
2249:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2250:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2251:   const PetscInt   NqJ       = TJ->Np;
2252:   const PetscInt   NbJ       = TJ->Nb;
2253:   const PetscInt   NcJ       = TJ->Nc;
2254:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2255:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2256:   const PetscInt   Ns = isHybridI ? 1 : 2;
2257:   const PetscInt   Nt = isHybridJ ? 1 : 2;
2258:   PetscInt         f, fc, g, gc, df, dg, s, t;
2259:   PetscErrorCode   ierr;

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

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

2275:       tmpBasisJ[gidx] = basisJ[gidx];
2276:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2277:     }
2278:   }
2279:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2280:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2281:   for (s = 0; s < Ns; ++s) {
2282:     for (f = 0; f < NbI; ++f) {
2283:       for (fc = 0; fc < NcI; ++fc) {
2284:         const PetscInt sc   = NcI*s+fc;  /* components from each side of the surface */
2285:         const PetscInt fidx = f*NcI+fc;  /* Test function basis index */
2286:         const PetscInt i    = offsetI+NbI*s+f; /* Element matrix row */
2287:         for (t = 0; t < Nt; ++t) {
2288:           for (g = 0; g < NbJ; ++g) {
2289:             for (gc = 0; gc < NcJ; ++gc) {
2290:               const PetscInt tc   = NcJ*t+gc;  /* components from each side of the surface */
2291:               const PetscInt gidx = g*NcJ+gc;  /* Trial function basis index */
2292:               const PetscInt j    = offsetJ+NbJ*t+g; /* Element matrix column */
2293:               const PetscInt fOff = eOffset+i*totDim+j;

2295:               elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2296:               for (df = 0; df < dE; ++df) {
2297:                 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2298:                 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2299:                 for (dg = 0; dg < dE; ++dg) {
2300:                   elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2301:                 }
2302:               }
2303:             }
2304:           }
2305:         }
2306:       }
2307:     }
2308:   }
2309:   return(0);
2310: }

2312: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2313: {
2314:   PetscDualSpace  dsp;
2315:   DM              dm;
2316:   PetscQuadrature quadDef;
2317:   PetscInt        dim, cdim, Nq;
2318:   PetscErrorCode  ierr;

2321:   PetscFEGetDualSpace(fe, &dsp);
2322:   PetscDualSpaceGetDM(dsp, &dm);
2323:   DMGetDimension(dm, &dim);
2324:   DMGetCoordinateDim(dm, &cdim);
2325:   PetscFEGetQuadrature(fe, &quadDef);
2326:   quad = quad ? quad : quadDef;
2327:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2328:   PetscMalloc1(Nq*cdim,      &cgeom->v);
2329:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2330:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2331:   PetscMalloc1(Nq,           &cgeom->detJ);
2332:   cgeom->dim       = dim;
2333:   cgeom->dimEmbed  = cdim;
2334:   cgeom->numCells  = 1;
2335:   cgeom->numPoints = Nq;
2336:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2337:   return(0);
2338: }

2340: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2341: {

2345:   PetscFree(cgeom->v);
2346:   PetscFree(cgeom->J);
2347:   PetscFree(cgeom->invJ);
2348:   PetscFree(cgeom->detJ);
2349:   return(0);
2350: }