Actual source code: fe.c

petsc-3.15.0 2021-04-05
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
791: - k   - The highest derivative we need to tabulate, very often 1

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

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

801:   Level: intermediate

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

814:   PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
815:   if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);}
816:   if (fem->T && k > fem->T->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->T->K);
817:   *T = fem->T;
818:   return(0);
819: }

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

824:   Not collective

826:   Input Parameter:
827: + fem - The PetscFE object
828: - k   - The highest derivative we need to tabulate, very often 1

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

833:   Note:
834: $ 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
835: $ 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
836: $ 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

838:   Level: intermediate

840: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
841: @*/
842: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
843: {
844:   PetscErrorCode   ierr;

849:   if (!fem->Tf) {
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:       PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
874:       PetscFree(facePoints);
875:     }
876:   }
877:   if (fem->Tf && k > fem->Tf->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->Tf->K);
878:   *Tf = fem->Tf;
879:   return(0);
880: }

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

885:   Not collective

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

890:   Output Parameters:
891: . Tc - The basis function values at face centroid points

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

896:   Level: intermediate

898: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
899: @*/
900: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
901: {
902:   PetscErrorCode   ierr;

907:   if (!fem->Tc) {
908:     PetscDualSpace  sp;
909:     DM              dm;
910:     const PetscInt *cone;
911:     PetscReal      *centroids;
912:     PetscInt        dim, numFaces, f;

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

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

931:   Not collective

933:   Input Parameters:
934: + fem     - The PetscFE object
935: . nrepl   - The number of replicas
936: . npoints - The number of tabulation points in a replica
937: . points  - The tabulation point coordinates
938: - K       - The number of derivatives calculated

940:   Output Parameter:
941: . T - The basis function values and derivatives at tabulation points

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

948:   Level: intermediate

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

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

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

993:   Not collective

995:   Input Parameters:
996: + fem     - The PetscFE object
997: . npoints - The number of tabulation points
998: . points  - The tabulation point coordinates
999: . K       - The number of derivatives calculated
1000: - T       - An existing tabulation object with enough allocated space

1002:   Output Parameter:
1003: . T - The basis function values and derivatives at tabulation points

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

1010:   Level: intermediate

1012: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1013: @*/
1014: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1015: {

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

1030:     PetscFEGetDualSpace(fem, &Q);
1031:     PetscDualSpaceGetDM(Q, &dm);
1032:     DMGetDimension(dm, &cdim);
1033:     PetscDualSpaceGetDimension(Q, &Nb);
1034:     PetscFEGetNumComponents(fem, &Nc);
1035:     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);
1036:     if (T->Nb   != Nb)              SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1037:     if (T->Nc   != Nc)              SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1038:     if (T->cdim != cdim)            SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1039:   }
1040:   T->Nr = 1;
1041:   T->Np = npoints;
1042:   (*fem->ops->createtabulation)(fem, npoints, points, K, T);
1043:   return(0);
1044: }

1046: /*@C
1047:   PetscTabulationDestroy - Frees memory from the associated tabulation.

1049:   Not collective

1051:   Input Parameter:
1052: . T - The tabulation

1054:   Level: intermediate

1056: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1057: @*/
1058: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1059: {
1060:   PetscInt       k;

1065:   if (!T || !(*T)) return(0);
1066:   for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);}
1067:   PetscFree((*T)->T);
1068:   PetscFree(*T);
1069:   *T = NULL;
1070:   return(0);
1071: }

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

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

1136: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1137: {
1138:   PetscInt       hStart, hEnd;
1139:   PetscDualSpace dsp;
1140:   DM             dm;

1146:   *trFE = NULL;
1147:   PetscFEGetDualSpace(fe,&dsp);
1148:   PetscDualSpaceGetDM(dsp,&dm);
1149:   DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1150:   if (hEnd <= hStart) return(0);
1151:   PetscFECreatePointTrace(fe,hStart,trFE);
1152:   return(0);
1153: }


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

1159:   Not collective

1161:   Input Parameter:
1162: . fe - The PetscFE

1164:   Output Parameter:
1165: . dim - The dimension

1167:   Level: intermediate

1169: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1170: @*/
1171: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1172: {

1178:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1179:   return(0);
1180: }

1182: /*@C
1183:   PetscFEPushforward - Map the reference element function to real space

1185:   Input Parameters:
1186: + fe     - The PetscFE
1187: . fegeom - The cell geometry
1188: . Nv     - The number of function values
1189: - vals   - The function values

1191:   Output Parameter:
1192: . vals   - The transformed function values

1194:   Level: advanced

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

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

1200: .seealso: PetscDualSpacePushforward()
1201: @*/
1202: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1203: {

1207:   PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1208:   return(0);
1209: }

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

