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

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

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

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:   }
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);

335:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*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;
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: /*@

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: @*/
636: {
641:   return(0);
642: }

644: /*@

647:   Not collective

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

653:   Level: intermediate

655: .seealso: PetscFECreate()
656: @*/
658: {
659:   PetscInt       Nc, qNc;

664:   if (q == fem->quadrature) return(0);
665:   PetscFEGetNumComponents(fem, &Nc);
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);
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: @*/
692: {
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: @*/
714: {
715:   PetscInt       Nc, qNc;

720:   PetscFEGetNumComponents(fem, &Nc);
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);
726:   PetscObjectReference((PetscObject) q);
727:   return(0);
728: }

730: /*@

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

742: @*/
743: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
744: {
746:   PetscErrorCode  ierr;

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;
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);
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.

1670: .seealso: PetscFECreateDefault()
1671: @*/
1672: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1673: {
1674:   PetscSpace      P, subP;
1675:   PetscDualSpace  Q, 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);
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);
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

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;
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);
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);
1784:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1785:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &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: {
1817:   DM              K;
1818:   PetscSpace      P;
1819:   PetscDualSpace  Q;
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);
1861:   PetscOptionsEnd();
1862:   quadPointsPerEdge = PetscMax(qorder + 1,1);
1863:   if (isSimplex) {
1866:   } else {
1869:   }
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: {
1903:   DM              K;
1904:   PetscSpace      P;
1905:   PetscDualSpace  Q;
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) {
1947:   } else {
1950:   }
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]);
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]);
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);
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);
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);
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);
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);
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);
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: }

2313: {
2314:   PetscDualSpace  dsp;
2315:   DM              dm;
2317:   PetscInt        dim, cdim, Nq;
2318:   PetscErrorCode  ierr;

2321:   PetscFEGetDualSpace(fe, &dsp);
2322:   PetscDualSpaceGetDM(dsp, &dm);
2323:   DMGetDimension(dm, &dim);
2324:   DMGetCoordinateDim(dm, &cdim);
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: }