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: /*@
59: PetscFERegister - Adds a new `PetscFEType`
61: Not Collective, No Fortran Support
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()`
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: /*@
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 (basic|opencl|composite|vector) - Sets the `PetscFEType`
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: /*@
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: /*@
160: PetscFEViewFromOptions - View 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, pass `NULL` to use the options prefix of `A`
167: - name - command line option name
169: Options Database Key:
170: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
172: Level: intermediate
174: Note:
175: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
176: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
178: .seealso: `PetscFE`, `PetscFEView()`, `PetscObjectViewFromOptions()`, `PetscFECreate()`, `PetscOptionsCreateViewer()`
179: @*/
180: PetscErrorCode PetscFEViewFromOptions(PetscFE A, PeOp PetscObject obj, const char name[])
181: {
182: PetscFunctionBegin;
184: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@
189: PetscFEView - Views a `PetscFE`
191: Collective
193: Input Parameters:
194: + fem - the `PetscFE` object to view
195: - viewer - the viewer
197: Level: beginner
199: .seealso: `PetscFE`, `PetscViewer`, `PetscFEDestroy()`, `PetscFEViewFromOptions()`
200: @*/
201: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
202: {
203: PetscBool isascii;
205: PetscFunctionBegin;
208: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fem), &viewer));
209: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer));
210: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
211: PetscTryTypeMethod(fem, view, viewer);
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: PetscFESetFromOptions - sets parameters in a `PetscFE` from the options database
218: Collective
220: Input Parameter:
221: . fem - the `PetscFE` object to set options for
223: Options Database Keys:
224: + -petscfe_num_blocks nblocks - the number of cell blocks to integrate concurrently
225: - -petscfe_num_batches nbatches - the number of cell batches to integrate serially
227: Level: intermediate
229: .seealso: `PetscFE`, `PetscFEView()`
230: @*/
231: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
232: {
233: const char *defaultType;
234: char name[256];
235: PetscBool flg;
237: PetscFunctionBegin;
239: if (!((PetscObject)fem)->type_name) defaultType = PETSCFEBASIC;
240: else defaultType = ((PetscObject)fem)->type_name;
241: if (!PetscFERegisterAllCalled) PetscCall(PetscFERegisterAll());
243: PetscObjectOptionsBegin((PetscObject)fem);
244: PetscCall(PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, sizeof(name), &flg));
245: if (flg) PetscCall(PetscFESetType(fem, name));
246: else if (!((PetscObject)fem)->type_name) PetscCall(PetscFESetType(fem, defaultType));
247: PetscCall(PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL, 1));
248: PetscCall(PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL, 1));
249: PetscTryTypeMethod(fem, setfromoptions, PetscOptionsObject);
250: /* process any options handlers added with PetscObjectAddOptionsHandler() */
251: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fem, PetscOptionsObject));
252: PetscOptionsEnd();
253: PetscCall(PetscFEViewFromOptions(fem, NULL, "-petscfe_view"));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@
258: PetscFESetUp - Construct data structures for the `PetscFE` after the `PetscFEType` has been set
260: Collective
262: Input Parameter:
263: . fem - the `PetscFE` object to setup
265: Level: intermediate
267: .seealso: `PetscFE`, `PetscFEView()`, `PetscFEDestroy()`
268: @*/
269: PetscErrorCode PetscFESetUp(PetscFE fem)
270: {
271: PetscFunctionBegin;
273: if (fem->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
274: PetscCall(PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0));
275: fem->setupcalled = PETSC_TRUE;
276: PetscTryTypeMethod(fem, setup);
277: PetscCall(PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@
282: PetscFEDestroy - Destroys a `PetscFE` object
284: Collective
286: Input Parameter:
287: . fem - the `PetscFE` object to destroy
289: Level: beginner
291: .seealso: `PetscFE`, `PetscFEView()`
292: @*/
293: PetscErrorCode PetscFEDestroy(PetscFE *fem)
294: {
295: PetscFunctionBegin;
296: if (!*fem) PetscFunctionReturn(PETSC_SUCCESS);
299: if (--((PetscObject)*fem)->refct > 0) {
300: *fem = NULL;
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
303: ((PetscObject)*fem)->refct = 0;
305: if ((*fem)->subspaces) {
306: PetscInt dim;
308: PetscCall(PetscDualSpaceGetDimension((*fem)->dualSpace, &dim));
309: for (PetscInt d = 0; d < dim; ++d) PetscCall(PetscFEDestroy(&(*fem)->subspaces[d]));
310: }
311: PetscCall(PetscFree((*fem)->subspaces));
312: PetscCall(PetscFree((*fem)->invV));
313: PetscCall(PetscTabulationDestroy(&(*fem)->T));
314: PetscCall(PetscTabulationDestroy(&(*fem)->Tf));
315: PetscCall(PetscTabulationDestroy(&(*fem)->Tc));
316: PetscCall(PetscSpaceDestroy(&(*fem)->basisSpace));
317: PetscCall(PetscDualSpaceDestroy(&(*fem)->dualSpace));
318: PetscCall(PetscQuadratureDestroy(&(*fem)->quadrature));
319: PetscCall(PetscQuadratureDestroy(&(*fem)->faceQuadrature));
320: #if PetscDefined(HAVE_LIBCEED)
321: PetscCallCEED(CeedBasisDestroy(&(*fem)->ceedBasis));
322: PetscCallCEED(CeedDestroy(&(*fem)->ceed));
323: #endif
325: PetscTryTypeMethod(*fem, destroy);
326: PetscCall(PetscHeaderDestroy(fem));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*@
331: PetscFECreate - Creates an empty `PetscFE` object. The type can then be set with `PetscFESetType()`.
333: Collective
335: Input Parameter:
336: . comm - The communicator for the `PetscFE` object
338: Output Parameter:
339: . fem - The `PetscFE` object
341: Level: beginner
343: .seealso: `PetscFE`, `PetscFEType`, `PetscFESetType()`, `PetscFECreateDefault()`, `PETSCFEGALERKIN`
344: @*/
345: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
346: {
347: PetscFE f;
349: PetscFunctionBegin;
350: PetscAssertPointer(fem, 2);
351: PetscCall(PetscCitationsRegister(FECitation, &FEcite));
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, pass `NULL` if not needed
484: . numBlocks - The number of blocks in a batch, pass `NULL` if not needed
485: . batchSize - The number of elements in a batch, pass `NULL` if not needed
486: - numBatches - The number of batches in a chunk, pass `NULL` if not needed
488: Level: intermediate
490: .seealso: `PetscFE`, `PetscFECreate()`, `PetscFESetTileSizes()`
491: @*/
492: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PeOp PetscInt *blockSize, PeOp PetscInt *numBlocks, PeOp PetscInt *batchSize, PeOp 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: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetDualSpace()`
543: @*/
544: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
545: {
546: PetscFunctionBegin;
549: PetscCall(PetscSpaceDestroy(&fem->basisSpace));
550: fem->basisSpace = sp;
551: PetscCall(PetscObjectReference((PetscObject)fem->basisSpace));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@
556: PetscFEGetDualSpace - Returns the `PetscDualSpace` used to define the inner product for a `PetscFE`
558: Not Collective
560: Input Parameter:
561: . fem - The `PetscFE` object
563: Output Parameter:
564: . sp - The `PetscDualSpace` object
566: Level: intermediate
568: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
569: @*/
570: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
571: {
572: PetscFunctionBegin;
574: PetscAssertPointer(sp, 2);
575: *sp = fem->dualSpace;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: /*@
580: PetscFESetDualSpace - Sets the `PetscDualSpace` used to define the inner product
582: Not Collective
584: Input Parameters:
585: + fem - The `PetscFE` object
586: - sp - The `PetscDualSpace` object
588: Level: intermediate
590: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`, `PetscFESetBasisSpace()`
591: @*/
592: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
593: {
594: PetscFunctionBegin;
597: PetscCall(PetscDualSpaceDestroy(&fem->dualSpace));
598: fem->dualSpace = sp;
599: PetscCall(PetscObjectReference((PetscObject)fem->dualSpace));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /*@
604: PetscFEGetQuadrature - Returns the `PetscQuadrature` used to calculate inner products
606: Not Collective
608: Input Parameter:
609: . fem - The `PetscFE` object
611: Output Parameter:
612: . q - The `PetscQuadrature` object
614: Level: intermediate
616: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`
617: @*/
618: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
619: {
620: PetscFunctionBegin;
622: PetscAssertPointer(q, 2);
623: *q = fem->quadrature;
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: /*@
628: PetscFESetQuadrature - Sets the `PetscQuadrature` used to calculate inner products
630: Not Collective
632: Input Parameters:
633: + fem - The `PetscFE` object
634: - q - The `PetscQuadrature` object
636: Level: intermediate
638: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFEGetFaceQuadrature()`
639: @*/
640: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
641: {
642: PetscInt Nc, qNc;
644: PetscFunctionBegin;
646: if (q == fem->quadrature) PetscFunctionReturn(PETSC_SUCCESS);
647: PetscCall(PetscFEGetNumComponents(fem, &Nc));
648: PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
649: 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);
650: PetscCall(PetscTabulationDestroy(&fem->T));
651: PetscCall(PetscTabulationDestroy(&fem->Tc));
652: PetscCall(PetscObjectReference((PetscObject)q));
653: PetscCall(PetscQuadratureDestroy(&fem->quadrature));
654: fem->quadrature = q;
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*@
659: PetscFEGetFaceQuadrature - Returns the `PetscQuadrature` used to calculate inner products on faces
661: Not Collective
663: Input Parameter:
664: . fem - The `PetscFE` object
666: Output Parameter:
667: . q - The `PetscQuadrature` object
669: Level: intermediate
671: Developer Notes:
672: There is a special face quadrature but not edge, likely this API would benefit from a refactorization
674: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
675: @*/
676: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
677: {
678: PetscFunctionBegin;
680: PetscAssertPointer(q, 2);
681: *q = fem->faceQuadrature;
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@
686: PetscFESetFaceQuadrature - Sets the `PetscQuadrature` used to calculate inner products on faces
688: Not Collective
690: Input Parameters:
691: + fem - The `PetscFE` object
692: - q - The `PetscQuadrature` object
694: Level: intermediate
696: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`
697: @*/
698: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
699: {
700: PetscInt Nc, qNc;
702: PetscFunctionBegin;
704: if (q == fem->faceQuadrature) PetscFunctionReturn(PETSC_SUCCESS);
705: PetscCall(PetscFEGetNumComponents(fem, &Nc));
706: PetscCall(PetscQuadratureGetNumComponents(q, &qNc));
707: 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);
708: PetscCall(PetscTabulationDestroy(&fem->Tf));
709: PetscCall(PetscObjectReference((PetscObject)q));
710: PetscCall(PetscQuadratureDestroy(&fem->faceQuadrature));
711: fem->faceQuadrature = q;
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: /*@
716: PetscFECopyQuadrature - Copy both volumetric and surface quadrature to a new `PetscFE`
718: Not Collective
720: Input Parameters:
721: + sfe - The `PetscFE` source for the quadratures
722: - tfe - The `PetscFE` target for the quadratures
724: Level: intermediate
726: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`, `PetscFECreate()`, `PetscFESetQuadrature()`, `PetscFESetFaceQuadrature()`
727: @*/
728: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
729: {
730: PetscQuadrature q;
732: PetscFunctionBegin;
735: PetscCall(PetscFEGetQuadrature(sfe, &q));
736: PetscCall(PetscFESetQuadrature(tfe, q));
737: PetscCall(PetscFEGetFaceQuadrature(sfe, &q));
738: PetscCall(PetscFESetFaceQuadrature(tfe, q));
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: /*@
743: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
745: Not Collective
747: Input Parameter:
748: . fem - The `PetscFE` object
750: Output Parameter:
751: . numDof - Array of length `dim` with the number of dofs in each dimension
753: Level: intermediate
755: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscFECreate()`
756: @*/
757: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt *numDof[])
758: {
759: PetscFunctionBegin;
761: PetscAssertPointer(numDof, 2);
762: PetscCall(PetscDualSpaceGetNumDof(fem->dualSpace, numDof));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: /*@
767: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
769: Not Collective
771: Input Parameters:
772: + fem - The `PetscFE` object
773: - k - The highest derivative we need to tabulate, very often 1
775: Output Parameter:
776: . T - The basis function values and derivatives at quadrature points
778: Level: intermediate
780: Note:
781: .vb
782: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
783: 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
784: T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the Hessian value at point p for basis function i, component c, in directions d and e
785: .ve
787: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
788: @*/
789: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
790: {
791: PetscInt npoints;
792: const PetscReal *points;
794: PetscFunctionBegin;
796: PetscAssertPointer(T, 3);
797: PetscCall(PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL));
798: if (!fem->T) PetscCall(PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T));
799: 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);
800: *T = fem->T;
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*@
805: PetscFEExpandFaceQuadrature - Expand a face quadrature into a cell quadrature by mapping the face
806: quadrature points and weights through each face of the cell reference geometry.
808: Not Collective
810: Input Parameters:
811: + fe - the `PetscFE` object whose cell geometry defines the faces
812: - fq - the face quadrature to expand
814: Output Parameter:
815: . efq - the expanded quadrature covering all faces of the cell
817: Level: developer
819: .seealso: `PetscFE`, `PetscQuadrature`, `PetscFECreateFaceQuadrature()`, `PetscFEGetQuadrature()`
820: @*/
821: PetscErrorCode PetscFEExpandFaceQuadrature(PetscFE fe, PetscQuadrature fq, PetscQuadrature *efq)
822: {
823: DM dm;
824: PetscDualSpace sp;
825: const PetscInt *faces;
826: const PetscReal *points, *weights;
827: DMPolytopeType ct;
828: PetscReal *facePoints, *faceWeights;
829: PetscInt dim, cStart, Nf, Nc, Np, order;
831: PetscFunctionBegin;
832: PetscCall(PetscFEGetDualSpace(fe, &sp));
833: PetscCall(PetscDualSpaceGetDM(sp, &dm));
834: PetscCall(DMGetDimension(dm, &dim));
835: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
836: PetscCall(DMPlexGetConeSize(dm, cStart, &Nf));
837: PetscCall(DMPlexGetCone(dm, cStart, &faces));
838: PetscCall(PetscQuadratureGetData(fq, NULL, &Nc, &Np, &points, &weights));
839: PetscCall(PetscMalloc1(Nf * Np * dim, &facePoints));
840: PetscCall(PetscMalloc1(Nf * Np * Nc, &faceWeights));
841: for (PetscInt f = 0; f < Nf; ++f) {
842: const PetscReal xi0[3] = {-1., -1., -1.};
843: PetscReal v0[3], J[9], detJ;
845: PetscCall(DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ));
846: for (PetscInt q = 0; q < Np; ++q) {
847: CoordinatesRefToReal(dim, dim - 1, xi0, v0, J, &points[q * (dim - 1)], &facePoints[(f * Np + q) * dim]);
848: for (PetscInt c = 0; c < Nc; ++c) faceWeights[(f * Np + q) * Nc + c] = weights[q * Nc + c];
849: }
850: }
851: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)fq), efq));
852: PetscCall(PetscQuadratureGetCellType(fq, &ct));
853: PetscCall(PetscQuadratureSetCellType(*efq, ct));
854: PetscCall(PetscQuadratureGetOrder(fq, &order));
855: PetscCall(PetscQuadratureSetOrder(*efq, order));
856: PetscCall(PetscQuadratureSetData(*efq, dim, Nc, Nf * Np, facePoints, faceWeights));
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
863: Not Collective
865: Input Parameters:
866: + fem - The `PetscFE` object
867: - k - The highest derivative we need to tabulate, very often 1
869: Output Parameter:
870: . Tf - The basis function values and derivatives at face quadrature points
872: Level: intermediate
874: Note:
875: .vb
876: 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
877: 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
878: T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the Hessian value at point f,q for basis function i, component c, in directions d and e
879: .ve
881: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
882: @*/
883: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
884: {
885: PetscFunctionBegin;
887: PetscAssertPointer(Tf, 3);
888: if (!fem->Tf) {
889: PetscQuadrature fq;
891: PetscCall(PetscFEGetFaceQuadrature(fem, &fq));
892: if (fq) {
893: PetscQuadrature efq;
894: const PetscReal *facePoints;
895: PetscInt Np, eNp;
897: PetscCall(PetscFEExpandFaceQuadrature(fem, fq, &efq));
898: PetscCall(PetscQuadratureGetData(fq, NULL, NULL, &Np, NULL, NULL));
899: PetscCall(PetscQuadratureGetData(efq, NULL, NULL, &eNp, &facePoints, NULL));
900: if (PetscDefined(USE_DEBUG)) {
901: PetscDualSpace sp;
902: DM dm;
903: PetscInt cStart, Nf;
905: PetscCall(PetscFEGetDualSpace(fem, &sp));
906: PetscCall(PetscDualSpaceGetDM(sp, &dm));
907: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
908: PetscCall(DMPlexGetConeSize(dm, cStart, &Nf));
909: PetscCheck(Nf == eNp / Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of faces %" PetscInt_FMT " != %" PetscInt_FMT " number of quadrature replicas", Nf, eNp / Np);
910: }
911: PetscCall(PetscFECreateTabulation(fem, eNp / Np, Np, facePoints, k, &fem->Tf));
912: PetscCall(PetscQuadratureDestroy(&efq));
913: }
914: }
915: 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);
916: *Tf = fem->Tf;
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@
921: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
923: Not Collective
925: Input Parameter:
926: . fem - The `PetscFE` object
928: Output Parameter:
929: . Tc - The basis function values at face centroid points
931: Level: intermediate
933: Note:
934: .vb
935: T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
936: .ve
938: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscTabulation`, `PetscFEGetFaceTabulation()`, `PetscFEGetCellTabulation()`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`
939: @*/
940: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
941: {
942: PetscFunctionBegin;
944: PetscAssertPointer(Tc, 2);
945: if (!fem->Tc) {
946: PetscDualSpace sp;
947: DM dm;
948: const PetscInt *cone;
949: PetscReal *centroids;
950: PetscInt dim, numFaces, f;
952: PetscCall(PetscFEGetDualSpace(fem, &sp));
953: PetscCall(PetscDualSpaceGetDM(sp, &dm));
954: PetscCall(DMGetDimension(dm, &dim));
955: PetscCall(DMPlexGetConeSize(dm, 0, &numFaces));
956: PetscCall(DMPlexGetCone(dm, 0, &cone));
957: PetscCall(PetscMalloc1(numFaces * dim, ¢roids));
958: for (f = 0; f < numFaces; ++f) PetscCall(DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f * dim], NULL));
959: PetscCall(PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc));
960: PetscCall(PetscFree(centroids));
961: }
962: *Tc = fem->Tc;
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: /*@
967: PetscFECreateTabulation - Creates a `PetscTabulation` object to hold the basis functions, and perhaps derivatives, at the points provided.
969: Not Collective
971: Input Parameters:
972: + fem - The `PetscFE` object
973: . nrepl - The number of replicas
974: . npoints - The number of tabulation points in a replica
975: . points - The tabulation point coordinates
976: - K - The number of derivatives calculated
978: Output Parameter:
979: . T - The `PetscTabulation` to hold the basis function values and derivatives at tabulation points
981: Level: intermediate
983: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`, `PetscFEComputeTabulation()`
984: @*/
985: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
986: {
987: DM dm;
988: PetscDualSpace Q;
989: PetscInt Nb; /* Dimension of FE space P */
990: PetscInt Nc; /* Field components */
991: PetscInt cdim; /* Reference coordinate dimension */
993: PetscFunctionBegin;
994: if (!npoints || !fem->dualSpace || K < 0) {
995: *T = NULL;
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: PetscAssertPointer(points, 4);
1000: PetscAssertPointer(T, 6);
1001: PetscCall(PetscFEGetDualSpace(fem, &Q));
1002: PetscCall(PetscDualSpaceGetDM(Q, &dm));
1003: PetscCall(DMGetDimension(dm, &cdim));
1004: PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
1005: PetscCall(PetscFEGetNumComponents(fem, &Nc));
1006: {
1007: PetscSpace sp;
1008: PetscInt Nv;
1010: PetscCall(PetscFEGetBasisSpace(fem, &sp));
1011: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
1012: PetscCheck(cdim == Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Dual space mesh dim %" PetscInt_FMT " != %" PetscInt_FMT " number of space variables", cdim, Nv);
1013: }
1014: PetscCall(PetscMalloc1(1, T));
1015: (*T)->K = !cdim ? 0 : K;
1016: (*T)->Nr = nrepl;
1017: (*T)->Np = npoints;
1018: (*T)->Nb = Nb;
1019: (*T)->Nc = Nc;
1020: (*T)->cdim = cdim;
1021: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1022: for (PetscInt k = 0; k <= (*T)->K; ++k) PetscCall(PetscCalloc1(nrepl * npoints * Nb * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1023: PetscUseTypeMethod(fem, computetabulation, nrepl * npoints, points, K, *T);
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: /*@
1028: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1030: Not Collective
1032: Input Parameters:
1033: + fem - The `PetscFE` object
1034: . npoints - The number of tabulation points
1035: . points - The tabulation point coordinates
1036: . K - The number of derivatives calculated
1037: - T - An existing tabulation object with enough allocated space, created with `PetscFECreateTabulation()`
1039: Output Parameter:
1040: . T - The basis function values and derivatives at tabulation points
1042: Level: intermediate
1044: Note:
1045: .vb
1046: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1047: 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
1048: T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the Hessian value at point p for basis function i, component c, in directions d and e
1049: .ve
1051: .seealso: `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscTabulationDestroy()`, `PetscFECreateTabulation()`
1052: @*/
1053: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1054: {
1055: PetscFunctionBeginHot;
1056: if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(PETSC_SUCCESS);
1058: PetscAssertPointer(points, 3);
1059: PetscAssertPointer(T, 5);
1060: if (PetscDefined(USE_DEBUG)) {
1061: DM dm;
1062: PetscDualSpace Q;
1063: PetscInt Nb; /* Dimension of FE space P */
1064: PetscInt Nc; /* Field components */
1065: PetscInt cdim; /* Reference coordinate dimension */
1067: PetscCall(PetscFEGetDualSpace(fem, &Q));
1068: PetscCall(PetscDualSpaceGetDM(Q, &dm));
1069: PetscCall(DMGetDimension(dm, &cdim));
1070: PetscCall(PetscDualSpaceGetDimension(Q, &Nb));
1071: PetscCall(PetscFEGetNumComponents(fem, &Nc));
1072: 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);
1073: PetscCheck(T->Nb == Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %" PetscInt_FMT " must match requested Nb %" PetscInt_FMT, T->Nb, Nb);
1074: PetscCheck(T->Nc == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %" PetscInt_FMT " must match requested Nc %" PetscInt_FMT, T->Nc, Nc);
1075: PetscCheck(T->cdim == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %" PetscInt_FMT " must match requested cdim %" PetscInt_FMT, T->cdim, cdim);
1076: }
1077: T->Nr = 1;
1078: T->Np = npoints;
1079: PetscUseTypeMethod(fem, computetabulation, npoints, points, K, T);
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: PetscTabulationDestroy - Frees memory from the associated tabulation.
1086: Not Collective
1088: Input Parameter:
1089: . T - The tabulation
1091: Level: intermediate
1093: .seealso: `PetscTabulation`, `PetscFECreateTabulation()`, `PetscFEGetCellTabulation()`
1094: @*/
1095: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1096: {
1097: PetscFunctionBegin;
1098: PetscAssertPointer(T, 1);
1099: if (!T || !*T) PetscFunctionReturn(PETSC_SUCCESS);
1100: for (PetscInt k = 0; k <= (*T)->K; ++k) PetscCall(PetscFree((*T)->T[k]));
1101: PetscCall(PetscFree((*T)->T));
1102: PetscCall(PetscFree(*T));
1103: *T = NULL;
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: static PetscErrorCode PetscFECreatePointTraceDefault_Internal(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1108: {
1109: PetscSpace bsp, bsubsp;
1110: PetscDualSpace dsp, dsubsp;
1111: PetscInt dim, depth, numComp, i, j, coneSize, order;
1112: DM dm;
1113: DMLabel label;
1114: PetscReal *xi, *v, *J, detJ;
1115: const char *name;
1116: PetscQuadrature origin, fullQuad, subQuad;
1118: PetscFunctionBegin;
1119: PetscCall(PetscFEGetBasisSpace(fe, &bsp));
1120: PetscCall(PetscFEGetDualSpace(fe, &dsp));
1121: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1122: PetscCall(DMGetDimension(dm, &dim));
1123: PetscCall(DMPlexGetDepthLabel(dm, &label));
1124: PetscCall(DMLabelGetValue(label, refPoint, &depth));
1125: PetscCall(PetscCalloc1(depth, &xi));
1126: PetscCall(PetscMalloc1(dim, &v));
1127: PetscCall(PetscMalloc1(dim * dim, &J));
1128: for (i = 0; i < depth; i++) xi[i] = 0.;
1129: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &origin));
1130: PetscCall(PetscQuadratureSetData(origin, depth, 0, 1, xi, NULL));
1131: PetscCall(DMPlexComputeCellGeometryFEM(dm, refPoint, origin, v, J, NULL, &detJ));
1132: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1133: for (i = 1; i < dim; i++) {
1134: for (j = 0; j < depth; j++) J[i * depth + j] = J[i * dim + j];
1135: }
1136: PetscCall(PetscQuadratureDestroy(&origin));
1137: PetscCall(PetscDualSpaceGetPointSubspace(dsp, refPoint, &dsubsp));
1138: PetscCall(PetscSpaceCreateSubspace(bsp, dsubsp, v, J, NULL, NULL, PETSC_OWN_POINTER, &bsubsp));
1139: PetscCall(PetscSpaceSetUp(bsubsp));
1140: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), trFE));
1141: PetscCall(PetscFESetType(*trFE, PETSCFEBASIC));
1142: PetscCall(PetscFEGetNumComponents(fe, &numComp));
1143: PetscCall(PetscFESetNumComponents(*trFE, numComp));
1144: PetscCall(PetscFESetBasisSpace(*trFE, bsubsp));
1145: PetscCall(PetscFESetDualSpace(*trFE, dsubsp));
1146: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1147: if (name) PetscCall(PetscFESetName(*trFE, name));
1148: PetscCall(PetscFEGetQuadrature(fe, &fullQuad));
1149: PetscCall(PetscQuadratureGetOrder(fullQuad, &order));
1150: PetscCall(DMPlexGetConeSize(dm, refPoint, &coneSize));
1151: if (coneSize == 2 * depth) PetscCall(PetscDTGaussTensorQuadrature(depth, 1, (order + 2) / 2, -1., 1., &subQuad));
1152: else PetscCall(PetscDTSimplexQuadrature(depth, order, PETSCDTSIMPLEXQUAD_DEFAULT, &subQuad));
1153: PetscCall(PetscFESetQuadrature(*trFE, subQuad));
1154: PetscCall(PetscFESetUp(*trFE));
1155: PetscCall(PetscQuadratureDestroy(&subQuad));
1156: PetscCall(PetscSpaceDestroy(&bsubsp));
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1161: {
1162: PetscFunctionBegin;
1164: PetscAssertPointer(trFE, 3);
1165: if (fe->ops->createpointtrace) PetscUseTypeMethod(fe, createpointtrace, refPoint, trFE);
1166: else PetscCall(PetscFECreatePointTraceDefault_Internal(fe, refPoint, trFE));
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: /*@
1171: PetscFECreateHeightTrace - Create the trace `PetscFE` for the first mesh point of the given height stratum.
1173: Not Collective
1175: Input Parameters:
1176: + fe - the `PetscFE` object
1177: - height - the height of the stratum whose first point is used to construct the trace element
1179: Output Parameter:
1180: . trFE - the trace `PetscFE`, or `NULL` if the requested height stratum is empty
1182: Level: developer
1184: .seealso: `PetscFE`, `PetscFECreatePointTrace()`, `PetscFEGetHeightSubspace()`, `DMPlexGetHeightStratum()`
1185: @*/
1186: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1187: {
1188: PetscInt hStart, hEnd;
1189: PetscDualSpace dsp;
1190: DM dm;
1192: PetscFunctionBegin;
1194: PetscAssertPointer(trFE, 3);
1195: *trFE = NULL;
1196: PetscCall(PetscFEGetDualSpace(fe, &dsp));
1197: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
1198: PetscCall(DMPlexGetHeightStratum(dm, height, &hStart, &hEnd));
1199: if (hEnd <= hStart) PetscFunctionReturn(PETSC_SUCCESS);
1200: PetscCall(PetscFECreatePointTrace(fe, hStart, trFE));
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1207: Not Collective
1209: Input Parameter:
1210: . fem - The `PetscFE`
1212: Output Parameter:
1213: . dim - The dimension
1215: Level: intermediate
1217: .seealso: `PetscFE`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
1218: @*/
1219: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1220: {
1221: PetscFunctionBegin;
1223: PetscAssertPointer(dim, 2);
1224: PetscTryTypeMethod(fem, getdimension, dim);
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: /*@
1229: PetscFEPushforward - Map the reference element function to real space
1231: Input Parameters:
1232: + fe - The `PetscFE`
1233: . fegeom - The cell geometry
1234: . Nv - The number of function values
1235: - vals - The function values
1237: Output Parameter:
1238: . vals - The transformed function values
1240: Level: advanced
1242: Notes:
1243: This just forwards the call onto `PetscDualSpacePushforward()`.
1245: It only handles transformations when the embedding dimension of the geometry in `fegeom` is the same as the reference dimension.
1247: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscDualSpacePushforward()`
1248: @*/
1249: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1250: {
1251: PetscFunctionBeginHot;
1252: PetscCall(PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1259: Input Parameters:
1260: + fe - The `PetscFE`
1261: . fegeom - The cell geometry
1262: . Nv - The number of function gradient values
1263: - vals - The function gradient values
1265: Output Parameter:
1266: . vals - The transformed function gradient values
1268: Level: advanced
1270: Notes:
1271: This just forwards the call onto `PetscDualSpacePushforwardGradient()`.
1273: It only handles transformations when the embedding dimension of the geometry in `fegeom` is the same as the reference dimension.
1275: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardGradient()`, `PetscDualSpacePushforward()`
1276: @*/
1277: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1278: {
1279: PetscFunctionBeginHot;
1280: PetscCall(PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1284: /*@
1285: PetscFEPushforwardHessian - Map the reference element function Hessian to real space
1287: Input Parameters:
1288: + fe - The `PetscFE`
1289: . fegeom - The cell geometry
1290: . Nv - The number of function Hessian values
1291: - vals - The function Hessian values
1293: Output Parameter:
1294: . vals - The transformed function Hessian values
1296: Level: advanced
1298: Notes:
1299: This just forwards the call onto `PetscDualSpacePushforwardHessian()`.
1301: It only handles transformations when the embedding dimension of the geometry in `fegeom` is the same as the reference dimension.
1303: Developer Note:
1304: It is unclear why all these one line convenience routines are desirable
1306: .seealso: `PetscFE`, `PetscFEGeom`, `PetscDualSpace`, `PetscFEPushforward()`, `PetscDualSpacePushforwardHessian()`, `PetscDualSpacePushforward()`
1307: @*/
1308: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1309: {
1310: PetscFunctionBeginHot;
1311: PetscCall(PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals));
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: /*
1316: Purpose: Compute element vector for chunk of elements
1318: Input:
1319: Sizes:
1320: Ne: number of elements
1321: Nf: number of fields
1322: PetscFE
1323: dim: spatial dimension
1324: Nb: number of basis functions
1325: Nc: number of field components
1326: PetscQuadrature
1327: Nq: number of quadrature points
1329: Geometry:
1330: PetscFEGeom[Ne] possibly *Nq
1331: PetscReal v0s[dim]
1332: PetscReal n[dim]
1333: PetscReal jacobians[dim*dim]
1334: PetscReal jacobianInverses[dim*dim]
1335: PetscReal jacobianDeterminants
1336: FEM:
1337: PetscFE
1338: PetscQuadrature
1339: PetscReal quadPoints[Nq*dim]
1340: PetscReal quadWeights[Nq]
1341: PetscReal basis[Nq*Nb*Nc]
1342: PetscReal basisDer[Nq*Nb*Nc*dim]
1343: PetscScalar coefficients[Ne*Nb*Nc]
1344: PetscScalar elemVec[Ne*Nb*Nc]
1346: Problem:
1347: PetscInt f: the active field
1348: f0, f1
1350: Work Space:
1351: PetscFE
1352: PetscScalar f0[Nq*dim];
1353: PetscScalar f1[Nq*dim*dim];
1354: PetscScalar u[Nc];
1355: PetscScalar gradU[Nc*dim];
1356: PetscReal x[dim];
1357: PetscScalar realSpaceDer[dim];
1359: Purpose: Compute element vector for N_cb batches of elements
1361: Input:
1362: Sizes:
1363: N_cb: Number of serial cell batches
1365: Geometry:
1366: PetscReal v0s[Ne*dim]
1367: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1368: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1369: PetscReal jacobianDeterminants[Ne] possibly *Nq
1370: FEM:
1371: static PetscReal quadPoints[Nq*dim]
1372: static PetscReal quadWeights[Nq]
1373: static PetscReal basis[Nq*Nb*Nc]
1374: static PetscReal basisDer[Nq*Nb*Nc*dim]
1375: PetscScalar coefficients[Ne*Nb*Nc]
1376: PetscScalar elemVec[Ne*Nb*Nc]
1378: ex62.c:
1379: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1380: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1381: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1382: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1384: ex52.c:
1385: 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)
1386: 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)
1388: ex52_integrateElement.cu
1389: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1391: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1392: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1393: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1395: ex52_integrateElementOpenCL.c:
1396: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1397: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1398: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1400: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1401: */
1403: /*@
1404: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1406: Not Collective
1408: Input Parameters:
1409: + prob - The `PetscDS` specifying the discretizations and continuum functions
1410: . field - The field being integrated
1411: . Ne - The number of elements in the chunk
1412: . cgeom - The cell geometry for each cell in the chunk
1413: . coefficients - The array of FEM basis coefficients for the elements
1414: . probAux - The `PetscDS` specifying the auxiliary discretizations
1415: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1417: Output Parameter:
1418: . integral - the integral for this field
1420: Level: intermediate
1422: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrateBd()`
1423: @*/
1424: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1425: {
1426: PetscFE fe;
1428: PetscFunctionBegin;
1430: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1431: if (fe->ops->integrate) PetscCall((*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral));
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1435: /*@
1436: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1438: Not Collective
1440: Input Parameters:
1441: + prob - The `PetscDS` specifying the discretizations and continuum functions
1442: . field - The field being integrated
1443: . obj_func - The function to be integrated
1444: . Ne - The number of elements in the chunk
1445: . geom - The face geometry for each face in the chunk
1446: . coefficients - The array of FEM basis coefficients for the elements
1447: . probAux - The `PetscDS` specifying the auxiliary discretizations
1448: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1450: Output Parameter:
1451: . integral - the integral for this field
1453: Level: intermediate
1455: .seealso: `PetscFE`, `PetscDS`, `PetscFEIntegrateResidual()`, `PetscFEIntegrate()`
1456: @*/
1457: 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[])
1458: {
1459: PetscFE fe;
1461: PetscFunctionBegin;
1463: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
1464: if (fe->ops->integratebd) PetscCall((*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral));
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: /*@
1469: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1471: Not Collective
1473: Input Parameters:
1474: + ds - The `PetscDS` specifying the discretizations and continuum functions
1475: . key - The (label+value, field) being integrated
1476: . Ne - The number of elements in the chunk
1477: . cgeom - The cell geometry for each cell in the chunk
1478: . coefficients - The array of FEM basis coefficients for the elements
1479: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1480: . probAux - The `PetscDS` specifying the auxiliary discretizations
1481: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1482: - t - The time
1484: Output Parameter:
1485: . elemVec - the element residual vectors from each element
1487: Level: intermediate
1489: Note:
1490: .vb
1491: Loop over batch of elements (e):
1492: Loop over quadrature points (q):
1493: Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1494: Call f_0 and f_1
1495: Loop over element vector entries (f,fc --> i):
1496: elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1497: .ve
1499: .seealso: `PetscFEIntegrateBdResidual()`
1500: @*/
1501: 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[])
1502: {
1503: PetscFE fe;
1505: PetscFunctionBeginHot;
1507: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1508: if (fe->ops->integrateresidual) PetscCall((*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /*@
1513: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1515: Not Collective
1517: Input Parameters:
1518: + ds - The `PetscDS` specifying the discretizations and continuum functions
1519: . wf - The PetscWeakForm object holding the pointwise functions
1520: . key - The (label+value, field) being integrated
1521: . Ne - The number of elements in the chunk
1522: . fgeom - The face geometry for each cell in the chunk
1523: . coefficients - The array of FEM basis coefficients for the elements
1524: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1525: . probAux - The `PetscDS` specifying the auxiliary discretizations
1526: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1527: - t - The time
1529: Output Parameter:
1530: . elemVec - the element residual vectors from each element
1532: Level: intermediate
1534: .seealso: `PetscFEIntegrateResidual()`
1535: @*/
1536: 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[])
1537: {
1538: PetscFE fe;
1540: PetscFunctionBegin;
1542: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1543: if (fe->ops->integratebdresidual) PetscCall((*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1544: PetscFunctionReturn(PETSC_SUCCESS);
1545: }
1547: /*@
1548: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1550: Not Collective
1552: Input Parameters:
1553: + ds - The `PetscDS` specifying the discretizations and continuum functions
1554: . dsIn - The `PetscDS` specifying the discretizations and continuum functions for input
1555: . key - The (label+value, field) being integrated
1556: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1557: . Ne - The number of elements in the chunk
1558: . fgeom - The face geometry for each cell in the chunk
1559: . cgeom - The cell geometry for each neighbor cell in the chunk
1560: . coefficients - The array of FEM basis coefficients for the elements
1561: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1562: . probAux - The `PetscDS` specifying the auxiliary discretizations
1563: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1564: - t - The time
1566: Output Parameter:
1567: . elemVec - the element residual vectors from each element
1569: Level: developer
1571: .seealso: `PetscFEIntegrateResidual()`
1572: @*/
1573: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1574: {
1575: PetscFE fe;
1577: PetscFunctionBegin;
1580: PetscCall(PetscDSGetDiscretization(ds, key.field, (PetscObject *)&fe));
1581: if (fe->ops->integratehybridresidual) PetscCall((*fe->ops->integratehybridresidual)(ds, dsIn, key, s, Ne, fgeom, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec));
1582: PetscFunctionReturn(PETSC_SUCCESS);
1583: }
1585: /*@
1586: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1588: Not Collective
1590: Input Parameters:
1591: + rds - The `PetscDS` specifying the row discretizations and continuum functions
1592: . cds - The `PetscDS` specifying the column discretizations
1593: . jtype - The type of matrix pointwise functions that should be used
1594: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1595: . Ne - The number of elements in the chunk
1596: . cgeom - The cell geometry for each cell in the chunk
1597: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1598: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1599: . dsAux - The `PetscDS` specifying the auxiliary discretizations
1600: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1601: . t - The time
1602: - u_tshift - A multiplier for the $dF/du_t$ term (as opposed to the $dF/du$ term)
1604: Output Parameter:
1605: . elemMat - the element matrices for the Jacobian from each element
1607: Level: intermediate
1609: Note:
1610: .vb
1611: Loop over batch of elements (e):
1612: Loop over element matrix entries (f,fc,g,gc --> i,j):
1613: Loop over quadrature points (q):
1614: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1615: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1616: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1617: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1618: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1619: .ve
1621: .seealso: `PetscFEIntegrateResidual()`
1622: @*/
1623: PetscErrorCode PetscFEIntegrateJacobian(PetscDS rds, PetscDS cds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1624: {
1625: PetscFE fe;
1626: PetscInt Nf;
1628: PetscFunctionBegin;
1631: PetscCall(PetscDSGetNumFields(rds, &Nf));
1632: PetscCall(PetscDSGetDiscretization(rds, key.field / Nf, (PetscObject *)&fe));
1633: if (fe->ops->integratejacobian) PetscCall((*fe->ops->integratejacobian)(rds, cds, jtype, key, Ne, cgeom, coefficients, coefficients_t, dsAux, coefficientsAux, t, u_tshift, elemMat));
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1637: /*@
1638: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1640: Not Collective
1642: Input Parameters:
1643: + ds - The `PetscDS` specifying the discretizations and continuum functions
1644: . wf - The PetscWeakForm holding the pointwise functions
1645: . jtype - The type of matrix pointwise functions that should be used
1646: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1647: . Ne - The number of elements in the chunk
1648: . fgeom - The face geometry for each cell in the chunk
1649: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1650: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1651: . probAux - The `PetscDS` specifying the auxiliary discretizations
1652: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1653: . t - The time
1654: - u_tshift - A multiplier for the $dF/du_t$ term (as opposed to the $dF/du$ term)
1656: Output Parameter:
1657: . elemMat - the element matrices for the Jacobian from each element
1659: Level: intermediate
1661: Note:
1662: .vb
1663: Loop over batch of elements (e):
1664: Loop over element matrix entries (f,fc,g,gc --> i,j):
1665: Loop over quadrature points (q):
1666: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1667: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1668: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1669: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1670: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1671: .ve
1673: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1674: @*/
1675: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFEJacobianType jtype, 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[])
1676: {
1677: PetscFE fe;
1678: PetscInt Nf;
1680: PetscFunctionBegin;
1682: PetscCall(PetscDSGetNumFields(ds, &Nf));
1683: PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1684: if (fe->ops->integratebdjacobian) PetscCall((*fe->ops->integratebdjacobian)(ds, wf, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@
1689: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1691: Not Collective
1693: Input Parameters:
1694: + ds - The `PetscDS` specifying the discretizations and continuum functions for the output
1695: . dsIn - The `PetscDS` specifying the discretizations and continuum functions for the input
1696: . jtype - The type of matrix pointwise functions that should be used
1697: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1698: . s - The side of the cell being integrated, 0 for negative and 1 for positive
1699: . Ne - The number of elements in the chunk
1700: . fgeom - The face geometry for each cell in the chunk
1701: . cgeom - The cell geometry for each neighbor cell in the chunk
1702: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1703: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1704: . probAux - The `PetscDS` specifying the auxiliary discretizations
1705: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1706: . t - The time
1707: - u_tshift - A multiplier for the $dF/du_t$ term (as opposed to the $dF/du$ term)
1709: Output Parameter:
1710: . elemMat - the element matrices for the Jacobian from each element
1712: Level: developer
1714: Note:
1715: .vb
1716: Loop over batch of elements (e):
1717: Loop over element matrix entries (f,fc,g,gc --> i,j):
1718: Loop over quadrature points (q):
1719: Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1720: elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1721: + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1722: + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1723: + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1724: .ve
1726: .seealso: `PetscFEIntegrateJacobian()`, `PetscFEIntegrateResidual()`
1727: @*/
1728: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscDS dsIn, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1729: {
1730: PetscFE fe;
1731: PetscInt Nf;
1733: PetscFunctionBegin;
1735: PetscCall(PetscDSGetNumFields(ds, &Nf));
1736: PetscCall(PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *)&fe));
1737: if (fe->ops->integratehybridjacobian) PetscCall((*fe->ops->integratehybridjacobian)(ds, dsIn, jtype, key, s, Ne, fgeom, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat));
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: /*@
1742: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1744: Input Parameters:
1745: + fe - The finite element space
1746: - height - The height of the `DMPLEX` point
1748: Output Parameter:
1749: . subfe - The subspace of this `PetscFE` space
1751: Level: advanced
1753: Note:
1754: For example, if we want the subspace of this space for a face, we would choose height = 1.
1756: .seealso: `PetscFECreateDefault()`
1757: @*/
1758: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1759: {
1760: PetscSpace P, subP;
1761: PetscDualSpace Q, subQ;
1762: PetscQuadrature subq;
1763: PetscInt dim, Nc;
1765: PetscFunctionBegin;
1767: PetscAssertPointer(subfe, 3);
1768: if (height == 0) {
1769: *subfe = fe;
1770: PetscFunctionReturn(PETSC_SUCCESS);
1771: }
1772: PetscCall(PetscFEGetBasisSpace(fe, &P));
1773: PetscCall(PetscFEGetDualSpace(fe, &Q));
1774: PetscCall(PetscFEGetNumComponents(fe, &Nc));
1775: PetscCall(PetscFEGetFaceQuadrature(fe, &subq));
1776: PetscCall(PetscDualSpaceGetDimension(Q, &dim));
1777: 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);
1778: if (!fe->subspaces) PetscCall(PetscCalloc1(dim, &fe->subspaces));
1779: if (height <= dim) {
1780: if (!fe->subspaces[height - 1]) {
1781: PetscFE sub = NULL;
1782: const char *name;
1784: PetscCall(PetscSpaceGetHeightSubspace(P, height, &subP));
1785: PetscCall(PetscDualSpaceGetHeightSubspace(Q, height, &subQ));
1786: if (subQ) {
1787: PetscCall(PetscObjectReference((PetscObject)subP));
1788: PetscCall(PetscObjectReference((PetscObject)subQ));
1789: PetscCall(PetscObjectReference((PetscObject)subq));
1790: PetscCall(PetscFECreateFromSpaces(subP, subQ, subq, NULL, &sub));
1791: }
1792: if (sub) {
1793: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
1794: if (name) PetscCall(PetscFESetName(sub, name));
1795: }
1796: fe->subspaces[height - 1] = sub;
1797: }
1798: *subfe = fe->subspaces[height - 1];
1799: } else {
1800: *subfe = NULL;
1801: }
1802: PetscFunctionReturn(PETSC_SUCCESS);
1803: }
1805: /*@
1806: PetscFERefine - Create a "refined" `PetscFE` object that refines the reference cell into
1807: smaller copies.
1809: Collective
1811: Input Parameter:
1812: . fe - The initial `PetscFE`
1814: Output Parameter:
1815: . feRef - The refined `PetscFE`
1817: Level: advanced
1819: Notes:
1820: This is typically used to generate a preconditioner for a higher order method from a lower order method on a
1821: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1822: interpolation between regularly refined meshes.
1824: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1825: @*/
1826: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1827: {
1828: PetscSpace P, Pref;
1829: PetscDualSpace Q, Qref;
1830: DM K, Kref;
1831: PetscQuadrature q, qref;
1832: const PetscReal *v0, *jac;
1833: PetscInt numComp, numSubelements;
1834: PetscInt cStart, cEnd, c;
1835: PetscDualSpace *cellSpaces;
1837: PetscFunctionBegin;
1838: PetscCall(PetscFEGetBasisSpace(fe, &P));
1839: PetscCall(PetscFEGetDualSpace(fe, &Q));
1840: PetscCall(PetscFEGetQuadrature(fe, &q));
1841: PetscCall(PetscDualSpaceGetDM(Q, &K));
1842: /* Create space */
1843: PetscCall(PetscObjectReference((PetscObject)P));
1844: Pref = P;
1845: /* Create dual space */
1846: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1847: PetscCall(PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED));
1848: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fe), &Kref));
1849: PetscCall(DMGetCoordinatesLocalSetUp(Kref));
1850: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1851: PetscCall(DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd));
1852: PetscCall(PetscMalloc1(cEnd - cStart, &cellSpaces));
1853: /* TODO: fix for non-uniform refinement */
1854: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1855: PetscCall(PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces));
1856: PetscCall(PetscFree(cellSpaces));
1857: PetscCall(DMDestroy(&Kref));
1858: PetscCall(PetscDualSpaceSetUp(Qref));
1859: /* Create element */
1860: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)fe), feRef));
1861: PetscCall(PetscFESetType(*feRef, PETSCFECOMPOSITE));
1862: PetscCall(PetscFESetBasisSpace(*feRef, Pref));
1863: PetscCall(PetscFESetDualSpace(*feRef, Qref));
1864: PetscCall(PetscFEGetNumComponents(fe, &numComp));
1865: PetscCall(PetscFESetNumComponents(*feRef, numComp));
1866: PetscCall(PetscFESetUp(*feRef));
1867: PetscCall(PetscSpaceDestroy(&Pref));
1868: PetscCall(PetscDualSpaceDestroy(&Qref));
1869: /* Create quadrature */
1870: PetscCall(PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL));
1871: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1872: PetscCall(PetscFESetQuadrature(*feRef, qref));
1873: PetscCall(PetscQuadratureDestroy(&qref));
1874: PetscFunctionReturn(PETSC_SUCCESS);
1875: }
1877: static PetscErrorCode PetscFESetDefaultName_Private(PetscFE fe)
1878: {
1879: PetscSpace P;
1880: PetscDualSpace Q;
1881: DM K;
1882: DMPolytopeType ct;
1883: PetscInt degree;
1884: char name[64];
1886: PetscFunctionBegin;
1887: PetscCall(PetscFEGetBasisSpace(fe, &P));
1888: PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
1889: PetscCall(PetscFEGetDualSpace(fe, &Q));
1890: PetscCall(PetscDualSpaceGetDM(Q, &K));
1891: PetscCall(DMPlexGetCellType(K, 0, &ct));
1892: switch (ct) {
1893: case DM_POLYTOPE_SEGMENT:
1894: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1895: case DM_POLYTOPE_QUADRILATERAL:
1896: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1897: case DM_POLYTOPE_HEXAHEDRON:
1898: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1899: PetscCall(PetscSNPrintf(name, sizeof(name), "Q%" PetscInt_FMT, degree));
1900: break;
1901: case DM_POLYTOPE_TRIANGLE:
1902: case DM_POLYTOPE_TETRAHEDRON:
1903: PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT, degree));
1904: break;
1905: case DM_POLYTOPE_TRI_PRISM:
1906: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1907: PetscCall(PetscSNPrintf(name, sizeof(name), "P%" PetscInt_FMT "xQ%" PetscInt_FMT, degree, degree));
1908: break;
1909: default:
1910: PetscCall(PetscSNPrintf(name, sizeof(name), "FE"));
1911: }
1912: PetscCall(PetscFESetName(fe, name));
1913: PetscFunctionReturn(PETSC_SUCCESS);
1914: }
1916: /*@
1917: PetscFECreateFromSpaces - Create a `PetscFE` from the basis and dual spaces
1919: Collective
1921: Input Parameters:
1922: + P - The basis space
1923: . Q - The dual space
1924: . q - The cell quadrature
1925: - fq - The face quadrature
1927: Output Parameter:
1928: . fem - The `PetscFE` object
1930: Level: beginner
1932: Note:
1933: 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,
1934: the caller must use `PetscObjectReference()` before this call.
1936: .seealso: `PetscFE`, `PetscSpace`, `PetscDualSpace`, `PetscQuadrature`,
1937: `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
1938: @*/
1939: PetscErrorCode PetscFECreateFromSpaces(PetscSpace P, PetscDualSpace Q, PetscQuadrature q, PetscQuadrature fq, PetscFE *fem)
1940: {
1941: PetscInt Nc;
1942: PetscInt p_Ns = -1, p_Nc = -1, q_Ns = -1, q_Nc = -1;
1943: PetscBool p_is_uniform_sum = PETSC_FALSE, p_interleave_basis = PETSC_FALSE, p_interleave_components = PETSC_FALSE;
1944: PetscBool q_is_uniform_sum = PETSC_FALSE, q_interleave_basis = PETSC_FALSE, q_interleave_components = PETSC_FALSE;
1945: const char *prefix;
1947: PetscFunctionBegin;
1948: PetscCall(PetscObjectTypeCompare((PetscObject)P, PETSCSPACESUM, &p_is_uniform_sum));
1949: if (p_is_uniform_sum) {
1950: PetscSpace subsp_0 = NULL;
1951: PetscCall(PetscSpaceSumGetNumSubspaces(P, &p_Ns));
1952: PetscCall(PetscSpaceGetNumComponents(P, &p_Nc));
1953: PetscCall(PetscSpaceSumGetConcatenate(P, &p_is_uniform_sum));
1954: PetscCall(PetscSpaceSumGetInterleave(P, &p_interleave_basis, &p_interleave_components));
1955: for (PetscInt s = 0; s < p_Ns; s++) {
1956: PetscSpace subsp;
1958: PetscCall(PetscSpaceSumGetSubspace(P, s, &subsp));
1959: if (!s) {
1960: subsp_0 = subsp;
1961: } else if (subsp != subsp_0) {
1962: p_is_uniform_sum = PETSC_FALSE;
1963: }
1964: }
1965: }
1966: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &q_is_uniform_sum));
1967: if (q_is_uniform_sum) {
1968: PetscDualSpace subsp_0 = NULL;
1969: PetscCall(PetscDualSpaceSumGetNumSubspaces(Q, &q_Ns));
1970: PetscCall(PetscDualSpaceGetNumComponents(Q, &q_Nc));
1971: PetscCall(PetscDualSpaceSumGetConcatenate(Q, &q_is_uniform_sum));
1972: PetscCall(PetscDualSpaceSumGetInterleave(Q, &q_interleave_basis, &q_interleave_components));
1973: for (PetscInt s = 0; s < q_Ns; s++) {
1974: PetscDualSpace subsp;
1976: PetscCall(PetscDualSpaceSumGetSubspace(Q, s, &subsp));
1977: if (!s) {
1978: subsp_0 = subsp;
1979: } else if (subsp != subsp_0) {
1980: q_is_uniform_sum = PETSC_FALSE;
1981: }
1982: }
1983: }
1984: 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)) {
1985: PetscSpace scalar_space;
1986: PetscDualSpace scalar_dspace;
1987: PetscFE scalar_fe;
1989: PetscCall(PetscSpaceSumGetSubspace(P, 0, &scalar_space));
1990: PetscCall(PetscDualSpaceSumGetSubspace(Q, 0, &scalar_dspace));
1991: PetscCall(PetscObjectReference((PetscObject)scalar_space));
1992: PetscCall(PetscObjectReference((PetscObject)scalar_dspace));
1993: PetscCall(PetscObjectReference((PetscObject)q));
1994: PetscCall(PetscObjectReference((PetscObject)fq));
1995: PetscCall(PetscFECreateFromSpaces(scalar_space, scalar_dspace, q, fq, &scalar_fe));
1996: PetscCall(PetscFECreateVector(scalar_fe, p_Ns, p_interleave_basis, p_interleave_components, fem));
1997: PetscCall(PetscFEDestroy(&scalar_fe));
1998: } else {
1999: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)P), fem));
2000: PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
2001: }
2002: PetscCall(PetscSpaceGetNumComponents(P, &Nc));
2003: PetscCall(PetscFESetNumComponents(*fem, Nc));
2004: PetscCall(PetscFESetBasisSpace(*fem, P));
2005: PetscCall(PetscFESetDualSpace(*fem, Q));
2006: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)P, &prefix));
2007: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
2008: PetscCall(PetscFESetUp(*fem));
2009: PetscCall(PetscSpaceDestroy(&P));
2010: PetscCall(PetscDualSpaceDestroy(&Q));
2011: PetscCall(PetscFESetQuadrature(*fem, q));
2012: PetscCall(PetscFESetFaceQuadrature(*fem, fq));
2013: PetscCall(PetscQuadratureDestroy(&q));
2014: PetscCall(PetscQuadratureDestroy(&fq));
2015: PetscCall(PetscFESetDefaultName_Private(*fem));
2016: PetscFunctionReturn(PETSC_SUCCESS);
2017: }
2019: static PetscErrorCode PetscFECreate_Internal(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt degree, PetscInt qorder, PetscBool setFromOptions, PetscFE *fem)
2020: {
2021: DM K;
2022: PetscSpace P;
2023: PetscDualSpace Q;
2024: PetscQuadrature q, fq;
2025: PetscBool tensor;
2026: PetscDTSimplexQuadratureType qtype = PETSCDTSIMPLEXQUAD_DEFAULT;
2028: PetscFunctionBegin;
2029: if (prefix) PetscAssertPointer(prefix, 5);
2030: PetscAssertPointer(fem, 9);
2031: switch (ct) {
2032: case DM_POLYTOPE_SEGMENT:
2033: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2034: case DM_POLYTOPE_QUADRILATERAL:
2035: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2036: case DM_POLYTOPE_HEXAHEDRON:
2037: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2038: tensor = PETSC_TRUE;
2039: break;
2040: default:
2041: tensor = PETSC_FALSE;
2042: }
2043: /* Create space */
2044: PetscCall(PetscSpaceCreate(comm, &P));
2045: PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
2046: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
2047: PetscCall(PetscSpacePolynomialSetTensor(P, tensor));
2048: PetscCall(PetscSpaceSetNumComponents(P, Nc));
2049: PetscCall(PetscSpaceSetNumVariables(P, dim));
2050: if (degree >= 0) {
2051: PetscCall(PetscSpaceSetDegree(P, degree, PETSC_DETERMINE));
2052: if (ct == DM_POLYTOPE_TRI_PRISM || ct == DM_POLYTOPE_TRI_PRISM_TENSOR) {
2053: PetscSpace Pend, Pside;
2055: PetscCall(PetscSpaceSetNumComponents(P, 1));
2056: PetscCall(PetscSpaceCreate(comm, &Pend));
2057: PetscCall(PetscSpaceSetType(Pend, PETSCSPACEPOLYNOMIAL));
2058: PetscCall(PetscSpacePolynomialSetTensor(Pend, PETSC_FALSE));
2059: PetscCall(PetscSpaceSetNumComponents(Pend, 1));
2060: PetscCall(PetscSpaceSetNumVariables(Pend, dim - 1));
2061: PetscCall(PetscSpaceSetDegree(Pend, degree, PETSC_DETERMINE));
2062: PetscCall(PetscSpaceCreate(comm, &Pside));
2063: PetscCall(PetscSpaceSetType(Pside, PETSCSPACEPOLYNOMIAL));
2064: PetscCall(PetscSpacePolynomialSetTensor(Pside, PETSC_FALSE));
2065: PetscCall(PetscSpaceSetNumComponents(Pside, 1));
2066: PetscCall(PetscSpaceSetNumVariables(Pside, 1));
2067: PetscCall(PetscSpaceSetDegree(Pside, degree, PETSC_DETERMINE));
2068: PetscCall(PetscSpaceSetType(P, PETSCSPACETENSOR));
2069: PetscCall(PetscSpaceTensorSetNumSubspaces(P, 2));
2070: PetscCall(PetscSpaceTensorSetSubspace(P, 0, Pend));
2071: PetscCall(PetscSpaceTensorSetSubspace(P, 1, Pside));
2072: PetscCall(PetscSpaceDestroy(&Pend));
2073: PetscCall(PetscSpaceDestroy(&Pside));
2075: if (Nc > 1) {
2076: PetscSpace scalar_P = P;
2078: PetscCall(PetscSpaceCreate(comm, &P));
2079: PetscCall(PetscSpaceSetNumVariables(P, dim));
2080: PetscCall(PetscSpaceSetNumComponents(P, Nc));
2081: PetscCall(PetscSpaceSetType(P, PETSCSPACESUM));
2082: PetscCall(PetscSpaceSumSetNumSubspaces(P, Nc));
2083: PetscCall(PetscSpaceSumSetConcatenate(P, PETSC_TRUE));
2084: PetscCall(PetscSpaceSumSetInterleave(P, PETSC_TRUE, PETSC_FALSE));
2085: for (PetscInt i = 0; i < Nc; i++) PetscCall(PetscSpaceSumSetSubspace(P, i, scalar_P));
2086: PetscCall(PetscSpaceDestroy(&scalar_P));
2087: }
2088: }
2089: }
2090: if (setFromOptions) PetscCall(PetscSpaceSetFromOptions(P));
2091: PetscCall(PetscSpaceSetUp(P));
2092: PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
2093: PetscCall(PetscSpacePolynomialGetTensor(P, &tensor));
2094: PetscCall(PetscSpaceGetNumComponents(P, &Nc));
2095: /* Create dual space */
2096: PetscCall(PetscDualSpaceCreate(comm, &Q));
2097: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
2098: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
2099: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
2100: PetscCall(PetscDualSpaceSetDM(Q, K));
2101: PetscCall(DMDestroy(&K));
2102: PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
2103: PetscCall(PetscDualSpaceSetOrder(Q, degree));
2104: PetscCall(PetscDualSpaceLagrangeSetTensor(Q, (tensor || (ct == DM_POLYTOPE_TRI_PRISM)) ? PETSC_TRUE : PETSC_FALSE));
2105: if (setFromOptions) PetscCall(PetscDualSpaceSetFromOptions(Q));
2106: PetscCall(PetscDualSpaceSetUp(Q));
2108: qorder = qorder >= 0 ? qorder : degree;
2109: if (setFromOptions) {
2110: PetscObjectOptionsBegin((PetscObject)P);
2111: PetscCall(PetscOptionsBoundedInt("-petscfe_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "PetscFECreateDefault", qorder, &qorder, NULL, 0));
2112: PetscCall(PetscOptionsEnum("-petscfe_default_quadrature_type", "Simplex quadrature type", "PetscDTSimplexQuadratureType", PetscDTSimplexQuadratureTypes, (PetscEnum)qtype, (PetscEnum *)&qtype, NULL));
2113: PetscOptionsEnd();
2114: }
2115: PetscCall(PetscDTCreateQuadratureByCell(ct, qorder, qtype, &q, &fq));
2116: /* Create finite element */
2117: PetscCall(PetscFECreateFromSpaces(P, Q, q, fq, fem));
2118: if (setFromOptions) PetscCall(PetscFESetFromOptions(*fem));
2119: PetscFunctionReturn(PETSC_SUCCESS);
2120: }
2122: /*@
2123: PetscFECreateDefault - Create a `PetscFE` for basic FEM computation
2125: Collective
2127: Input Parameters:
2128: + comm - The MPI comm
2129: . dim - The spatial dimension
2130: . Nc - The number of components
2131: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2132: . prefix - The options prefix, or `NULL`
2133: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2135: Output Parameter:
2136: . fem - The `PetscFE` object
2138: Level: beginner
2140: Notes:
2141: Preferred usage is `PetscFECreateByCell()`
2143: Each subobject is SetFromOption() during creation, so that the object may be customized from the command line, using the prefix specified above.
2144: See the links below for the particular options available.
2146: .seealso: `PetscFE`, `PetscFECreateLagrange()`, `PetscFECreateByCell()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`, `PetscFESetFromOptions()`,
2147: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2148: @*/
2149: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
2150: {
2151: PetscFunctionBegin;
2152: PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2153: PetscFunctionReturn(PETSC_SUCCESS);
2154: }
2156: /*@
2157: PetscFECreateByCell - Create a `PetscFE` for basic FEM computation
2159: Collective
2161: Input Parameters:
2162: + comm - The MPI comm
2163: . dim - The spatial dimension
2164: . Nc - The number of components
2165: . ct - The celltype of the reference cell
2166: . prefix - The options prefix, or `NULL`
2167: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2169: Output Parameter:
2170: . fem - The `PetscFE` object
2172: Level: beginner
2174: Note:
2175: 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.
2177: Developer Notes:
2178: This should be called `PetscFECreateDefaultByCell()` since it is the extension/replacement for `PetscFECreateDefault()`
2180: Since this generalizes/replaces `PetscFECreateDefault()` for different `DMPolytopeType` its name should be `PetscFECreateDefaultByPolytopeType()`
2182: .seealso: `PetscFE`, `PetscFECreateDefault()`, `PetscFECreateLagrange()`, `PetscSpaceSetFromOptions()`, `PetscDualSpaceSetFromOptions()`,
2183: `PetscFESetFromOptions()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`, `DMPolytopeType`
2184: @*/
2185: PetscErrorCode PetscFECreateByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, const char prefix[], PetscInt qorder, PetscFE *fem)
2186: {
2187: PetscFunctionBegin;
2188: PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, prefix, PETSC_DECIDE, qorder, PETSC_TRUE, fem));
2189: PetscFunctionReturn(PETSC_SUCCESS);
2190: }
2192: /*@
2193: PetscFECreateLagrange - Create a `PetscFE` for the basic Lagrange space of degree `k`
2195: Collective
2197: Input Parameters:
2198: + comm - The MPI comm
2199: . dim - The spatial dimension
2200: . Nc - The number of components
2201: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
2202: . k - The degree of the space
2203: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2205: Output Parameter:
2206: . fem - The `PetscFE` object
2208: Level: beginner
2210: Notes:
2211: Preferred usage is `PetscFECreateLagrangeByCell()`
2213: 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`.
2215: .seealso: `PetscFE`, `PetscFECreateLagrangeByCell()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2216: @*/
2217: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
2218: {
2219: PetscFunctionBegin;
2220: PetscCall(PetscFECreate_Internal(comm, dim, Nc, DMPolytopeTypeSimpleShape(dim, isSimplex), NULL, k, qorder, PETSC_FALSE, fem));
2221: PetscFunctionReturn(PETSC_SUCCESS);
2222: }
2224: /*@
2225: PetscFECreateLagrangeByCell - Create a `PetscFE` for the basic Lagrange space of degree `k`
2227: Collective
2229: Input Parameters:
2230: + comm - The MPI comm
2231: . dim - The spatial dimension
2232: . Nc - The number of components
2233: . ct - The celltype of the reference cell
2234: . k - The degree of the space
2235: - qorder - The quadrature order or `PETSC_DETERMINE` to use `PetscSpace` polynomial degree
2237: Output Parameter:
2238: . fem - The `PetscFE` object
2240: Level: beginner
2242: Note:
2243: 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`.
2245: Developer Note:
2246: Since this generalizes/replaces `PetscFECreateLagrange()` for different `DMPolytopeType` its name should be `PetscFECreateLagrangeByPolytopeType()`
2248: .seealso: `PetscFE`, `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`,
2249: `DMPolytopeType`
2250: @*/
2251: PetscErrorCode PetscFECreateLagrangeByCell(MPI_Comm comm, PetscInt dim, PetscInt Nc, DMPolytopeType ct, PetscInt k, PetscInt qorder, PetscFE *fem)
2252: {
2253: PetscFunctionBegin;
2254: PetscCall(PetscFECreate_Internal(comm, dim, Nc, ct, NULL, k, qorder, PETSC_FALSE, fem));
2255: PetscFunctionReturn(PETSC_SUCCESS);
2256: }
2258: /*@
2259: PetscFELimitDegree - Copy a `PetscFE` but limit the degree to be in the given range
2261: Collective
2263: Input Parameters:
2264: + fe - The `PetscFE`
2265: . minDegree - The minimum degree, or `PETSC_DETERMINE` for no limit
2266: - maxDegree - The maximum degree, or `PETSC_DETERMINE` for no limit
2268: Output Parameter:
2269: . newfe - The `PetscFE` object
2271: Level: advanced
2273: Note:
2274: This currently only works for Lagrange elements.
2276: .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2277: @*/
2278: PetscErrorCode PetscFELimitDegree(PetscFE fe, PetscInt minDegree, PetscInt maxDegree, PetscFE *newfe)
2279: {
2280: PetscDualSpace Q;
2281: PetscBool islag, issum;
2282: PetscInt oldk = 0, k;
2284: PetscFunctionBegin;
2285: PetscCall(PetscFEGetDualSpace(fe, &Q));
2286: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACELAGRANGE, &islag));
2287: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &issum));
2288: if (islag) {
2289: PetscCall(PetscDualSpaceGetOrder(Q, &oldk));
2290: } else if (issum) {
2291: PetscDualSpace subQ;
2293: PetscCall(PetscDualSpaceSumGetSubspace(Q, 0, &subQ));
2294: PetscCall(PetscDualSpaceGetOrder(subQ, &oldk));
2295: } else {
2296: PetscCall(PetscObjectReference((PetscObject)fe));
2297: *newfe = fe;
2298: PetscFunctionReturn(PETSC_SUCCESS);
2299: }
2300: k = oldk;
2301: if (minDegree >= 0) k = PetscMax(k, minDegree);
2302: if (maxDegree >= 0) k = PetscMin(k, maxDegree);
2303: if (k != oldk) {
2304: DM K;
2305: PetscSpace P;
2306: PetscQuadrature q;
2307: DMPolytopeType ct;
2308: PetscInt dim, Nc;
2310: PetscCall(PetscFEGetBasisSpace(fe, &P));
2311: PetscCall(PetscSpaceGetNumVariables(P, &dim));
2312: PetscCall(PetscSpaceGetNumComponents(P, &Nc));
2313: PetscCall(PetscDualSpaceGetDM(Q, &K));
2314: PetscCall(DMPlexGetCellType(K, 0, &ct));
2315: PetscCall(PetscFECreateLagrangeByCell(PetscObjectComm((PetscObject)fe), dim, Nc, ct, k, PETSC_DETERMINE, newfe));
2316: PetscCall(PetscFEGetQuadrature(fe, &q));
2317: PetscCall(PetscFESetQuadrature(*newfe, q));
2318: } else {
2319: PetscCall(PetscObjectReference((PetscObject)fe));
2320: *newfe = fe;
2321: }
2322: PetscFunctionReturn(PETSC_SUCCESS);
2323: }
2325: /*@
2326: PetscFECreateBrokenElement - Create a discontinuous version of the input `PetscFE`
2328: Collective
2330: Input Parameters:
2331: . cgfe - The continuous `PetscFE` object
2333: Output Parameter:
2334: . dgfe - The discontinuous `PetscFE` object
2336: Level: advanced
2338: Note:
2339: This only works for Lagrange elements.
2341: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`, `PetscFECreateLagrange()`, `PetscFECreateLagrangeByCell()`, `PetscDualSpaceLagrangeSetContinuity()`
2342: @*/
2343: PetscErrorCode PetscFECreateBrokenElement(PetscFE cgfe, PetscFE *dgfe)
2344: {
2345: PetscSpace P;
2346: PetscDualSpace Q, dgQ;
2347: PetscQuadrature q, fq;
2348: PetscBool is_lagrange, is_sum;
2350: PetscFunctionBegin;
2351: PetscCall(PetscFEGetBasisSpace(cgfe, &P));
2352: PetscCall(PetscObjectReference((PetscObject)P));
2353: PetscCall(PetscFEGetDualSpace(cgfe, &Q));
2354: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACELAGRANGE, &is_lagrange));
2355: PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &is_sum));
2356: PetscCheck(is_lagrange || is_sum, PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only create broken elements of Lagrange elements");
2357: PetscCall(PetscDualSpaceDuplicate(Q, &dgQ));
2358: PetscCall(PetscDualSpaceLagrangeSetContinuity(dgQ, PETSC_FALSE));
2359: PetscCall(PetscDualSpaceSetUp(dgQ));
2360: PetscCall(PetscFEGetQuadrature(cgfe, &q));
2361: PetscCall(PetscObjectReference((PetscObject)q));
2362: PetscCall(PetscFEGetFaceQuadrature(cgfe, &fq));
2363: PetscCall(PetscObjectReference((PetscObject)fq));
2364: PetscCall(PetscFECreateFromSpaces(P, dgQ, q, fq, dgfe));
2365: PetscFunctionReturn(PETSC_SUCCESS);
2366: }
2368: /*@
2369: PetscFESetName - Names the `PetscFE` and its subobjects
2371: Not Collective
2373: Input Parameters:
2374: + fe - The `PetscFE`
2375: - name - The name
2377: Level: intermediate
2379: .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()`
2380: @*/
2381: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2382: {
2383: PetscSpace P;
2384: PetscDualSpace Q;
2386: PetscFunctionBegin;
2387: PetscCall(PetscFEGetBasisSpace(fe, &P));
2388: PetscCall(PetscFEGetDualSpace(fe, &Q));
2389: PetscCall(PetscObjectSetName((PetscObject)fe, name));
2390: PetscCall(PetscObjectSetName((PetscObject)P, name));
2391: PetscCall(PetscObjectSetName((PetscObject)Q, name));
2392: PetscFunctionReturn(PETSC_SUCCESS);
2393: }
2395: 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[])
2396: {
2397: PetscInt dOffset = 0, fOffset = 0, f, g;
2399: for (f = 0; f < Nf; ++f) {
2400: 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);
2401: 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);
2402: PetscFE fe;
2403: const PetscInt k = ds->jetDegree[f];
2404: const PetscInt cdim = T[f]->cdim;
2405: const PetscInt dE = fegeom->dimEmbed;
2406: const PetscInt Nq = T[f]->Np;
2407: const PetscInt Nbf = T[f]->Nb;
2408: const PetscInt Ncf = T[f]->Nc;
2409: const PetscReal *Bq = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf];
2410: const PetscReal *Dq = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim];
2411: const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL;
2412: PetscInt hOffset = 0, b, c, d;
2414: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
2415: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2416: for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2417: for (b = 0; b < Nbf; ++b) {
2418: for (c = 0; c < Ncf; ++c) {
2419: const PetscInt cidx = b * Ncf + c;
2421: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2422: for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b];
2423: }
2424: }
2425: if (k > 1) {
2426: for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * dE;
2427: for (d = 0; d < dE * dE * Ncf; ++d) u_x[hOffset + fOffset * dE * dE + d] = 0.0;
2428: for (b = 0; b < Nbf; ++b) {
2429: for (c = 0; c < Ncf; ++c) {
2430: const PetscInt cidx = b * Ncf + c;
2432: for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * dE * dE + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b];
2433: }
2434: }
2435: PetscCall(PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * dE * dE]));
2436: }
2437: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset]));
2438: PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE]));
2439: if (u_t) {
2440: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2441: for (b = 0; b < Nbf; ++b) {
2442: for (c = 0; c < Ncf; ++c) {
2443: const PetscInt cidx = b * Ncf + c;
2445: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2446: }
2447: }
2448: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2449: }
2450: fOffset += Ncf;
2451: dOffset += Nbf;
2452: }
2453: return PETSC_SUCCESS;
2454: }
2456: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt rc, PetscInt qc, PetscTabulation Tab[], const PetscInt rf[], const PetscInt qf[], PetscTabulation Tabf[], PetscFEGeom *fegeom, PetscFEGeom *fegeomNbr, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2457: {
2458: PetscInt dOffset = 0, fOffset = 0, f;
2460: /* f is the field number in the DS */
2461: for (f = 0; f < Nf; ++f) {
2462: PetscBool isCohesive;
2463: PetscInt Ns;
2465: if (!Tab[f]) continue;
2466: PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
2467: Ns = isCohesive ? 1 : 2;
2468: {
2469: PetscTabulation T = isCohesive ? Tab[f] : Tabf[f];
2470: PetscFE fe = (PetscFE)ds->disc[f];
2471: const PetscInt dEt = T->cdim;
2472: const PetscInt dE = fegeom->dimEmbed;
2473: const PetscInt Nq = T->Np;
2474: const PetscInt Nbf = T->Nb;
2475: const PetscInt Ncf = T->Nc;
2477: for (PetscInt s = 0; s < Ns; ++s) {
2478: const PetscInt r = isCohesive ? rc : rf[s];
2479: const PetscInt q = isCohesive ? qc : qf[s];
2480: const PetscReal *Bq = &T->T[0][(r * Nq + q) * Nbf * Ncf];
2481: const PetscReal *Dq = &T->T[1][(r * Nq + q) * Nbf * Ncf * dEt];
2482: PetscInt b, c, d;
2484: 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);
2485: 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);
2486: for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0;
2487: for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0;
2488: for (b = 0; b < Nbf; ++b) {
2489: for (c = 0; c < Ncf; ++c) {
2490: const PetscInt cidx = b * Ncf + c;
2492: u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b];
2493: for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b];
2494: }
2495: }
2496: PetscCall(PetscFEPushforward(fe, isCohesive ? fegeom : &fegeomNbr[s], 1, &u[fOffset]));
2497: PetscCall(PetscFEPushforwardGradient(fe, isCohesive ? fegeom : &fegeomNbr[s], 1, &u_x[fOffset * dE]));
2498: if (u_t) {
2499: for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0;
2500: for (b = 0; b < Nbf; ++b) {
2501: for (c = 0; c < Ncf; ++c) {
2502: const PetscInt cidx = b * Ncf + c;
2504: u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b];
2505: }
2506: }
2507: PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]));
2508: }
2509: fOffset += Ncf;
2510: dOffset += Nbf;
2511: }
2512: }
2513: }
2514: return PETSC_SUCCESS;
2515: }
2517: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2518: {
2519: PetscFE fe;
2520: PetscTabulation Tc;
2521: PetscInt b, c;
2523: if (!prob) return PETSC_SUCCESS;
2524: PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe));
2525: PetscCall(PetscFEGetFaceCentroidTabulation(fe, &Tc));
2526: {
2527: const PetscReal *faceBasis = Tc->T[0];
2528: const PetscInt Nb = Tc->Nb;
2529: const PetscInt Nc = Tc->Nc;
2531: for (c = 0; c < Nc; ++c) u[c] = 0.0;
2532: for (b = 0; b < Nb; ++b) {
2533: for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c];
2534: }
2535: }
2536: return PETSC_SUCCESS;
2537: }
2539: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2540: {
2541: PetscFEGeom pgeom;
2542: const PetscInt dEt = T->cdim;
2543: const PetscInt dE = fegeom->dimEmbed;
2544: const PetscInt Nq = T->Np;
2545: const PetscInt Nb = T->Nb;
2546: const PetscInt Nc = T->Nc;
2547: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2548: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt];
2549: PetscInt q, b, c, d;
2551: for (q = 0; q < Nq; ++q) {
2552: for (b = 0; b < Nb; ++b) {
2553: for (c = 0; c < Nc; ++c) {
2554: const PetscInt bcidx = b * Nc + c;
2556: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2557: for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d];
2558: for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0;
2559: }
2560: }
2561: PetscCall(PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom));
2562: PetscCall(PetscFEPushforward(fe, &pgeom, Nb, tmpBasis));
2563: PetscCall(PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer));
2564: for (b = 0; b < Nb; ++b) {
2565: for (c = 0; c < Nc; ++c) {
2566: const PetscInt bcidx = b * Nc + c;
2567: const PetscInt qcidx = q * Nc + c;
2569: elemVec[b] += tmpBasis[bcidx] * f0[qcidx];
2570: for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2571: }
2572: }
2573: }
2574: return PETSC_SUCCESS;
2575: }
2577: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt side, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2578: {
2579: const PetscInt dE = T->cdim;
2580: const PetscInt Nq = T->Np;
2581: const PetscInt Nb = T->Nb;
2582: const PetscInt Nc = T->Nc;
2583: const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc];
2584: const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE];
2586: for (PetscInt q = 0; q < Nq; ++q) {
2587: for (PetscInt b = 0; b < Nb; ++b) {
2588: for (PetscInt c = 0; c < Nc; ++c) {
2589: const PetscInt bcidx = b * Nc + c;
2591: tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx];
2592: for (PetscInt d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d];
2593: }
2594: }
2595: PetscCall(PetscFEPushforward(fe, fegeom, Nb, tmpBasis));
2596: // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2597: // PetscCall(PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer));
2598: if (side == 2) {
2599: // Integrating over whole cohesive cell, so insert for both sides
2600: for (PetscInt s = 0; s < 2; ++s) {
2601: for (PetscInt b = 0; b < Nb; ++b) {
2602: for (PetscInt c = 0; c < Nc; ++c) {
2603: const PetscInt bcidx = b * Nc + c;
2604: const PetscInt qcidx = (q * 2 + s) * Nc + c;
2606: elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx];
2607: for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2608: }
2609: }
2610: }
2611: } else {
2612: // Integrating over endcaps of cohesive cell, so insert for correct side
2613: for (PetscInt b = 0; b < Nb; ++b) {
2614: for (PetscInt c = 0; c < Nc; ++c) {
2615: const PetscInt bcidx = b * Nc + c;
2616: const PetscInt qcidx = q * Nc + c;
2618: elemVec[Nb * side + b] += tmpBasis[bcidx] * f0[qcidx];
2619: for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * side + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d];
2620: }
2621: }
2622: }
2623: }
2624: return PETSC_SUCCESS;
2625: }
2627: #define petsc_elemmat_kernel_g1(_NbI, _NcI, _NbJ, _NcJ, _dE) \
2628: do { \
2629: for (PetscInt fc = 0; fc < (_NcI); ++fc) { \
2630: for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \
2631: const PetscScalar *G = g1 + (fc * (_NcJ) + gc) * _dE; \
2632: for (PetscInt f = 0; f < (_NbI); ++f) { \
2633: const PetscScalar tBIv = tmpBasisI[f * (_NcI) + fc]; \
2634: for (PetscInt g = 0; g < (_NbJ); ++g) { \
2635: const PetscScalar *tBDJ = tmpBasisDerJ + (g * (_NcJ) + gc) * (_dE); \
2636: PetscScalar s = 0.0; \
2637: for (PetscInt df = 0; df < _dE; ++df) s += G[df] * tBDJ[df]; \
2638: elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s * tBIv; \
2639: } \
2640: } \
2641: } \
2642: } \
2643: } while (0)
2645: #define petsc_elemmat_kernel_g2(_NbI, _NcI, _NbJ, _NcJ, _dE) \
2646: do { \
2647: for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \
2648: for (PetscInt fc = 0; fc < (_NcI); ++fc) { \
2649: const PetscScalar *G = g2 + (fc * (_NcJ) + gc) * _dE; \
2650: for (PetscInt g = 0; g < (_NbJ); ++g) { \
2651: const PetscScalar tBJv = tmpBasisJ[g * (_NcJ) + gc]; \
2652: for (PetscInt f = 0; f < (_NbI); ++f) { \
2653: const PetscScalar *tBDI = tmpBasisDerI + (f * (_NcI) + fc) * (_dE); \
2654: PetscScalar s = 0.0; \
2655: for (PetscInt df = 0; df < _dE; ++df) s += tBDI[df] * G[df]; \
2656: elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s * tBJv; \
2657: } \
2658: } \
2659: } \
2660: } \
2661: } while (0)
2663: #define petsc_elemmat_kernel_g3(_NbI, _NcI, _NbJ, _NcJ, _dE) \
2664: do { \
2665: for (PetscInt fc = 0; fc < (_NcI); ++fc) { \
2666: for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \
2667: const PetscScalar *G = g3 + (fc * (_NcJ) + gc) * (_dE) * (_dE); \
2668: for (PetscInt f = 0; f < (_NbI); ++f) { \
2669: const PetscScalar *tBDI = tmpBasisDerI + (f * (_NcI) + fc) * (_dE); \
2670: for (PetscInt g = 0; g < (_NbJ); ++g) { \
2671: PetscScalar s = 0.0; \
2672: const PetscScalar *tBDJ = tmpBasisDerJ + (g * (_NcJ) + gc) * (_dE); \
2673: for (PetscInt df = 0; df < (_dE); ++df) { \
2674: for (PetscInt dg = 0; dg < (_dE); ++dg) s += tBDI[df] * G[df * (_dE) + dg] * tBDJ[dg]; \
2675: } \
2676: elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s; \
2677: } \
2678: } \
2679: } \
2680: } \
2681: } while (0)
2683: 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 totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2684: {
2685: const PetscInt cdim = TI->cdim;
2686: const PetscInt dE = fegeom->dimEmbed;
2687: const PetscInt NqI = TI->Np;
2688: const PetscInt NbI = TI->Nb;
2689: const PetscInt NcI = TI->Nc;
2690: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2691: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * cdim];
2692: const PetscInt NqJ = TJ->Np;
2693: const PetscInt NbJ = TJ->Nb;
2694: const PetscInt NcJ = TJ->Nc;
2695: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2696: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * cdim];
2698: for (PetscInt f = 0; f < NbI; ++f) {
2699: for (PetscInt fc = 0; fc < NcI; ++fc) {
2700: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2702: tmpBasisI[fidx] = basisI[fidx];
2703: for (PetscInt df = 0; df < cdim; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * cdim + df];
2704: }
2705: }
2706: PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2707: PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2708: if (feI != feJ) {
2709: for (PetscInt g = 0; g < NbJ; ++g) {
2710: for (PetscInt gc = 0; gc < NcJ; ++gc) {
2711: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2713: tmpBasisJ[gidx] = basisJ[gidx];
2714: for (PetscInt dg = 0; dg < cdim; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * cdim + dg];
2715: }
2716: }
2717: PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2718: PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2719: } else {
2720: tmpBasisJ = tmpBasisI;
2721: tmpBasisDerJ = tmpBasisDerI;
2722: }
2723: if (PetscUnlikely(g0)) {
2724: for (PetscInt f = 0; f < NbI; ++f) {
2725: const PetscInt i = offsetI + f; /* Element matrix row */
2727: for (PetscInt fc = 0; fc < NcI; ++fc) {
2728: const PetscScalar bI = tmpBasisI[f * NcI + fc]; /* Test function basis value */
2730: for (PetscInt g = 0; g < NbJ; ++g) {
2731: const PetscInt j = offsetJ + g; /* Element matrix column */
2732: const PetscInt fOff = i * totDim + j;
2734: for (PetscInt gc = 0; gc < NcJ; ++gc) elemMat[fOff] += bI * g0[fc * NcJ + gc] * tmpBasisJ[g * NcJ + gc];
2735: }
2736: }
2737: }
2738: }
2739: if (PetscUnlikely(g1)) {
2740: #if 1
2741: if (dE == 2) {
2742: petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, 2);
2743: } else if (dE == 3) {
2744: petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, 3);
2745: } else {
2746: petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, dE);
2747: }
2748: #else
2749: for (PetscInt f = 0; f < NbI; ++f) {
2750: const PetscInt i = offsetI + f; /* Element matrix row */
2752: for (PetscInt fc = 0; fc < NcI; ++fc) {
2753: const PetscScalar bI = tmpBasisI[f * NcI + fc]; /* Test function basis value */
2755: for (PetscInt g = 0; g < NbJ; ++g) {
2756: const PetscInt j = offsetJ + g; /* Element matrix column */
2757: const PetscInt fOff = i * totDim + j;
2759: for (PetscInt gc = 0; gc < NcJ; ++gc) {
2760: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2762: for (PetscInt df = 0; df < dE; ++df) elemMat[fOff] += bI * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2763: }
2764: }
2765: }
2766: }
2767: #endif
2768: }
2769: if (PetscUnlikely(g2)) {
2770: #if 1
2771: if (dE == 2) {
2772: petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, 2);
2773: } else if (dE == 3) {
2774: petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, 3);
2775: } else {
2776: petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, dE);
2777: }
2778: #else
2779: for (PetscInt g = 0; g < NbJ; ++g) {
2780: const PetscInt j = offsetJ + g; /* Element matrix column */
2782: for (PetscInt gc = 0; gc < NcJ; ++gc) {
2783: const PetscScalar bJ = tmpBasisJ[g * NcJ + gc]; /* Trial function basis value */
2785: for (PetscInt f = 0; f < NbI; ++f) {
2786: const PetscInt i = offsetI + f; /* Element matrix row */
2787: const PetscInt fOff = i * totDim + j;
2789: for (PetscInt fc = 0; fc < NcI; ++fc) {
2790: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2792: for (PetscInt df = 0; df < dE; ++df) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * bJ;
2793: }
2794: }
2795: }
2796: }
2797: #endif
2798: }
2799: if (PetscUnlikely(g3)) {
2800: #if 1
2801: if (dE == 2) {
2802: petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, 2);
2803: } else if (dE == 3) {
2804: petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, 3);
2805: } else {
2806: petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, dE);
2807: }
2808: #else
2809: for (PetscInt f = 0; f < NbI; ++f) {
2810: const PetscInt i = offsetI + f; /* Element matrix row */
2812: for (PetscInt fc = 0; fc < NcI; ++fc) {
2813: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2815: for (PetscInt g = 0; g < NbJ; ++g) {
2816: const PetscInt j = offsetJ + g; /* Element matrix column */
2817: const PetscInt fOff = i * totDim + j;
2819: for (PetscInt gc = 0; gc < NcJ; ++gc) {
2820: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2822: for (PetscInt df = 0; df < dE; ++df) {
2823: for (PetscInt dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2824: }
2825: }
2826: }
2827: }
2828: }
2829: #endif
2830: }
2831: return PETSC_SUCCESS;
2832: }
2834: #undef petsc_elemmat_kernel_g1
2835: #undef petsc_elemmat_kernel_g2
2836: #undef petsc_elemmat_kernel_g3
2838: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt s, PetscInt t, 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[])
2839: {
2840: const PetscInt dE = TI->cdim;
2841: const PetscInt NqI = TI->Np;
2842: const PetscInt NbI = TI->Nb;
2843: const PetscInt NcI = TI->Nc;
2844: const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI];
2845: const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE];
2846: const PetscInt NqJ = TJ->Np;
2847: const PetscInt NbJ = TJ->Nb;
2848: const PetscInt NcJ = TJ->Nc;
2849: const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ];
2850: const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE];
2851: const PetscInt so = isHybridI ? 0 : s;
2852: const PetscInt to = isHybridJ ? 0 : t;
2853: PetscInt f, fc, g, gc, df, dg;
2855: for (f = 0; f < NbI; ++f) {
2856: for (fc = 0; fc < NcI; ++fc) {
2857: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2859: tmpBasisI[fidx] = basisI[fidx];
2860: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df];
2861: }
2862: }
2863: PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI));
2864: PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI));
2865: for (g = 0; g < NbJ; ++g) {
2866: for (gc = 0; gc < NcJ; ++gc) {
2867: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2869: tmpBasisJ[gidx] = basisJ[gidx];
2870: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg];
2871: }
2872: }
2873: PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ));
2874: // TODO This is currently broken since we do not pull the geometry down to the lower dimension
2875: // PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ));
2876: for (f = 0; f < NbI; ++f) {
2877: for (fc = 0; fc < NcI; ++fc) {
2878: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
2879: const PetscInt i = offsetI + NbI * so + f; /* Element matrix row */
2880: for (g = 0; g < NbJ; ++g) {
2881: for (gc = 0; gc < NcJ; ++gc) {
2882: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
2883: const PetscInt j = offsetJ + NbJ * to + g; /* Element matrix column */
2884: const PetscInt fOff = eOffset + i * totDim + j;
2886: elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx];
2887: for (df = 0; df < dE; ++df) {
2888: elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df];
2889: elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx];
2890: for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg];
2891: }
2892: }
2893: }
2894: }
2895: }
2896: return PETSC_SUCCESS;
2897: }
2899: /*@
2900: PetscFECreateCellGeometry - Populates the arrays in a `PetscFEGeom` for a single reference cell of a `PetscFE`.
2902: Not Collective
2904: Input Parameters:
2905: + fe - the `PetscFE` whose dual-space `DM` provides the reference cell
2906: - quad - the quadrature at which to evaluate the geometry, or `NULL` to use the `PetscFE`'s own quadrature
2908: Output Parameter:
2909: . cgeom - the `PetscFEGeom` populated with reference-cell coordinates, Jacobians, inverse Jacobians, and their determinants
2911: Level: developer
2913: Notes:
2914: This does not create `cgeom`, it allocates the arrays within one
2916: Free the storage with `PetscFEDestroyCellGeometry()`.
2918: .seealso: `PetscFE`, `PetscFEGeom`, `PetscFEDestroyCellGeometry()`, `PetscFEGetQuadrature()`, `DMPlexComputeCellGeometryFEM()`
2919: @*/
2920: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2921: {
2922: PetscDualSpace dsp;
2923: DM dm;
2924: PetscQuadrature quadDef;
2925: PetscInt dim, cdim, Nq;
2927: PetscFunctionBegin;
2928: PetscCall(PetscFEGetDualSpace(fe, &dsp));
2929: PetscCall(PetscDualSpaceGetDM(dsp, &dm));
2930: PetscCall(DMGetDimension(dm, &dim));
2931: PetscCall(DMGetCoordinateDim(dm, &cdim));
2932: PetscCall(PetscFEGetQuadrature(fe, &quadDef));
2933: quad = quad ? quad : quadDef;
2934: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
2935: PetscCall(PetscMalloc1(Nq * cdim, &cgeom->v));
2936: PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->J));
2937: PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ));
2938: PetscCall(PetscMalloc1(Nq, &cgeom->detJ));
2939: cgeom->dim = dim;
2940: cgeom->dimEmbed = cdim;
2941: cgeom->numCells = 1;
2942: cgeom->numPoints = Nq;
2943: PetscCall(DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ));
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: /*@
2948: PetscFEDestroyCellGeometry - Free the arrays inside a `PetscFEGeom` allocated by `PetscFECreateCellGeometry()`.
2950: Not Collective
2952: Input Parameters:
2953: + fe - the `PetscFE` (unused, kept for API symmetry with `PetscFECreateCellGeometry()`)
2954: - cgeom - the `PetscFEGeom` whose owned arrays should be freed
2956: Level: developer
2958: .seealso: `PetscFE`, `PetscFEGeom`, `PetscFECreateCellGeometry()`
2959: @*/
2960: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2961: {
2962: PetscFunctionBegin;
2963: PetscCall(PetscFree(cgeom->v));
2964: PetscCall(PetscFree(cgeom->J));
2965: PetscCall(PetscFree(cgeom->invJ));
2966: PetscCall(PetscFree(cgeom->detJ));
2967: PetscFunctionReturn(PETSC_SUCCESS);
2968: }
2970: #if 0
2971: PetscErrorCode PetscFEUpdateElementMat_Internal_SparseIndices(PetscTabulation TI, PetscTabulation TJ, PetscInt dimEmbed, const PetscInt g0[], const PetscInt g1[], const PetscInt g2[], const PetscInt g3[], PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscInt *n_g0, PetscInt **g0_idxs_out, PetscInt *n_g1, PetscInt **g1_idxs_out, PetscInt *n_g2, PetscInt **g2_idxs_out, PetscInt *n_g3, PetscInt **g3_idxs_out)
2972: {
2973: const PetscInt dE = dimEmbed;
2974: const PetscInt NbI = TI->Nb;
2975: const PetscInt NcI = TI->Nc;
2976: const PetscInt NbJ = TJ->Nb;
2977: const PetscInt NcJ = TJ->Nc;
2978: PetscBool has_g0 = g0 ? PETSC_TRUE : PETSC_FALSE;
2979: PetscBool has_g1 = g1 ? PETSC_TRUE : PETSC_FALSE;
2980: PetscBool has_g2 = g2 ? PETSC_TRUE : PETSC_FALSE;
2981: PetscBool has_g3 = g3 ? PETSC_TRUE : PETSC_FALSE;
2982: PetscInt *g0_idxs = NULL, *g1_idxs = NULL, *g2_idxs = NULL, *g3_idxs = NULL;
2983: PetscInt g0_i, g1_i, g2_i, g3_i;
2985: PetscFunctionBegin;
2986: g0_i = g1_i = g2_i = g3_i = 0;
2987: if (has_g0)
2988: for (PetscInt i = 0; i < NcI * NcJ; i++)
2989: if (g0[i]) g0_i += NbI * NbJ;
2990: if (has_g1)
2991: for (PetscInt i = 0; i < NcI * NcJ * dE; i++)
2992: if (g1[i]) g1_i += NbI * NbJ;
2993: if (has_g2)
2994: for (PetscInt i = 0; i < NcI * NcJ * dE; i++)
2995: if (g2[i]) g2_i += NbI * NbJ;
2996: if (has_g3)
2997: for (PetscInt i = 0; i < NcI * NcJ * dE * dE; i++)
2998: if (g3[i]) g3_i += NbI * NbJ;
2999: if (g0_i == NbI * NbJ * NcI * NcJ) g0_i = 0;
3000: if (g1_i == NbI * NbJ * NcI * NcJ * dE) g1_i = 0;
3001: if (g2_i == NbI * NbJ * NcI * NcJ * dE) g2_i = 0;
3002: if (g3_i == NbI * NbJ * NcI * NcJ * dE * dE) g3_i = 0;
3003: has_g0 = g0_i ? PETSC_TRUE : PETSC_FALSE;
3004: has_g1 = g1_i ? PETSC_TRUE : PETSC_FALSE;
3005: has_g2 = g2_i ? PETSC_TRUE : PETSC_FALSE;
3006: has_g3 = g3_i ? PETSC_TRUE : PETSC_FALSE;
3007: if (has_g0) PetscCall(PetscMalloc1(4 * g0_i, &g0_idxs));
3008: if (has_g1) PetscCall(PetscMalloc1(4 * g1_i, &g1_idxs));
3009: if (has_g2) PetscCall(PetscMalloc1(4 * g2_i, &g2_idxs));
3010: if (has_g3) PetscCall(PetscMalloc1(4 * g3_i, &g3_idxs));
3011: g0_i = g1_i = g2_i = g3_i = 0;
3013: for (PetscInt f = 0; f < NbI; ++f) {
3014: const PetscInt i = offsetI + f; /* Element matrix row */
3015: for (PetscInt fc = 0; fc < NcI; ++fc) {
3016: const PetscInt fidx = f * NcI + fc; /* Test function basis index */
3018: for (PetscInt g = 0; g < NbJ; ++g) {
3019: const PetscInt j = offsetJ + g; /* Element matrix column */
3020: const PetscInt fOff = i * totDim + j;
3021: for (PetscInt gc = 0; gc < NcJ; ++gc) {
3022: const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */
3024: if (has_g0) {
3025: if (g0[fc * NcJ + gc]) {
3026: g0_idxs[4 * g0_i + 0] = fidx;
3027: g0_idxs[4 * g0_i + 1] = fc * NcJ + gc;
3028: g0_idxs[4 * g0_i + 2] = gidx;
3029: g0_idxs[4 * g0_i + 3] = fOff;
3030: g0_i++;
3031: }
3032: }
3034: for (PetscInt df = 0; df < dE; ++df) {
3035: if (has_g1) {
3036: if (g1[(fc * NcJ + gc) * dE + df]) {
3037: g1_idxs[4 * g1_i + 0] = fidx;
3038: g1_idxs[4 * g1_i + 1] = (fc * NcJ + gc) * dE + df;
3039: g1_idxs[4 * g1_i + 2] = gidx * dE + df;
3040: g1_idxs[4 * g1_i + 3] = fOff;
3041: g1_i++;
3042: }
3043: }
3044: if (has_g2) {
3045: if (g2[(fc * NcJ + gc) * dE + df]) {
3046: g2_idxs[4 * g2_i + 0] = fidx * dE + df;
3047: g2_idxs[4 * g2_i + 1] = (fc * NcJ + gc) * dE + df;
3048: g2_idxs[4 * g2_i + 2] = gidx;
3049: g2_idxs[4 * g2_i + 3] = fOff;
3050: g2_i++;
3051: }
3052: }
3053: if (has_g3) {
3054: for (PetscInt dg = 0; dg < dE; ++dg) {
3055: if (g3[((fc * NcJ + gc) * dE + df) * dE + dg]) {
3056: g3_idxs[4 * g3_i + 0] = fidx * dE + df;
3057: g3_idxs[4 * g3_i + 1] = ((fc * NcJ + gc) * dE + df) * dE + dg;
3058: g3_idxs[4 * g3_i + 2] = gidx * dE + dg;
3059: g3_idxs[4 * g3_i + 3] = fOff;
3060: g3_i++;
3061: }
3062: }
3063: }
3064: }
3065: }
3066: }
3067: }
3068: }
3069: *n_g0 = g0_i;
3070: *n_g1 = g1_i;
3071: *n_g2 = g2_i;
3072: *n_g3 = g3_i;
3074: *g0_idxs_out = g0_idxs;
3075: *g1_idxs_out = g1_idxs;
3076: *g2_idxs_out = g2_idxs;
3077: *g3_idxs_out = g3_idxs;
3078: PetscFunctionReturn(PETSC_SUCCESS);
3079: }
3081: //example HOW TO USE
3082: for (PetscInt i = 0; i < g0_sparse_n; i++) {
3083: PetscInt bM = g0_sparse_idxs[4 * i + 0];
3084: PetscInt bN = g0_sparse_idxs[4 * i + 1];
3085: PetscInt bK = g0_sparse_idxs[4 * i + 2];
3086: PetscInt bO = g0_sparse_idxs[4 * i + 3];
3087: elemMat[bO] += tmpBasisI[bM] * g0[bN] * tmpBasisJ[bK];
3088: }
3089: #endif