1214:   Input Parameters:
1215: + fe     - The PetscFE
1216: . fegeom - The cell geometry
1217: . Nv     - The number of function gradient values
1218: - vals   - The function gradient values

1220:   Output Parameter:
1221: . vals   - The transformed function gradient values

1223:   Level: advanced

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

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

1229: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1230: @*/
1231: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1232: {

1236:   PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1237:   return(0);
1238: }

1240: /*@C
1241:   PetscFEPushforwardHessian - Map the reference element function Hessian to real space

1243:   Input Parameters:
1244: + fe     - The PetscFE
1245: . fegeom - The cell geometry
1246: . Nv     - The number of function Hessian values
1247: - vals   - The function Hessian values

1249:   Output Parameter:
1250: . vals   - The transformed function Hessian values

1252:   Level: advanced

1254:   Note: This just forwards the call onto PetscDualSpacePushforwardHessian().

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

1258: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward()
1259: @*/
1260: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1261: {

1265:   PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1266:   return(0);
1267: }

1269: /*
1270: Purpose: Compute element vector for chunk of elements

1272: Input:
1273:   Sizes:
1274:      Ne:  number of elements
1275:      Nf:  number of fields
1276:      PetscFE
1277:        dim: spatial dimension
1278:        Nb:  number of basis functions
1279:        Nc:  number of field components
1280:        PetscQuadrature
1281:          Nq:  number of quadrature points

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

1300:   Problem:
1301:      PetscInt f: the active field
1302:      f0, f1

1304:   Work Space:
1305:      PetscFE
1306:        PetscScalar f0[Nq*dim];
1307:        PetscScalar f1[Nq*dim*dim];
1308:        PetscScalar u[Nc];
1309:        PetscScalar gradU[Nc*dim];
1310:        PetscReal   x[dim];
1311:        PetscScalar realSpaceDer[dim];

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

1315: Input:
1316:   Sizes:
1317:      N_cb: Number of serial cell batches

1319:   Geometry:
1320:      PetscReal v0s[Ne*dim]
1321:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
1322:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1323:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
1324:   FEM:
1325:      static PetscReal   quadPoints[Nq*dim]
1326:      static PetscReal   quadWeights[Nq]
1327:      static PetscReal   basis[Nq*Nb*Nc]
1328:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
1329:      PetscScalar coefficients[Ne*Nb*Nc]
1330:      PetscScalar elemVec[Ne*Nb*Nc]

1332: ex62.c:
1333:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1334:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1335:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1336:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

1338: ex52.c:
1339:   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)
1340:   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)

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

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

1349: ex52_integrateElementOpenCL.c:
1350: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1351:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1352:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

1360:   Not collective

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

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

1374:   Level: intermediate

1376: .seealso: PetscFEIntegrateResidual()
1377: @*/
1378: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1379:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1380: {
1381:   PetscFE        fe;

1386:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1387:   if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1388:   return(0);
1389: }

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

1394:   Not collective

1396:   Input Parameters:
1397: + prob         - The PetscDS specifying the discretizations and continuum functions
1398: . field        - The field being integrated
1399: . obj_func     - The function to be integrated
1400: . Ne           - The number of elements in the chunk
1401: . fgeom        - The face geometry for each face in the chunk
1402: . coefficients - The array of FEM basis coefficients for the elements
1403: . probAux      - The PetscDS specifying the auxiliary discretizations
1404: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

1406:   Output Parameter:
1407: . integral     - the integral for this field

1409:   Level: intermediate

1411: .seealso: PetscFEIntegrateResidual()
1412: @*/
1413: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1414:                                   void (*obj_func)(PetscInt, PetscInt, PetscInt,
1415:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1416:                                                    const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1417:                                                    PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1418:                                   PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1419: {
1420:   PetscFE        fe;

1425:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1426:   if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1427:   return(0);
1428: }

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

1433:   Not collective

1435:   Input Parameters:
1436: + ds           - The PetscDS specifying the discretizations and continuum functions
1437: . key          - The (label+value, field) being integrated
1438: . Ne           - The number of elements in the chunk
1439: . cgeom        - The cell geometry for each cell in the chunk
1440: . coefficients - The array of FEM basis coefficients for the elements
1441: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1442: . probAux      - The PetscDS specifying the auxiliary discretizations
1443: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1444: - t            - The time

1446:   Output Parameter:
1447: . elemVec      - the element residual vectors from each element

1449:   Note:
1450: $ Loop over batch of elements (e):
1451: $   Loop over quadrature points (q):
1452: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1453: $     Call f_0 and f_1
1454: $   Loop over element vector entries (f,fc --> i):
1455: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

1457:   Level: intermediate

1459: .seealso: PetscFEIntegrateResidual()
1460: @*/
1461: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1462:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1463: {
1464:   PetscFE        fe;

1469:   PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);
1470:   if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1471:   return(0);
1472: }

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

