Actual source code: fe.c
petsc-3.12.5 2020-03-29
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: PetscFunctionList PetscFEList = NULL;
54: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
56: /*@C
57: PetscFERegister - Adds a new PetscFE implementation
59: Not Collective
61: Input Parameters:
62: + name - The name of a new user-defined creation routine
63: - create_func - The creation routine itself
65: Notes:
66: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
68: Sample usage:
69: .vb
70: PetscFERegister("my_fe", MyPetscFECreate);
71: .ve
73: Then, your PetscFE type can be chosen with the procedural interface via
74: .vb
75: PetscFECreate(MPI_Comm, PetscFE *);
76: PetscFESetType(PetscFE, "my_fe");
77: .ve
78: or at runtime via the option
79: .vb
80: -petscfe_type my_fe
81: .ve
83: Level: advanced
85: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
87: @*/
88: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
89: {
93: PetscFunctionListAdd(&PetscFEList, sname, function);
94: return(0);
95: }
97: /*@C
98: PetscFESetType - Builds a particular PetscFE
100: Collective on fem
102: Input Parameters:
103: + fem - The PetscFE object
104: - name - The kind of FEM space
106: Options Database Key:
107: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
109: Level: intermediate
111: .seealso: PetscFEGetType(), PetscFECreate()
112: @*/
113: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
114: {
115: PetscErrorCode (*r)(PetscFE);
116: PetscBool match;
121: PetscObjectTypeCompare((PetscObject) fem, name, &match);
122: if (match) return(0);
124: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
125: PetscFunctionListFind(PetscFEList, name, &r);
126: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
128: if (fem->ops->destroy) {
129: (*fem->ops->destroy)(fem);
130: fem->ops->destroy = NULL;
131: }
132: (*r)(fem);
133: PetscObjectChangeTypeName((PetscObject) fem, name);
134: return(0);
135: }
137: /*@C
138: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
140: Not Collective
142: Input Parameter:
143: . fem - The PetscFE
145: Output Parameter:
146: . name - The PetscFE type name
148: Level: intermediate
150: .seealso: PetscFESetType(), PetscFECreate()
151: @*/
152: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
153: {
159: if (!PetscFERegisterAllCalled) {
160: PetscFERegisterAll();
161: }
162: *name = ((PetscObject) fem)->type_name;
163: return(0);
164: }
166: /*@C
167: PetscFEView - Views a PetscFE
169: Collective on fem
171: Input Parameter:
172: + fem - the PetscFE object to view
173: - viewer - the viewer
175: Level: beginner
177: .seealso PetscFEDestroy()
178: @*/
179: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
180: {
181: PetscBool iascii;
187: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
188: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
189: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
190: if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
191: return(0);
192: }
194: /*@
195: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
197: Collective on fem
199: Input Parameter:
200: . fem - the PetscFE object to set options for
202: Options Database:
203: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
204: - -petscfe_num_batches - the number of cell batches to integrate serially
206: Level: intermediate
208: .seealso PetscFEView()
209: @*/
210: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
211: {
212: const char *defaultType;
213: char name[256];
214: PetscBool flg;
219: if (!((PetscObject) fem)->type_name) {
220: defaultType = PETSCFEBASIC;
221: } else {
222: defaultType = ((PetscObject) fem)->type_name;
223: }
224: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
226: PetscObjectOptionsBegin((PetscObject) fem);
227: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
228: if (flg) {
229: PetscFESetType(fem, name);
230: } else if (!((PetscObject) fem)->type_name) {
231: PetscFESetType(fem, defaultType);
232: }
233: PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
234: PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
235: if (fem->ops->setfromoptions) {
236: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
237: }
238: /* process any options handlers added with PetscObjectAddOptionsHandler() */
239: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
240: PetscOptionsEnd();
241: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
242: return(0);
243: }
245: /*@C
246: PetscFESetUp - Construct data structures for the PetscFE
248: Collective on fem
250: Input Parameter:
251: . fem - the PetscFE object to setup
253: Level: intermediate
255: .seealso PetscFEView(), PetscFEDestroy()
256: @*/
257: PetscErrorCode PetscFESetUp(PetscFE fem)
258: {
263: if (fem->setupcalled) return(0);
264: fem->setupcalled = PETSC_TRUE;
265: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
266: return(0);
267: }
269: /*@
270: PetscFEDestroy - Destroys a PetscFE object
272: Collective on fem
274: Input Parameter:
275: . fem - the PetscFE object to destroy
277: Level: beginner
279: .seealso PetscFEView()
280: @*/
281: PetscErrorCode PetscFEDestroy(PetscFE *fem)
282: {
286: if (!*fem) return(0);
289: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
290: ((PetscObject) (*fem))->refct = 0;
292: if ((*fem)->subspaces) {
293: PetscInt dim, d;
295: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
296: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
297: }
298: PetscFree((*fem)->subspaces);
299: PetscFree((*fem)->invV);
300: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
301: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
302: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
303: PetscSpaceDestroy(&(*fem)->basisSpace);
304: PetscDualSpaceDestroy(&(*fem)->dualSpace);
305: PetscQuadratureDestroy(&(*fem)->quadrature);
306: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
308: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
309: PetscHeaderDestroy(fem);
310: return(0);
311: }
313: /*@
314: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
316: Collective
318: Input Parameter:
319: . comm - The communicator for the PetscFE object
321: Output Parameter:
322: . fem - The PetscFE object
324: Level: beginner
326: .seealso: PetscFESetType(), PETSCFEGALERKIN
327: @*/
328: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
329: {
330: PetscFE f;
335: PetscCitationsRegister(FECitation,&FEcite);
336: *fem = NULL;
337: PetscFEInitializePackage();
339: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
341: f->basisSpace = NULL;
342: f->dualSpace = NULL;
343: f->numComponents = 1;
344: f->subspaces = NULL;
345: f->invV = NULL;
346: f->B = NULL;
347: f->D = NULL;
348: f->H = NULL;
349: f->Bf = NULL;
350: f->Df = NULL;
351: f->Hf = NULL;
352: PetscArrayzero(&f->quadrature, 1);
353: PetscArrayzero(&f->faceQuadrature, 1);
354: f->blockSize = 0;
355: f->numBlocks = 1;
356: f->batchSize = 0;
357: f->numBatches = 1;
359: *fem = f;
360: return(0);
361: }
363: /*@
364: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
366: Not collective
368: Input Parameter:
369: . fem - The PetscFE object
371: Output Parameter:
372: . dim - The spatial dimension
374: Level: intermediate
376: .seealso: PetscFECreate()
377: @*/
378: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
379: {
380: DM dm;
386: PetscDualSpaceGetDM(fem->dualSpace, &dm);
387: DMGetDimension(dm, dim);
388: return(0);
389: }
391: /*@
392: PetscFESetNumComponents - Sets the number of components in the element
394: Not collective
396: Input Parameters:
397: + fem - The PetscFE object
398: - comp - The number of field components
400: Level: intermediate
402: .seealso: PetscFECreate()
403: @*/
404: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
405: {
408: fem->numComponents = comp;
409: return(0);
410: }
412: /*@
413: PetscFEGetNumComponents - Returns the number of components in the element
415: Not collective
417: Input Parameter:
418: . fem - The PetscFE object
420: Output Parameter:
421: . comp - The number of field components
423: Level: intermediate
425: .seealso: PetscFECreate()
426: @*/
427: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
428: {
432: *comp = fem->numComponents;
433: return(0);
434: }
436: /*@
437: PetscFESetTileSizes - Sets the tile sizes for evaluation
439: Not collective
441: Input Parameters:
442: + fem - The PetscFE object
443: . blockSize - The number of elements in a block
444: . numBlocks - The number of blocks in a batch
445: . batchSize - The number of elements in a batch
446: - numBatches - The number of batches in a chunk
448: Level: intermediate
450: .seealso: PetscFECreate()
451: @*/
452: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
453: {
456: fem->blockSize = blockSize;
457: fem->numBlocks = numBlocks;
458: fem->batchSize = batchSize;
459: fem->numBatches = numBatches;
460: return(0);
461: }
463: /*@
464: PetscFEGetTileSizes - Returns the tile sizes for evaluation
466: Not collective
468: Input Parameter:
469: . fem - The PetscFE object
471: Output Parameters:
472: + blockSize - The number of elements in a block
473: . numBlocks - The number of blocks in a batch
474: . batchSize - The number of elements in a batch
475: - numBatches - The number of batches in a chunk
477: Level: intermediate
479: .seealso: PetscFECreate()
480: @*/
481: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
482: {
489: if (blockSize) *blockSize = fem->blockSize;
490: if (numBlocks) *numBlocks = fem->numBlocks;
491: if (batchSize) *batchSize = fem->batchSize;
492: if (numBatches) *numBatches = fem->numBatches;
493: return(0);
494: }
496: /*@
497: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
499: Not collective
501: Input Parameter:
502: . fem - The PetscFE object
504: Output Parameter:
505: . sp - The PetscSpace object
507: Level: intermediate
509: .seealso: PetscFECreate()
510: @*/
511: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
512: {
516: *sp = fem->basisSpace;
517: return(0);
518: }
520: /*@
521: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
523: Not collective
525: Input Parameters:
526: + fem - The PetscFE object
527: - sp - The PetscSpace object
529: Level: intermediate
531: .seealso: PetscFECreate()
532: @*/
533: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
534: {
540: PetscSpaceDestroy(&fem->basisSpace);
541: fem->basisSpace = sp;
542: PetscObjectReference((PetscObject) fem->basisSpace);
543: return(0);
544: }
546: /*@
547: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
549: Not collective
551: Input Parameter:
552: . fem - The PetscFE object
554: Output Parameter:
555: . sp - The PetscDualSpace object
557: Level: intermediate
559: .seealso: PetscFECreate()
560: @*/
561: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
562: {
566: *sp = fem->dualSpace;
567: return(0);
568: }
570: /*@
571: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
573: Not collective
575: Input Parameters:
576: + fem - The PetscFE object
577: - sp - The PetscDualSpace object
579: Level: intermediate
581: .seealso: PetscFECreate()
582: @*/
583: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
584: {
590: PetscDualSpaceDestroy(&fem->dualSpace);
591: fem->dualSpace = sp;
592: PetscObjectReference((PetscObject) fem->dualSpace);
593: return(0);
594: }
596: /*@
597: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
599: Not collective
601: Input Parameter:
602: . fem - The PetscFE object
604: Output Parameter:
605: . q - The PetscQuadrature object
607: Level: intermediate
609: .seealso: PetscFECreate()
610: @*/
611: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
612: {
616: *q = fem->quadrature;
617: return(0);
618: }
620: /*@
621: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
623: Not collective
625: Input Parameters:
626: + fem - The PetscFE object
627: - q - The PetscQuadrature object
629: Level: intermediate
631: .seealso: PetscFECreate()
632: @*/
633: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
634: {
635: PetscInt Nc, qNc;
640: PetscFEGetNumComponents(fem, &Nc);
641: PetscQuadratureGetNumComponents(q, &qNc);
642: if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
643: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
644: PetscQuadratureDestroy(&fem->quadrature);
645: fem->quadrature = q;
646: PetscObjectReference((PetscObject) q);
647: return(0);
648: }
650: /*@
651: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
653: Not collective
655: Input Parameter:
656: . fem - The PetscFE object
658: Output Parameter:
659: . q - The PetscQuadrature object
661: Level: intermediate
663: .seealso: PetscFECreate()
664: @*/
665: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
666: {
670: *q = fem->faceQuadrature;
671: return(0);
672: }
674: /*@
675: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
677: Not collective
679: Input Parameters:
680: + fem - The PetscFE object
681: - q - The PetscQuadrature object
683: Level: intermediate
685: .seealso: PetscFECreate()
686: @*/
687: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
688: {
693: PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
694: PetscQuadratureDestroy(&fem->faceQuadrature);
695: fem->faceQuadrature = q;
696: PetscObjectReference((PetscObject) q);
697: return(0);
698: }
700: /*@
701: PetscFECopyQuadrature - Copy both volumetric and surface quadrature
703: Not collective
705: Input Parameters:
706: + sfe - The PetscFE source for the quadratures
707: - tfe - The PetscFE target for the quadratures
709: Level: intermediate
711: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
712: @*/
713: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
714: {
715: PetscQuadrature q;
716: PetscErrorCode ierr;
721: PetscFEGetQuadrature(sfe, &q);
722: PetscFESetQuadrature(tfe, q);
723: PetscFEGetFaceQuadrature(sfe, &q);
724: PetscFESetFaceQuadrature(tfe, q);
725: return(0);
726: }
728: /*@C
729: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
731: Not collective
733: Input Parameter:
734: . fem - The PetscFE object
736: Output Parameter:
737: . numDof - Array with the number of dofs per dimension
739: Level: intermediate
741: .seealso: PetscFECreate()
742: @*/
743: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
744: {
750: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
751: return(0);
752: }
754: /*@C
755: PetscFEGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points
757: Not collective
759: Input Parameter:
760: . fem - The PetscFE object
762: Output Parameters:
763: + B - The basis function values at quadrature points
764: . D - The basis function derivatives at quadrature points
765: - H - The basis function second derivatives at quadrature points
767: Note:
768: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
769: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
770: $ H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
772: Level: intermediate
774: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
775: @*/
776: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
777: {
778: PetscInt npoints;
779: const PetscReal *points;
780: PetscErrorCode ierr;
787: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
788: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
789: if (B) *B = fem->B;
790: if (D) *D = fem->D;
791: if (H) *H = fem->H;
792: return(0);
793: }
795: /*@C
796: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points
798: Not collective
800: Input Parameter:
801: . fem - The PetscFE object
803: Output Parameters:
804: + B - The basis function values at face quadrature points
805: . D - The basis function derivatives at face quadrature points
806: - H - The basis function second derivatives at face quadrature points
808: Note:
809: $ Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
810: $ 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
811: $ Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e
813: Level: intermediate
815: .seealso: PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
816: @*/
817: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
818: {
819: PetscErrorCode ierr;
826: if (!fem->Bf) {
827: const PetscReal xi0[3] = {-1., -1., -1.};
828: PetscReal v0[3], J[9], detJ;
829: PetscQuadrature fq;
830: PetscDualSpace sp;
831: DM dm;
832: const PetscInt *faces;
833: PetscInt dim, numFaces, f, npoints, q;
834: const PetscReal *points;
835: PetscReal *facePoints;
837: PetscFEGetDualSpace(fem, &sp);
838: PetscDualSpaceGetDM(sp, &dm);
839: DMGetDimension(dm, &dim);
840: DMPlexGetConeSize(dm, 0, &numFaces);
841: DMPlexGetCone(dm, 0, &faces);
842: PetscFEGetFaceQuadrature(fem, &fq);
843: if (fq) {
844: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
845: PetscMalloc1(numFaces*npoints*dim, &facePoints);
846: for (f = 0; f < numFaces; ++f) {
847: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
848: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
849: }
850: PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
851: PetscFree(facePoints);
852: }
853: }
854: if (Bf) *Bf = fem->Bf;
855: if (Df) *Df = fem->Df;
856: if (Hf) *Hf = fem->Hf;
857: return(0);
858: }
860: /*@C
861: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face centroid points
863: Not collective
865: Input Parameter:
866: . fem - The PetscFE object
868: Output Parameters:
869: + B - The basis function values at face centroid points
870: . D - The basis function derivatives at face centroid points
871: - H - The basis function second derivatives at face centroid points
873: Note:
874: $ Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
875: $ Df[((f*pdim + i)*Nc + c)*dim + d] is the derivative value at point f for basis function i, component c, in direction d
876: $ Hf[(((f*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f for basis function i, component c, in directions d and e
878: Level: intermediate
880: .seealso: PetscFEGetFaceTabulation(), PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation()
881: @*/
882: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
883: {
884: PetscErrorCode ierr;
889: if (!fem->F) {
890: PetscDualSpace sp;
891: DM dm;
892: const PetscInt *cone;
893: PetscReal *centroids;
894: PetscInt dim, numFaces, f;
896: PetscFEGetDualSpace(fem, &sp);
897: PetscDualSpaceGetDM(sp, &dm);
898: DMGetDimension(dm, &dim);
899: DMPlexGetConeSize(dm, 0, &numFaces);
900: DMPlexGetCone(dm, 0, &cone);
901: PetscMalloc1(numFaces*dim, ¢roids);
902: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
903: PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
904: PetscFree(centroids);
905: }
906: *F = fem->F;
907: return(0);
908: }
910: /*@C
911: PetscFEGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
913: Not collective
915: Input Parameters:
916: + fem - The PetscFE object
917: . npoints - The number of tabulation points
918: - points - The tabulation point coordinates
920: Output Parameters:
921: + B - The basis function values at tabulation points
922: . D - The basis function derivatives at tabulation points
923: - H - The basis function second derivatives at tabulation points
925: Note:
926: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
927: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
928: $ H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
930: Level: intermediate
932: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
933: @*/
934: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
935: {
936: DM dm;
937: PetscInt pdim; /* Dimension of FE space P */
938: PetscInt dim; /* Spatial dimension */
939: PetscInt comp; /* Field components */
940: PetscErrorCode ierr;
943: if (!npoints) {
944: if (B) *B = NULL;
945: if (D) *D = NULL;
946: if (H) *H = NULL;
947: return(0);
948: }
954: PetscDualSpaceGetDM(fem->dualSpace, &dm);
955: DMGetDimension(dm, &dim);
956: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
957: PetscFEGetNumComponents(fem, &comp);
958: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
959: if (!dim) {
960: if (D) *D = NULL;
961: if (H) *H = NULL;
962: } else {
963: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
964: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
965: }
966: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
967: return(0);
968: }
970: /*@C
971: PetscFERestoreTabulation - Frees memory from the associated tabulation.
973: Not collective
975: Input Parameters:
976: + fem - The PetscFE object
977: . npoints - The number of tabulation points
978: . points - The tabulation point coordinates
979: . B - The basis function values at tabulation points
980: . D - The basis function derivatives at tabulation points
981: - H - The basis function second derivatives at tabulation points
983: Note:
984: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
985: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
986: $ H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
988: Level: intermediate
990: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation()
991: @*/
992: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
993: {
994: DM dm;
999: PetscDualSpaceGetDM(fem->dualSpace, &dm);
1000: if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
1001: if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
1002: if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
1003: return(0);
1004: }
1006: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1007: {
1008: PetscSpace bsp, bsubsp;
1009: PetscDualSpace dsp, dsubsp;
1010: PetscInt dim, depth, numComp, i, j, coneSize, order;
1011: PetscFEType type;
1012: DM dm;
1013: DMLabel label;
1014: PetscReal *xi, *v, *J, detJ;
1015: const char *name;
1016: PetscQuadrature origin, fullQuad, subQuad;
1022: PetscFEGetBasisSpace(fe,&bsp);
1023: PetscFEGetDualSpace(fe,&dsp);
1024: PetscDualSpaceGetDM(dsp,&dm);
1025: DMGetDimension(dm,&dim);
1026: DMPlexGetDepthLabel(dm,&label);
1027: DMLabelGetValue(label,refPoint,&depth);
1028: PetscCalloc1(depth,&xi);
1029: PetscMalloc1(dim,&v);
1030: PetscMalloc1(dim*dim,&J);
1031: for (i = 0; i < depth; i++) xi[i] = 0.;
1032: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1033: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1034: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1035: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1036: for (i = 1; i < dim; i++) {
1037: for (j = 0; j < depth; j++) {
1038: J[i * depth + j] = J[i * dim + j];
1039: }
1040: }
1041: PetscQuadratureDestroy(&origin);
1042: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1043: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1044: PetscSpaceSetUp(bsubsp);
1045: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1046: PetscFEGetType(fe,&type);
1047: PetscFESetType(*trFE,type);
1048: PetscFEGetNumComponents(fe,&numComp);
1049: PetscFESetNumComponents(*trFE,numComp);
1050: PetscFESetBasisSpace(*trFE,bsubsp);
1051: PetscFESetDualSpace(*trFE,dsubsp);
1052: PetscObjectGetName((PetscObject) fe, &name);
1053: if (name) {PetscFESetName(*trFE, name);}
1054: PetscFEGetQuadrature(fe,&fullQuad);
1055: PetscQuadratureGetOrder(fullQuad,&order);
1056: DMPlexGetConeSize(dm,refPoint,&coneSize);
1057: if (coneSize == 2 * depth) {
1058: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1059: } else {
1060: PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1061: }
1062: PetscFESetQuadrature(*trFE,subQuad);
1063: PetscFESetUp(*trFE);
1064: PetscQuadratureDestroy(&subQuad);
1065: PetscSpaceDestroy(&bsubsp);
1066: return(0);
1067: }
1069: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1070: {
1071: PetscInt hStart, hEnd;
1072: PetscDualSpace dsp;
1073: DM dm;
1079: *trFE = NULL;
1080: PetscFEGetDualSpace(fe,&dsp);
1081: PetscDualSpaceGetDM(dsp,&dm);
1082: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1083: if (hEnd <= hStart) return(0);
1084: PetscFECreatePointTrace(fe,hStart,trFE);
1085: return(0);
1086: }
1089: /*@
1090: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1092: Not collective
1094: Input Parameter:
1095: . fe - The PetscFE
1097: Output Parameter:
1098: . dim - The dimension
1100: Level: intermediate
1102: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1103: @*/
1104: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1105: {
1111: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1112: return(0);
1113: }
1115: /*@C
1116: PetscFEPushforward - Map the reference element function to real space
1118: Input Parameters:
1119: + fe - The PetscFE
1120: . fegeom - The cell geometry
1121: . Nv - The number of function values
1122: - vals - The function values
1124: Output Parameter:
1125: . vals - The transformed function values
1127: Level: advanced
1129: Note: This just forwards the call onto PetscDualSpacePushforward().
1131: .seealso: PetscDualSpacePushforward()
1132: @*/
1133: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1134: {
1138: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1139: return(0);
1140: }
1142: /*@C
1143: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1145: Input Parameters:
1146: + fe - The PetscFE
1147: . fegeom - The cell geometry
1148: . Nv - The number of function gradient values
1149: - vals - The function gradient values
1151: Output Parameter:
1152: . vals - The transformed function gradient values
1154: Level: advanced
1156: Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1158: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1159: @*/
1160: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1161: {
1165: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1166: return(0);
1167: }
1169: /*
1170: Purpose: Compute element vector for chunk of elements
1172: Input:
1173: Sizes:
1174: Ne: number of elements
1175: Nf: number of fields
1176: PetscFE
1177: dim: spatial dimension
1178: Nb: number of basis functions
1179: Nc: number of field components
1180: PetscQuadrature
1181: Nq: number of quadrature points
1183: Geometry:
1184: PetscFEGeom[Ne] possibly *Nq
1185: PetscReal v0s[dim]
1186: PetscReal n[dim]
1187: PetscReal jacobians[dim*dim]
1188: PetscReal jacobianInverses[dim*dim]
1189: PetscReal jacobianDeterminants
1190: FEM:
1191: PetscFE
1192: PetscQuadrature
1193: PetscReal quadPoints[Nq*dim]
1194: PetscReal quadWeights[Nq]
1195: PetscReal basis[Nq*Nb*Nc]
1196: PetscReal basisDer[Nq*Nb*Nc*dim]
1197: PetscScalar coefficients[Ne*Nb*Nc]
1198: PetscScalar elemVec[Ne*Nb*Nc]
1200: Problem:
1201: PetscInt f: the active field
1202: f0, f1
1204: Work Space:
1205: PetscFE
1206: PetscScalar f0[Nq*dim];
1207: PetscScalar f1[Nq*dim*dim];
1208: PetscScalar u[Nc];
1209: PetscScalar gradU[Nc*dim];
1210: PetscReal x[dim];
1211: PetscScalar realSpaceDer[dim];
1213: Purpose: Compute element vector for N_cb batches of elements
1215: Input:
1216: Sizes:
1217: N_cb: Number of serial cell batches
1219: Geometry:
1220: PetscReal v0s[Ne*dim]
1221: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1222: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1223: PetscReal jacobianDeterminants[Ne] possibly *Nq
1224: FEM:
1225: static PetscReal quadPoints[Nq*dim]
1226: static PetscReal quadWeights[Nq]
1227: static PetscReal basis[Nq*Nb*Nc]
1228: static PetscReal basisDer[Nq*Nb*Nc*dim]
1229: PetscScalar coefficients[Ne*Nb*Nc]
1230: PetscScalar elemVec[Ne*Nb*Nc]
1232: ex62.c:
1233: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1234: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1235: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1236: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1238: ex52.c:
1239: 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)
1240: 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)
1242: ex52_integrateElement.cu
1243: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1245: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1246: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1247: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1249: ex52_integrateElementOpenCL.c:
1250: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1251: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1252: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1254: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1255: */
1257: /*@C
1258: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1260: Not collective
1262: Input Parameters:
1263: + fem - The PetscFE object for the field being integrated
1264: . prob - The PetscDS specifying the discretizations and continuum functions
1265: . field - The field being integrated
1266: . Ne - The number of elements in the chunk
1267: . cgeom - The cell geometry for each cell in the chunk
1268: . coefficients - The array of FEM basis coefficients for the elements
1269: . probAux - The PetscDS specifying the auxiliary discretizations
1270: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1272: Output Parameter
1273: . integral - the integral for this field
1275: Level: intermediate
1277: .seealso: PetscFEIntegrateResidual()
1278: @*/
1279: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1280: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1281: {
1282: PetscFE fe;
1287: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1288: if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1289: return(0);
1290: }
1292: /*@C
1293: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1295: Not collective
1297: Input Parameters:
1298: + fem - The PetscFE object for the field being integrated
1299: . prob - The PetscDS specifying the discretizations and continuum functions
1300: . field - The field being integrated
1301: . obj_func - The function to be integrated
1302: . Ne - The number of elements in the chunk
1303: . fgeom - The face geometry for each face in the chunk
1304: . coefficients - The array of FEM basis coefficients for the elements
1305: . probAux - The PetscDS specifying the auxiliary discretizations
1306: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1308: Output Parameter
1309: . integral - the integral for this field
1311: Level: intermediate
1313: .seealso: PetscFEIntegrateResidual()
1314: @*/
1315: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1316: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1317: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1318: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1319: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1320: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1321: {
1322: PetscFE fe;
1327: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1328: if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1329: return(0);
1330: }
1332: /*@C
1333: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1335: Not collective
1337: Input Parameters:
1338: + fem - The PetscFE object for the field being integrated
1339: . prob - The PetscDS specifying the discretizations and continuum functions
1340: . field - The field being integrated
1341: . Ne - The number of elements in the chunk
1342: . cgeom - The cell geometry for each cell in the chunk
1343: . coefficients - The array of FEM basis coefficients for the elements
1344: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1345: . probAux - The PetscDS specifying the auxiliary discretizations
1346: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1347: - t - The time
1349: Output Parameter
1350: . elemVec - the element residual vectors from each element
1352: Note:
1353: $ Loop over batch of elements (e):
1354: $ Loop over quadrature points (q):
1355: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1356: $ Call f_0 and f_1
1357: $ Loop over element vector entries (f,fc --> i):
1358: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1360: Level: intermediate
1362: .seealso: PetscFEIntegrateResidual()
1363: @*/
1364: PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1365: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1366: {
1367: PetscFE fe;
1372: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1373: if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1374: return(0);
1375: }
1377: /*@C
1378: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1380: Not collective
1382: Input Parameters:
1383: + fem - The PetscFE object for the field being integrated
1384: . prob - The PetscDS specifying the discretizations and continuum functions
1385: . field - The field being integrated
1386: . Ne - The number of elements in the chunk
1387: . fgeom - The face geometry for each cell in the chunk
1388: . coefficients - The array of FEM basis coefficients for the elements
1389: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1390: . probAux - The PetscDS specifying the auxiliary discretizations
1391: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1392: - t - The time
1394: Output Parameter
1395: . elemVec - the element residual vectors from each element
1397: Level: intermediate
1399: .seealso: PetscFEIntegrateResidual()
1400: @*/
1401: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1402: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1403: {
1404: PetscFE fe;
1409: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1410: if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1411: return(0);
1412: }
1414: /*@C
1415: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1417: Not collective
1419: Input Parameters:
1420: + fem - The PetscFE object for the field being integrated
1421: . prob - The PetscDS specifying the discretizations and continuum functions
1422: . jtype - The type of matrix pointwise functions that should be used
1423: . fieldI - The test field being integrated
1424: . fieldJ - The basis field being integrated
1425: . Ne - The number of elements in the chunk
1426: . cgeom - The cell geometry for each cell in the chunk
1427: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1428: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1429: . probAux - The PetscDS specifying the auxiliary discretizations
1430: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1431: . t - The time
1432: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1434: Output Parameter
1435: . elemMat - the element matrices for the Jacobian from each element
1437: Note:
1438: $ Loop over batch of elements (e):
1439: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1440: $ Loop over quadrature points (q):
1441: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1442: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1443: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1444: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1445: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1446: Level: intermediate
1448: .seealso: PetscFEIntegrateResidual()
1449: @*/
1450: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1451: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1452: {
1453: PetscFE fe;
1458: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1459: if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1460: return(0);
1461: }
1463: /*@C
1464: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1466: Not collective
1468: Input Parameters:
1469: . prob - The PetscDS specifying the discretizations and continuum functions
1470: . fieldI - The test field being integrated
1471: . fieldJ - The basis field being integrated
1472: . Ne - The number of elements in the chunk
1473: . fgeom - The face geometry for each cell in the chunk
1474: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1475: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1476: . probAux - The PetscDS specifying the auxiliary discretizations
1477: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1478: . t - The time
1479: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1481: Output Parameter
1482: . elemMat - the element matrices for the Jacobian from each element
1484: Note:
1485: $ Loop over batch of elements (e):
1486: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1487: $ Loop over quadrature points (q):
1488: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1489: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1490: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1491: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1492: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1493: Level: intermediate
1495: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1496: @*/
1497: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1498: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1499: {
1500: PetscFE fe;
1505: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1506: if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1507: return(0);
1508: }
1510: /*@
1511: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1513: Input Parameters:
1514: + fe - The finite element space
1515: - height - The height of the Plex point
1517: Output Parameter:
1518: . subfe - The subspace of this FE space
1520: Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1522: Level: advanced
1524: .seealso: PetscFECreateDefault()
1525: @*/
1526: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1527: {
1528: PetscSpace P, subP;
1529: PetscDualSpace Q, subQ;
1530: PetscQuadrature subq;
1531: PetscFEType fetype;
1532: PetscInt dim, Nc;
1533: PetscErrorCode ierr;
1538: if (height == 0) {
1539: *subfe = fe;
1540: return(0);
1541: }
1542: PetscFEGetBasisSpace(fe, &P);
1543: PetscFEGetDualSpace(fe, &Q);
1544: PetscFEGetNumComponents(fe, &Nc);
1545: PetscFEGetFaceQuadrature(fe, &subq);
1546: PetscDualSpaceGetDimension(Q, &dim);
1547: if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);}
1548: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1549: if (height <= dim) {
1550: if (!fe->subspaces[height-1]) {
1551: PetscFE sub;
1552: const char *name;
1554: PetscSpaceGetHeightSubspace(P, height, &subP);
1555: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1556: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1557: PetscObjectGetName((PetscObject) fe, &name);
1558: PetscObjectSetName((PetscObject) sub, name);
1559: PetscFEGetType(fe, &fetype);
1560: PetscFESetType(sub, fetype);
1561: PetscFESetBasisSpace(sub, subP);
1562: PetscFESetDualSpace(sub, subQ);
1563: PetscFESetNumComponents(sub, Nc);
1564: PetscFESetUp(sub);
1565: PetscFESetQuadrature(sub, subq);
1566: fe->subspaces[height-1] = sub;
1567: }
1568: *subfe = fe->subspaces[height-1];
1569: } else {
1570: *subfe = NULL;
1571: }
1572: return(0);
1573: }
1575: /*@
1576: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1577: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1578: sparsity). It is also used to create an interpolation between regularly refined meshes.
1580: Collective on fem
1582: Input Parameter:
1583: . fe - The initial PetscFE
1585: Output Parameter:
1586: . feRef - The refined PetscFE
1588: Level: advanced
1590: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1591: @*/
1592: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1593: {
1594: PetscSpace P, Pref;
1595: PetscDualSpace Q, Qref;
1596: DM K, Kref;
1597: PetscQuadrature q, qref;
1598: const PetscReal *v0, *jac;
1599: PetscInt numComp, numSubelements;
1600: PetscErrorCode ierr;
1603: PetscFEGetBasisSpace(fe, &P);
1604: PetscFEGetDualSpace(fe, &Q);
1605: PetscFEGetQuadrature(fe, &q);
1606: PetscDualSpaceGetDM(Q, &K);
1607: /* Create space */
1608: PetscObjectReference((PetscObject) P);
1609: Pref = P;
1610: /* Create dual space */
1611: PetscDualSpaceDuplicate(Q, &Qref);
1612: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1613: PetscDualSpaceSetDM(Qref, Kref);
1614: DMDestroy(&Kref);
1615: PetscDualSpaceSetUp(Qref);
1616: /* Create element */
1617: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1618: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1619: PetscFESetBasisSpace(*feRef, Pref);
1620: PetscFESetDualSpace(*feRef, Qref);
1621: PetscFEGetNumComponents(fe, &numComp);
1622: PetscFESetNumComponents(*feRef, numComp);
1623: PetscFESetUp(*feRef);
1624: PetscSpaceDestroy(&Pref);
1625: PetscDualSpaceDestroy(&Qref);
1626: /* Create quadrature */
1627: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1628: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1629: PetscFESetQuadrature(*feRef, qref);
1630: PetscQuadratureDestroy(&qref);
1631: return(0);
1632: }
1634: /*@C
1635: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1637: Collective
1639: Input Parameters:
1640: + comm - The MPI comm
1641: . dim - The spatial dimension
1642: . Nc - The number of components
1643: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1644: . prefix - The options prefix, or NULL
1645: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1647: Output Parameter:
1648: . fem - The PetscFE object
1650: Level: beginner
1652: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1653: @*/
1654: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1655: {
1656: PetscQuadrature q, fq;
1657: DM K;
1658: PetscSpace P;
1659: PetscDualSpace Q;
1660: PetscInt order, quadPointsPerEdge;
1661: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1662: PetscErrorCode ierr;
1665: /* Create space */
1666: PetscSpaceCreate(comm, &P);
1667: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1668: PetscSpacePolynomialSetTensor(P, tensor);
1669: PetscSpaceSetNumComponents(P, Nc);
1670: PetscSpaceSetNumVariables(P, dim);
1671: PetscSpaceSetFromOptions(P);
1672: PetscSpaceSetUp(P);
1673: PetscSpaceGetDegree(P, &order, NULL);
1674: PetscSpacePolynomialGetTensor(P, &tensor);
1675: /* Create dual space */
1676: PetscDualSpaceCreate(comm, &Q);
1677: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1678: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1679: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1680: PetscDualSpaceSetDM(Q, K);
1681: DMDestroy(&K);
1682: PetscDualSpaceSetNumComponents(Q, Nc);
1683: PetscDualSpaceSetOrder(Q, order);
1684: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1685: PetscDualSpaceSetFromOptions(Q);
1686: PetscDualSpaceSetUp(Q);
1687: /* Create element */
1688: PetscFECreate(comm, fem);
1689: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1690: PetscFESetBasisSpace(*fem, P);
1691: PetscFESetDualSpace(*fem, Q);
1692: PetscFESetNumComponents(*fem, Nc);
1693: PetscFESetFromOptions(*fem);
1694: PetscFESetUp(*fem);
1695: PetscSpaceDestroy(&P);
1696: PetscDualSpaceDestroy(&Q);
1697: /* Create quadrature (with specified order if given) */
1698: qorder = qorder >= 0 ? qorder : order;
1699: PetscObjectOptionsBegin((PetscObject)*fem);
1700: PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1701: PetscOptionsEnd();
1702: quadPointsPerEdge = PetscMax(qorder + 1,1);
1703: if (isSimplex) {
1704: PetscDTGaussJacobiQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1705: PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1706: } else {
1707: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1708: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1709: }
1710: PetscFESetQuadrature(*fem, q);
1711: PetscFESetFaceQuadrature(*fem, fq);
1712: PetscQuadratureDestroy(&q);
1713: PetscQuadratureDestroy(&fq);
1714: return(0);
1715: }
1717: /*@C
1718: PetscFESetName - Names the FE and its subobjects
1720: Not collective
1722: Input Parameters:
1723: + fe - The PetscFE
1724: - name - The name
1726: Level: intermediate
1728: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1729: @*/
1730: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1731: {
1732: PetscSpace P;
1733: PetscDualSpace Q;
1737: PetscFEGetBasisSpace(fe, &P);
1738: PetscFEGetDualSpace(fe, &Q);
1739: PetscObjectSetName((PetscObject) fe, name);
1740: PetscObjectSetName((PetscObject) P, name);
1741: PetscObjectSetName((PetscObject) Q, name);
1742: return(0);
1743: }
1745: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1746: {
1747: PetscInt dOffset = 0, fOffset = 0, f;
1750: for (f = 0; f < Nf; ++f) {
1751: PetscFE fe;
1752: const PetscInt Nbf = Nb[f], Ncf = Nc[f];
1753: const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
1754: const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
1755: PetscInt b, c, d;
1757: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
1758: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
1759: for (d = 0; d < dim*Ncf; ++d) u_x[fOffset*dim+d] = 0.0;
1760: for (b = 0; b < Nbf; ++b) {
1761: for (c = 0; c < Ncf; ++c) {
1762: const PetscInt cidx = b*Ncf+c;
1764: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
1765: for (d = 0; d < dim; ++d) u_x[(fOffset+c)*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
1766: }
1767: }
1768: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
1769: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*dim]);
1770: if (u_t) {
1771: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
1772: for (b = 0; b < Nbf; ++b) {
1773: for (c = 0; c < Ncf; ++c) {
1774: const PetscInt cidx = b*Ncf+c;
1776: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
1777: }
1778: }
1779: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
1780: }
1781: fOffset += Ncf;
1782: dOffset += Nbf;
1783: }
1784: return 0;
1785: }
1787: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
1788: {
1789: PetscFE fe;
1790: PetscReal *faceBasis;
1791: PetscInt Nb, Nc, b, c;
1794: if (!prob) return 0;
1795: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1796: PetscFEGetDimension(fe, &Nb);
1797: PetscFEGetNumComponents(fe, &Nc);
1798: PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
1799: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
1800: for (b = 0; b < Nb; ++b) {
1801: for (c = 0; c < Nc; ++c) {
1802: const PetscInt cidx = b*Nc+c;
1804: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
1805: }
1806: }
1807: return 0;
1808: }
1810: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
1811: {
1812: PetscInt q, b, c, d;
1815: for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
1816: for (q = 0; q < Nq; ++q) {
1817: for (b = 0; b < Nb; ++b) {
1818: for (c = 0; c < Nc; ++c) {
1819: const PetscInt bcidx = b*Nc+c;
1821: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
1822: for (d = 0; d < dim; ++d) tmpBasisDer[bcidx*dim+d] = basisDer[q*Nb*Nc*dim+bcidx*dim+d];
1823: }
1824: }
1825: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
1826: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
1827: for (b = 0; b < Nb; ++b) {
1828: for (c = 0; c < Nc; ++c) {
1829: const PetscInt bcidx = b*Nc+c;
1830: const PetscInt qcidx = q*Nc+c;
1832: elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
1833: for (d = 0; d < dim; ++d) elemVec[b] += tmpBasisDer[bcidx*dim+d]*f1[qcidx*dim+d];
1834: }
1835: }
1836: }
1837: return(0);
1838: }
1840: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt dim, PetscInt NbI, PetscInt NcI, const PetscReal basisI[], const PetscReal basisDerI[], PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscInt NbJ, PetscInt NcJ, const PetscReal basisJ[], const PetscReal basisDerJ[], 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[])
1841: {
1842: PetscInt f, fc, g, gc, df, dg;
1845: for (f = 0; f < NbI; ++f) {
1846: for (fc = 0; fc < NcI; ++fc) {
1847: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
1849: tmpBasisI[fidx] = basisI[fidx];
1850: for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dim+df] = basisDerI[fidx*dim+df];
1851: }
1852: }
1853: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
1854: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
1855: for (g = 0; g < NbJ; ++g) {
1856: for (gc = 0; gc < NcJ; ++gc) {
1857: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
1859: tmpBasisJ[gidx] = basisJ[gidx];
1860: for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dim+dg] = basisDerJ[gidx*dim+dg];
1861: }
1862: }
1863: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
1864: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
1865: for (f = 0; f < NbI; ++f) {
1866: for (fc = 0; fc < NcI; ++fc) {
1867: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
1868: const PetscInt i = offsetI+f; /* Element matrix row */
1869: for (g = 0; g < NbJ; ++g) {
1870: for (gc = 0; gc < NcJ; ++gc) {
1871: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
1872: const PetscInt j = offsetJ+g; /* Element matrix column */
1873: const PetscInt fOff = eOffset+i*totDim+j;
1875: elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
1876: for (df = 0; df < dim; ++df) {
1877: elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dim+df]*tmpBasisDerJ[gidx*dim+df];
1878: elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g2[(fc*NcJ+gc)*dim+df]*tmpBasisJ[gidx];
1879: for (dg = 0; dg < dim; ++dg) {
1880: elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g3[((fc*NcJ+gc)*dim+df)*dim+dg]*tmpBasisDerJ[gidx*dim+dg];
1881: }
1882: }
1883: }
1884: }
1885: }
1886: }
1887: return(0);
1888: }
1890: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
1891: {
1892: PetscDualSpace dsp;
1893: DM dm;
1894: PetscQuadrature quadDef;
1895: PetscInt dim, cdim, Nq;
1896: PetscErrorCode ierr;
1899: PetscFEGetDualSpace(fe, &dsp);
1900: PetscDualSpaceGetDM(dsp, &dm);
1901: DMGetDimension(dm, &dim);
1902: DMGetCoordinateDim(dm, &cdim);
1903: PetscFEGetQuadrature(fe, &quadDef);
1904: quad = quad ? quad : quadDef;
1905: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
1906: PetscMalloc1(Nq*cdim, &cgeom->v);
1907: PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
1908: PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
1909: PetscMalloc1(Nq, &cgeom->detJ);
1910: cgeom->dim = dim;
1911: cgeom->dimEmbed = cdim;
1912: cgeom->numCells = 1;
1913: cgeom->numPoints = Nq;
1914: DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
1915: return(0);
1916: }
1918: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
1919: {
1923: PetscFree(cgeom->v);
1924: PetscFree(cgeom->J);
1925: PetscFree(cgeom->invJ);
1926: PetscFree(cgeom->detJ);
1927: return(0);
1928: }