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 `PetscFEType`
61: Not Collective
63: Input Parameters:
64: + sname - The name of a new user-defined creation routine
65: - function - The creation routine
67: Example Usage:
68: .vb
69: PetscFERegister("my_fe", MyPetscFECreate);
70: .ve
72: Then, your PetscFE type can be chosen with the procedural interface via
73: .vb
74: PetscFECreate(MPI_Comm, PetscFE *);
75: PetscFESetType(PetscFE, "my_fe");
76: .ve
77: or at runtime via the option
78: .vb
79: -petscfe_type my_fe
80: .ve
82: Level: advanced
84: Note:
85: `PetscFERegister()` may be called multiple times to add several user-defined `PetscFE`s
87: .seealso: `PetscFE`, `PetscFEType`, `PetscFERegisterAll()`, `PetscFERegisterDestroy()`
88: @*/
89: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
90: {
91: PetscFunctionBegin;
92: PetscCall(PetscFunctionListAdd(&PetscFEList, sname, function));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@C
97: PetscFESetType - Builds a particular `PetscFE`
99: Collective
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: `PetscFEType`, `PetscFE`, `PetscFEGetType()`, `PetscFECreate()`
111: @*/
112: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
113: {
114: PetscErrorCode (*r)(PetscFE);
115: PetscBool match;
117: PetscFunctionBegin;
119: PetscCall(PetscObjectTypeCompare((PetscObject)fem, name, &match));
120: if (match) PetscFunctionReturn(PETSC_SUCCESS);
122: if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
123: PetscCall(PetscFunctionListFind(PetscFEList, name, &r));
124: PetscCheck(r, PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
126: PetscTryTypeMethod(fem, destroy);
127: fem->ops->destroy = NULL;
129: PetscCall((*r)(fem));
130: PetscCall(PetscObjectChangeTypeName((PetscObject)fem, name));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@C
135: PetscFEGetType - Gets the `PetscFEType` (as a string) from the `PetscFE` object.
137: Not Collective
139: Input Parameter:
140: . fem - The `PetscFE`
142: Output Parameter:
143: . name - The `PetscFEType` name
145: Level: intermediate
147: .seealso: `PetscFEType`, `PetscFE`, `PetscFESetType()`, `PetscFECreate()`
148: @*/
149: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
150: {
151: PetscFunctionBegin;
153: PetscAssertPointer(name, 2);
154: if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
155: *name = ((PetscObject)fem)->type_name;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@C
160: PetscFEViewFromOptions - View from a `PetscFE` based on values in the options database
162: Collective
164: Input Parameters:
165: + A - the `PetscFE` object
166: . obj - Optional object that provides the options prefix
167: - name - command line option name
169: Level: intermediate
171: .seealso: `PetscFE`, `PetscFEView()`, `PetscObjectViewFromOptions()`, `PetscFECreate()`
172: @*/
173: PetscErrorCode PetscFEViewFromOptions(PetscFE A, PetscObject obj, const char name[])
174: {
175: PetscFunctionBegin;
177: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@C
182: PetscFEView - Views a `PetscFE`
184: Collective
186: Input Parameters:
187: + fem - the `PetscFE` object to view
188: - viewer - the viewer
190: Level: beginner
192: .seealso: `PetscFE`, `PetscViewer`, `PetscFEDestroy()`, `PetscFEViewFromOptions()`
193: @*/
194: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
195: {
196: PetscBool iascii;
198: PetscFunctionBegin;
201: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fem), &viewer));
202: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer));
203: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
204: PetscTryTypeMethod(fem, view, viewer);
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: PetscFESetFromOptions - sets parameters in a `PetscFE` from the options database
211: Collective
213: Input Parameter:
214: . fem - the `PetscFE` object to set options for
216: Options Database Keys:
217: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
218: - -petscfe_num_batches - the number of cell batches to integrate serially
220: Level: intermediate
222: .seealso: `PetscFEV`, `PetscFEView()`
223: @*/
224: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
225: {
226: const char *defaultType;
227: char name[256];
228: PetscBool flg;
230: PetscFunctionBegin;
232: if (!((PetscObject)fem)->type_name) {
233: defaultType = PETSCFEBASIC;
234: } else {
235: defaultType = ((PetscObject)fem)->type_name;
236: }
237: if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
239: PetscObjectOptionsBegin((PetscObject)fem);
240: PetscCall(PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg));
241: if (flg) {
242: PetscCall(PetscFESetType(fem, name));
243: } else if (!((PetscObject)fem)->type_name) {
244: PetscCall(PetscFESetType(fem, defaultType));
245: }
246: PetscCall(PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL, 1));
247: PetscCall(PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL, 1));
248: PetscTryTypeMethod(fem, setfromoptions, PetscOptionsObject);
249: /* process any options handlers added with PetscObjectAddOptionsHandler() */
250: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fem, PetscOptionsObject));
251: PetscOptionsEnd();
252: PetscCall(PetscFEViewFromOptions(fem, NULL, "-petscfe_view"));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@C
257: PetscFESetUp - Construct data structures for the `PetscFE` after the `PetscFEType` has been set
259: Collective
261: Input Parameter:
262: . fem - the `PetscFE` object to setup
264: Level: intermediate
266: .seealso: `PetscFE`, `PetscFEView()`, `PetscFEDestroy()`
267: @*/
268: PetscErrorCode PetscFESetUp(PetscFE fem)
269: {
270: PetscFunctionBegin;
272: if (fem->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
273: PetscCall(PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0));
274: fem->setupcalled = PETSC_TRUE;
275: PetscTryTypeMethod(fem, setup);
276: PetscCall(PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@
281: PetscFEDestroy - Destroys a `PetscFE` object
283: Collective
285: Input Parameter:
286: . fem - the `PetscFE` object to destroy
288: Level: beginner
290: .seealso: `PetscFE`, `PetscFEView()`
291: @*/
292: PetscErrorCode PetscFEDestroy(PetscFE *fem)
293: {
294: PetscFunctionBegin;
295: if (!*fem) PetscFunctionReturn(PETSC_SUCCESS);
298: if (--((PetscObject)(*fem))->refct > 0) {
299: *fem = NULL;
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
302: ((PetscObject)(*fem))->refct = 0;
304: if ((*fem)->subspaces) {
305: PetscInt dim, d;
307: PetscCall(PetscDualSpaceGetDimension((*fem)->dualSpace, &dim));
308: for (d = 0; d < dim; ++d) PetscCall(PetscFEDestroy(&(*fem)->subspaces[d]));
309: }
310: PetscCall(PetscFree((*fem)->subspaces));
311: PetscCall(PetscFree((*fem)->invV));
312: PetscCall(PetscTabulationDestroy(&(*fem)->T));
313: PetscCall(PetscTabulationDestroy(&(*fem)->Tf));
314: PetscCall(PetscTabulationDestroy(&(*fem)->Tc));
315: PetscCall(PetscSpaceDestroy(&(*fem)->basisSpace));
316: PetscCall(PetscDualSpaceDestroy(&(*fem)->dualSpace));
317: PetscCall(PetscQuadratureDestroy(&(*fem)->quadrature));
318: PetscCall(PetscQuadratureDestroy(&(*fem)->faceQuadrature));
319: #ifdef PETSC_HAVE_LIBCEED
320: PetscCallCEED(CeedBasisDestroy(&(*fem)->ceedBasis));
321: PetscCallCEED(CeedDestroy(&(*fem)->ceed));
322: #endif
324: PetscTryTypeMethod((*fem), destroy);
325: PetscCall(PetscHeaderDestroy(fem));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*@
330: PetscFECreate - Creates an empty `PetscFE` object. The type can then be set with `PetscFESetType()`.
332: Collective
334: Input Parameter:
335: . comm - The communicator for the `PetscFE` object
337: Output Parameter:
338: . fem - The `PetscFE` object
340: Level: beginner
342: .seealso: `PetscFE`, `PetscFEType`, `PetscFESetType()`, `PetscFECreateDefault()`, `PETSCFEGALERKIN`
343: @*/
344: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
345: {
346: PetscFE f;
348: PetscFunctionBegin;
349: PetscAssertPointer(fem, 2);
350: PetscCall(PetscCitationsRegister(FECitation, &FEcite));
351: *fem = NULL;
352: PetscCall(PetscFEInitializePackage());
354: PetscCall(PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView));
356: f->basisSpace = NULL;
357: f->dualSpace = NULL;
358: f->numComponents = 1;
359: f->subspaces = NULL;
360: f->invV = NULL;
361: f->T = NULL;
362: f->Tf = NULL;
363: f->Tc = NULL;
364: PetscCall(PetscArrayzero(&f->quadrature, 1));
365: PetscCall(PetscArrayzero(&f->faceQuadrature, 1));
366: f->blockSize = 0;
367: f->numBlocks = 1;
368: f->batchSize = 0;
369: f->numBatches = 1;
371: *fem = f;
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /*@
376: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
378: Not Collective
380: Input Parameter:
381: . fem - The `PetscFE` object
383: Output Parameter:
384: . dim - The spatial dimension
386: Level: intermediate
388: .seealso: `PetscFE`, `PetscFECreate()`
389: @*/
390: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
391: {
392: DM dm;
394: PetscFunctionBegin;
396: PetscAssertPointer(dim, 2);
397: PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
398: PetscCall(DMGetDimension(dm, dim));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@
403: PetscFESetNumComponents - Sets the number of field components in the element
405: Not Collective
407: Input Parameters:
408: + fem - The `PetscFE` object
409: - comp - The number of field components
411: Level: intermediate
413: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`, `PetscFEGetNumComponents()`
414: @*/
415: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
416: {
417: PetscFunctionBegin;
419: fem->numComponents = comp;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /*@
424: PetscFEGetNumComponents - Returns the number of components in the element
426: Not Collective
428: Input Parameter:
429: . fem - The `PetscFE` object
431: Output Parameter:
432: . comp - The number of field components
434: Level: intermediate
436: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetSpatialDimension()`
437: @*/
438: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
439: {
440: PetscFunctionBegin;
442: PetscAssertPointer(comp, 2);
443: *comp = fem->numComponents;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*@
448: PetscFESetTileSizes - Sets the tile sizes for evaluation
450: Not Collective
452: Input Parameters:
453: + fem - The `PetscFE` object
454: . blockSize - The number of elements in a block
455: . numBlocks - The number of blocks in a batch
456: . batchSize - The number of elements in a batch
457: - numBatches - The number of batches in a chunk
459: Level: intermediate
461: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFEGetTileSizes()`
462: @*/
463: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
464: {
465: PetscFunctionBegin;
467: fem->blockSize = blockSize;
468: fem->numBlocks = numBlocks;
469: fem->batchSize = batchSize;
470: fem->numBatches = numBatches;
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /*@
475: PetscFEGetTileSizes - Returns the tile sizes for evaluation
477: Not Collective
479: Input Parameter:
480: . fem - The `PetscFE` object
482: Output Parameters:
483: + blockSize - The number of elements in a block
484: . numBlocks - The number of blocks in a batch
485: . batchSize - The number of elements in a batch
486: - numBatches - The number of batches in a chunk
488: Level: intermediate
490: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFESetTileSizes()`
491: @*/
492: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
493: {
494: PetscFunctionBegin;
496: if (blockSize) PetscAssertPointer(blockSize, 2);
497: if (numBlocks) PetscAssertPointer(numBlocks, 3);
498: if (batchSize) PetscAssertPointer(batchSize, 4);
499: if (numBatches) PetscAssertPointer(numBatches, 5);
500: if (blockSize) *blockSize = fem->blockSize;
501: if (numBlocks) *numBlocks = fem->numBlocks;
502: if (batchSize) *batchSize = fem->batchSize;
503: if (numBatches) *numBatches = fem->numBatches;
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@
508: PetscFEGetBasisSpace - Returns the `PetscSpace` used for the approximation of the solution for the `PetscFE`
510: Not Collective
512: Input Parameter:
513: . fem - The `PetscFE` object
515: Output Parameter:
516: . sp - The `PetscSpace` object
518: Level: intermediate
520: .seealso: `PetscFE`, `PetscSpace`, `PetscFECreate()`
521: @*/
522: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
523: {
524: PetscFunctionBegin;
526: PetscAssertPointer(sp, 2);
527: *sp = fem->basisSpace;
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: /*@
532: PetscFESetBasisSpace - Sets the `PetscSpace` used for the approximation of the solution
534: Not Collective
536: Input Parameters:
537: + fem - The `PetscFE` object
538: - sp - The `PetscSpace` object
540: Level: intermediate
542: Developer Notes:
543: There is `PetscFESetBasisSpace()` but the `PetscFESetDualSpace()`, likely the Basis is unneeded in the function name
545: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetDualSpace()`
546: @*/
547: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
548: {
549: PetscFunctionBegin;
552: PetscCall(PetscSpaceDestroy(&fem->basisSpace));
553: fem->basisSpace = sp;
554: PetscCall(PetscObjectReference((PetscObject)fem->basisSpace));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: PetscFEGetDualSpace - Returns the `PetscDualSpace` used to define the inner product for a `PetscFE`
561: Not Collective
563: Input Parameter:
564: . fem - The `PetscFE` object
566: Output Parameter:
567: . sp - The `PetscDualSpace` object
569: Level: intermediate
571: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
572: @*/
573: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
574: {
575: PetscFunctionBegin;
577: PetscAssertPointer(sp, 2);
578: *sp = fem->dualSpace;
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@
583: PetscFESetDualSpace - Sets the `PetscDualSpace` used to define the inner product
585: Not Collective
587: Input Parameters:
588: + fem - The `PetscFE` object
589: - sp - The `PetscDualSpace` object
591: Level: intermediate
593: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetBasisSpace()`
594: @*/
595: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
596: {
597: PetscFunctionBegin;
600: PetscCall(PetscDualSpaceDestroy(&fem->dualSpace));
601: fem->dualSpace = sp;
602: PetscCall(PetscObjectReference((PetscObject)fem->dualSpace));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@
607: PetscFEGetQuadrature - Returns the `PetscQuadrature` used to calculate inner products
609: Not Collective
611: Input Parameter:
612: . fem - The `PetscFE` object
614: Output Parameter:
615: . q - The `PetscQuadrature` object
617: Level: intermediate
619: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`
620: @*/
621: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
622: {
623: PetscFunctionBegin;
625: PetscAssertPointer(q, 2);
626: *q = fem->quadrature;
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: /*@
631: PetscFESetQuadrature - Sets the `PetscQuadrature` used to calculate inner products
633: Not Collective
635: Input Parameters:
636: + fem - The `PetscFE` object
637: - q - The `PetscQuadrature` object
639: Level: intermediate
641: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFEGetFaceQuadrature()`
642: @*/
643: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
644: {
645: PetscInt Nc, qNc;
647: PetscFunctionBegin;
649: if (q == fem->quadrature) PetscFunctionReturn(PETSC_SUCCESS);
650: PetscCall(PetscFEGetNumComponents(fem, &Nc));
651: PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
652: PetscCheck(!(qNc != 1) || !(Nc != qNc), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_SIZ, "FE components %" PetscInt_FMT " != Quadrature components %" PetscInt_FMT " and non-scalar quadrature", Nc, qNc);
653: PetscCall(PetscTabulationDestroy(&fem->T));
654: PetscCall(PetscTabulationDestroy(&fem->Tc));
655: PetscCall(PetscObjectReference((PetscObject)q));
656: PetscCall(PetscQuadratureDestroy(&fem->quadrature));
657: fem->quadrature = q;
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*@
662: PetscFEGetFaceQuadrature - Returns the `PetscQuadrature` used to calculate inner products on faces
664: Not Collective
666: Input Parameter:
667: . fem - The `PetscFE` object
669: Output Parameter:
670: . q - The `PetscQuadrature` object
672: Level: intermediate
674: Developer Notes:
675: There is a special face quadrature but not edge, likely this API would benefit from a refactorization
677: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
678: @*/
679: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
680: {
681: PetscFunctionBegin;
683: PetscAssertPointer(q, 2);
684: *q = fem->faceQuadrature;
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: /*@
689: PetscFESetFaceQuadrature - Sets the `PetscQuadrature` used to calculate inner products on faces
691: Not Collective
693: Input Parameters:
694: + fem - The `PetscFE` object
695: - q - The `PetscQuadrature` object
697: Level: intermediate
699: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`
700: @*/
701: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
702: {
703: PetscInt Nc, qNc;
705: PetscFunctionBegin;
707: if (q == fem->faceQuadrature) PetscFunctionReturn(PETSC_SUCCESS);
708: PetscCall(PetscFEGetNumComponents(fem, &Nc));
709: PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
710: PetscCheck(!(qNc != 1) || !(Nc != qNc), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_SIZ, "FE components %" PetscInt_FMT " != Quadrature components %" PetscInt_FMT " and non-scalar quadrature", Nc, qNc);
711: PetscCall(PetscTabulationDestroy(&fem->Tf));
712: PetscCall(PetscObjectReference((PetscObject)q));
713: PetscCall(PetscQuadratureDestroy(&fem->faceQuadrature));
714: fem->faceQuadrature = q;
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: /*@
719: PetscFECopyQuadrature - Copy both volumetric and surface quadrature to a new `PetscFE`
721: Not Collective
723: Input Parameters:
724: + sfe - The `PetscFE` source for the quadratures
725: - tfe - The `PetscFE` target for the quadratures
727: Level: intermediate
729: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
730: @*/
731: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
732: {
733: PetscQuadrature q;
735: PetscFunctionBegin;
738: PetscCall(PetscFEGetQuadrature(sfe, &q));
739: PetscCall(PetscFESetQuadrature(tfe, q));
740: PetscCall(PetscFEGetFaceQuadrature(sfe, &q));
741: PetscCall(PetscFESetFaceQuadrature(tfe, q));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: /*@C
746: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
748: Not Collective
750: Input Parameter:
751: . fem - The `PetscFE` object
753: Output Parameter:
754: . numDof - Array with the number of dofs per dimension
756: Level: intermediate
758: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
759: @*/
760: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
761: {
762: PetscFunctionBegin;
764: PetscAssertPointer(numDof, 2);
765: PetscCall(PetscDualSpaceGetNumDof(fem->dualSpace, numDof));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /*@C
770: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
772: Not Collective
774: Input Parameters:
775: + fem - The `PetscFE` object
776: - k - The highest derivative we need to tabulate, very often 1
778: Output Parameter:
779: . T - The basis function values and derivatives at quadrature points
781: Level: intermediate
783: Note:
784: .vb
785: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
786: 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
787: 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
788: .ve
790: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
791: @*/
792: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
793: {
794: PetscInt npoints;
795: const PetscReal *points;
797: PetscFunctionBegin;
799: PetscAssertPointer(T, 3);
800: PetscCall(PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL));
801: if (!fem->T) PetscCall(PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T));
802: PetscCheck(!fem->T || k <= fem->T->K || (!fem->T->cdim && !fem->T->K), PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %" PetscInt_FMT " derivatives, but only tabulated %" PetscInt_FMT, k, fem->T->K);
803: *T = fem->T;
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*@C
808: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
810: Not Collective
812: Input Parameters:
813: + fem - The `PetscFE` object
814: - k - The highest derivative we need to tabulate, very often 1
816: Output Parameter:
817: . Tf - The basis function values and derivatives at face quadrature points
819: Level: intermediate
821: Note:
822: .vb
823: 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
824: 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
825: 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
826: .ve
828: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
829: @*/
830: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
831: {
832: PetscFunctionBegin;
834: PetscAssertPointer(Tf, 3);
835: if (!fem->Tf) {
836: const PetscReal xi0[3] = {-1., -1., -1.};
837: PetscReal v0[3], J[9], detJ;
838: PetscQuadrature fq;
839: PetscDualSpace sp;
840: DM dm;
841: const PetscInt *faces;
842: PetscInt dim, numFaces, f, npoints, q;
843: const PetscReal *points;
844: PetscReal *facePoints;
846: PetscCall(PetscFEGetDualSpace(fem, &sp));
847: PetscCall(PetscDualSpaceGetDM(sp, &dm));
848: PetscCall(DMGetDimension(dm, &dim));
849: PetscCall(DMPlexGetConeSize(dm, 0, &numFaces));
850: PetscCall(DMPlexGetCone(dm, 0, &faces));
851: PetscCall(PetscFEGetFaceQuadrature(fem, &fq));
852: if (fq) {
853: PetscCall(PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL));
854: PetscCall(PetscMalloc1(numFaces * npoints * dim, &facePoints));
855: for (f = 0; f < numFaces; ++f) {
856: PetscCall(DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ));
857: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * npoints + q) * dim]);
858: }
859: PetscCall(PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf));
860: PetscCall(PetscFree(facePoints));
861: }
862: }
863: PetscCheck(!fem->Tf || k <= fem->Tf->K, PetscObjectComm((PetscObject)fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %" PetscInt_FMT " derivatives, but only tabulated %" PetscInt_FMT, k, fem->Tf->K);
864: *Tf = fem->Tf;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@C
869: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
871: Not Collective
873: Input Parameter:
874: . fem - The `PetscFE` object
876: Output Parameter:
877: . Tc - The basis function values at face centroid points
879: Level: intermediate
881: Note:
882: .vb
883: T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
884: .ve
886: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
887: @*/
888: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
889: {
890: PetscFunctionBegin;
892: PetscAssertPointer(Tc, 2);
893: if (!fem->Tc) {
894: PetscDualSpace sp;
895: DM dm;
896: const PetscInt *cone;
897: PetscReal *centroids;
898: PetscInt dim, numFaces, f;
900: PetscCall(PetscFEGetDualSpace(fem, &sp));
901: PetscCall(PetscDualSpaceGetDM(sp, &dm));
902: PetscCall(DMGetDimension(dm, &dim));
903: PetscCall(DMPlexGetConeSize(dm, 0, &numFaces));
904: PetscCall(DMPlexGetCone(dm, 0, &cone));
905: PetscCall(PetscMalloc1(numFaces * dim, ¢roids));
906: for (f = 0; f < numFaces; ++f) PetscCall(DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f * dim], NULL));
907: PetscCall(PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc));
908: PetscCall(PetscFree(centroids));
909: }
910: *Tc = fem->Tc;
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: /*@C
915: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
917: Not Collective
919: Input Parameters:
920: + fem - The `PetscFE` object
921: . nrepl - The number of replicas
922: . npoints - The number of tabulation points in a replica
923: . points - The tabulation point coordinates
924: - K - The number of derivatives calculated
926: Output Parameter:
927: . T - The basis function values and derivatives at tabulation points
929: Level: intermediate
931: Note:
932: .vb
933: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
934: 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
935: T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis
936: T->function i, component c, in directions d and e
937: .ve
939: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
940: @*/
941: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
942: {
943: DM dm;
944: PetscDualSpace Q;
945: PetscInt Nb; /* Dimension of FE space P */
946: PetscInt Nc; /* Field components */
947: PetscInt cdim; /* Reference coordinate dimension */
948: PetscInt k;
950: PetscFunctionBegin;
951: if (!npoints || !fem->dualSpace || K < 0) {
952: *T = NULL;
953: PetscFunctionReturn(PETSC_SUCCESS);
954: }
956: PetscAssertPointer(points, 4);
957: PetscAssertPointer(T, 6);
958: PetscCall(PetscFEGetDualSpace(fem, &Q));
959: PetscCall(PetscDualSpaceGetDM(Q, &dm));
960: PetscCall(DMGetDimension(dm, &cdim));
961: PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
962: PetscCall(PetscFEGetNumComponents(fem, &Nc));
963: PetscCall(PetscMalloc1(1, T));
964: (*T)->K = !cdim ? 0 : K;
965: (*T)->Nr = nrepl;
966: (*T)->Np = npoints;
967: (*T)->Nb = Nb;
968: (*T)->Nc = Nc;
969: (*T)->cdim = cdim;
970: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
971: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscCalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
972: PetscUseTypeMethod(fem, createtabulation, nrepl * npoints, points, K, *T);
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@C
977: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
979: Not Collective
981: Input Parameters:
982: + fem - The `PetscFE` object
983: . npoints - The number of tabulation points
984: . points - The tabulation point coordinates
985: . K - The number of derivatives calculated
986: - T - An existing tabulation object with enough allocated space
988: Output Parameter:
989: . T - The basis function values and derivatives at tabulation points
991: Level: intermediate
993: Note:
994: .vb
995: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
996: 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
997: 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
998: .ve
1000: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`
1001: @*/
1002: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1003: {
1004: PetscFunctionBeginHot;
1005: if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(PETSC_SUCCESS);
1007: PetscAssertPointer(points, 3);
1008: PetscAssertPointer(T, 5);
1009: if (PetscDefined(USE_DEBUG)) {
1010: DM dm;
1011: PetscDualSpace Q;
1012: PetscInt Nb; /* Dimension of FE space P */
1013: PetscInt Nc; /* Field components */
1014: PetscInt cdim; /* Reference coordinate dimension */
1016: PetscCall(PetscFEGetDualSpace(fem, &Q));
1017: PetscCall(PetscDualSpaceGetDM(Q, &dm));
1018: PetscCall(DMGetDimension(dm, &cdim));
1019: PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
1020: PetscCall(PetscFEGetNumComponents(fem, &Nc));
1021: PetscCheck(T->K == (!cdim ? 0 : K), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %" PetscInt_FMT " must match requested K %" PetscInt_FMT, T->K, !cdim ? 0 : K);
1022: PetscCheck(T->Nb == Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %" PetscInt_FMT " must match requested Nb %" PetscInt_FMT, T->Nb, Nb);
1023: PetscCheck(T->Nc == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %" PetscInt_FMT " must match requested Nc %" PetscInt_FMT, T->Nc, Nc);
1024: PetscCheck(T->cdim == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %" PetscInt_FMT " must match requested cdim %" PetscInt_FMT, T->cdim, cdim);
1025: }
1026: T->Nr = 1;
1027: T->Np = npoints;
1028: PetscUseTypeMethod(fem, createtabulation, npoints, points, K, T);
1029: PetscFunctionReturn(PETSC_SUCCESS);
1030: }
1032: /*@C
1033: PetscTabulationDestroy - Frees memory from the associated tabulation.
1035: Not Collective
1037: Input Parameter:
1038: . T - The tabulation
1040: Level: intermediate
1042: .seealso: `PetscTabulation`, `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
1043: @*/
1044: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1045: {
1046: PetscInt k;
1048: PetscFunctionBegin;
1049: PetscAssertPointer(T, 1);
1050: if (!T || !(*T)) PetscFunctionReturn(PETSC_SUCCESS);
1051: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscFree((*T)->T[k]));
1052: PetscCall(PetscFree((*T)->T));
1053: PetscCall(PetscFree(*T));
1054: *T = NULL;
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: static PetscErrorCode PetscFECreatePointTraceDefault_Internal(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1059: {
1060: PetscSpace bsp, bsubsp;
1061: PetscDualSpace dsp, dsubsp;
1062: PetscInt dim, depth, numComp, i, j, coneSize, order;
1063: DM dm;
1064: DMLabel label;
1065: PetscReal *xi, *v, *J, detJ;
1066: const char *name;
1067: PetscQuadrature origin, fullQuad, subQuad;
1069: PetscFunctionBegin;
1070: PetscCall(PetscFEGetBasisSpace(fe, &bsp));
1071: PetscCall(PetscFEGetDualSpace(fe, &dsp));
1072: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1073: PetscCall(DMGetDimension(dm, &dim));
1074: PetscCall(DMPlexGetDepthLabel(dm, &label));
1075: PetscCall(DMLabelGetValue(label, refPoint, &depth));
1076: PetscCall(PetscCalloc1(depth, &xi));
1077: PetscCall(PetscMalloc1(dim, &v));
1078: PetscCall(PetscMalloc1(dim * dim, &J));
1079: for (i = 0; i < depth; i++) xi[i] = 0.;
1080: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &origin));
1081: PetscCall(PetscQuadratureSetData(origin, depth, 0, 1, xi, NULL));
1082: PetscCall(DMPlexComputeCellGeometryFEM(dm, refPoint, origin, v, J, NULL, &detJ));
1083: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1084: for (i = 1; i < dim; i++) {
1085: for (j = 0; j < depth; j++) J[i * depth + j] = J[i * dim + j];
1086: }
1087: PetscCall(PetscQuadratureDestroy(&origin));
1088: PetscCall(PetscDualSpaceGetPointSubspace(dsp, refPoint, &dsubsp));
1089: PetscCall(PetscSpaceCreateSubspace(bsp, dsubsp, v, J, NULL, NULL, PETSC_OWN_POINTER, &bsubsp));
1090: PetscCall(PetscSpaceSetUp(bsubsp));
1091: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), trFE));
1092: PetscCall(PetscFESetType(*trFE, PETSCFEBASIC));
1093: PetscCall(PetscFEGetNumComponents(fe, &numComp));
1094: PetscCall(PetscFESetNumComponents(*trFE, numComp));
1095: PetscCall(PetscFESetBasisSpace(*trFE, bsubsp));
1096: PetscCall(PetscFESetDualSpace(*trFE, dsubsp));
1097: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1098: if (name) PetscCall(PetscFESetName(*trFE, name));
1099: PetscCall(PetscFEGetQuadrature(fe, &fullQuad));
1100: PetscCall(PetscQuadratureGetOrder(fullQuad, &order));
1101: PetscCall(DMPlexGetConeSize(dm, refPoint, &coneSize));
1102: if (coneSize == 2 * depth) PetscCall(PetscDTGaussTensorQuadrature(depth, 1, (order + 2) / 2, -1., 1., &subQuad));
1103: else PetscCall(PetscDTSimplexQuadrature(depth, order, PETSCDTSIMPLEXQUAD_DEFAULT, &subQuad));
1104: PetscCall(PetscFESetQuadrature(*trFE, subQuad));
1105: PetscCall(PetscFESetUp(*trFE));
1106: PetscCall(PetscQuadratureDestroy(&subQuad));
1107: PetscCall(PetscSpaceDestroy(&bsubsp));
1108: PetscFunctionReturn(PETSC_SUCCESS);
1109: }
1111: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1112: {
1113: PetscErrorCode (*createpointtrace)(PetscFE, PetscInt, PetscFE *);
1115: PetscFunctionBegin;
1117: PetscAssertPointer(trFE, 3);
1118: createpointtrace = fe->ops->createpointtrace;
1119: if (createpointtrace) {
1120: PetscCall((*createpointtrace)(fe, refPoint, trFE));
1121: } else {
1122: PetscCall(PetscFECreatePointTraceDefault_Internal(fe, refPoint, trFE));
1123: }
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1128: {
1129: PetscInt hStart, hEnd;
1130: PetscDualSpace dsp;
1131: DM dm;
1133: PetscFunctionBegin;
1135: PetscAssertPointer(trFE, 3);
1136: *trFE = NULL;
1137: PetscCall(PetscFEGetDualSpace(fe, &dsp));
1138: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1139: PetscCall(DMPlexGetHeightStratum(dm, height, &hStart, &hEnd));
1140: if (hEnd <= hStart) PetscFunctionReturn(PETSC_SUCCESS);
1141: PetscCall(PetscFECreatePointTrace(fe, hStart, trFE));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: /*@
1146: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1148: Not Collective
1150: Input Parameter:
1151: . fem - The `PetscFE`
1153: Output Parameter:
1154: . dim - The dimension
1156: Level: intermediate
1158: .seealso: `PetscFE`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1159: @*/
1160: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1161: {
1162: PetscFunctionBegin;
1164: PetscAssertPointer(dim, 2);
1165: PetscTryTypeMethod(fem, getdimension, dim);
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: /*@C
1170: PetscFEPushforward - Map the reference element function to real space
1172: Input Parameters:
1173: + fe - The `PetscFE`
1174: . fegeom - The cell geometry
1175: . Nv - The number of function values
1176: - vals - The function values
1178: Output Parameter:
1179: . vals - The transformed function values
1181: Level: advanced
1183: Notes:
1184: This just forwards the call onto `PetscDualSpacePushforward()`.
1186: It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1188: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscDualSpacePushforward()`
1189: @*/
1190: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1191: {
1192: PetscFunctionBeginHot;
1193: PetscCall(PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: /*@C
1198: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1200: Input Parameters:
1201: + fe - The `PetscFE`
1202: . fegeom - The cell geometry
1203: . Nv - The number of function gradient values
1204: - vals - The function gradient values
1206: Output Parameter:
1207: . vals - The transformed function gradient values
1209: Level: advanced
1211: Notes:
1212: This just forwards the call onto `PetscDualSpacePushforwardGradient()`.
1214: It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1216: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1217: @*/
1218: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1219: {
1220: PetscFunctionBeginHot;
1221: PetscCall(PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: /*@C
1226: PetscFEPushforwardHessian - Map the reference element function Hessian to real space
1228: Input Parameters:
1229: + fe - The `PetscFE`
1230: . fegeom - The cell geometry
1231: . Nv - The number of function Hessian values
1232: - vals - The function Hessian values
1234: Output Parameter:
1235: . vals - The transformed function Hessian values
1237: Level: advanced
1239: Notes:
1240: This just forwards the call onto `PetscDualSpacePushforwardHessian()`.
1242: It only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1244: Developer Notes:
1245: It is unclear why all these one line convenience routines are desirable
1247: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1248: @*/
1249: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1250: {
1251: PetscFunctionBeginHot;
1252: PetscCall(PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*
1257: Purpose: Compute element vector for chunk of elements
1259: Input:
1260: Sizes:
1261: Ne: number of elements
1262: Nf: number of fields
1263: PetscFE
1264: dim: spatial dimension
1265: Nb: number of basis functions
1266: Nc: number of field components
1267: PetscQuadrature
1268: Nq: number of quadrature points
1270: Geometry:
1271: PetscFEGeom[Ne] possibly *Nq
1272: PetscReal v0s[dim]
1273: PetscReal n[dim]
1274: PetscReal jacobians[dim*dim]
1275: PetscReal jacobianInverses[dim*dim]
1276: PetscReal jacobianDeterminants
1277: FEM:
1278: PetscFE
1279: PetscQuadrature
1280: PetscReal quadPoints[Nq*dim]
1281: PetscReal quadWeights[Nq]
1282: PetscReal basis[Nq*Nb*Nc]
1283: PetscReal basisDer[Nq*Nb*Nc*dim]
1284: PetscScalar coefficients[Ne*Nb*Nc]
1285: PetscScalar elemVec[Ne*Nb*Nc]
1287: Problem:
1288: PetscInt f: the active field
1289: f0, f1
1291: Work Space:
1292: PetscFE
1293: PetscScalar f0[Nq*dim];
1294: PetscScalar f1[Nq*dim*dim];
1295: PetscScalar u[Nc];
1296: PetscScalar gradU[Nc*dim];
1297: PetscReal x[dim];
1298: PetscScalar realSpaceDer[dim];
1300: Purpose: Compute element vector for N_cb batches of elements
1302: Input:
1303: Sizes:
1304: N_cb: Number of serial cell batches
1306: Geometry:
1307: PetscReal v0s[Ne*dim]
1308: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1309: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1310: PetscReal jacobianDeterminants[Ne] possibly *Nq
1311: FEM:
1312: static PetscReal quadPoints[Nq*dim]
1313: static PetscReal quadWeights[Nq]
1314: static PetscReal basis[Nq*Nb*Nc]
1315: static PetscReal basisDer[Nq*Nb*Nc*dim]
1316: PetscScalar coefficients[Ne*Nb*Nc]
1317: PetscScalar elemVec[Ne*Nb*Nc]
1319: ex62.c:
1320: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1321: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1322: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1323: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1325: ex52.c:
1326: 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)
1327: 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)
1329: ex52_integrateElement.cu
1330: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1332: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1333: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1334: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1336: ex52_integrateElementOpenCL.c:
1337: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1338: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1339: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1341: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1342: */
1344: /*@C
1345: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1347: Not Collective
1349: Input Parameters:
1350: + prob - The `PetscDS` specifying the discretizations and continuum functions
1351: . field - The field being integrated
1352: . Ne - The number of elements in the chunk
1353: . cgeom - The cell geometry for each cell in the chunk
1354: . coefficients - The array of FEM basis coefficients for the elements
1355: . probAux - The `PetscDS` specifying the auxiliary discretizations
1356: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1358: Output Parameter:
1359: . integral - the integral for this field
1361: Level: intermediate
1363: Developer Notes:
1364: The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.
1366: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrateBd()`
1367: @*/
1368: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1369: {
1370: PetscFE fe;
1372: PetscFunctionBegin;
1374: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1375: if (fe->ops->integrate) PetscCall((*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral));
1376: PetscFunctionReturn(PETSC_SUCCESS);
1377: }
1379: /*@C
1380: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1382: Not Collective
1384: Input Parameters:
1385: + prob - The `PetscDS` specifying the discretizations and continuum functions
1386: . field - The field being integrated
1387: . obj_func - The function to be integrated
1388: . Ne - The number of elements in the chunk
1389: . geom - The face geometry for each face in the chunk
1390: . coefficients - The array of FEM basis coefficients for the elements
1391: . probAux - The `PetscDS` specifying the auxiliary discretizations
1392: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1394: Output Parameter:
1395: . integral - the integral for this field
1397: Level: intermediate
1399: Developer Notes:
1400: The function name begins with `PetscFE` and yet the first argument is `PetscDS` and it has no `PetscFE` arguments.
1402: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrate()`
1403: @*/
1404: 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[])
1405: {
1406: PetscFE fe;
1408: PetscFunctionBegin;
1410: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1411: if (fe->ops->integratebd) PetscCall((*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral));
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: /*@C
1416: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1418: Not Collective
1420: Input Parameters:
1421: + ds - The `PetscDS` specifying the discretizations and continuum functions
1422: . key - The (label+value, field) being integrated
1423: . Ne - The number of elements in the chunk
1424: . cgeom - The cell geometry for each cell in the chunk
1425: . coefficients - The array of FEM basis coefficients for the elements
1426: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1427: . probAux - The `PetscDS` specifying the auxiliary discretizations
1428: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1429: - t - The time
1431: Output Parameter:
1432: . elemVec - the element residual vectors from each element
1434: Level: intermediate
1436: Note:
1437: .vb
1438: Loop over batch of elements (e):
1439: Loop over quadrature points (q):
1440: Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1441: Call f_0 and f_1
1442: Loop over element vector entries (f,fc --> i):
1443: elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1444: .ve
1446: .seealso: `PetscFEIntegrateBdResidual()`
1447: @*/
1448: 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[])
1449: {
1450: PetscFE fe;
1452: PetscFunctionBeginHot;
1454: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1455: if (fe->ops->integrateresidual) PetscCall((*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1456: PetscFunctionReturn(PETSC_SUCCESS);
1457: }
1459: /*@C
1460: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1462: Not Collective
1464: Input Parameters:
1465: + ds - The `PetscDS` specifying the discretizations and continuum functions
1466: . wf - The PetscWeakForm object holding the pointwise functions
1467: . key - The (label+value, field) being integrated
1468: . Ne - The number of elements in the chunk
1469: . fgeom - The face geometry for each cell in the chunk
1470: . coefficients - The array of FEM basis coefficients for the elements
1471: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1472: . probAux - The `PetscDS` specifying the auxiliary discretizations
1473: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1474: - t - The time
1476: Output Parameter:
1477: . elemVec - the element residual vectors from each element
1479: Level: intermediate
1481: .seealso: `PetscFEIntegrateResidual()`
1482: @*/
1483: 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[])
1484: {
1485: PetscFE fe;
1487: PetscFunctionBegin;
1489: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1490: if (fe->ops->integratebdresidual) PetscCall((*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1491: PetscFunctionReturn(PETSC_SUCCESS);
1492: }
1494: /*@C
1495: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1497: Not Collective
1499: Input Parameters:
1500: + ds - The `PetscDS` specifying the discretizations and continuum functions
1501: . dsIn - The `PetscDS` specifying the discretizations and continuum functions for input
1502: . key - The (label+value, field) being integrated
1503: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1504: . Ne - The number of elements in the chunk
1505: . fgeom - The face geometry for each cell in the chunk
1506: . coefficients - The array of FEM basis coefficients for the elements
1507: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1508: . probAux - The `PetscDS` specifying the auxiliary discretizations
1509: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1510: - t - The time
1512: Output Parameter:
1513: . elemVec - the element residual vectors from each element
1515: Level: developer
1517: .seealso: `PetscFEIntegrateResidual()`
1518: @*/
1519: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1520: {
1521: PetscFE fe;
1523: PetscFunctionBegin;
1526: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1527: if (fe->ops->integratehybridresidual) PetscCall((*fe->ops->integratehybridresidual)(ds, dsIn, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1528: PetscFunctionReturn(PETSC_SUCCESS);
1529: }
1531: /*@C
1532: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1534: Not Collective
1536: Input Parameters:
1537: + ds - The `PetscDS` specifying the discretizations and continuum functions
1538: . jtype - The type of matrix pointwise functions that should be used
1539: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1540: . Ne - The number of elements in the chunk
1541: . cgeom - The cell geometry for each cell in the chunk
1542: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1543: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1544: . probAux - The `PetscDS` specifying the auxiliary discretizations
1545: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1546: . t - The time
1547: - u_tshift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1549: Output Parameter:
1550: . elemMat - the element matrices for the Jacobian from each element
1552: Level: intermediate
1554: Note:
1555: .vb
1556: Loop over batch of elements (e):
1557: Loop over element matrix entries (f,fc,g,gc --> i,j):
1558: Loop over quadrature points (q):
1559: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1560: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1561: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1562: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1563: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1564: .ve
1566: .seealso: `PetscFEIntegrateResidual()`
1567: @*/
1568: 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[])
1569: {
1570: PetscFE fe;
1571: PetscInt Nf;
1573: PetscFunctionBegin;
1575: PetscCall(PetscDSGetNumFields(ds, &Nf));
1576: PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1577: if (fe->ops->integratejacobian) PetscCall((*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1578: PetscFunctionReturn(PETSC_SUCCESS);
1579: }
1581: /*@C
1582: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1584: Not Collective
1586: Input Parameters:
1587: + ds - The `PetscDS` specifying the discretizations and continuum functions
1588: . wf - The PetscWeakForm holding the pointwise functions
1589: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1590: . Ne - The number of elements in the chunk
1591: . fgeom - The face geometry for each cell in the chunk
1592: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1593: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1594: . probAux - The `PetscDS` specifying the auxiliary discretizations
1595: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1596: . t - The time
1597: - u_tshift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1599: Output Parameter:
1600: . elemMat - the element matrices for the Jacobian from each element
1602: Level: intermediate
1604: Note:
1605: .vb
1606: Loop over batch of elements (e):
1607: Loop over element matrix entries (f,fc,g,gc --> i,j):
1608: Loop over quadrature points (q):
1609: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1610: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1611: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1612: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1613: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1614: .ve
1616: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1617: @*/
1618: 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[])
1619: {
1620: PetscFE fe;
1621: PetscInt Nf;
1623: PetscFunctionBegin;
1625: PetscCall(PetscDSGetNumFields(ds, &Nf));
1626: PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1627: if (fe->ops->integratebdjacobian) PetscCall((*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: /*@C
1632: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1634: Not Collective
1636: Input Parameters:
1637: + ds - The `PetscDS` specifying the discretizations and continuum functions for the output
1638: . dsIn - The `PetscDS` specifying the discretizations and continuum functions for the input
1639: . jtype - The type of matrix pointwise functions that should be used
1640: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1641: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1642: . Ne - The number of elements in the chunk
1643: . fgeom - The face geometry for each cell in the chunk
1644: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1645: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1646: . probAux - The `PetscDS` specifying the auxiliary discretizations
1647: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1648: . t - The time
1649: - u_tshift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1651: Output Parameter:
1652: . elemMat - the element matrices for the Jacobian from each element
1654: Level: developer
1656: Note:
1657: .vb
1658: Loop over batch of elements (e):
1659: Loop over element matrix entries (f,fc,g,gc --> i,j):
1660: Loop over quadrature points (q):
1661: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1662: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1663: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1664: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1665: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1666: .ve
1668: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1669: @*/
1670: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscDS dsIn, 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[])
1671: {
1672: PetscFE fe;
1673: PetscInt Nf;
1675: PetscFunctionBegin;
1677: PetscCall(PetscDSGetNumFields(ds, &Nf));
1678: PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1679: if (fe->ops->integratehybridjacobian) PetscCall((*fe->ops->integratehybridjacobian)(ds, dsIn, jtype, key, s, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: /*@
1684: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1686: Input Parameters:
1687: + fe - The finite element space
1688: - height - The height of the `DMPLEX` point
1690: Output Parameter:
1691: . subfe - The subspace of this `PetscFE` space
1693: Level: advanced
1695: Note:
1696: For example, if we want the subspace of this space for a face, we would choose height = 1.
1698: .seealso: `PetscFECreateDefault()`
1699: @*/
1700: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1701: {
1702: PetscSpace P, subP;
1703: PetscDualSpace Q, subQ;
1704: PetscQuadrature subq;
1705: PetscInt dim, Nc;
1707: PetscFunctionBegin;
1709: PetscAssertPointer(subfe, 3);
1710: if (height == 0) {
1711: *subfe = fe;
1712: PetscFunctionReturn(PETSC_SUCCESS);
1713: }
1714: PetscCall(PetscFEGetBasisSpace(fe, &P));
1715: PetscCall(PetscFEGetDualSpace(fe, &Q));
1716: PetscCall(PetscFEGetNumComponents(fe, &Nc));
1717: PetscCall(PetscFEGetFaceQuadrature(fe, &subq));
1718: PetscCall(PetscDualSpaceGetDimension(Q, &dim));
1719: PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
1720: if (!fe->subspaces) PetscCall(PetscCalloc1(dim, &fe->subspaces));
1721: if (height <= dim) {
1722: if (!fe->subspaces[height - 1]) {
1723: PetscFE sub = NULL;
1724: const char *name;
1726: PetscCall(PetscSpaceGetHeightSubspace(P, height, &subP));
1727: PetscCall(PetscDualSpaceGetHeightSubspace(Q, height, &subQ));
1728: if (subQ) {
1729: PetscCall(PetscObjectReference((PetscObject)subP));
1730: PetscCall(PetscObjectReference((PetscObject)subQ));
1731: PetscCall(PetscObjectReference((PetscObject)subq));
1732: PetscCall(PetscFECreateFromSpaces(subP, subQ, subq, NULL, &sub));
1733: }
1734: if (sub) {
1735: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1736: if (name) PetscCall(PetscFESetName(sub, name));
1737: }
1738: fe->subspaces[height - 1] = sub;
1739: }
1740: *subfe = fe->subspaces[height - 1];
1741: } else {
1742: *subfe = NULL;
1743: }
1744: PetscFunctionReturn(PETSC_SUCCESS);
1745: }
1747: /*@
1748: PetscFERefine - Create a "refined" `PetscFE` object that refines the reference cell into
1749: smaller copies.
1751: Collective
1753: Input Parameter:
1754: . fe - The initial `PetscFE`
1756: Output Parameter:
1757: . feRef - The refined `PetscFE`
1759: Level: advanced
1761: Notes:
1762: This is typically used to generate a preconditioner for a higher order method from a lower order method on a
1763: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1764: interpolation between regularly refined meshes.
1766: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1767: @*/
1768: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1769: {
1770: PetscSpace P, Pref;
1771: PetscDualSpace Q, Qref;
1772: DM K, Kref;
1773: PetscQuadrature q, qref;
1774: const PetscReal *v0, *jac;
1775: PetscInt numComp, numSubelements;
1776: PetscInt cStart, cEnd, c;
1777: PetscDualSpace *cellSpaces;
1779: PetscFunctionBegin;
1780: PetscCall(PetscFEGetBasisSpace(fe, &P));
1781: PetscCall(PetscFEGetDualSpace(fe, &Q));
1782: PetscCall(PetscFEGetQuadrature(fe, &q));
1783: PetscCall(PetscDualSpaceGetDM(Q, &K));
1784: /* Create space */
1785: PetscCall(PetscObjectReference((PetscObject)P));
1786: Pref = P;
1787: /* Create dual space */
1788: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1789: PetscCall(PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED));
1790: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref));
1791: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1792: PetscCall(DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd));
1793: PetscCall(PetscMalloc1(cEnd - cStart, &cellSpaces));
1794: /* TODO: fix for non-uniform refinement */
1795: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1796: PetscCall(PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces));
1797: PetscCall(PetscFree(cellSpaces));
1798: PetscCall(DMDestroy(&Kref));
1799: PetscCall(PetscDualSpaceSetUp(Qref));
1800: /* Create element */
1801: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), feRef));
1802: PetscCall(PetscFESetType(*feRef, PETSCFECOMPOSITE));
1803: PetscCall(PetscFESetBasisSpace(*feRef, Pref));
1804: PetscCall(PetscFESetDualSpace(*feRef, Qref));
1805: PetscCall(PetscFEGetNumComponents(fe, &numComp));
1806: PetscCall(PetscFESetNumComponents(*feRef, numComp));
1807: PetscCall(PetscFESetUp(*feRef));
1808: PetscCall(PetscSpaceDestroy(&Pref));
1809: PetscCall(PetscDualSpaceDestroy(&Qref));
1810: /* Create quadrature */
1811: PetscCall(PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL));
1812: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1813: PetscCall(PetscFESetQuadrature(*feRef, qref));
1814: PetscCall(PetscQuadratureDestroy(&qref));
1815: PetscFunctionReturn(PETSC_SUCCESS);
1816: }
1818: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1819: {
1820: PetscSpace P;
1821: PetscDualSpace Q;
1822: DM K;
1823: DMPolytopeType ct;
1824: PetscInt degree;
1825: char name[64];
1827: PetscFunctionBegin;
1828: PetscCall(PetscFEGetBasisSpace(fe, &P));
1829: PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
1830: PetscCall(PetscFEGetDualSpace(fe, &Q));
1831: PetscCall(PetscDualSpaceGetDM(Q, &K));
1832: PetscCall(DMPlexGetCellType(K, 0, &ct));
1833: switch (ct) {
1834: case DM_POLYTOPE_SEGMENT:
1835: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1836: case DM_POLYTOPE_QUADRILATERAL:
1837: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1838: case DM_POLYTOPE_HEXAHEDRON:
1839: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1840: PetscCall(PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree));
1841: break;
1842: case DM_POLYTOPE_TRIANGLE:
1843: case DM_POLYTOPE_TETRAHEDRON:
1844: PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree));
1845: break;
1846: case DM_POLYTOPE_TRI_PRISM:
1847: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1848: PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree));
1849: break;
1850: default:
1851: PetscCall(PetscSNPrintf(name, sizeof(name), "FE"));
1852: }
1853: PetscCall(PetscFESetName(fe, name));
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: /*@
1858: PetscFECreateFromSpaces - Create a `PetscFE` from the basis and dual spaces
1860: Collective
1862: Input Parameters:
1863: + P - The basis space
1864: . Q - The dual space
1865: . q - The cell quadrature
1866: - fq - The face quadrature
1868: Output Parameter:
1869: . fem - The `PetscFE` object
1871: Level: beginner
1873: Note:
1874: 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.
1876: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`,
1877: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1878: @*/
1879: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1880: {
1881: PetscInt Nc;
1882: PetscInt p_Ns = -1, p_Nc = -1, q_Ns = -1, q_Nc = -1;
1883: PetscBool p_is_uniform_sum = PETSC_FALSE, p_interleave_basis = PETSC_FALSE, p_interleave_components = PETSC_FALSE;
1884: PetscBool q_is_uniform_sum = PETSC_FALSE, q_interleave_basis = PETSC_FALSE, q_interleave_components = PETSC_FALSE;
1885: const char *prefix;
1887: PetscFunctionBegin;
1888: PetscCall(PetscObjectTypeCompare((PetscObject)P, PETSCSPACESUM, &p_is_uniform_sum));
1889: if (p_is_uniform_sum) {
1890: PetscSpace subsp_0 = NULL;
1891: PetscCall(PetscSpaceSumGetNumSubspaces(P, &p_Ns));
1892: PetscCall(PetscSpaceGetNumComponents(P, &p_Nc));
1893: PetscCall(PetscSpaceSumGetConcatenate(P, &p_is_uniform_sum));
1894: PetscCall(PetscSpaceSumGetInterleave(P, &p_interleave_basis, &p_interleave_components));
1895: for (PetscInt s = 0; s < p_Ns; s++) {
1896: PetscSpace subsp;
1898: PetscCall(PetscSpaceSumGetSubspace(P, s, &subsp));
1899: if (!s) {
1900: subsp_0 = subsp;
1901: } else if (subsp != subsp_0) {
1902: p_is_uniform_sum = PETSC_FALSE;
1903: }
1904: }
1905: }
1906: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &q_is_uniform_sum));
1907: if (q_is_uniform_sum) {
1908: PetscDualSpace subsp_0 = NULL;
1909: PetscCall(PetscDualSpaceSumGetNumSubspaces(Q, &q_Ns));
1910: PetscCall(PetscDualSpaceGetNumComponents(Q, &q_Nc));
1911: PetscCall(PetscDualSpaceSumGetConcatenate(Q, &q_is_uniform_sum));
1912: PetscCall(PetscDualSpaceSumGetInterleave(Q, &q_interleave_basis, &q_interleave_components));
1913: for (PetscInt s = 0; s < q_Ns; s++) {
1914: PetscDualSpace subsp;
1916: PetscCall(PetscDualSpaceSumGetSubspace(Q, s, &subsp));
1917: if (!s) {
1918: subsp_0 = subsp;
1919: } else if (subsp != subsp_0) {
1920: q_is_uniform_sum = PETSC_FALSE;
1921: }
1922: }
1923: }
1924: if (p_is_uniform_sum && q_is_uniform_sum && (p_interleave_basis == q_interleave_basis) && (p_interleave_components == q_interleave_components) && (p_Ns == q_Ns) && (p_Nc == q_Nc)) {
1925: PetscSpace scalar_space;
1926: PetscDualSpace scalar_dspace;
1927: PetscFE scalar_fe;
1929: PetscCall(PetscSpaceSumGetSubspace(P, 0, &scalar_space));
1930: PetscCall(PetscDualSpaceSumGetSubspace(Q, 0, &scalar_dspace));
1931: PetscCall(PetscObjectReference((PetscObject)scalar_space));
1932: PetscCall(PetscObjectReference((PetscObject)scalar_dspace));
1933: PetscCall(PetscObjectReference((PetscObject)q));
1934: PetscCall(PetscObjectReference((PetscObject)fq));
1935: PetscCall(PetscFECreateFromSpaces(scalar_space, scalar_dspace, q, fq, &scalar_fe));
1936: PetscCall(PetscFECreateVector(scalar_fe, p_Ns, p_interleave_basis, p_interleave_components, fem));
1937: PetscCall(PetscFEDestroy(&scalar_fe));
1938: } else {
1939: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)P), fem));
1940: PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
1941: }
1942: PetscCall(PetscSpaceGetNumComponents(P, &Nc));
1943: PetscCall(PetscFESetNumComponents(*fem, Nc));
1944: PetscCall(PetscFESetBasisSpace(*fem, P));
1945: PetscCall(PetscFESetDualSpace(*fem, Q));
1946: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)P, &prefix));
1947: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
1948: PetscCall(PetscFESetUp(*fem));
1949: PetscCall(PetscSpaceDestroy(&P));
1950: PetscCall(PetscDualSpaceDestroy(&Q));
1951: PetscCall(PetscFESetQuadrature(*fem, q));
1952: PetscCall(PetscFESetFaceQuadrature(*fem, fq));
1953: PetscCall(PetscQuadratureDestroy(&q));
1954: PetscCall(PetscQuadratureDestroy(&fq));
1955: PetscCall(PetscFESetDefaultName_Private(*fem));
1956: PetscFunctionReturn(PETSC_SUCCESS);
1957: }
1959: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
1960: {
1961: DM K;
1962: PetscSpace P;
1963: PetscDualSpace Q;
1964: PetscQuadrature q, fq;
1965: PetscBool tensor;
1967: PetscFunctionBegin;
1968: if (prefix) PetscAssertPointer(prefix, 5);
1969: PetscAssertPointer(fem, 9);
1970: switch (ct) {
1971: case DM_POLYTOPE_SEGMENT:
1972: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1973: case DM_POLYTOPE_QUADRILATERAL:
1974: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1975: case DM_POLYTOPE_HEXAHEDRON:
1976: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1977: tensor = PETSC_TRUE;
1978: break;
1979: default:
1980: tensor = PETSC_FALSE;
1981: }
1982: /* Create space */
1983: PetscCall(PetscSpaceCreate(comm, &P));
1984: PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
1985: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
1986: PetscCall(PetscSpacePolynomialSetTensor(P, tensor));
1987: PetscCall(PetscSpaceSetNumComponents(P, Nc));
1988: PetscCall(PetscSpaceSetNumVariables(P, dim));
1989: if (degree >= 0) {
1990: PetscCall(PetscSpaceSetDegree(P, degree, PETSC_DETERMINE));
1991: if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
1992: PetscSpace Pend, Pside;
1994: PetscCall(PetscSpaceSetNumComponents(P, 1));
1995: PetscCall(PetscSpaceCreate(comm, &Pend));
1996: PetscCall(PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL));
1997: PetscCall(PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE));
1998: PetscCall(PetscSpaceSetNumComponents(Pend, 1));
1999: PetscCall(PetscSpaceSetNumVariables(Pend, dim - 1));
2000: PetscCall(PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE));
2001: PetscCall(PetscSpaceCreate(comm, &Pside));
2002: PetscCall(PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL));
2003: PetscCall(PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE));
2004: PetscCall(PetscSpaceSetNumComponents(Pside, 1));
2005: PetscCall(PetscSpaceSetNumVariables(Pside, 1));
2006: PetscCall(PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE));
2007: PetscCall(PetscSpaceSetType(P, PETSCSPACETENSOR));
2008: PetscCall(PetscSpaceTensorSetNumSubspaces(P, 2));
2009: PetscCall(PetscSpaceTensorSetSubspace(P, 0, Pend));
2010: PetscCall(PetscSpaceTensorSetSubspace(P, 1, Pside));
2011: PetscCall(PetscSpaceDestroy(&Pend));
2012: PetscCall(PetscSpaceDestroy(&Pside));
2014: if (Nc > 1) {
2015: PetscSpace scalar_P = P;
2017: PetscCall(PetscSpaceCreate(comm, &P));
2018: PetscCall(PetscSpaceSetNumVariables(P, dim));
2019: PetscCall(PetscSpaceSetNumComponents(P, Nc));
2020: PetscCall(PetscSpaceSetType(P, PETSCSPACESUM));
2021: PetscCall(PetscSpaceSumSetNumSubspaces(P, Nc));
2022: PetscCall(PetscSpaceSumSetConcatenate(P, PETSC_TRUE));
2023: PetscCall(PetscSpaceSumSetInterleave(P, PETSC_TRUE, PETSC_FALSE));
2024: for (PetscInt i = 0; i < Nc; i++) PetscCall(PetscSpaceSumSetSubspace(P, i, scalar_P));
2025: PetscCall(PetscSpaceDestroy(&scalar_P));
2026: }
2027: }
2028: }
2029: if (setFromOptions) PetscCall(PetscSpaceSetFromOptions(P));
2030: PetscCall(PetscSpaceSetUp(P));
2031: PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
2032: PetscCall(PetscSpacePolynomialGetTensor(P, &tensor));
2033: PetscCall(PetscSpaceGetNumComponents(P, &Nc));
2034: /* Create dual space */
2035: PetscCall(PetscDualSpaceCreate(comm, &Q));
2036: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
2037: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
2038: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
2039: PetscCall(PetscDualSpaceSetDM(Q, K));
2040: PetscCall(DMDestroy(&K));
2041: PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
2042: PetscCall(PetscDualSpaceSetOrder(Q, degree));
2043: PetscCall(PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE));
2044: if (setFromOptions) PetscCall(PetscDualSpaceSetFromOptions(Q));
2045: PetscCall(PetscDualSpaceSetUp(Q));
2046: /* Create quadrature */
2047: qorder = qorder >= 0 ? qorder : degree;
2048: if (setFromOptions) {
2049: PetscObjectOptionsBegin((PetscObject)P);
2050: PetscCall(PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0));
2051: PetscOptionsEnd();
2052: }
2053: PetscCall(PetscDTCreateDefaultQuadrature(ct, qorder, &q, &fq));
2054: /* Create finite element */
2055: PetscCall(PetscFECreateFromSpaces(P, Q, q, fq, fem));
2056: if (setFromOptions) PetscCall(PetscFESetFromOptions(*fem));
2057: PetscFunctionReturn(PETSC_SUCCESS);
2058: }
2060: /*@C
2061: PetscFECreateDefault - Create a `PetscFE` for basic FEM computation
2063: Collective
2065: Input Parameters:
2066: + comm - The MPI comm
2067: . dim - The spatial dimension
2068: . Nc - The number of components
2069: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2070: . prefix - The options prefix, or `NULL`
2071: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2073: Output Parameter:
2074: . fem - The `PetscFE` object
2076: Level: beginner
2078: Note:
2079: 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.
2081: .seealso: `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2082: @*/
2083: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
2084: {
2085: PetscFunctionBegin;
2086: PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2087: PetscFunctionReturn(PETSC_SUCCESS);
2088: }
2090: /*@C
2091: PetscFECreateByCell - Create a `PetscFE` for basic FEM computation
2093: Collective
2095: Input Parameters:
2096: + comm - The MPI comm
2097: . dim - The spatial dimension
2098: . Nc - The number of components
2099: . ct - The celltype of the reference cell
2100: . prefix - The options prefix, or `NULL`
2101: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2103: Output Parameter:
2104: . fem - The `PetscFE` object
2106: Level: beginner
2108: Note:
2109: 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.
2111: .seealso: `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2112: @*/
2113: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
2114: {
2115: PetscFunctionBegin;
2116: PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }
2120: /*@
2121: PetscFECreateLagrange - Create a `PetscFE` for the basic Lagrange space of degree k
2123: Collective
2125: Input Parameters:
2126: + comm - The MPI comm
2127: . dim - The spatial dimension
2128: . Nc - The number of components
2129: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2130: . k - The degree k of the space
2131: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2133: Output Parameter:
2134: . fem - The `PetscFE` object
2136: Level: beginner
2138: Note:
2139: 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.
2141: .seealso: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2142: @*/
2143: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2144: {
2145: PetscFunctionBegin;
2146: PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem));
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: /*@
2151: PetscFECreateLagrangeByCell - Create a `PetscFE` for the basic Lagrange space of degree k
2153: Collective
2155: Input Parameters:
2156: + comm - The MPI comm
2157: . dim - The spatial dimension
2158: . Nc - The number of components
2159: . ct - The celltype of the reference cell
2160: . k - The degree k of the space
2161: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2163: Output Parameter:
2164: . fem - The `PetscFE` object
2166: Level: beginner
2168: Note:
2169: 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.
2171: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2172: @*/
2173: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2174: {
2175: PetscFunctionBegin;
2176: PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem));
2177: PetscFunctionReturn(PETSC_SUCCESS);
2178: }
2180: /*@C
2181: PetscFESetName - Names the `PetscFE` and its subobjects
2183: Not Collective
2185: Input Parameters:
2186: + fe - The `PetscFE`
2187: - name - The name
2189: Level: intermediate
2191: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2192: @*/
2193: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2194: {
2195: PetscSpace P;
2196: PetscDualSpace Q;
2198: PetscFunctionBegin;
2199: PetscCall(PetscFEGetBasisSpace(fe, &P));
2200: PetscCall(PetscFEGetDualSpace(fe, &Q));
2201: PetscCall(PetscObjectSetName((PetscObject)fe, name));
2202: PetscCall(PetscObjectSetName((PetscObject)P, name));
2203: PetscCall(PetscObjectSetName((PetscObject)Q, name));
2204: PetscFunctionReturn(PETSC_SUCCESS);
2205: }
2207: 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[])
2208: {
2209: PetscInt dOffset = 0, fOffset = 0, f, g;
2211: for (f = 0; f < Nf; ++f) {
2212: PetscCheck(r < T[f]->Nr, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", r, T[f]->Nr);
2213: PetscCheck(q < T[f]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", q, T[f]->Np);
2214: PetscFE fe;
2215: const PetscInt k = ds->jetDegree[f];
2216: const PetscInt cdim = T[f]->cdim;
2217: const PetscInt dE = fegeom->dimEmbed;
2218: const PetscInt Nq = T[f]->Np;
2219: const PetscInt Nbf = T[f]->Nb;
2220: const PetscInt Ncf = T[f]->Nc;
2221: const PetscReal *Bq = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2222: const PetscReal *Dq = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2223: const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2224: PetscInt hOffset = 0, b, c, d;
2226: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
2227: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2228: for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2229: for (b = 0; b < Nbf; ++b) {
2230: for (c = 0; c < Ncf; ++c) {
2231: const PetscInt cidx = b * Ncf + c;
2233: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2234: for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2235: }
2236: }
2237: if (k > 1) {
2238: for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * dE;
2239: for (d = 0; d < dE * dE * Ncf; ++d) u_x[hOffset + fOffset * dE * dE + d] = 0.0;
2240: for (b = 0; b < Nbf; ++b) {
2241: for (c = 0; c < Ncf; ++c) {
2242: const PetscInt cidx = b * Ncf + c;
2244: for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * dE * dE + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2245: }
2246: }
2247: PetscCall(PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * dE * dE]));
2248: }
2249: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset]));
2250: PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]));
2251: if (u_t) {
2252: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2253: for (b = 0; b < Nbf; ++b) {
2254: for (c = 0; c < Ncf; ++c) {
2255: const PetscInt cidx = b * Ncf + c;
2257: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2258: }
2259: }
2260: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2261: }
2262: fOffset += Ncf;
2263: dOffset += Nbf;
2264: }
2265: return PETSC_SUCCESS;
2266: }
2268: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt rc, PetscInt qc, PetscTabulation Tab[], const PetscInt rf[], const PetscInt qf[], PetscTabulation Tabf[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2269: {
2270: PetscInt dOffset = 0, fOffset = 0, f, g;
2272: /* f is the field number in the DS, g is the field number in u[] */
2273: for (f = 0, g = 0; f < Nf; ++f) {
2274: PetscBool isCohesive;
2275: PetscInt Ns, s;
2277: if (!Tab[f]) continue;
2278: PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
2279: Ns = isCohesive ? 1 : 2;
2280: {
2281: PetscTabulation T = isCohesive ? Tab[f] : Tabf[f];
2282: PetscFE fe = (PetscFE)ds->disc[f];
2283: const PetscInt dEt = T->cdim;
2284: const PetscInt dE = fegeom->dimEmbed;
2285: const PetscInt Nq = T->Np;
2286: const PetscInt Nbf = T->Nb;
2287: const PetscInt Ncf = T->Nc;
2289: for (s = 0; s < Ns; ++s, ++g) {
2290: const PetscInt r = isCohesive ? rc : rf[s];
2291: const PetscInt q = isCohesive ? qc : qf[s];
2292: const PetscReal *Bq = &T->T[0][(r * Nq + q) * Nbf * Ncf];
2293: const PetscReal *Dq = &T->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2294: PetscInt b, c, d;
2296: PetscCheck(r < T->Nr, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " Side %" PetscInt_FMT " Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", f, s, r, T->Nr);
2297: PetscCheck(q < T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " Side %" PetscInt_FMT " Point number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", f, s, q, T->Np);
2298: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2299: for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2300: for (b = 0; b < Nbf; ++b) {
2301: for (c = 0; c < Ncf; ++c) {
2302: const PetscInt cidx = b * Ncf + c;
2304: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2305: for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2306: }
2307: }
2308: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset]));
2309: PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]));
2310: if (u_t) {
2311: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2312: for (b = 0; b < Nbf; ++b) {
2313: for (c = 0; c < Ncf; ++c) {
2314: const PetscInt cidx = b * Ncf + c;
2316: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2317: }
2318: }
2319: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2320: }
2321: fOffset += Ncf;
2322: dOffset += Nbf;
2323: }
2324: }
2325: }
2326: return PETSC_SUCCESS;
2327: }
2329: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2330: {
2331: PetscFE fe;
2332: PetscTabulation Tc;
2333: PetscInt b, c;
2335: if (!prob) return PETSC_SUCCESS;
2336: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
2337: PetscCall(PetscFEGetFaceCentroidTabulation(fe, &Tc));
2338: {
2339: const PetscReal *faceBasis = Tc->T[0];
2340: const PetscInt Nb = Tc->Nb;
2341: const PetscInt Nc = Tc->Nc;
2343: for (c = 0; c < Nc; ++c) u[c] = 0.0;
2344: for (b = 0; b < Nb; ++b) {
2345: for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2346: }
2347: }
2348: return PETSC_SUCCESS;
2349: }
2351: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2352: {
2353: PetscFEGeom pgeom;
2354: const PetscInt dEt = T->cdim;
2355: const PetscInt dE = fegeom->dimEmbed;
2356: const PetscInt Nq = T->Np;
2357: const PetscInt Nb = T->Nb;
2358: const PetscInt Nc = T->Nc;
2359: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2360: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt];
2361: PetscInt q, b, c, d;
2363: for (q = 0; q < Nq; ++q) {
2364: for (b = 0; b < Nb; ++b) {
2365: for (c = 0; c < Nc; ++c) {
2366: const PetscInt bcidx = b * Nc + c;
2368: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2369: for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d];
2370: for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0;
2371: }
2372: }
2373: PetscCall(PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom));
2374: PetscCall(PetscFEPushforward(fe, &pgeom, Nb, tmpBasis));
2375: PetscCall(PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer));
2376: for (b = 0; b < Nb; ++b) {
2377: for (c = 0; c < Nc; ++c) {
2378: const PetscInt bcidx = b * Nc + c;
2379: const PetscInt qcidx = q * Nc + c;
2381: elemVec[b] += tmpBasis[bcidx] * f0[qcidx];
2382: for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2383: }
2384: }
2385: }
2386: return PETSC_SUCCESS;
2387: }
2389: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt s, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2390: {
2391: const PetscInt dE = T->cdim;
2392: const PetscInt Nq = T->Np;
2393: const PetscInt Nb = T->Nb;
2394: const PetscInt Nc = T->Nc;
2395: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2396: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];
2397: PetscInt q, b, c, d;
2399: for (q = 0; q < Nq; ++q) {
2400: for (b = 0; b < Nb; ++b) {
2401: for (c = 0; c < Nc; ++c) {
2402: const PetscInt bcidx = b * Nc + c;
2404: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2405: for (d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2406: }
2407: }
2408: PetscCall(PetscFEPushforward(fe, fegeom, Nb, tmpBasis));
2409: // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2410: // PetscCall(PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer));
2411: for (b = 0; b < Nb; ++b) {
2412: for (c = 0; c < Nc; ++c) {
2413: const PetscInt bcidx = b * Nc + c;
2414: const PetscInt qcidx = q * Nc + c;
2416: elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2417: for (d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2418: }
2419: }
2420: }
2421: return PETSC_SUCCESS;
2422: }
2424: 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[])
2425: {
2426: const PetscInt cdim = TI->cdim;
2427: const PetscInt dE = fegeom->dimEmbed;
2428: const PetscInt NqI = TI->Np;
2429: const PetscInt NbI = TI->Nb;
2430: const PetscInt NcI = TI->Nc;
2431: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2432: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * cdim];
2433: const PetscInt NqJ = TJ->Np;
2434: const PetscInt NbJ = TJ->Nb;
2435: const PetscInt NcJ = TJ->Nc;
2436: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2437: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * cdim];
2438: PetscInt f, fc, g, gc, df, dg;
2440: for (f = 0; f < NbI; ++f) {
2441: for (fc = 0; fc < NcI; ++fc) {
2442: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2444: tmpBasisI[fidx] = basisI[fidx];
2445: for (df = 0; df < cdim; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * cdim + df];
2446: }
2447: }
2448: PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2449: PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2450: for (g = 0; g < NbJ; ++g) {
2451: for (gc = 0; gc < NcJ; ++gc) {
2452: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2454: tmpBasisJ[gidx] = basisJ[gidx];
2455: for (dg = 0; dg < cdim; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * cdim + dg];
2456: }
2457: }
2458: PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2459: PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2460: for (f = 0; f < NbI; ++f) {
2461: for (fc = 0; fc < NcI; ++fc) {
2462: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2463: const PetscInt i = offsetI + f; /* Element matrix row */
2464: for (g = 0; g < NbJ; ++g) {
2465: for (gc = 0; gc < NcJ; ++gc) {
2466: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2467: const PetscInt j = offsetJ + g; /* Element matrix column */
2468: const PetscInt fOff = eOffset + i * totDim + j;
2470: elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2471: for (df = 0; df < dE; ++df) {
2472: elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2473: elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2474: for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2475: }
2476: }
2477: }
2478: }
2479: }
2480: return PETSC_SUCCESS;
2481: }
2483: 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[])
2484: {
2485: const PetscInt dE = TI->cdim;
2486: const PetscInt NqI = TI->Np;
2487: const PetscInt NbI = TI->Nb;
2488: const PetscInt NcI = TI->Nc;
2489: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2490: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2491: const PetscInt NqJ = TJ->Np;
2492: const PetscInt NbJ = TJ->Nb;
2493: const PetscInt NcJ = TJ->Nc;
2494: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2495: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2496: const PetscInt so = isHybridI ? 0 : s;
2497: const PetscInt to = isHybridJ ? 0 : s;
2498: PetscInt f, fc, g, gc, df, dg;
2500: for (f = 0; f < NbI; ++f) {
2501: for (fc = 0; fc < NcI; ++fc) {
2502: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2504: tmpBasisI[fidx] = basisI[fidx];
2505: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2506: }
2507: }
2508: PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2509: PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2510: for (g = 0; g < NbJ; ++g) {
2511: for (gc = 0; gc < NcJ; ++gc) {
2512: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2514: tmpBasisJ[gidx] = basisJ[gidx];
2515: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2516: }
2517: }
2518: PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2519: // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2520: // PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2521: for (f = 0; f < NbI; ++f) {
2522: for (fc = 0; fc < NcI; ++fc) {
2523: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2524: const PetscInt i = offsetI + NbI * so + f; /* Element matrix row */
2525: for (g = 0; g < NbJ; ++g) {
2526: for (gc = 0; gc < NcJ; ++gc) {
2527: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2528: const PetscInt j = offsetJ + NbJ * to + g; /* Element matrix column */
2529: const PetscInt fOff = eOffset + i * totDim + j;
2531: elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2532: for (df = 0; df < dE; ++df) {
2533: elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2534: elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2535: for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2536: }
2537: }
2538: }
2539: }
2540: }
2541: return PETSC_SUCCESS;
2542: }
2544: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2545: {
2546: PetscDualSpace dsp;
2547: DM dm;
2548: PetscQuadrature quadDef;
2549: PetscInt dim, cdim, Nq;
2551: PetscFunctionBegin;
2552: PetscCall(PetscFEGetDualSpace(fe, &dsp));
2553: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
2554: PetscCall(DMGetDimension(dm, &dim));
2555: PetscCall(DMGetCoordinateDim(dm, &cdim));
2556: PetscCall(PetscFEGetQuadrature(fe, &quadDef));
2557: quad = quad ? quad : quadDef;
2558: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
2559: PetscCall(PetscMalloc1(Nq * cdim, &cgeom->v));
2560: PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->J));
2561: PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ));
2562: PetscCall(PetscMalloc1(Nq, &cgeom->detJ));
2563: cgeom->dim = dim;
2564: cgeom->dimEmbed = cdim;
2565: cgeom->numCells = 1;
2566: cgeom->numPoints = Nq;
2567: PetscCall(DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ));
2568: PetscFunctionReturn(PETSC_SUCCESS);
2569: }
2571: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2572: {
2573: PetscFunctionBegin;
2574: PetscCall(PetscFree(cgeom->v));
2575: PetscCall(PetscFree(cgeom->J));
2576: PetscCall(PetscFree(cgeom->invJ));
2577: PetscCall(PetscFree(cgeom->detJ));
2578: PetscFunctionReturn(PETSC_SUCCESS);
2579: }