1477:   Not collective

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

1490:   Output Parameter:
1491: . elemVec      - the element residual vectors from each element

1493:   Level: intermediate

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

1505:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1506:   if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1507:   return(0);
1508: }

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

1513:   Not collective

1515:   Input Parameters:
1516: + prob         - The PetscDS specifying the discretizations and continuum functions
1517: . key          - The (label+value, field) being integrated
1518: . Ne           - The number of elements in the chunk
1519: . fgeom        - The face geometry for each cell in the chunk
1520: . coefficients - The array of FEM basis coefficients for the elements
1521: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1522: . probAux      - The PetscDS specifying the auxiliary discretizations
1523: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1524: - t            - The time

1526:   Output Parameter
1527: . elemVec      - the element residual vectors from each element

1529:   Level: developer

1531: .seealso: PetscFEIntegrateResidual()
1532: @*/
1533: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1534:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1535: {
1536:   PetscFE        fe;

1541:   PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe);
1542:   if (fe->ops->integratehybridresidual) {(*fe->ops->integratehybridresidual)(prob, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1543:   return(0);
1544: }

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

1549:   Not collective

1551:   Input Parameters:
1552: + ds           - The PetscDS specifying the discretizations and continuum functions
1553: . jtype        - The type of matrix pointwise functions that should be used
1554: . key          - The (label+value, fieldI*Nf + fieldJ) being integrated
1555: . Ne           - The number of elements in the chunk
1556: . cgeom        - The cell geometry for each cell in the chunk
1557: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1558: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1559: . probAux      - The PetscDS specifying the auxiliary discretizations
1560: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1561: . t            - The time
1562: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1564:   Output Parameter:
1565: . elemMat      - the element matrices for the Jacobian from each element

1567:   Note:
1568: $ Loop over batch of elements (e):
1569: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1570: $     Loop over quadrature points (q):
1571: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1572: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1573: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1574: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1575: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1576:   Level: intermediate

1578: .seealso: PetscFEIntegrateResidual()
1579: @*/
1580: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1581:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1582: {
1583:   PetscFE        fe;
1584:   PetscInt       Nf;

1589:   PetscDSGetNumFields(ds, &Nf);
1590:   PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1591:   if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1592:   return(0);
1593: }

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

1598:   Not collective

1600:   Input Parameters:
1601: + prob         - The PetscDS specifying the discretizations and continuum functions
1602: . fieldI       - The test field being integrated
1603: . fieldJ       - The basis field being integrated
1604: . Ne           - The number of elements in the chunk
1605: . fgeom        - The face geometry for each cell in the chunk
1606: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1607: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1608: . probAux      - The PetscDS specifying the auxiliary discretizations
1609: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1610: . t            - The time
1611: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1613:   Output Parameter:
1614: . elemMat              - the element matrices for the Jacobian from each element

1616:   Note:
1617: $ Loop over batch of elements (e):
1618: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1619: $     Loop over quadrature points (q):
1620: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1621: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1622: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1623: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1624: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1625:   Level: intermediate

1627: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1628: @*/
1629: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1630:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1631: {
1632:   PetscFE        fe;

1637:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1638:   if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1639:   return(0);
1640: }

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

1645:   Not collective

1647:   Input Parameters:
1648: + prob         - The PetscDS specifying the discretizations and continuum functions
1649: . jtype        - The type of matrix pointwise functions that should be used
1650: . fieldI       - The test field being integrated
1651: . fieldJ       - The basis field being integrated
1652: . Ne           - The number of elements in the chunk
1653: . fgeom        - The face geometry for each cell in the chunk
1654: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1655: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1656: . probAux      - The PetscDS specifying the auxiliary discretizations
1657: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1658: . t            - The time
1659: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

1661:   Output Parameter
1662: . elemMat              - the element matrices for the Jacobian from each element

1664:   Note:
1665: $ Loop over batch of elements (e):
1666: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
1667: $     Loop over quadrature points (q):
1668: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1669: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1670: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1671: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1672: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1673:   Level: developer

1675: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1676: @*/
1677: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1678:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1679: {
1680:   PetscFE        fe;

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

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

1693:   Input Parameters:
1694: + fe     - The finite element space
1695: - height - The height of the Plex point

1697:   Output Parameter:
1698: . subfe  - The subspace of this FE space

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

1702:   Level: advanced

1704: .seealso: PetscFECreateDefault()
1705: @*/
1706: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1707: {
1708:   PetscSpace      P, subP;
1709:   PetscDualSpace  Q, subQ;
1710:   PetscQuadrature subq;
1711:   PetscFEType     fetype;
1712:   PetscInt        dim, Nc;
1713:   PetscErrorCode  ierr;

1718:   if (height == 0) {
1719:     *subfe = fe;
1720:     return(0);
1721:   }
1722:   PetscFEGetBasisSpace(fe, &P);
1723:   PetscFEGetDualSpace(fe, &Q);
1724:   PetscFEGetNumComponents(fe, &Nc);
1725:   PetscFEGetFaceQuadrature(fe, &subq);
1726:   PetscDualSpaceGetDimension(Q, &dim);
1727:   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);
1728:   if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1729:   if (height <= dim) {
1730:     if (!fe->subspaces[height-1]) {
1731:       PetscFE     sub = NULL;
1732:       const char *name;

1734:       PetscSpaceGetHeightSubspace(P, height, &subP);
1735:       PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1736:       if (subQ) {
1737:         PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1738:         PetscObjectGetName((PetscObject) fe,  &name);
1739:         PetscObjectSetName((PetscObject) sub,  name);
1740:         PetscFEGetType(fe, &fetype);
1741:         PetscFESetType(sub, fetype);
1742:         PetscFESetBasisSpace(sub, subP);
1743:         PetscFESetDualSpace(sub, subQ);
1744:         PetscFESetNumComponents(sub, Nc);
1745:         PetscFESetUp(sub);
1746:         PetscFESetQuadrature(sub, subq);
1747:       }
1748:       fe->subspaces[height-1] = sub;
1749:     }
1750:     *subfe = fe->subspaces[height-1];
1751:   } else {
1752:     *subfe = NULL;
1753:   }
1754:   return(0);
1755: }

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

1762:   Collective on fem

1764:   Input Parameter:
1765: . fe - The initial PetscFE

1767:   Output Parameter:
1768: . feRef - The refined PetscFE

1770:   Level: advanced

1772: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1773: @*/
1774: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1775: {
1776:   PetscSpace       P, Pref;
1777:   PetscDualSpace   Q, Qref;
1778:   DM               K, Kref;
1779:   PetscQuadrature  q, qref;
1780:   const PetscReal *v0, *jac;
1781:   PetscInt         numComp, numSubelements;
1782:   PetscInt         cStart, cEnd, c;
1783:   PetscDualSpace  *cellSpaces;
1784:   PetscErrorCode   ierr;

1787:   PetscFEGetBasisSpace(fe, &P);
1788:   PetscFEGetDualSpace(fe, &Q);
1789:   PetscFEGetQuadrature(fe, &q);
1790:   PetscDualSpaceGetDM(Q, &K);
1791:   /* Create space */
1792:   PetscObjectReference((PetscObject) P);
1793:   Pref = P;
1794:   /* Create dual space */
1795:   PetscDualSpaceDuplicate(Q, &Qref);
1796:   PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1797:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1798:   PetscDualSpaceSetDM(Qref, Kref);
1799:   DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1800:   PetscMalloc1(cEnd - cStart, &cellSpaces);
1801:   /* TODO: fix for non-uniform refinement */
1802:   for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1803:   PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1804:   PetscFree(cellSpaces);
1805:   DMDestroy(&Kref);
1806:   PetscDualSpaceSetUp(Qref);
1807:   /* Create element */
1808:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1809:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
1810:   PetscFESetBasisSpace(*feRef, Pref);
1811:   PetscFESetDualSpace(*feRef, Qref);
1812:   PetscFEGetNumComponents(fe,    &numComp);
1813:   PetscFESetNumComponents(*feRef, numComp);
1814:   PetscFESetUp(*feRef);
1815:   PetscSpaceDestroy(&Pref);
1816:   PetscDualSpaceDestroy(&Qref);
1817:   /* Create quadrature */
1818:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1819:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1820:   PetscFESetQuadrature(*feRef, qref);
1821:   PetscQuadratureDestroy(&qref);
1822:   return(0);
1823: }

1825: /*@C
1826:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

1828:   Collective

1830:   Input Parameters:
1831: + comm      - The MPI comm
1832: . dim       - The spatial dimension
1833: . Nc        - The number of components
1834: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1835: . prefix    - The options prefix, or NULL
1836: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1838:   Output Parameter:
1839: . fem - The PetscFE object

1841:   Note:
1842:   Each subobject is SetFromOption() during creation, so that the object may be customized from the command line, using the prefix specified above. See the links below for the particular options available.

1844:   Level: beginner

1846: .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1847: @*/
1848: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1849: {
1850:   PetscQuadrature q, fq;
1851:   DM              K;
1852:   PetscSpace      P;
1853:   PetscDualSpace  Q;
1854:   PetscInt        order, quadPointsPerEdge;
1855:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1856:   PetscErrorCode  ierr;

1859:   /* Create space */
1860:   PetscSpaceCreate(comm, &P);
1861:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1862:   PetscSpacePolynomialSetTensor(P, tensor);
1863:   PetscSpaceSetNumComponents(P, Nc);
1864:   PetscSpaceSetNumVariables(P, dim);
1865:   PetscSpaceSetFromOptions(P);
1866:   PetscSpaceSetUp(P);
1867:   PetscSpaceGetDegree(P, &order, NULL);
1868:   PetscSpacePolynomialGetTensor(P, &tensor);
1869:   /* Create dual space */
1870:   PetscDualSpaceCreate(comm, &Q);
1871:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1872:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1873:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1874:   PetscDualSpaceSetDM(Q, K);
1875:   DMDestroy(&K);
1876:   PetscDualSpaceSetNumComponents(Q, Nc);
1877:   PetscDualSpaceSetOrder(Q, order);
1878:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1879:   PetscDualSpaceSetFromOptions(Q);
1880:   PetscDualSpaceSetUp(Q);
1881:   /* Create element */
1882:   PetscFECreate(comm, fem);
1883:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1884:   PetscFESetBasisSpace(*fem, P);
1885:   PetscFESetDualSpace(*fem, Q);
1886:   PetscFESetNumComponents(*fem, Nc);
1887:   PetscFESetFromOptions(*fem);
1888:   PetscFESetUp(*fem);
1889:   PetscSpaceDestroy(&P);
1890:   PetscDualSpaceDestroy(&Q);
1891:   /* Create quadrature (with specified order if given) */
1892:   qorder = qorder >= 0 ? qorder : order;
1893:   PetscObjectOptionsBegin((PetscObject)*fem);
1894:   PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1895:   PetscOptionsEnd();
1896:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1897:   if (isSimplex) {
1898:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1899:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1900:   } else {
1901:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1902:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1903:   }
1904:   PetscFESetQuadrature(*fem, q);
1905:   PetscFESetFaceQuadrature(*fem, fq);
1906:   PetscQuadratureDestroy(&q);
1907:   PetscQuadratureDestroy(&fq);
1908:   return(0);
1909: }

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

1914:   Collective

1916:   Input Parameters:
1917: + comm      - The MPI comm
1918: . dim       - The spatial dimension
1919: . Nc        - The number of components
1920: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1921: . k         - The degree k of the space
1922: - qorder    - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree

1924:   Output Parameter:
1925: . fem       - The PetscFE object

1927:   Level: beginner

1929:   Notes:
1930:   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.

1932: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1933: @*/
1934: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1935: {
1936:   PetscQuadrature q, fq;
1937:   DM              K;
1938:   PetscSpace      P;
1939:   PetscDualSpace  Q;
1940:   PetscInt        quadPointsPerEdge;
1941:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1942:   char            name[64];
1943:   PetscErrorCode  ierr;

1946:   /* Create space */
1947:   PetscSpaceCreate(comm, &P);
1948:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1949:   PetscSpacePolynomialSetTensor(P, tensor);
1950:   PetscSpaceSetNumComponents(P, Nc);
1951:   PetscSpaceSetNumVariables(P, dim);
1952:   PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1953:   PetscSpaceSetUp(P);
1954:   /* Create dual space */
1955:   PetscDualSpaceCreate(comm, &Q);
1956:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1957:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1958:   PetscDualSpaceSetDM(Q, K);
1959:   DMDestroy(&K);
1960:   PetscDualSpaceSetNumComponents(Q, Nc);
1961:   PetscDualSpaceSetOrder(Q, k);
1962:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
1963:   PetscDualSpaceSetUp(Q);
1964:   /* Create finite element */
1965:   PetscFECreate(comm, fem);
1966:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1967:   PetscObjectSetName((PetscObject) *fem, name);
1968:   PetscFESetType(*fem, PETSCFEBASIC);
1969:   PetscFESetBasisSpace(*fem, P);
1970:   PetscFESetDualSpace(*fem, Q);
1971:   PetscFESetNumComponents(*fem, Nc);
1972:   PetscFESetUp(*fem);
1973:   PetscSpaceDestroy(&P);
1974:   PetscDualSpaceDestroy(&Q);
1975:   /* Create quadrature (with specified order if given) */
1976:   qorder = qorder >= 0 ? qorder : k;
1977:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1978:   if (isSimplex) {
1979:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1980:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1981:   } else {
1982:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
1983:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1984:   }
1985:   PetscFESetQuadrature(*fem, q);
1986:   PetscFESetFaceQuadrature(*fem, fq);
1987:   PetscQuadratureDestroy(&q);
1988:   PetscQuadratureDestroy(&fq);
1989:   /* Set finite element name */
1990:   PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1991:   PetscFESetName(*fem, name);
1992:   return(0);
1993: }

1995: /*@C
1996:   PetscFESetName - Names the FE and its subobjects

1998:   Not collective

2000:   Input Parameters:
2001: + fe   - The PetscFE
2002: - name - The name

2004:   Level: intermediate

2006: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
2007: @*/
2008: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2009: {
2010:   PetscSpace     P;
2011:   PetscDualSpace Q;

2015:   PetscFEGetBasisSpace(fe, &P);
2016:   PetscFEGetDualSpace(fe, &Q);
2017:   PetscObjectSetName((PetscObject) fe, name);
2018:   PetscObjectSetName((PetscObject) P,  name);
2019:   PetscObjectSetName((PetscObject) Q,  name);
2020:   return(0);
2021: }

2023: 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[])
2024: {
2025:   PetscInt       dOffset = 0, fOffset = 0, f, g;

2028:   for (f = 0; f < Nf; ++f) {
2029:     PetscFE          fe;
2030:     const PetscInt   k    = ds->jetDegree[f];
2031:     const PetscInt   cdim = T[f]->cdim;
2032:     const PetscInt   Nq   = T[f]->Np;
2033:     const PetscInt   Nbf  = T[f]->Nb;
2034:     const PetscInt   Ncf  = T[f]->Nc;
2035:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2036:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2037:     const PetscReal *Hq   = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL;
2038:     PetscInt         hOffset = 0, b, c, d;

2040:     PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2041:     for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2042:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2043:     for (b = 0; b < Nbf; ++b) {
2044:       for (c = 0; c < Ncf; ++c) {
2045:         const PetscInt cidx = b*Ncf+c;

2047:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2048:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2049:       }
2050:     }
2051:     if (k > 1) {
2052:       for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim;
2053:       for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0;
2054:       for (b = 0; b < Nbf; ++b) {
2055:         for (c = 0; c < Ncf; ++c) {
2056:           const PetscInt cidx = b*Ncf+c;

2058:           for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b];
2059:         }
2060:       }
2061:       PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);
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:   return 0;
2080: }

