Actual source code: fe.c
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: {
92: PetscFunctionListAdd(&PetscFEList, sname, function);
93: return 0;
94: }
96: /*@C
97: PetscFESetType - Builds a particular PetscFE
99: Collective on fem
101: Input Parameters:
102: + fem - The PetscFE object
103: - name - The kind of FEM space
105: Options Database Key:
106: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
108: Level: intermediate
110: .seealso: `PetscFEGetType()`, `PetscFECreate()`
111: @*/
112: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
113: {
114: PetscErrorCode (*r)(PetscFE);
115: PetscBool match;
118: PetscObjectTypeCompare((PetscObject)fem, name, &match);
119: if (match) return 0;
121: if (!PetscFERegisterAllCalled) PetscFERegisterAll();
122: PetscFunctionListFind(PetscFEList, name, &r);
125: PetscTryTypeMethod(fem, destroy);
126: fem->ops->destroy = NULL;
128: (*r)(fem);
129: PetscObjectChangeTypeName((PetscObject)fem, name);
130: return 0;
131: }
133: /*@C
134: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
136: Not Collective
138: Input Parameter:
139: . fem - The PetscFE
141: Output Parameter:
142: . name - The PetscFE type name
144: Level: intermediate
146: .seealso: `PetscFESetType()`, `PetscFECreate()`
147: @*/
148: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
149: {
152: if (!PetscFERegisterAllCalled) PetscFERegisterAll();
153: *name = ((PetscObject)fem)->type_name;
154: return 0;
155: }
157: /*@C
158: PetscFEViewFromOptions - View from Options
160: Collective on PetscFE
162: Input Parameters:
163: + A - the PetscFE object
164: . obj - Optional object
165: - name - command line option
167: Level: intermediate
168: .seealso: `PetscFE()`, `PetscFEView()`, `PetscObjectViewFromOptions()`, `PetscFECreate()`
169: @*/
170: PetscErrorCode PetscFEViewFromOptions(PetscFE A, PetscObject obj, const char name[])
171: {
173: PetscObjectViewFromOptions((PetscObject)A, obj, name);
174: return 0;
175: }
177: /*@C
178: PetscFEView - Views a PetscFE
180: Collective on fem
182: Input Parameters:
183: + fem - the PetscFE object to view
184: - viewer - the viewer
186: Level: beginner
188: .seealso `PetscFEDestroy()`
189: @*/
190: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
191: {
192: PetscBool iascii;
196: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fem), &viewer);
197: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
198: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
199: PetscTryTypeMethod(fem, view, viewer);
200: return 0;
201: }
203: /*@
204: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
206: Collective on fem
208: Input Parameter:
209: . fem - the PetscFE object to set options for
211: Options Database:
212: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
213: - -petscfe_num_batches - the number of cell batches to integrate serially
215: Level: intermediate
217: .seealso `PetscFEView()`
218: @*/
219: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
220: {
221: const char *defaultType;
222: char name[256];
223: PetscBool flg;
226: if (!((PetscObject)fem)->type_name) {
227: defaultType = PETSCFEBASIC;
228: } else {
229: defaultType = ((PetscObject)fem)->type_name;
230: }
231: if (!PetscFERegisterAllCalled) PetscFERegisterAll();
233: PetscObjectOptionsBegin((PetscObject)fem);
234: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
235: if (flg) {
236: PetscFESetType(fem, name);
237: } else if (!((PetscObject)fem)->type_name) {
238: PetscFESetType(fem, defaultType);
239: }
240: PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL, 1);
241: PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL, 1);
242: PetscTryTypeMethod(fem, setfromoptions, PetscOptionsObject);
243: /* process any options handlers added with PetscObjectAddOptionsHandler() */
244: PetscObjectProcessOptionsHandlers((PetscObject)fem, PetscOptionsObject);
245: PetscOptionsEnd();
246: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
247: return 0;
248: }
250: /*@C
251: PetscFESetUp - Construct data structures for the PetscFE
253: Collective on fem
255: Input Parameter:
256: . fem - the PetscFE object to setup
258: Level: intermediate
260: .seealso `PetscFEView()`, `PetscFEDestroy()`
261: @*/
262: PetscErrorCode PetscFESetUp(PetscFE fem)
263: {
265: if (fem->setupcalled) return 0;
266: PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
267: fem->setupcalled = PETSC_TRUE;
268: PetscTryTypeMethod(fem, setup);
269: PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
270: return 0;
271: }
273: /*@
274: PetscFEDestroy - Destroys a PetscFE object
276: Collective on fem
278: Input Parameter:
279: . fem - the PetscFE object to destroy
281: Level: beginner
283: .seealso `PetscFEView()`
284: @*/
285: PetscErrorCode PetscFEDestroy(PetscFE *fem)
286: {
287: if (!*fem) return 0;
290: if (--((PetscObject)(*fem))->refct > 0) {
291: *fem = NULL;
292: return 0;
293: }
294: ((PetscObject)(*fem))->refct = 0;
296: if ((*fem)->subspaces) {
297: PetscInt dim, d;
299: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
300: for (d = 0; d < dim; ++d) PetscFEDestroy(&(*fem)->subspaces[d]);
301: }
302: PetscFree((*fem)->subspaces);
303: PetscFree((*fem)->invV);
304: PetscTabulationDestroy(&(*fem)->T);
305: PetscTabulationDestroy(&(*fem)->Tf);
306: PetscTabulationDestroy(&(*fem)->Tc);
307: PetscSpaceDestroy(&(*fem)->basisSpace);
308: PetscDualSpaceDestroy(&(*fem)->dualSpace);
309: PetscQuadratureDestroy(&(*fem)->quadrature);
310: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
311: #ifdef PETSC_HAVE_LIBCEED
312: CeedBasisDestroy(&(*fem)->ceedBasis);
313: CeedDestroy(&(*fem)->ceed);
314: #endif
316: PetscTryTypeMethod((*fem), destroy);
317: PetscHeaderDestroy(fem);
318: return 0;
319: }
321: /*@
322: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
324: Collective
326: Input Parameter:
327: . comm - The communicator for the PetscFE object
329: Output Parameter:
330: . fem - The PetscFE object
332: Level: beginner
334: .seealso: `PetscFESetType()`, `PETSCFEGALERKIN`
335: @*/
336: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
337: {
338: PetscFE f;
341: PetscCitationsRegister(FECitation, &FEcite);
342: *fem = NULL;
343: PetscFEInitializePackage();
345: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
347: f->basisSpace = NULL;
348: f->dualSpace = NULL;
349: f->numComponents = 1;
350: f->subspaces = NULL;
351: f->invV = NULL;
352: f->T = NULL;
353: f->Tf = NULL;
354: f->Tc = NULL;
355: PetscArrayzero(&f->quadrature, 1);
356: PetscArrayzero(&f->faceQuadrature, 1);
357: f->blockSize = 0;
358: f->numBlocks = 1;
359: f->batchSize = 0;
360: f->numBatches = 1;
362: *fem = f;
363: return 0;
364: }
366: /*@
367: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
369: Not collective
371: Input Parameter:
372: . fem - The PetscFE object
374: Output Parameter:
375: . dim - The spatial dimension
377: Level: intermediate
379: .seealso: `PetscFECreate()`
380: @*/
381: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
382: {
383: DM dm;
387: PetscDualSpaceGetDM(fem->dualSpace, &dm);
388: DMGetDimension(dm, dim);
389: return 0;
390: }
392: /*@
393: PetscFESetNumComponents - Sets the number of components in the element
395: Not collective
397: Input Parameters:
398: + fem - The PetscFE object
399: - comp - The number of field components
401: Level: intermediate
403: .seealso: `PetscFECreate()`
404: @*/
405: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
406: {
408: fem->numComponents = comp;
409: return 0;
410: }
412: /*@
413: PetscFEGetNumComponents - Returns the number of components in the element
415: Not collective
417: Input Parameter:
418: . fem - The PetscFE object
420: Output Parameter:
421: . comp - The number of field components
423: Level: intermediate
425: .seealso: `PetscFECreate()`
426: @*/
427: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
428: {
431: *comp = fem->numComponents;
432: return 0;
433: }
435: /*@
436: PetscFESetTileSizes - Sets the tile sizes for evaluation
438: Not collective
440: Input Parameters:
441: + fem - The PetscFE object
442: . blockSize - The number of elements in a block
443: . numBlocks - The number of blocks in a batch
444: . batchSize - The number of elements in a batch
445: - numBatches - The number of batches in a chunk
447: Level: intermediate
449: .seealso: `PetscFECreate()`
450: @*/
451: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
452: {
454: fem->blockSize = blockSize;
455: fem->numBlocks = numBlocks;
456: fem->batchSize = batchSize;
457: fem->numBatches = numBatches;
458: return 0;
459: }
461: /*@
462: PetscFEGetTileSizes - Returns the tile sizes for evaluation
464: Not collective
466: Input Parameter:
467: . fem - The PetscFE object
469: Output Parameters:
470: + blockSize - The number of elements in a block
471: . numBlocks - The number of blocks in a batch
472: . batchSize - The number of elements in a batch
473: - numBatches - The number of batches in a chunk
475: Level: intermediate
477: .seealso: `PetscFECreate()`
478: @*/
479: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
480: {
486: if (blockSize) *blockSize = fem->blockSize;
487: if (numBlocks) *numBlocks = fem->numBlocks;
488: if (batchSize) *batchSize = fem->batchSize;
489: if (numBatches) *numBatches = fem->numBatches;
490: return 0;
491: }
493: /*@
494: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
496: Not collective
498: Input Parameter:
499: . fem - The PetscFE object
501: Output Parameter:
502: . sp - The PetscSpace object
504: Level: intermediate
506: .seealso: `PetscFECreate()`
507: @*/
508: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
509: {
512: *sp = fem->basisSpace;
513: return 0;
514: }
516: /*@
517: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
519: Not collective
521: Input Parameters:
522: + fem - The PetscFE object
523: - sp - The PetscSpace object
525: Level: intermediate
527: .seealso: `PetscFECreate()`
528: @*/
529: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
530: {
533: PetscSpaceDestroy(&fem->basisSpace);
534: fem->basisSpace = sp;
535: PetscObjectReference((PetscObject)fem->basisSpace);
536: return 0;
537: }
539: /*@
540: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
542: Not collective
544: Input Parameter:
545: . fem - The PetscFE object
547: Output Parameter:
548: . sp - The PetscDualSpace object
550: Level: intermediate
552: .seealso: `PetscFECreate()`
553: @*/
554: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
555: {
558: *sp = fem->dualSpace;
559: return 0;
560: }
562: /*@
563: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
565: Not collective
567: Input Parameters:
568: + fem - The PetscFE object
569: - sp - The PetscDualSpace object
571: Level: intermediate
573: .seealso: `PetscFECreate()`
574: @*/
575: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
576: {
579: PetscDualSpaceDestroy(&fem->dualSpace);
580: fem->dualSpace = sp;
581: PetscObjectReference((PetscObject)fem->dualSpace);
582: return 0;
583: }
585: /*@
586: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
588: Not collective
590: Input Parameter:
591: . fem - The PetscFE object
593: Output Parameter:
594: . q - The PetscQuadrature object
596: Level: intermediate
598: .seealso: `PetscFECreate()`
599: @*/
600: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
601: {
604: *q = fem->quadrature;
605: return 0;
606: }
608: /*@
609: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
611: Not collective
613: Input Parameters:
614: + fem - The PetscFE object
615: - q - The PetscQuadrature object
617: Level: intermediate
619: .seealso: `PetscFECreate()`
620: @*/
621: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
622: {
623: PetscInt Nc, qNc;
626: if (q == fem->quadrature) return 0;
627: PetscFEGetNumComponents(fem, &Nc);
628: PetscQuadratureGetNumComponents(q, &qNc);
630: PetscTabulationDestroy(&fem->T);
631: PetscTabulationDestroy(&fem->Tc);
632: PetscObjectReference((PetscObject)q);
633: PetscQuadratureDestroy(&fem->quadrature);
634: fem->quadrature = q;
635: return 0;
636: }
638: /*@
639: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
641: Not collective
643: Input Parameter:
644: . fem - The PetscFE object
646: Output Parameter:
647: . q - The PetscQuadrature object
649: Level: intermediate
651: .seealso: `PetscFECreate()`
652: @*/
653: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
654: {
657: *q = fem->faceQuadrature;
658: return 0;
659: }
661: /*@
662: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
664: Not collective
666: Input Parameters:
667: + fem - The PetscFE object
668: - q - The PetscQuadrature object
670: Level: intermediate
672: .seealso: `PetscFECreate()`
673: @*/
674: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
675: {
676: PetscInt Nc, qNc;
679: PetscFEGetNumComponents(fem, &Nc);
680: PetscQuadratureGetNumComponents(q, &qNc);
682: PetscTabulationDestroy(&fem->Tf);
683: PetscQuadratureDestroy(&fem->faceQuadrature);
684: fem->faceQuadrature = q;
685: PetscObjectReference((PetscObject)q);
686: return 0;
687: }
689: /*@
690: PetscFECopyQuadrature - Copy both volumetric and surface quadrature
692: Not collective
694: Input Parameters:
695: + sfe - The PetscFE source for the quadratures
696: - tfe - The PetscFE target for the quadratures
698: Level: intermediate
700: .seealso: `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
701: @*/
702: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
703: {
704: PetscQuadrature q;
708: PetscFEGetQuadrature(sfe, &q);
709: PetscFESetQuadrature(tfe, q);
710: PetscFEGetFaceQuadrature(sfe, &q);
711: PetscFESetFaceQuadrature(tfe, q);
712: return 0;
713: }
715: /*@C
716: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
718: Not collective
720: Input Parameter:
721: . fem - The PetscFE object
723: Output Parameter:
724: . numDof - Array with the number of dofs per dimension
726: Level: intermediate
728: .seealso: `PetscFECreate()`
729: @*/
730: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
731: {
734: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
735: return 0;
736: }
738: /*@C
739: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
741: Not collective
743: Input Parameters:
744: + fem - The PetscFE object
745: - k - The highest derivative we need to tabulate, very often 1
747: Output Parameter:
748: . T - The basis function values and derivatives at quadrature points
750: Note:
751: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
752: $ 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
753: $ 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
755: Level: intermediate
757: .seealso: `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
758: @*/
759: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
760: {
761: PetscInt npoints;
762: const PetscReal *points;
766: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
767: if (!fem->T) PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);
769: *T = fem->T;
770: return 0;
771: }
773: /*@C
774: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
776: Not collective
778: Input Parameters:
779: + fem - The PetscFE object
780: - k - The highest derivative we need to tabulate, very often 1
782: Output Parameters:
783: . Tf - The basis function values and derivatives at face quadrature points
785: Note:
786: $ 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
787: $ 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
788: $ 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
790: Level: intermediate
792: .seealso: `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
793: @*/
794: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
795: {
798: if (!fem->Tf) {
799: const PetscReal xi0[3] = {-1., -1., -1.};
800: PetscReal v0[3], J[9], detJ;
801: PetscQuadrature fq;
802: PetscDualSpace sp;
803: DM dm;
804: const PetscInt *faces;
805: PetscInt dim, numFaces, f, npoints, q;
806: const PetscReal *points;
807: PetscReal *facePoints;
809: PetscFEGetDualSpace(fem, &sp);
810: PetscDualSpaceGetDM(sp, &dm);
811: DMGetDimension(dm, &dim);
812: DMPlexGetConeSize(dm, 0, &numFaces);
813: DMPlexGetCone(dm, 0, &faces);
814: PetscFEGetFaceQuadrature(fem, &fq);
815: if (fq) {
816: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
817: PetscMalloc1(numFaces * npoints * dim, &facePoints);
818: for (f = 0; f < numFaces; ++f) {
819: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
820: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * npoints + q) * dim]);
821: }
822: PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
823: PetscFree(facePoints);
824: }
825: }
827: *Tf = fem->Tf;
828: return 0;
829: }
831: /*@C
832: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
834: Not collective
836: Input Parameter:
837: . fem - The PetscFE object
839: Output Parameters:
840: . Tc - The basis function values at face centroid points
842: Note:
843: $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
845: Level: intermediate
847: .seealso: `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
848: @*/
849: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
850: {
853: if (!fem->Tc) {
854: PetscDualSpace sp;
855: DM dm;
856: const PetscInt *cone;
857: PetscReal *centroids;
858: PetscInt dim, numFaces, f;
860: PetscFEGetDualSpace(fem, &sp);
861: PetscDualSpaceGetDM(sp, &dm);
862: DMGetDimension(dm, &dim);
863: DMPlexGetConeSize(dm, 0, &numFaces);
864: DMPlexGetCone(dm, 0, &cone);
865: PetscMalloc1(numFaces * dim, ¢roids);
866: for (f = 0; f < numFaces; ++f) DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f * dim], NULL);
867: PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
868: PetscFree(centroids);
869: }
870: *Tc = fem->Tc;
871: return 0;
872: }
874: /*@C
875: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
877: Not collective
879: Input Parameters:
880: + fem - The PetscFE object
881: . nrepl - The number of replicas
882: . npoints - The number of tabulation points in a replica
883: . points - The tabulation point coordinates
884: - K - The number of derivatives calculated
886: Output Parameter:
887: . T - The basis function values and derivatives at tabulation points
889: Note:
890: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
891: $ 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
892: $ 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
894: Level: intermediate
896: .seealso: `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
897: @*/
898: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
899: {
900: DM dm;
901: PetscDualSpace Q;
902: PetscInt Nb; /* Dimension of FE space P */
903: PetscInt Nc; /* Field components */
904: PetscInt cdim; /* Reference coordinate dimension */
905: PetscInt k;
907: if (!npoints || !fem->dualSpace || K < 0) {
908: *T = NULL;
909: return 0;
910: }
914: PetscFEGetDualSpace(fem, &Q);
915: PetscDualSpaceGetDM(Q, &dm);
916: DMGetDimension(dm, &cdim);
917: PetscDualSpaceGetDimension(Q, &Nb);
918: PetscFEGetNumComponents(fem, &Nc);
919: PetscMalloc1(1, T);
920: (*T)->K = !cdim ? 0 : K;
921: (*T)->Nr = nrepl;
922: (*T)->Np = npoints;
923: (*T)->Nb = Nb;
924: (*T)->Nc = Nc;
925: (*T)->cdim = cdim;
926: PetscMalloc1((*T)->K + 1, &(*T)->T);
927: for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
928: PetscUseTypeMethod(fem, createtabulation, nrepl * npoints, points, K, *T);
929: return 0;
930: }
932: /*@C
933: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
935: Not collective
937: Input Parameters:
938: + fem - The PetscFE object
939: . npoints - The number of tabulation points
940: . points - The tabulation point coordinates
941: . K - The number of derivatives calculated
942: - T - An existing tabulation object with enough allocated space
944: Output Parameter:
945: . T - The basis function values and derivatives at tabulation points
947: Note:
948: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
949: $ 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
950: $ 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
952: Level: intermediate
954: .seealso: `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
955: @*/
956: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
957: {
959: if (!npoints || !fem->dualSpace || K < 0) return 0;
963: if (PetscDefined(USE_DEBUG)) {
964: DM dm;
965: PetscDualSpace Q;
966: PetscInt Nb; /* Dimension of FE space P */
967: PetscInt Nc; /* Field components */
968: PetscInt cdim; /* Reference coordinate dimension */
970: PetscFEGetDualSpace(fem, &Q);
971: PetscDualSpaceGetDM(Q, &dm);
972: DMGetDimension(dm, &cdim);
973: PetscDualSpaceGetDimension(Q, &Nb);
974: PetscFEGetNumComponents(fem, &Nc);
979: }
980: T->Nr = 1;
981: T->Np = npoints;
982: PetscUseTypeMethod(fem, createtabulation, npoints, points, K, T);
983: return 0;
984: }
986: /*@C
987: PetscTabulationDestroy - Frees memory from the associated tabulation.
989: Not collective
991: Input Parameter:
992: . T - The tabulation
994: Level: intermediate
996: .seealso: `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
997: @*/
998: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
999: {
1000: PetscInt k;
1003: if (!T || !(*T)) return 0;
1004: for (k = 0; k <= (*T)->K; ++k) PetscFree((*T)->T[k]);
1005: PetscFree((*T)->T);
1006: PetscFree(*T);
1007: *T = NULL;
1008: return 0;
1009: }
1011: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1012: {
1013: PetscSpace bsp, bsubsp;
1014: PetscDualSpace dsp, dsubsp;
1015: PetscInt dim, depth, numComp, i, j, coneSize, order;
1016: PetscFEType type;
1017: DM dm;
1018: DMLabel label;
1019: PetscReal *xi, *v, *J, detJ;
1020: const char *name;
1021: PetscQuadrature origin, fullQuad, subQuad;
1025: PetscFEGetBasisSpace(fe, &bsp);
1026: PetscFEGetDualSpace(fe, &dsp);
1027: PetscDualSpaceGetDM(dsp, &dm);
1028: DMGetDimension(dm, &dim);
1029: DMPlexGetDepthLabel(dm, &label);
1030: DMLabelGetValue(label, refPoint, &depth);
1031: PetscCalloc1(depth, &xi);
1032: PetscMalloc1(dim, &v);
1033: PetscMalloc1(dim * dim, &J);
1034: for (i = 0; i < depth; i++) xi[i] = 0.;
1035: PetscQuadratureCreate(PETSC_COMM_SELF, &origin);
1036: PetscQuadratureSetData(origin, depth, 0, 1, xi, NULL);
1037: DMPlexComputeCellGeometryFEM(dm, refPoint, origin, v, J, NULL, &detJ);
1038: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1039: for (i = 1; i < dim; i++) {
1040: for (j = 0; j < depth; j++) J[i * depth + j] = J[i * dim + j];
1041: }
1042: PetscQuadratureDestroy(&origin);
1043: PetscDualSpaceGetPointSubspace(dsp, refPoint, &dsubsp);
1044: PetscSpaceCreateSubspace(bsp, dsubsp, v, J, NULL, NULL, PETSC_OWN_POINTER, &bsubsp);
1045: PetscSpaceSetUp(bsubsp);
1046: PetscFECreate(PetscObjectComm((PetscObject)fe), trFE);
1047: PetscFEGetType(fe, &type);
1048: PetscFESetType(*trFE, type);
1049: PetscFEGetNumComponents(fe, &numComp);
1050: PetscFESetNumComponents(*trFE, numComp);
1051: PetscFESetBasisSpace(*trFE, bsubsp);
1052: PetscFESetDualSpace(*trFE, dsubsp);
1053: PetscObjectGetName((PetscObject)fe, &name);
1054: if (name) PetscFESetName(*trFE, name);
1055: PetscFEGetQuadrature(fe, &fullQuad);
1056: PetscQuadratureGetOrder(fullQuad, &order);
1057: DMPlexGetConeSize(dm, refPoint, &coneSize);
1058: if (coneSize == 2 * depth) PetscDTGaussTensorQuadrature(depth, 1, (order + 1) / 2, -1., 1., &subQuad);
1059: else PetscDTStroudConicalQuadrature(depth, 1, (order + 1) / 2, -1., 1., &subQuad);
1060: PetscFESetQuadrature(*trFE, subQuad);
1061: PetscFESetUp(*trFE);
1062: PetscQuadratureDestroy(&subQuad);
1063: PetscSpaceDestroy(&bsubsp);
1064: return 0;
1065: }
1067: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1068: {
1069: PetscInt hStart, hEnd;
1070: PetscDualSpace dsp;
1071: DM dm;
1075: *trFE = NULL;
1076: PetscFEGetDualSpace(fe, &dsp);
1077: PetscDualSpaceGetDM(dsp, &dm);
1078: DMPlexGetHeightStratum(dm, height, &hStart, &hEnd);
1079: if (hEnd <= hStart) return 0;
1080: PetscFECreatePointTrace(fe, hStart, trFE);
1081: return 0;
1082: }
1084: /*@
1085: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1087: Not collective
1089: Input Parameter:
1090: . fe - The PetscFE
1092: Output Parameter:
1093: . dim - The dimension
1095: Level: intermediate
1097: .seealso: `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1098: @*/
1099: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1100: {
1103: PetscTryTypeMethod(fem, getdimension, dim);
1104: return 0;
1105: }
1107: /*@C
1108: PetscFEPushforward - Map the reference element function to real space
1110: Input Parameters:
1111: + fe - The PetscFE
1112: . fegeom - The cell geometry
1113: . Nv - The number of function values
1114: - vals - The function values
1116: Output Parameter:
1117: . vals - The transformed function values
1119: Level: advanced
1121: Note: This just forwards the call onto PetscDualSpacePushforward().
1123: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1125: .seealso: `PetscDualSpacePushforward()`
1126: @*/
1127: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1128: {
1130: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1131: return 0;
1132: }
1134: /*@C
1135: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1137: Input Parameters:
1138: + fe - The PetscFE
1139: . fegeom - The cell geometry
1140: . Nv - The number of function gradient values
1141: - vals - The function gradient values
1143: Output Parameter:
1144: . vals - The transformed function gradient values
1146: Level: advanced
1148: Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1150: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1152: .seealso: `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1153: @*/
1154: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1155: {
1157: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1158: return 0;
1159: }
1161: /*@C
1162: PetscFEPushforwardHessian - Map the reference element function Hessian to real space
1164: Input Parameters:
1165: + fe - The PetscFE
1166: . fegeom - The cell geometry
1167: . Nv - The number of function Hessian values
1168: - vals - The function Hessian values
1170: Output Parameter:
1171: . vals - The transformed function Hessian values
1173: Level: advanced
1175: Note: This just forwards the call onto PetscDualSpacePushforwardHessian().
1177: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1179: .seealso: `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1180: @*/
1181: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1182: {
1184: PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1185: return 0;
1186: }
1188: /*
1189: Purpose: Compute element vector for chunk of elements
1191: Input:
1192: Sizes:
1193: Ne: number of elements
1194: Nf: number of fields
1195: PetscFE
1196: dim: spatial dimension
1197: Nb: number of basis functions
1198: Nc: number of field components
1199: PetscQuadrature
1200: Nq: number of quadrature points
1202: Geometry:
1203: PetscFEGeom[Ne] possibly *Nq
1204: PetscReal v0s[dim]
1205: PetscReal n[dim]
1206: PetscReal jacobians[dim*dim]
1207: PetscReal jacobianInverses[dim*dim]
1208: PetscReal jacobianDeterminants
1209: FEM:
1210: PetscFE
1211: PetscQuadrature
1212: PetscReal quadPoints[Nq*dim]
1213: PetscReal quadWeights[Nq]
1214: PetscReal basis[Nq*Nb*Nc]
1215: PetscReal basisDer[Nq*Nb*Nc*dim]
1216: PetscScalar coefficients[Ne*Nb*Nc]
1217: PetscScalar elemVec[Ne*Nb*Nc]
1219: Problem:
1220: PetscInt f: the active field
1221: f0, f1
1223: Work Space:
1224: PetscFE
1225: PetscScalar f0[Nq*dim];
1226: PetscScalar f1[Nq*dim*dim];
1227: PetscScalar u[Nc];
1228: PetscScalar gradU[Nc*dim];
1229: PetscReal x[dim];
1230: PetscScalar realSpaceDer[dim];
1232: Purpose: Compute element vector for N_cb batches of elements
1234: Input:
1235: Sizes:
1236: N_cb: Number of serial cell batches
1238: Geometry:
1239: PetscReal v0s[Ne*dim]
1240: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1241: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1242: PetscReal jacobianDeterminants[Ne] possibly *Nq
1243: FEM:
1244: static PetscReal quadPoints[Nq*dim]
1245: static PetscReal quadWeights[Nq]
1246: static PetscReal basis[Nq*Nb*Nc]
1247: static PetscReal basisDer[Nq*Nb*Nc*dim]
1248: PetscScalar coefficients[Ne*Nb*Nc]
1249: PetscScalar elemVec[Ne*Nb*Nc]
1251: ex62.c:
1252: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1253: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1254: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1255: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1257: ex52.c:
1258: 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)
1259: 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)
1261: ex52_integrateElement.cu
1262: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1264: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1265: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1266: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1268: ex52_integrateElementOpenCL.c:
1269: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1270: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1271: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1273: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1274: */
1276: /*@C
1277: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1279: Not collective
1281: Input Parameters:
1282: + prob - The PetscDS specifying the discretizations and continuum functions
1283: . field - The field being integrated
1284: . Ne - The number of elements in the chunk
1285: . cgeom - The cell geometry for each cell in the chunk
1286: . coefficients - The array of FEM basis coefficients for the elements
1287: . probAux - The PetscDS specifying the auxiliary discretizations
1288: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1290: Output Parameter:
1291: . integral - the integral for this field
1293: Level: intermediate
1295: .seealso: `PetscFEIntegrateResidual()`
1296: @*/
1297: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1298: {
1299: PetscFE fe;
1302: PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1303: if (fe->ops->integrate) (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);
1304: return 0;
1305: }
1307: /*@C
1308: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1310: Not collective
1312: Input Parameters:
1313: + prob - The PetscDS specifying the discretizations and continuum functions
1314: . field - The field being integrated
1315: . obj_func - The function to be integrated
1316: . Ne - The number of elements in the chunk
1317: . fgeom - The face geometry for each face in the chunk
1318: . coefficients - The array of FEM basis coefficients for the elements
1319: . probAux - The PetscDS specifying the auxiliary discretizations
1320: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1322: Output Parameter:
1323: . integral - the integral for this field
1325: Level: intermediate
1327: .seealso: `PetscFEIntegrateResidual()`
1328: @*/
1329: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, void (*obj_func)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1330: {
1331: PetscFE fe;
1334: PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
1335: if (fe->ops->integratebd) (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);
1336: return 0;
1337: }
1339: /*@C
1340: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1342: Not collective
1344: Input Parameters:
1345: + ds - The PetscDS specifying the discretizations and continuum functions
1346: . key - The (label+value, field) being integrated
1347: . Ne - The number of elements in the chunk
1348: . cgeom - The cell geometry for each cell in the chunk
1349: . coefficients - The array of FEM basis coefficients for the elements
1350: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1351: . probAux - The PetscDS specifying the auxiliary discretizations
1352: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1353: - t - The time
1355: Output Parameter:
1356: . elemVec - the element residual vectors from each element
1358: Note:
1359: $ Loop over batch of elements (e):
1360: $ Loop over quadrature points (q):
1361: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1362: $ Call f_0 and f_1
1363: $ Loop over element vector entries (f,fc --> i):
1364: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1366: Level: intermediate
1368: .seealso: `PetscFEIntegrateResidual()`
1369: @*/
1370: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1371: {
1372: PetscFE fe;
1376: PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1377: if (fe->ops->integrateresidual) (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1378: return 0;
1379: }
1381: /*@C
1382: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1384: Not collective
1386: Input Parameters:
1387: + ds - The PetscDS specifying the discretizations and continuum functions
1388: . wf - The PetscWeakForm object holding the pointwise functions
1389: . key - The (label+value, field) being integrated
1390: . Ne - The number of elements in the chunk
1391: . fgeom - The face geometry for each cell in the chunk
1392: . coefficients - The array of FEM basis coefficients for the elements
1393: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1394: . probAux - The PetscDS specifying the auxiliary discretizations
1395: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1396: - t - The time
1398: Output Parameter:
1399: . elemVec - the element residual vectors from each element
1401: Level: intermediate
1403: .seealso: `PetscFEIntegrateResidual()`
1404: @*/
1405: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1406: {
1407: PetscFE fe;
1410: PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe);
1411: if (fe->ops->integratebdresidual) (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1412: return 0;
1413: }
1415: /*@C
1416: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1418: Not collective
1420: Input Parameters:
1421: + prob - The PetscDS specifying the discretizations and continuum functions
1422: . key - The (label+value, field) being integrated
1423: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1424: . Ne - The number of elements in the chunk
1425: . fgeom - The face geometry for each cell in the chunk
1426: . coefficients - The array of FEM basis coefficients for the elements
1427: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1428: . probAux - The PetscDS specifying the auxiliary discretizations
1429: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1430: - t - The time
1432: Output Parameter
1433: . elemVec - the element residual vectors from each element
1435: Level: developer
1437: .seealso: `PetscFEIntegrateResidual()`
1438: @*/
1439: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1440: {
1441: PetscFE fe;
1444: PetscDSGetDiscretization(prob, key.field, (PetscObject *)&fe);
1445: if (fe->ops->integratehybridresidual) (*fe->ops->integratehybridresidual)(prob, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
1446: return 0;
1447: }
1449: /*@C
1450: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1452: Not collective
1454: Input Parameters:
1455: + ds - The PetscDS specifying the discretizations and continuum functions
1456: . jtype - The type of matrix pointwise functions that should be used
1457: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1458: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1459: . Ne - The number of elements in the chunk
1460: . cgeom - The cell geometry for each cell in the chunk
1461: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1462: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1463: . probAux - The PetscDS specifying the auxiliary discretizations
1464: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1465: . t - The time
1466: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1468: Output Parameter:
1469: . elemMat - the element matrices for the Jacobian from each element
1471: Note:
1472: $ Loop over batch of elements (e):
1473: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1474: $ Loop over quadrature points (q):
1475: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1476: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1477: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1478: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1479: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1480: Level: intermediate
1482: .seealso: `PetscFEIntegrateResidual()`
1483: @*/
1484: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1485: {
1486: PetscFE fe;
1487: PetscInt Nf;
1490: PetscDSGetNumFields(ds, &Nf);
1491: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1492: if (fe->ops->integratejacobian) (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1493: return 0;
1494: }
1496: /*@C
1497: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1499: Not collective
1501: Input Parameters:
1502: + ds - The PetscDS specifying the discretizations and continuum functions
1503: . wf - The PetscWeakForm holding the pointwise functions
1504: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1505: . Ne - The number of elements in the chunk
1506: . fgeom - The face geometry for each cell in the chunk
1507: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1508: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1509: . probAux - The PetscDS specifying the auxiliary discretizations
1510: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1511: . t - The time
1512: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1514: Output Parameter:
1515: . elemMat - the element matrices for the Jacobian from each element
1517: Note:
1518: $ Loop over batch of elements (e):
1519: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1520: $ Loop over quadrature points (q):
1521: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1522: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1523: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1524: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1525: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1526: Level: intermediate
1528: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1529: @*/
1530: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1531: {
1532: PetscFE fe;
1533: PetscInt Nf;
1536: PetscDSGetNumFields(ds, &Nf);
1537: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1538: if (fe->ops->integratebdjacobian) (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1539: return 0;
1540: }
1542: /*@C
1543: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1545: Not collective
1547: Input Parameters:
1548: + ds - The PetscDS specifying the discretizations and continuum functions
1549: . jtype - The type of matrix pointwise functions that should be used
1550: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1551: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1552: . Ne - The number of elements in the chunk
1553: . fgeom - The face geometry for each cell in the chunk
1554: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1555: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1556: . probAux - The PetscDS specifying the auxiliary discretizations
1557: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1558: . t - The time
1559: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1561: Output Parameter
1562: . elemMat - the element matrices for the Jacobian from each element
1564: Note:
1565: $ Loop over batch of elements (e):
1566: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1567: $ Loop over quadrature points (q):
1568: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1569: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1570: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1571: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1572: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1573: Level: developer
1575: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1576: @*/
1577: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1578: {
1579: PetscFE fe;
1580: PetscInt Nf;
1583: PetscDSGetNumFields(ds, &Nf);
1584: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe);
1585: if (fe->ops->integratehybridjacobian) (*fe->ops->integratehybridjacobian)(ds, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);
1586: return 0;
1587: }
1589: /*@
1590: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1592: Input Parameters:
1593: + fe - The finite element space
1594: - height - The height of the Plex point
1596: Output Parameter:
1597: . subfe - The subspace of this FE space
1599: Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1601: Level: advanced
1603: .seealso: `PetscFECreateDefault()`
1604: @*/
1605: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1606: {
1607: PetscSpace P, subP;
1608: PetscDualSpace Q, subQ;
1609: PetscQuadrature subq;
1610: PetscFEType fetype;
1611: PetscInt dim, Nc;
1615: if (height == 0) {
1616: *subfe = fe;
1617: return 0;
1618: }
1619: PetscFEGetBasisSpace(fe, &P);
1620: PetscFEGetDualSpace(fe, &Q);
1621: PetscFEGetNumComponents(fe, &Nc);
1622: PetscFEGetFaceQuadrature(fe, &subq);
1623: PetscDualSpaceGetDimension(Q, &dim);
1625: if (!fe->subspaces) PetscCalloc1(dim, &fe->subspaces);
1626: if (height <= dim) {
1627: if (!fe->subspaces[height - 1]) {
1628: PetscFE sub = NULL;
1629: const char *name;
1631: PetscSpaceGetHeightSubspace(P, height, &subP);
1632: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1633: if (subQ) {
1634: PetscFECreate(PetscObjectComm((PetscObject)fe), &sub);
1635: PetscObjectGetName((PetscObject)fe, &name);
1636: PetscObjectSetName((PetscObject)sub, name);
1637: PetscFEGetType(fe, &fetype);
1638: PetscFESetType(sub, fetype);
1639: PetscFESetBasisSpace(sub, subP);
1640: PetscFESetDualSpace(sub, subQ);
1641: PetscFESetNumComponents(sub, Nc);
1642: PetscFESetUp(sub);
1643: PetscFESetQuadrature(sub, subq);
1644: }
1645: fe->subspaces[height - 1] = sub;
1646: }
1647: *subfe = fe->subspaces[height - 1];
1648: } else {
1649: *subfe = NULL;
1650: }
1651: return 0;
1652: }
1654: /*@
1655: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1656: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1657: sparsity). It is also used to create an interpolation between regularly refined meshes.
1659: Collective on fem
1661: Input Parameter:
1662: . fe - The initial PetscFE
1664: Output Parameter:
1665: . feRef - The refined PetscFE
1667: Level: advanced
1669: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1670: @*/
1671: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1672: {
1673: PetscSpace P, Pref;
1674: PetscDualSpace Q, Qref;
1675: DM K, Kref;
1676: PetscQuadrature q, qref;
1677: const PetscReal *v0, *jac;
1678: PetscInt numComp, numSubelements;
1679: PetscInt cStart, cEnd, c;
1680: PetscDualSpace *cellSpaces;
1682: PetscFEGetBasisSpace(fe, &P);
1683: PetscFEGetDualSpace(fe, &Q);
1684: PetscFEGetQuadrature(fe, &q);
1685: PetscDualSpaceGetDM(Q, &K);
1686: /* Create space */
1687: PetscObjectReference((PetscObject)P);
1688: Pref = P;
1689: /* Create dual space */
1690: PetscDualSpaceDuplicate(Q, &Qref);
1691: PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1692: DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref);
1693: PetscDualSpaceSetDM(Qref, Kref);
1694: DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1695: PetscMalloc1(cEnd - cStart, &cellSpaces);
1696: /* TODO: fix for non-uniform refinement */
1697: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1698: PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1699: PetscFree(cellSpaces);
1700: DMDestroy(&Kref);
1701: PetscDualSpaceSetUp(Qref);
1702: /* Create element */
1703: PetscFECreate(PetscObjectComm((PetscObject)fe), feRef);
1704: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1705: PetscFESetBasisSpace(*feRef, Pref);
1706: PetscFESetDualSpace(*feRef, Qref);
1707: PetscFEGetNumComponents(fe, &numComp);
1708: PetscFESetNumComponents(*feRef, numComp);
1709: PetscFESetUp(*feRef);
1710: PetscSpaceDestroy(&Pref);
1711: PetscDualSpaceDestroy(&Qref);
1712: /* Create quadrature */
1713: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1714: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1715: PetscFESetQuadrature(*feRef, qref);
1716: PetscQuadratureDestroy(&qref);
1717: return 0;
1718: }
1720: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1721: {
1722: PetscSpace P;
1723: PetscDualSpace Q;
1724: DM K;
1725: DMPolytopeType ct;
1726: PetscInt degree;
1727: char name[64];
1729: PetscFEGetBasisSpace(fe, &P);
1730: PetscSpaceGetDegree(P, °ree, NULL);
1731: PetscFEGetDualSpace(fe, &Q);
1732: PetscDualSpaceGetDM(Q, &K);
1733: DMPlexGetCellType(K, 0, &ct);
1734: switch (ct) {
1735: case DM_POLYTOPE_SEGMENT:
1736: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1737: case DM_POLYTOPE_QUADRILATERAL:
1738: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1739: case DM_POLYTOPE_HEXAHEDRON:
1740: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1741: PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree);
1742: break;
1743: case DM_POLYTOPE_TRIANGLE:
1744: case DM_POLYTOPE_TETRAHEDRON:
1745: PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree);
1746: break;
1747: case DM_POLYTOPE_TRI_PRISM:
1748: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1749: PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree);
1750: break;
1751: default:
1752: PetscSNPrintf(name, sizeof(name), "FE");
1753: }
1754: PetscFESetName(fe, name);
1755: return 0;
1756: }
1758: static PetscErrorCode PetscFECreateDefaultQuadrature_Private(PetscInt dim, DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq)
1759: {
1760: const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1);
1762: switch (ct) {
1763: case DM_POLYTOPE_SEGMENT:
1764: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1765: case DM_POLYTOPE_QUADRILATERAL:
1766: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1767: case DM_POLYTOPE_HEXAHEDRON:
1768: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1769: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1770: PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1771: break;
1772: case DM_POLYTOPE_TRIANGLE:
1773: case DM_POLYTOPE_TETRAHEDRON:
1774: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q);
1775: PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1776: break;
1777: case DM_POLYTOPE_TRI_PRISM:
1778: case DM_POLYTOPE_TRI_PRISM_TENSOR: {
1779: PetscQuadrature q1, q2;
1781: PetscDTStroudConicalQuadrature(2, 1, quadPointsPerEdge, -1.0, 1.0, &q1);
1782: PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2);
1783: PetscDTTensorQuadratureCreate(q1, q2, q);
1784: PetscQuadratureDestroy(&q1);
1785: PetscQuadratureDestroy(&q2);
1786: }
1787: PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq);
1788: /* TODO Need separate quadratures for each face */
1789: break;
1790: default:
1791: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
1792: }
1793: return 0;
1794: }
1796: /*@
1797: PetscFECreateFromSpaces - Create a PetscFE from the basis and dual spaces
1799: Collective
1801: Input Parameters:
1802: + P - The basis space
1803: . Q - The dual space
1804: . q - The cell quadrature
1805: - fq - The face quadrature
1807: Output Parameter:
1808: . fem - The PetscFE object
1810: Note:
1811: The PetscFE takes ownership of these spaces by calling destroy on each. They should not be used after this call, and for borrowed references from `PetscFEGetSpace()` and the like, the caller must use `PetscObjectReference` before this call.
1813: Level: beginner
1815: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1816: @*/
1817: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1818: {
1819: PetscInt Nc;
1820: const char *prefix;
1822: PetscFECreate(PetscObjectComm((PetscObject)P), fem);
1823: PetscObjectGetOptionsPrefix((PetscObject)P, &prefix);
1824: PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix);
1825: PetscFESetType(*fem, PETSCFEBASIC);
1826: PetscFESetBasisSpace(*fem, P);
1827: PetscFESetDualSpace(*fem, Q);
1828: PetscSpaceGetNumComponents(P, &Nc);
1829: PetscFESetNumComponents(*fem, Nc);
1830: PetscFESetUp(*fem);
1831: PetscSpaceDestroy(&P);
1832: PetscDualSpaceDestroy(&Q);
1833: PetscFESetQuadrature(*fem, q);
1834: PetscFESetFaceQuadrature(*fem, fq);
1835: PetscQuadratureDestroy(&q);
1836: PetscQuadratureDestroy(&fq);
1837: PetscFESetDefaultName_Private(*fem);
1838: return 0;
1839: }
1841: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1842: {
1843: DM K;
1844: PetscSpace P;
1845: PetscDualSpace Q;
1846: PetscQuadrature q, fq;
1847: PetscBool tensor;
1851: switch (ct) {
1852: case DM_POLYTOPE_SEGMENT:
1853: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1854: case DM_POLYTOPE_QUADRILATERAL:
1855: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1856: case DM_POLYTOPE_HEXAHEDRON:
1857: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1858: tensor = PETSC_TRUE;
1859: break;
1860: default:
1861: tensor = PETSC_FALSE;
1862: }
1863: /* Create space */
1864: PetscSpaceCreate(comm, &P);
1865: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1866: PetscObjectSetOptionsPrefix((PetscObject)P, prefix);
1867: PetscSpacePolynomialSetTensor(P, tensor);
1868: PetscSpaceSetNumComponents(P, Nc);
1869: PetscSpaceSetNumVariables(P, dim);
1870: if (degree >= 0) {
1871: PetscSpaceSetDegree(P, degree, PETSC_DETERMINE);
1872: if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1873: PetscSpace Pend, Pside;
1875: PetscSpaceCreate(comm, &Pend);
1876: PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL);
1877: PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE);
1878: PetscSpaceSetNumComponents(Pend, Nc);
1879: PetscSpaceSetNumVariables(Pend, dim - 1);
1880: PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE);
1881: PetscSpaceCreate(comm, &Pside);
1882: PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL);
1883: PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE);
1884: PetscSpaceSetNumComponents(Pside, 1);
1885: PetscSpaceSetNumVariables(Pside, 1);
1886: PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE);
1887: PetscSpaceSetType(P, PETSCSPACETENSOR);
1888: PetscSpaceTensorSetNumSubspaces(P, 2);
1889: PetscSpaceTensorSetSubspace(P, 0, Pend);
1890: PetscSpaceTensorSetSubspace(P, 1, Pside);
1891: PetscSpaceDestroy(&Pend);
1892: PetscSpaceDestroy(&Pside);
1893: }
1894: }
1895: if (setFromOptions) PetscSpaceSetFromOptions(P);
1896: PetscSpaceSetUp(P);
1897: PetscSpaceGetDegree(P, °ree, NULL);
1898: PetscSpacePolynomialGetTensor(P, &tensor);
1899: PetscSpaceGetNumComponents(P, &Nc);
1900: /* Create dual space */
1901: PetscDualSpaceCreate(comm, &Q);
1902: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1903: PetscObjectSetOptionsPrefix((PetscObject)Q, prefix);
1904: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1905: PetscDualSpaceSetDM(Q, K);
1906: DMDestroy(&K);
1907: PetscDualSpaceSetNumComponents(Q, Nc);
1908: PetscDualSpaceSetOrder(Q, degree);
1909: /* TODO For some reason, we need a tensor dualspace with wedges */
1910: PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE);
1911: if (setFromOptions) PetscDualSpaceSetFromOptions(Q);
1912: PetscDualSpaceSetUp(Q);
1913: /* Create quadrature */
1914: qorder = qorder >= 0 ? qorder : degree;
1915: if (setFromOptions) {
1916: PetscObjectOptionsBegin((PetscObject)P);
1917: PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0);
1918: PetscOptionsEnd();
1919: }
1920: PetscFECreateDefaultQuadrature_Private(dim, ct, qorder, &q, &fq);
1921: /* Create finite element */
1922: PetscFECreateFromSpaces(P, Q, q, fq, fem);
1923: if (setFromOptions) PetscFESetFromOptions(*fem);
1924: return 0;
1925: }
1927: /*@C
1928: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1930: Collective
1932: Input Parameters:
1933: + comm - The MPI comm
1934: . dim - The spatial dimension
1935: . Nc - The number of components
1936: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1937: . prefix - The options prefix, or NULL
1938: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1940: Output Parameter:
1941: . fem - The PetscFE object
1943: Note:
1944: 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.
1946: Level: beginner
1948: .seealso: `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1949: @*/
1950: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1951: {
1952: PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1953: return 0;
1954: }
1956: /*@C
1957: PetscFECreateByCell - Create a PetscFE for basic FEM computation
1959: Collective
1961: Input Parameters:
1962: + comm - The MPI comm
1963: . dim - The spatial dimension
1964: . Nc - The number of components
1965: . ct - The celltype of the reference cell
1966: . prefix - The options prefix, or NULL
1967: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1969: Output Parameter:
1970: . fem - The PetscFE object
1972: Note:
1973: 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.
1975: Level: beginner
1977: .seealso: `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1978: @*/
1979: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
1980: {
1981: PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem);
1982: return 0;
1983: }
1985: /*@
1986: PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k
1988: Collective
1990: Input Parameters:
1991: + comm - The MPI comm
1992: . dim - The spatial dimension
1993: . Nc - The number of components
1994: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1995: . k - The degree k of the space
1996: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1998: Output Parameter:
1999: . fem - The PetscFE object
2001: Level: beginner
2003: Notes:
2004: 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.
2006: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2007: @*/
2008: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2009: {
2010: PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem);
2011: return 0;
2012: }
2014: /*@
2015: PetscFECreateLagrangeByCell - Create a PetscFE for the basic Lagrange space of degree k
2017: Collective
2019: Input Parameters:
2020: + comm - The MPI comm
2021: . dim - The spatial dimension
2022: . Nc - The number of components
2023: . ct - The celltype of the reference cell
2024: . k - The degree k of the space
2025: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
2027: Output Parameter:
2028: . fem - The PetscFE object
2030: Level: beginner
2032: Notes:
2033: 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.
2035: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2036: @*/
2037: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2038: {
2039: PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem);
2040: return 0;
2041: }
2043: /*@C
2044: PetscFESetName - Names the FE and its subobjects
2046: Not collective
2048: Input Parameters:
2049: + fe - The PetscFE
2050: - name - The name
2052: Level: intermediate
2054: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2055: @*/
2056: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2057: {
2058: PetscSpace P;
2059: PetscDualSpace Q;
2061: PetscFEGetBasisSpace(fe, &P);
2062: PetscFEGetDualSpace(fe, &Q);
2063: PetscObjectSetName((PetscObject)fe, name);
2064: PetscObjectSetName((PetscObject)P, name);
2065: PetscObjectSetName((PetscObject)Q, name);
2066: return 0;
2067: }
2069: 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[])
2070: {
2071: PetscInt dOffset = 0, fOffset = 0, f, g;
2073: for (f = 0; f < Nf; ++f) {
2074: PetscFE fe;
2075: const PetscInt k = ds->jetDegree[f];
2076: const PetscInt cdim = T[f]->cdim;
2077: const PetscInt Nq = T[f]->Np;
2078: const PetscInt Nbf = T[f]->Nb;
2079: const PetscInt Ncf = T[f]->Nc;
2080: const PetscReal *Bq = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2081: const PetscReal *Dq = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2082: const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2083: PetscInt hOffset = 0, b, c, d;
2085: PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
2086: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2087: for (d = 0; d < cdim * Ncf; ++d) u_x[fOffset * cdim + d] = 0.0;
2088: for (b = 0; b < Nbf; ++b) {
2089: for (c = 0; c < Ncf; ++c) {
2090: const PetscInt cidx = b * Ncf + c;
2092: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2093: for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * cdim + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2094: }
2095: }
2096: if (k > 1) {
2097: for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * cdim;
2098: for (d = 0; d < cdim * cdim * Ncf; ++d) u_x[hOffset + fOffset * cdim * cdim + d] = 0.0;
2099: for (b = 0; b < Nbf; ++b) {
2100: for (c = 0; c < Ncf; ++c) {
2101: const PetscInt cidx = b * Ncf + c;
2103: for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * cdim * cdim + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2104: }
2105: }
2106: PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * cdim * cdim]);
2107: }
2108: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2109: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * cdim]);
2110: if (u_t) {
2111: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2112: for (b = 0; b < Nbf; ++b) {
2113: for (c = 0; c < Ncf; ++c) {
2114: const PetscInt cidx = b * Ncf + c;
2116: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2117: }
2118: }
2119: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2120: }
2121: fOffset += Ncf;
2122: dOffset += Nbf;
2123: }
2124: return 0;
2125: }
2127: 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[])
2128: {
2129: PetscInt dOffset = 0, fOffset = 0, f, g;
2131: /* f is the field number in the DS, g is the field number in u[] */
2132: for (f = 0, g = 0; f < Nf; ++f) {
2133: PetscFE fe = (PetscFE)ds->disc[f];
2134: const PetscInt dEt = T[f]->cdim;
2135: const PetscInt dE = fegeom->dimEmbed;
2136: const PetscInt Nq = T[f]->Np;
2137: const PetscInt Nbf = T[f]->Nb;
2138: const PetscInt Ncf = T[f]->Nc;
2139: const PetscReal *Bq = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2140: const PetscReal *Dq = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2141: PetscBool isCohesive;
2142: PetscInt Ns, s;
2144: if (!T[f]) continue;
2145: PetscDSGetCohesive(ds, f, &isCohesive);
2146: Ns = isCohesive ? 1 : 2;
2147: for (s = 0; s < Ns; ++s, ++g) {
2148: PetscInt b, c, d;
2150: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2151: for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2152: for (b = 0; b < Nbf; ++b) {
2153: for (c = 0; c < Ncf; ++c) {
2154: const PetscInt cidx = b * Ncf + c;
2156: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2157: for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2158: }
2159: }
2160: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2161: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]);
2162: if (u_t) {
2163: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2164: for (b = 0; b < Nbf; ++b) {
2165: for (c = 0; c < Ncf; ++c) {
2166: const PetscInt cidx = b * Ncf + c;
2168: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2169: }
2170: }
2171: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2172: }
2173: fOffset += Ncf;
2174: dOffset += Nbf;
2175: }
2176: }
2177: return 0;
2178: }
2180: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2181: {
2182: PetscFE fe;
2183: PetscTabulation Tc;
2184: PetscInt b, c;
2186: if (!prob) return 0;
2187: PetscDSGetDiscretization(prob, field, (PetscObject *)&fe);
2188: PetscFEGetFaceCentroidTabulation(fe, &Tc);
2189: {
2190: const PetscReal *faceBasis = Tc->T[0];
2191: const PetscInt Nb = Tc->Nb;
2192: const PetscInt Nc = Tc->Nc;
2194: for (c = 0; c < Nc; ++c) u[c] = 0.0;
2195: for (b = 0; b < Nb; ++b) {
2196: for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2197: }
2198: }
2199: return 0;
2200: }
2202: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2203: {
2204: PetscFEGeom pgeom;
2205: const PetscInt dEt = T->cdim;
2206: const PetscInt dE = fegeom->dimEmbed;
2207: const PetscInt Nq = T->Np;
2208: const PetscInt Nb = T->Nb;
2209: const PetscInt Nc = T->Nc;
2210: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2211: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt];
2212: PetscInt q, b, c, d;
2214: for (q = 0; q < Nq; ++q) {
2215: for (b = 0; b < Nb; ++b) {
2216: for (c = 0; c < Nc; ++c) {
2217: const PetscInt bcidx = b * Nc + c;
2219: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2220: for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d];
2221: for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0;
2222: }
2223: }
2224: PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);
2225: PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);
2226: PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer);
2227: for (b = 0; b < Nb; ++b) {
2228: for (c = 0; c < Nc; ++c) {
2229: const PetscInt bcidx = b * Nc + c;
2230: const PetscInt qcidx = q * Nc + c;
2232: elemVec[b] += tmpBasis[bcidx] * f0[qcidx];
2233: for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2234: }
2235: }
2236: }
2237: return (0);
2238: }
2240: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2241: {
2242: const PetscInt dE = T->cdim;
2243: const PetscInt Nq = T->Np;
2244: const PetscInt Nb = T->Nb;
2245: const PetscInt Nc = T->Nc;
2246: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2247: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];
2248: PetscInt q, b, c, d;
2250: for (q = 0; q < Nq; ++q) {
2251: for (b = 0; b < Nb; ++b) {
2252: for (c = 0; c < Nc; ++c) {
2253: const PetscInt bcidx = b * Nc + c;
2255: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2256: for (d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2257: }
2258: }
2259: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2260: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2261: for (b = 0; b < Nb; ++b) {
2262: for (c = 0; c < Nc; ++c) {
2263: const PetscInt bcidx = b * Nc + c;
2264: const PetscInt qcidx = q * Nc + c;
2266: elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2267: for (d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2268: }
2269: }
2270: }
2271: return (0);
2272: }
2274: 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[])
2275: {
2276: const PetscInt dE = TI->cdim;
2277: const PetscInt NqI = TI->Np;
2278: const PetscInt NbI = TI->Nb;
2279: const PetscInt NcI = TI->Nc;
2280: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2281: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2282: const PetscInt NqJ = TJ->Np;
2283: const PetscInt NbJ = TJ->Nb;
2284: const PetscInt NcJ = TJ->Nc;
2285: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2286: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2287: PetscInt f, fc, g, gc, df, dg;
2289: for (f = 0; f < NbI; ++f) {
2290: for (fc = 0; fc < NcI; ++fc) {
2291: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2293: tmpBasisI[fidx] = basisI[fidx];
2294: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2295: }
2296: }
2297: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2298: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2299: for (g = 0; g < NbJ; ++g) {
2300: for (gc = 0; gc < NcJ; ++gc) {
2301: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2303: tmpBasisJ[gidx] = basisJ[gidx];
2304: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2305: }
2306: }
2307: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2308: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
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 */
2312: const PetscInt i = offsetI + f; /* Element matrix row */
2313: for (g = 0; g < NbJ; ++g) {
2314: for (gc = 0; gc < NcJ; ++gc) {
2315: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2316: const PetscInt j = offsetJ + g; /* Element matrix column */
2317: const PetscInt fOff = eOffset + i * totDim + j;
2319: elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2320: for (df = 0; df < dE; ++df) {
2321: elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2322: elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2323: for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2324: }
2325: }
2326: }
2327: }
2328: }
2329: return (0);
2330: }
2332: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt s, 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[])
2333: {
2334: const PetscInt dE = TI->cdim;
2335: const PetscInt NqI = TI->Np;
2336: const PetscInt NbI = TI->Nb;
2337: const PetscInt NcI = TI->Nc;
2338: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2339: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2340: const PetscInt NqJ = TJ->Np;
2341: const PetscInt NbJ = TJ->Nb;
2342: const PetscInt NcJ = TJ->Nc;
2343: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2344: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2345: const PetscInt so = isHybridI ? 0 : s;
2346: const PetscInt to = isHybridJ ? 0 : s;
2347: PetscInt f, fc, g, gc, df, dg;
2349: for (f = 0; f < NbI; ++f) {
2350: for (fc = 0; fc < NcI; ++fc) {
2351: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2353: tmpBasisI[fidx] = basisI[fidx];
2354: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2355: }
2356: }
2357: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2358: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2359: for (g = 0; g < NbJ; ++g) {
2360: for (gc = 0; gc < NcJ; ++gc) {
2361: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2363: tmpBasisJ[gidx] = basisJ[gidx];
2364: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2365: }
2366: }
2367: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2368: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2369: for (f = 0; f < NbI; ++f) {
2370: for (fc = 0; fc < NcI; ++fc) {
2371: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2372: const PetscInt i = offsetI + NbI * so + f; /* Element matrix row */
2373: for (g = 0; g < NbJ; ++g) {
2374: for (gc = 0; gc < NcJ; ++gc) {
2375: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2376: const PetscInt j = offsetJ + NbJ * to + g; /* Element matrix column */
2377: const PetscInt fOff = eOffset + i * totDim + j;
2379: elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2380: for (df = 0; df < dE; ++df) {
2381: elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2382: elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2383: for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2384: }
2385: }
2386: }
2387: }
2388: }
2389: return (0);
2390: }
2392: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2393: {
2394: PetscDualSpace dsp;
2395: DM dm;
2396: PetscQuadrature quadDef;
2397: PetscInt dim, cdim, Nq;
2399: PetscFEGetDualSpace(fe, &dsp);
2400: PetscDualSpaceGetDM(dsp, &dm);
2401: DMGetDimension(dm, &dim);
2402: DMGetCoordinateDim(dm, &cdim);
2403: PetscFEGetQuadrature(fe, &quadDef);
2404: quad = quad ? quad : quadDef;
2405: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2406: PetscMalloc1(Nq * cdim, &cgeom->v);
2407: PetscMalloc1(Nq * cdim * cdim, &cgeom->J);
2408: PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ);
2409: PetscMalloc1(Nq, &cgeom->detJ);
2410: cgeom->dim = dim;
2411: cgeom->dimEmbed = cdim;
2412: cgeom->numCells = 1;
2413: cgeom->numPoints = Nq;
2414: DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2415: return 0;
2416: }
2418: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2419: {
2420: PetscFree(cgeom->v);
2421: PetscFree(cgeom->J);
2422: PetscFree(cgeom->invJ);
2423: PetscFree(cgeom->detJ);
2424: return 0;
2425: }