2082: 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[])
2083: {
2084:   PetscInt       dOffset = 0, fOffset = 0, g;

2087:   for (g = 0; g < 2*Nf-1; ++g) {
2088:     if (!T[g/2]) continue;
2089:     {
2090:     PetscFE          fe;
2091:     const PetscInt   f   = g/2;
2092:     const PetscInt   cdim = T[f]->cdim;
2093:     const PetscInt   Nq   = T[f]->Np;
2094:     const PetscInt   Nbf  = T[f]->Nb;
2095:     const PetscInt   Ncf  = T[f]->Nc;
2096:     const PetscReal *Bq   = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2097:     const PetscReal *Dq   = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2098:     PetscInt         b, c, d;

2100:     fe = (PetscFE) ds->disc[f];
2101:     for (c = 0; c < Ncf; ++c)      u[fOffset+c] = 0.0;
2102:     for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2103:     for (b = 0; b < Nbf; ++b) {
2104:       for (c = 0; c < Ncf; ++c) {
2105:         const PetscInt cidx = b*Ncf+c;

2107:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2108:         for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2109:       }
2110:     }
2111:     PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2112:     PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2113:     if (u_t) {
2114:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2115:       for (b = 0; b < Nbf; ++b) {
2116:         for (c = 0; c < Ncf; ++c) {
2117:           const PetscInt cidx = b*Ncf+c;

2119:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2120:         }
2121:       }
2122:       PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2123:     }
2124:     fOffset += Ncf;
2125:     dOffset += Nbf;
2126:   }
2127:   }
2128:   return 0;
2129: }

2131: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2132: {
2133:   PetscFE         fe;
2134:   PetscTabulation Tc;
2135:   PetscInt        b, c;
2136:   PetscErrorCode  ierr;

2138:   if (!prob) return 0;
2139:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2140:   PetscFEGetFaceCentroidTabulation(fe, &Tc);
2141:   {
2142:     const PetscReal *faceBasis = Tc->T[0];
2143:     const PetscInt   Nb        = Tc->Nb;
2144:     const PetscInt   Nc        = Tc->Nc;

2146:     for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2147:     for (b = 0; b < Nb; ++b) {
2148:       for (c = 0; c < Nc; ++c) {
2149:         u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2150:       }
2151:     }
2152:   }
2153:   return 0;
2154: }

2156: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2157: {
2158:   const PetscInt   dE       = T->cdim; /* fegeom->dimEmbed */
2159:   const PetscInt   Nq       = T->Np;
2160:   const PetscInt   Nb       = T->Nb;
2161:   const PetscInt   Nc       = T->Nc;
2162:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2163:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2164:   PetscInt         q, b, c, d;
2165:   PetscErrorCode   ierr;

2167:   for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2168:   for (q = 0; q < Nq; ++q) {
2169:     for (b = 0; b < Nb; ++b) {
2170:       for (c = 0; c < Nc; ++c) {
2171:         const PetscInt bcidx = b*Nc+c;

2173:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2174:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2175:       }
2176:     }
2177:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2178:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2179:     for (b = 0; b < Nb; ++b) {
2180:       for (c = 0; c < Nc; ++c) {
2181:         const PetscInt bcidx = b*Nc+c;
2182:         const PetscInt qcidx = q*Nc+c;

2184:         elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2185:         for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2186:       }
2187:     }
2188:   }
2189:   return(0);
2190: }

2192: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2193: {
2194:   const PetscInt   dE       = T->cdim;
2195:   const PetscInt   Nq       = T->Np;
2196:   const PetscInt   Nb       = T->Nb;
2197:   const PetscInt   Nc       = T->Nc;
2198:   const PetscReal *basis    = &T->T[0][r*Nq*Nb*Nc];
2199:   const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2200:   PetscInt         q, b, c, d, s;
2201:   PetscErrorCode   ierr;

2203:   for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0;
2204:   for (q = 0; q < Nq; ++q) {
2205:     for (b = 0; b < Nb; ++b) {
2206:       for (c = 0; c < Nc; ++c) {
2207:         const PetscInt bcidx = b*Nc+c;

2209:         tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2210:         for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2211:       }
2212:     }
2213:     PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2214:     PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2215:     for (s = 0; s < 2; ++s) {
2216:       for (b = 0; b < Nb; ++b) {
2217:         for (c = 0; c < Nc; ++c) {
2218:           const PetscInt bcidx = b*Nc+c;
2219:           const PetscInt qcidx = (q*2+s)*Nc+c;

2221:           elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2222:           for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2223:         }
2224:       }
2225:     }
2226:   }
2227:   return(0);
2228: }

2230: 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[])
2231: {
2232:   const PetscInt   dE        = TI->cdim;
2233:   const PetscInt   NqI       = TI->Np;
2234:   const PetscInt   NbI       = TI->Nb;
2235:   const PetscInt   NcI       = TI->Nc;
2236:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2237:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2238:   const PetscInt   NqJ       = TJ->Np;
2239:   const PetscInt   NbJ       = TJ->Nb;
2240:   const PetscInt   NcJ       = TJ->Nc;
2241:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2242:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2243:   PetscInt         f, fc, g, gc, df, dg;
2244:   PetscErrorCode   ierr;

2246:   for (f = 0; f < NbI; ++f) {
2247:     for (fc = 0; fc < NcI; ++fc) {
2248:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2250:       tmpBasisI[fidx] = basisI[fidx];
2251:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2252:     }
2253:   }
2254:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2255:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2256:   for (g = 0; g < NbJ; ++g) {
2257:     for (gc = 0; gc < NcJ; ++gc) {
2258:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2260:       tmpBasisJ[gidx] = basisJ[gidx];
2261:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2262:     }
2263:   }
2264:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2265:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2266:   for (f = 0; f < NbI; ++f) {
2267:     for (fc = 0; fc < NcI; ++fc) {
2268:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2269:       const PetscInt i    = offsetI+f; /* Element matrix row */
2270:       for (g = 0; g < NbJ; ++g) {
2271:         for (gc = 0; gc < NcJ; ++gc) {
2272:           const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2273:           const PetscInt j    = offsetJ+g; /* Element matrix column */
2274:           const PetscInt fOff = eOffset+i*totDim+j;

2276:           elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2277:           for (df = 0; df < dE; ++df) {
2278:             elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2279:             elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2280:             for (dg = 0; dg < dE; ++dg) {
2281:               elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2282:             }
2283:           }
2284:         }
2285:       }
2286:     }
2287:   }
2288:   return(0);
2289: }

2291: 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[])
2292: {
2293:   const PetscInt   dE        = TI->cdim;
2294:   const PetscInt   NqI       = TI->Np;
2295:   const PetscInt   NbI       = TI->Nb;
2296:   const PetscInt   NcI       = TI->Nc;
2297:   const PetscReal *basisI    = &TI->T[0][(r*NqI+q)*NbI*NcI];
2298:   const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2299:   const PetscInt   NqJ       = TJ->Np;
2300:   const PetscInt   NbJ       = TJ->Nb;
2301:   const PetscInt   NcJ       = TJ->Nc;
2302:   const PetscReal *basisJ    = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2303:   const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2304:   const PetscInt   Ns = isHybridI ? 1 : 2;
2305:   const PetscInt   Nt = isHybridJ ? 1 : 2;
2306:   PetscInt         f, fc, g, gc, df, dg, s, t;
2307:   PetscErrorCode   ierr;

2309:   for (f = 0; f < NbI; ++f) {
2310:     for (fc = 0; fc < NcI; ++fc) {
2311:       const PetscInt fidx = f*NcI+fc; /* Test function basis index */

2313:       tmpBasisI[fidx] = basisI[fidx];
2314:       for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2315:     }
2316:   }
2317:   PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2318:   PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2319:   for (g = 0; g < NbJ; ++g) {
2320:     for (gc = 0; gc < NcJ; ++gc) {
2321:       const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */

2323:       tmpBasisJ[gidx] = basisJ[gidx];
2324:       for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2325:     }
2326:   }
2327:   PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2328:   PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2329:   for (s = 0; s < Ns; ++s) {
2330:     for (f = 0; f < NbI; ++f) {
2331:       for (fc = 0; fc < NcI; ++fc) {
2332:         const PetscInt sc   = NcI*s+fc;  /* components from each side of the surface */
2333:         const PetscInt fidx = f*NcI+fc;  /* Test function basis index */
2334:         const PetscInt i    = offsetI+NbI*s+f; /* Element matrix row */
2335:         for (t = 0; t < Nt; ++t) {
2336:           for (g = 0; g < NbJ; ++g) {
2337:             for (gc = 0; gc < NcJ; ++gc) {
2338:               const PetscInt tc   = NcJ*t+gc;  /* components from each side of the surface */
2339:               const PetscInt gidx = g*NcJ+gc;  /* Trial function basis index */
2340:               const PetscInt j    = offsetJ+NbJ*t+g; /* Element matrix column */
2341:               const PetscInt fOff = eOffset+i*totDim+j;

2343:               elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2344:               for (df = 0; df < dE; ++df) {
2345:                 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2346:                 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2347:                 for (dg = 0; dg < dE; ++dg) {
2348:                   elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2349:                 }
2350:               }
2351:             }
2352:           }
2353:         }
2354:       }
2355:     }
2356:   }
2357:   return(0);
2358: }

2360: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2361: {
2362:   PetscDualSpace  dsp;
2363:   DM              dm;
2364:   PetscQuadrature quadDef;
2365:   PetscInt        dim, cdim, Nq;
2366:   PetscErrorCode  ierr;

2369:   PetscFEGetDualSpace(fe, &dsp);
2370:   PetscDualSpaceGetDM(dsp, &dm);
2371:   DMGetDimension(dm, &dim);
2372:   DMGetCoordinateDim(dm, &cdim);
2373:   PetscFEGetQuadrature(fe, &quadDef);
2374:   quad = quad ? quad : quadDef;
2375:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2376:   PetscMalloc1(Nq*cdim,      &cgeom->v);
2377:   PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2378:   PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2379:   PetscMalloc1(Nq,           &cgeom->detJ);
2380:   cgeom->dim       = dim;
2381:   cgeom->dimEmbed  = cdim;
2382:   cgeom->numCells  = 1;
2383:   cgeom->numPoints = Nq;
2384:   DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2385:   return(0);
2386: }

2388: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2389: {

2393:   PetscFree(cgeom->v);
2394:   PetscFree(cgeom->J);
2395:   PetscFree(cgeom->invJ);
2396:   PetscFree(cgeom->detJ);
2397:   return(0);
2398: }