Actual source code: fe.c
petsc-3.14.6 2021-03-30
1: /* Basis Jet Tabulation
3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6: as a prime basis.
8: \psi_i = \sum_k \alpha_{ki} \phi_k
10: Our nodal basis is defined in terms of the dual basis $n_j$
12: n_j \cdot \psi_i = \delta_{ji}
14: and we may act on the first equation to obtain
16: n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17: \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18: I = V \alpha
20: so the coefficients of the nodal basis in the prime basis are
22: \alpha = V^{-1}
24: We will define the dual basis vectors $n_j$ using a quadrature rule.
26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28: be implemented exactly as in FIAT using functionals $L_j$.
30: I will have to count the degrees correctly for the Legendre product when we are on simplices.
32: We will have three objects:
33: - Space, P: this just need point evaluation I think
34: - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35: - FEM: This keeps {P, P', Q}
36: */
37: #include <petsc/private/petscfeimpl.h>
38: #include <petscdmplex.h>
40: PetscBool FEcite = PETSC_FALSE;
41: const char FECitation[] = "@article{kirby2004,\n"
42: " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
43: " journal = {ACM Transactions on Mathematical Software},\n"
44: " author = {Robert C. Kirby},\n"
45: " volume = {30},\n"
46: " number = {4},\n"
47: " pages = {502--516},\n"
48: " doi = {10.1145/1039813.1039820},\n"
49: " year = {2004}\n}\n";
51: PetscClassId PETSCFE_CLASSID = 0;
53: PetscLogEvent PETSCFE_SetUp;
55: PetscFunctionList PetscFEList = NULL;
56: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
58: /*@C
59: PetscFERegister - Adds a new PetscFE implementation
61: Not Collective
63: Input Parameters:
64: + name - The name of a new user-defined creation routine
65: - create_func - The creation routine itself
67: Notes:
68: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
70: Sample usage:
71: .vb
72: PetscFERegister("my_fe", MyPetscFECreate);
73: .ve
75: Then, your PetscFE type can be chosen with the procedural interface via
76: .vb
77: PetscFECreate(MPI_Comm, PetscFE *);
78: PetscFESetType(PetscFE, "my_fe");
79: .ve
80: or at runtime via the option
81: .vb
82: -petscfe_type my_fe
83: .ve
85: Level: advanced
87: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
89: @*/
90: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
91: {
95: PetscFunctionListAdd(&PetscFEList, sname, function);
96: return(0);
97: }
99: /*@C
100: PetscFESetType - Builds a particular PetscFE
102: Collective on fem
104: Input Parameters:
105: + fem - The PetscFE object
106: - name - The kind of FEM space
108: Options Database Key:
109: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
111: Level: intermediate
113: .seealso: PetscFEGetType(), PetscFECreate()
114: @*/
115: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
116: {
117: PetscErrorCode (*r)(PetscFE);
118: PetscBool match;
123: PetscObjectTypeCompare((PetscObject) fem, name, &match);
124: if (match) return(0);
126: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
127: PetscFunctionListFind(PetscFEList, name, &r);
128: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
130: if (fem->ops->destroy) {
131: (*fem->ops->destroy)(fem);
132: fem->ops->destroy = NULL;
133: }
134: (*r)(fem);
135: PetscObjectChangeTypeName((PetscObject) fem, name);
136: return(0);
137: }
139: /*@C
140: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
142: Not Collective
144: Input Parameter:
145: . fem - The PetscFE
147: Output Parameter:
148: . name - The PetscFE type name
150: Level: intermediate
152: .seealso: PetscFESetType(), PetscFECreate()
153: @*/
154: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
155: {
161: if (!PetscFERegisterAllCalled) {
162: PetscFERegisterAll();
163: }
164: *name = ((PetscObject) fem)->type_name;
165: return(0);
166: }
168: /*@C
169: PetscFEViewFromOptions - View from Options
171: Collective on PetscFE
173: Input Parameters:
174: + A - the PetscFE object
175: . obj - Optional object
176: - name - command line option
178: Level: intermediate
179: .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
180: @*/
181: PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
182: {
187: PetscObjectViewFromOptions((PetscObject)A,obj,name);
188: return(0);
189: }
191: /*@C
192: PetscFEView - Views a PetscFE
194: Collective on fem
196: Input Parameter:
197: + fem - the PetscFE object to view
198: - viewer - the viewer
200: Level: beginner
202: .seealso PetscFEDestroy()
203: @*/
204: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
205: {
206: PetscBool iascii;
212: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
213: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
214: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
215: if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
216: return(0);
217: }
219: /*@
220: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
222: Collective on fem
224: Input Parameter:
225: . fem - the PetscFE object to set options for
227: Options Database:
228: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
229: - -petscfe_num_batches - the number of cell batches to integrate serially
231: Level: intermediate
233: .seealso PetscFEView()
234: @*/
235: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
236: {
237: const char *defaultType;
238: char name[256];
239: PetscBool flg;
244: if (!((PetscObject) fem)->type_name) {
245: defaultType = PETSCFEBASIC;
246: } else {
247: defaultType = ((PetscObject) fem)->type_name;
248: }
249: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
251: PetscObjectOptionsBegin((PetscObject) fem);
252: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
253: if (flg) {
254: PetscFESetType(fem, name);
255: } else if (!((PetscObject) fem)->type_name) {
256: PetscFESetType(fem, defaultType);
257: }
258: PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
259: PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
260: if (fem->ops->setfromoptions) {
261: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
262: }
263: /* process any options handlers added with PetscObjectAddOptionsHandler() */
264: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
265: PetscOptionsEnd();
266: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
267: return(0);
268: }
270: /*@C
271: PetscFESetUp - Construct data structures for the PetscFE
273: Collective on fem
275: Input Parameter:
276: . fem - the PetscFE object to setup
278: Level: intermediate
280: .seealso PetscFEView(), PetscFEDestroy()
281: @*/
282: PetscErrorCode PetscFESetUp(PetscFE fem)
283: {
288: if (fem->setupcalled) return(0);
289: PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
290: fem->setupcalled = PETSC_TRUE;
291: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
292: PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
293: return(0);
294: }
296: /*@
297: PetscFEDestroy - Destroys a PetscFE object
299: Collective on fem
301: Input Parameter:
302: . fem - the PetscFE object to destroy
304: Level: beginner
306: .seealso PetscFEView()
307: @*/
308: PetscErrorCode PetscFEDestroy(PetscFE *fem)
309: {
313: if (!*fem) return(0);
316: if (--((PetscObject)(*fem))->refct > 0) {*fem = NULL; return(0);}
317: ((PetscObject) (*fem))->refct = 0;
319: if ((*fem)->subspaces) {
320: PetscInt dim, d;
322: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
323: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
324: }
325: PetscFree((*fem)->subspaces);
326: PetscFree((*fem)->invV);
327: PetscTabulationDestroy(&(*fem)->T);
328: PetscTabulationDestroy(&(*fem)->Tf);
329: PetscTabulationDestroy(&(*fem)->Tc);
330: PetscSpaceDestroy(&(*fem)->basisSpace);
331: PetscDualSpaceDestroy(&(*fem)->dualSpace);
332: PetscQuadratureDestroy(&(*fem)->quadrature);
333: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
335: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
336: PetscHeaderDestroy(fem);
337: return(0);
338: }
340: /*@
341: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
343: Collective
345: Input Parameter:
346: . comm - The communicator for the PetscFE object
348: Output Parameter:
349: . fem - The PetscFE object
351: Level: beginner
353: .seealso: PetscFESetType(), PETSCFEGALERKIN
354: @*/
355: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
356: {
357: PetscFE f;
362: PetscCitationsRegister(FECitation,&FEcite);
363: *fem = NULL;
364: PetscFEInitializePackage();
366: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
368: f->basisSpace = NULL;
369: f->dualSpace = NULL;
370: f->numComponents = 1;
371: f->subspaces = NULL;
372: f->invV = NULL;
373: f->T = NULL;
374: f->Tf = NULL;
375: f->Tc = NULL;
376: PetscArrayzero(&f->quadrature, 1);
377: PetscArrayzero(&f->faceQuadrature, 1);
378: f->blockSize = 0;
379: f->numBlocks = 1;
380: f->batchSize = 0;
381: f->numBatches = 1;
383: *fem = f;
384: return(0);
385: }
387: /*@
388: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
390: Not collective
392: Input Parameter:
393: . fem - The PetscFE object
395: Output Parameter:
396: . dim - The spatial dimension
398: Level: intermediate
400: .seealso: PetscFECreate()
401: @*/
402: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
403: {
404: DM dm;
410: PetscDualSpaceGetDM(fem->dualSpace, &dm);
411: DMGetDimension(dm, dim);
412: return(0);
413: }
415: /*@
416: PetscFESetNumComponents - Sets the number of components in the element
418: Not collective
420: Input Parameters:
421: + fem - The PetscFE object
422: - comp - The number of field components
424: Level: intermediate
426: .seealso: PetscFECreate()
427: @*/
428: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
429: {
432: fem->numComponents = comp;
433: return(0);
434: }
436: /*@
437: PetscFEGetNumComponents - Returns the number of components in the element
439: Not collective
441: Input Parameter:
442: . fem - The PetscFE object
444: Output Parameter:
445: . comp - The number of field components
447: Level: intermediate
449: .seealso: PetscFECreate()
450: @*/
451: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
452: {
456: *comp = fem->numComponents;
457: return(0);
458: }
460: /*@
461: PetscFESetTileSizes - Sets the tile sizes for evaluation
463: Not collective
465: Input Parameters:
466: + fem - The PetscFE object
467: . blockSize - The number of elements in a block
468: . numBlocks - The number of blocks in a batch
469: . batchSize - The number of elements in a batch
470: - numBatches - The number of batches in a chunk
472: Level: intermediate
474: .seealso: PetscFECreate()
475: @*/
476: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
477: {
480: fem->blockSize = blockSize;
481: fem->numBlocks = numBlocks;
482: fem->batchSize = batchSize;
483: fem->numBatches = numBatches;
484: return(0);
485: }
487: /*@
488: PetscFEGetTileSizes - Returns the tile sizes for evaluation
490: Not collective
492: Input Parameter:
493: . fem - The PetscFE object
495: Output Parameters:
496: + blockSize - The number of elements in a block
497: . numBlocks - The number of blocks in a batch
498: . batchSize - The number of elements in a batch
499: - numBatches - The number of batches in a chunk
501: Level: intermediate
503: .seealso: PetscFECreate()
504: @*/
505: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
506: {
513: if (blockSize) *blockSize = fem->blockSize;
514: if (numBlocks) *numBlocks = fem->numBlocks;
515: if (batchSize) *batchSize = fem->batchSize;
516: if (numBatches) *numBatches = fem->numBatches;
517: return(0);
518: }
520: /*@
521: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
523: Not collective
525: Input Parameter:
526: . fem - The PetscFE object
528: Output Parameter:
529: . sp - The PetscSpace object
531: Level: intermediate
533: .seealso: PetscFECreate()
534: @*/
535: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
536: {
540: *sp = fem->basisSpace;
541: return(0);
542: }
544: /*@
545: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
547: Not collective
549: Input Parameters:
550: + fem - The PetscFE object
551: - sp - The PetscSpace object
553: Level: intermediate
555: .seealso: PetscFECreate()
556: @*/
557: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
558: {
564: PetscSpaceDestroy(&fem->basisSpace);
565: fem->basisSpace = sp;
566: PetscObjectReference((PetscObject) fem->basisSpace);
567: return(0);
568: }
570: /*@
571: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
573: Not collective
575: Input Parameter:
576: . fem - The PetscFE object
578: Output Parameter:
579: . sp - The PetscDualSpace object
581: Level: intermediate
583: .seealso: PetscFECreate()
584: @*/
585: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
586: {
590: *sp = fem->dualSpace;
591: return(0);
592: }
594: /*@
595: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
597: Not collective
599: Input Parameters:
600: + fem - The PetscFE object
601: - sp - The PetscDualSpace object
603: Level: intermediate
605: .seealso: PetscFECreate()
606: @*/
607: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
608: {
614: PetscDualSpaceDestroy(&fem->dualSpace);
615: fem->dualSpace = sp;
616: PetscObjectReference((PetscObject) fem->dualSpace);
617: return(0);
618: }
620: /*@
621: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
623: Not collective
625: Input Parameter:
626: . fem - The PetscFE object
628: Output Parameter:
629: . q - The PetscQuadrature object
631: Level: intermediate
633: .seealso: PetscFECreate()
634: @*/
635: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
636: {
640: *q = fem->quadrature;
641: return(0);
642: }
644: /*@
645: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
647: Not collective
649: Input Parameters:
650: + fem - The PetscFE object
651: - q - The PetscQuadrature object
653: Level: intermediate
655: .seealso: PetscFECreate()
656: @*/
657: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
658: {
659: PetscInt Nc, qNc;
664: if (q == fem->quadrature) return(0);
665: PetscFEGetNumComponents(fem, &Nc);
666: PetscQuadratureGetNumComponents(q, &qNc);
667: if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
668: PetscTabulationDestroy(&fem->T);
669: PetscTabulationDestroy(&fem->Tc);
670: PetscObjectReference((PetscObject) q);
671: PetscQuadratureDestroy(&fem->quadrature);
672: fem->quadrature = q;
673: return(0);
674: }
676: /*@
677: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
679: Not collective
681: Input Parameter:
682: . fem - The PetscFE object
684: Output Parameter:
685: . q - The PetscQuadrature object
687: Level: intermediate
689: .seealso: PetscFECreate()
690: @*/
691: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
692: {
696: *q = fem->faceQuadrature;
697: return(0);
698: }
700: /*@
701: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
703: Not collective
705: Input Parameters:
706: + fem - The PetscFE object
707: - q - The PetscQuadrature object
709: Level: intermediate
711: .seealso: PetscFECreate()
712: @*/
713: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
714: {
715: PetscInt Nc, qNc;
720: PetscFEGetNumComponents(fem, &Nc);
721: PetscQuadratureGetNumComponents(q, &qNc);
722: if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
723: PetscTabulationDestroy(&fem->Tf);
724: PetscQuadratureDestroy(&fem->faceQuadrature);
725: fem->faceQuadrature = q;
726: PetscObjectReference((PetscObject) q);
727: return(0);
728: }
730: /*@
731: PetscFECopyQuadrature - Copy both volumetric and surface quadrature
733: Not collective
735: Input Parameters:
736: + sfe - The PetscFE source for the quadratures
737: - tfe - The PetscFE target for the quadratures
739: Level: intermediate
741: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
742: @*/
743: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
744: {
745: PetscQuadrature q;
746: PetscErrorCode ierr;
751: PetscFEGetQuadrature(sfe, &q);
752: PetscFESetQuadrature(tfe, q);
753: PetscFEGetFaceQuadrature(sfe, &q);
754: PetscFESetFaceQuadrature(tfe, q);
755: return(0);
756: }
758: /*@C
759: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
761: Not collective
763: Input Parameter:
764: . fem - The PetscFE object
766: Output Parameter:
767: . numDof - Array with the number of dofs per dimension
769: Level: intermediate
771: .seealso: PetscFECreate()
772: @*/
773: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
774: {
780: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
781: return(0);
782: }
784: /*@C
785: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
787: Not collective
789: Input Parameter:
790: . fem - The PetscFE object
792: Output Parameter:
793: . T - The basis function values and derivatives at quadrature points
795: Note:
796: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
797: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
798: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
800: Level: intermediate
802: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
803: @*/
804: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T)
805: {
806: PetscInt npoints;
807: const PetscReal *points;
808: PetscErrorCode ierr;
813: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
814: if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);}
815: *T = fem->T;
816: return(0);
817: }
819: /*@C
820: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
822: Not collective
824: Input Parameter:
825: . fem - The PetscFE object
827: Output Parameters:
828: . Tf - The basis function values and derviatives at face quadrature points
830: Note:
831: $ T->T[0] = Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
832: $ T->T[1] = Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d
833: $ T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e
835: Level: intermediate
837: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
838: @*/
839: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf)
840: {
841: PetscErrorCode ierr;
846: if (!fem->Tf) {
847: const PetscReal xi0[3] = {-1., -1., -1.};
848: PetscReal v0[3], J[9], detJ;
849: PetscQuadrature fq;
850: PetscDualSpace sp;
851: DM dm;
852: const PetscInt *faces;
853: PetscInt dim, numFaces, f, npoints, q;
854: const PetscReal *points;
855: PetscReal *facePoints;
857: PetscFEGetDualSpace(fem, &sp);
858: PetscDualSpaceGetDM(sp, &dm);
859: DMGetDimension(dm, &dim);
860: DMPlexGetConeSize(dm, 0, &numFaces);
861: DMPlexGetCone(dm, 0, &faces);
862: PetscFEGetFaceQuadrature(fem, &fq);
863: if (fq) {
864: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
865: PetscMalloc1(numFaces*npoints*dim, &facePoints);
866: for (f = 0; f < numFaces; ++f) {
867: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
868: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
869: }
870: PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);
871: PetscFree(facePoints);
872: }
873: }
874: *Tf = fem->Tf;
875: return(0);
876: }
878: /*@C
879: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
881: Not collective
883: Input Parameter:
884: . fem - The PetscFE object
886: Output Parameters:
887: . Tc - The basis function values at face centroid points
889: Note:
890: $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
892: Level: intermediate
894: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
895: @*/
896: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
897: {
898: PetscErrorCode ierr;
903: if (!fem->Tc) {
904: PetscDualSpace sp;
905: DM dm;
906: const PetscInt *cone;
907: PetscReal *centroids;
908: PetscInt dim, numFaces, f;
910: PetscFEGetDualSpace(fem, &sp);
911: PetscDualSpaceGetDM(sp, &dm);
912: DMGetDimension(dm, &dim);
913: DMPlexGetConeSize(dm, 0, &numFaces);
914: DMPlexGetCone(dm, 0, &cone);
915: PetscMalloc1(numFaces*dim, ¢roids);
916: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
917: PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
918: PetscFree(centroids);
919: }
920: *Tc = fem->Tc;
921: return(0);
922: }
924: /*@C
925: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
927: Not collective
929: Input Parameters:
930: + fem - The PetscFE object
931: . nrepl - The number of replicas
932: . npoints - The number of tabulation points in a replica
933: . points - The tabulation point coordinates
934: - K - The number of derivatives calculated
936: Output Parameter:
937: . T - The basis function values and derivatives at tabulation points
939: Note:
940: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
941: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
942: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
944: Level: intermediate
946: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
947: @*/
948: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
949: {
950: DM dm;
951: PetscDualSpace Q;
952: PetscInt Nb; /* Dimension of FE space P */
953: PetscInt Nc; /* Field components */
954: PetscInt cdim; /* Reference coordinate dimension */
955: PetscInt k;
956: PetscErrorCode ierr;
959: if (!npoints || !fem->dualSpace || K < 0) {
960: *T = NULL;
961: return(0);
962: }
966: PetscFEGetDualSpace(fem, &Q);
967: PetscDualSpaceGetDM(Q, &dm);
968: DMGetDimension(dm, &cdim);
969: PetscDualSpaceGetDimension(Q, &Nb);
970: PetscFEGetNumComponents(fem, &Nc);
971: PetscMalloc1(1, T);
972: (*T)->K = !cdim ? 0 : K;
973: (*T)->Nr = nrepl;
974: (*T)->Np = npoints;
975: (*T)->Nb = Nb;
976: (*T)->Nc = Nc;
977: (*T)->cdim = cdim;
978: PetscMalloc1((*T)->K+1, &(*T)->T);
979: for (k = 0; k <= (*T)->K; ++k) {
980: PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
981: }
982: (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
983: return(0);
984: }
986: /*@C
987: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
989: Not collective
991: Input Parameters:
992: + fem - The PetscFE object
993: . npoints - The number of tabulation points
994: . points - The tabulation point coordinates
995: . K - The number of derivatives calculated
996: - T - An existing tabulation object with enough allocated space
998: Output Parameter:
999: . T - The basis function values and derivatives at tabulation points
1001: Note:
1002: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1003: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1004: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1006: Level: intermediate
1008: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1009: @*/
1010: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1011: {
1015: if (!npoints || !fem->dualSpace || K < 0) return(0);
1019: if (PetscDefined(USE_DEBUG)) {
1020: DM dm;
1021: PetscDualSpace Q;
1022: PetscInt Nb; /* Dimension of FE space P */
1023: PetscInt Nc; /* Field components */
1024: PetscInt cdim; /* Reference coordinate dimension */
1026: PetscFEGetDualSpace(fem, &Q);
1027: PetscDualSpaceGetDM(Q, &dm);
1028: DMGetDimension(dm, &cdim);
1029: PetscDualSpaceGetDimension(Q, &Nb);
1030: PetscFEGetNumComponents(fem, &Nc);
1031: if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K);
1032: if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1033: if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1034: if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1035: }
1036: T->Nr = 1;
1037: T->Np = npoints;
1038: (*fem->ops->createtabulation)(fem, npoints, points, K, T);
1039: return(0);
1040: }
1042: /*@C
1043: PetscTabulationDestroy - Frees memory from the associated tabulation.
1045: Not collective
1047: Input Parameter:
1048: . T - The tabulation
1050: Level: intermediate
1052: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1053: @*/
1054: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1055: {
1056: PetscInt k;
1061: if (!T || !(*T)) return(0);
1062: for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);}
1063: PetscFree((*T)->T);
1064: PetscFree(*T);
1065: *T = NULL;
1066: return(0);
1067: }
1069: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1070: {
1071: PetscSpace bsp, bsubsp;
1072: PetscDualSpace dsp, dsubsp;
1073: PetscInt dim, depth, numComp, i, j, coneSize, order;
1074: PetscFEType type;
1075: DM dm;
1076: DMLabel label;
1077: PetscReal *xi, *v, *J, detJ;
1078: const char *name;
1079: PetscQuadrature origin, fullQuad, subQuad;
1085: PetscFEGetBasisSpace(fe,&bsp);
1086: PetscFEGetDualSpace(fe,&dsp);
1087: PetscDualSpaceGetDM(dsp,&dm);
1088: DMGetDimension(dm,&dim);
1089: DMPlexGetDepthLabel(dm,&label);
1090: DMLabelGetValue(label,refPoint,&depth);
1091: PetscCalloc1(depth,&xi);
1092: PetscMalloc1(dim,&v);
1093: PetscMalloc1(dim*dim,&J);
1094: for (i = 0; i < depth; i++) xi[i] = 0.;
1095: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1096: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1097: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1098: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1099: for (i = 1; i < dim; i++) {
1100: for (j = 0; j < depth; j++) {
1101: J[i * depth + j] = J[i * dim + j];
1102: }
1103: }
1104: PetscQuadratureDestroy(&origin);
1105: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1106: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1107: PetscSpaceSetUp(bsubsp);
1108: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1109: PetscFEGetType(fe,&type);
1110: PetscFESetType(*trFE,type);
1111: PetscFEGetNumComponents(fe,&numComp);
1112: PetscFESetNumComponents(*trFE,numComp);
1113: PetscFESetBasisSpace(*trFE,bsubsp);
1114: PetscFESetDualSpace(*trFE,dsubsp);
1115: PetscObjectGetName((PetscObject) fe, &name);
1116: if (name) {PetscFESetName(*trFE, name);}
1117: PetscFEGetQuadrature(fe,&fullQuad);
1118: PetscQuadratureGetOrder(fullQuad,&order);
1119: DMPlexGetConeSize(dm,refPoint,&coneSize);
1120: if (coneSize == 2 * depth) {
1121: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1122: } else {
1123: PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1124: }
1125: PetscFESetQuadrature(*trFE,subQuad);
1126: PetscFESetUp(*trFE);
1127: PetscQuadratureDestroy(&subQuad);
1128: PetscSpaceDestroy(&bsubsp);
1129: return(0);
1130: }
1132: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1133: {
1134: PetscInt hStart, hEnd;
1135: PetscDualSpace dsp;
1136: DM dm;
1142: *trFE = NULL;
1143: PetscFEGetDualSpace(fe,&dsp);
1144: PetscDualSpaceGetDM(dsp,&dm);
1145: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1146: if (hEnd <= hStart) return(0);
1147: PetscFECreatePointTrace(fe,hStart,trFE);
1148: return(0);
1149: }
1152: /*@
1153: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1155: Not collective
1157: Input Parameter:
1158: . fe - The PetscFE
1160: Output Parameter:
1161: . dim - The dimension
1163: Level: intermediate
1165: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1166: @*/
1167: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1168: {
1174: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1175: return(0);
1176: }
1178: /*@C
1179: PetscFEPushforward - Map the reference element function to real space
1181: Input Parameters:
1182: + fe - The PetscFE
1183: . fegeom - The cell geometry
1184: . Nv - The number of function values
1185: - vals - The function values
1187: Output Parameter:
1188: . vals - The transformed function values
1190: Level: advanced
1192: Note: This just forwards the call onto PetscDualSpacePushforward().
1194: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1196: .seealso: PetscDualSpacePushforward()
1197: @*/
1198: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1199: {
1203: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1204: return(0);
1205: }
1207: /*@C
1208: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1210: Input Parameters:
1211: + fe - The PetscFE
1212: . fegeom - The cell geometry
1213: . Nv - The number of function gradient values
1214: - vals - The function gradient values
1216: Output Parameter:
1217: . vals - The transformed function gradient values
1219: Level: advanced
1221: Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1223: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1225: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1226: @*/
1227: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1228: {
1232: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1233: return(0);
1234: }
1236: /*
1237: Purpose: Compute element vector for chunk of elements
1239: Input:
1240: Sizes:
1241: Ne: number of elements
1242: Nf: number of fields
1243: PetscFE
1244: dim: spatial dimension
1245: Nb: number of basis functions
1246: Nc: number of field components
1247: PetscQuadrature
1248: Nq: number of quadrature points
1250: Geometry:
1251: PetscFEGeom[Ne] possibly *Nq
1252: PetscReal v0s[dim]
1253: PetscReal n[dim]
1254: PetscReal jacobians[dim*dim]
1255: PetscReal jacobianInverses[dim*dim]
1256: PetscReal jacobianDeterminants
1257: FEM:
1258: PetscFE
1259: PetscQuadrature
1260: PetscReal quadPoints[Nq*dim]
1261: PetscReal quadWeights[Nq]
1262: PetscReal basis[Nq*Nb*Nc]
1263: PetscReal basisDer[Nq*Nb*Nc*dim]
1264: PetscScalar coefficients[Ne*Nb*Nc]
1265: PetscScalar elemVec[Ne*Nb*Nc]
1267: Problem:
1268: PetscInt f: the active field
1269: f0, f1
1271: Work Space:
1272: PetscFE
1273: PetscScalar f0[Nq*dim];
1274: PetscScalar f1[Nq*dim*dim];
1275: PetscScalar u[Nc];
1276: PetscScalar gradU[Nc*dim];
1277: PetscReal x[dim];
1278: PetscScalar realSpaceDer[dim];
1280: Purpose: Compute element vector for N_cb batches of elements
1282: Input:
1283: Sizes:
1284: N_cb: Number of serial cell batches
1286: Geometry:
1287: PetscReal v0s[Ne*dim]
1288: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1289: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1290: PetscReal jacobianDeterminants[Ne] possibly *Nq
1291: FEM:
1292: static PetscReal quadPoints[Nq*dim]
1293: static PetscReal quadWeights[Nq]
1294: static PetscReal basis[Nq*Nb*Nc]
1295: static PetscReal basisDer[Nq*Nb*Nc*dim]
1296: PetscScalar coefficients[Ne*Nb*Nc]
1297: PetscScalar elemVec[Ne*Nb*Nc]
1299: ex62.c:
1300: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1301: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1302: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1303: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1305: ex52.c:
1306: PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
1307: PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
1309: ex52_integrateElement.cu
1310: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1312: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1313: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1314: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1316: ex52_integrateElementOpenCL.c:
1317: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1318: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1319: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1321: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1322: */
1324: /*@C
1325: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1327: Not collective
1329: Input Parameters:
1330: + prob - The PetscDS specifying the discretizations and continuum functions
1331: . field - The field being integrated
1332: . Ne - The number of elements in the chunk
1333: . cgeom - The cell geometry for each cell in the chunk
1334: . coefficients - The array of FEM basis coefficients for the elements
1335: . probAux - The PetscDS specifying the auxiliary discretizations
1336: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1338: Output Parameter:
1339: . integral - the integral for this field
1341: Level: intermediate
1343: .seealso: PetscFEIntegrateResidual()
1344: @*/
1345: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1346: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1347: {
1348: PetscFE fe;
1353: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1354: if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1355: return(0);
1356: }
1358: /*@C
1359: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1361: Not collective
1363: Input Parameters:
1364: + prob - The PetscDS specifying the discretizations and continuum functions
1365: . field - The field being integrated
1366: . obj_func - The function to be integrated
1367: . Ne - The number of elements in the chunk
1368: . fgeom - The face geometry for each face in the chunk
1369: . coefficients - The array of FEM basis coefficients for the elements
1370: . probAux - The PetscDS specifying the auxiliary discretizations
1371: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1373: Output Parameter:
1374: . integral - the integral for this field
1376: Level: intermediate
1378: .seealso: PetscFEIntegrateResidual()
1379: @*/
1380: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1381: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1382: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1383: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1384: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1385: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1386: {
1387: PetscFE fe;
1392: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1393: if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1394: return(0);
1395: }
1397: /*@C
1398: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1400: Not collective
1402: Input Parameters:
1403: + prob - The PetscDS specifying the discretizations and continuum functions
1404: . field - The field being integrated
1405: . Ne - The number of elements in the chunk
1406: . cgeom - The cell geometry for each cell in the chunk
1407: . coefficients - The array of FEM basis coefficients for the elements
1408: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1409: . probAux - The PetscDS specifying the auxiliary discretizations
1410: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1411: - t - The time
1413: Output Parameter:
1414: . elemVec - the element residual vectors from each element
1416: Note:
1417: $ Loop over batch of elements (e):
1418: $ Loop over quadrature points (q):
1419: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1420: $ Call f_0 and f_1
1421: $ Loop over element vector entries (f,fc --> i):
1422: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1424: Level: intermediate
1426: .seealso: PetscFEIntegrateResidual()
1427: @*/
1428: PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1429: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1430: {
1431: PetscFE fe;
1436: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1437: if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1438: return(0);
1439: }
1441: /*@C
1442: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1444: Not collective
1446: Input Parameters:
1447: + prob - The PetscDS specifying the discretizations and continuum functions
1448: . field - The field being integrated
1449: . Ne - The number of elements in the chunk
1450: . fgeom - The face geometry for each cell in the chunk
1451: . coefficients - The array of FEM basis coefficients for the elements
1452: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1453: . probAux - The PetscDS specifying the auxiliary discretizations
1454: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1455: - t - The time
1457: Output Parameter:
1458: . elemVec - the element residual vectors from each element
1460: Level: intermediate
1462: .seealso: PetscFEIntegrateResidual()
1463: @*/
1464: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1465: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1466: {
1467: PetscFE fe;
1472: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1473: if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1474: return(0);
1475: }
1477: /*@C
1478: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1480: Not collective
1482: Input Parameters:
1483: + prob - The PetscDS specifying the discretizations and continuum functions
1484: . field - The field being integrated
1485: . Ne - The number of elements in the chunk
1486: . fgeom - The face geometry for each cell in the chunk
1487: . coefficients - The array of FEM basis coefficients for the elements
1488: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1489: . probAux - The PetscDS specifying the auxiliary discretizations
1490: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1491: - t - The time
1493: Output Parameter
1494: . elemVec - the element residual vectors from each element
1496: Level: developer
1498: .seealso: PetscFEIntegrateResidual()
1499: @*/
1500: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1501: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1502: {
1503: PetscFE fe;
1508: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1509: if (fe->ops->integratehybridresidual) {(*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1510: return(0);
1511: }
1513: /*@C
1514: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1516: Not collective
1518: Input Parameters:
1519: + prob - The PetscDS specifying the discretizations and continuum functions
1520: . jtype - The type of matrix pointwise functions that should be used
1521: . fieldI - The test field being integrated
1522: . fieldJ - The basis field being integrated
1523: . Ne - The number of elements in the chunk
1524: . cgeom - The cell geometry for each cell in the chunk
1525: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1526: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1527: . probAux - The PetscDS specifying the auxiliary discretizations
1528: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1529: . t - The time
1530: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1532: Output Parameter:
1533: . elemMat - the element matrices for the Jacobian from each element
1535: Note:
1536: $ Loop over batch of elements (e):
1537: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1538: $ Loop over quadrature points (q):
1539: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1540: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1541: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1542: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1543: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1544: Level: intermediate
1546: .seealso: PetscFEIntegrateResidual()
1547: @*/
1548: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1549: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1550: {
1551: PetscFE fe;
1556: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1557: if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1558: return(0);
1559: }
1561: /*@C
1562: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1564: Not collective
1566: Input Parameters:
1567: + prob - The PetscDS specifying the discretizations and continuum functions
1568: . fieldI - The test field being integrated
1569: . fieldJ - The basis field being integrated
1570: . Ne - The number of elements in the chunk
1571: . fgeom - The face geometry for each cell in the chunk
1572: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1573: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1574: . probAux - The PetscDS specifying the auxiliary discretizations
1575: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1576: . t - The time
1577: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1579: Output Parameter:
1580: . elemMat - the element matrices for the Jacobian from each element
1582: Note:
1583: $ Loop over batch of elements (e):
1584: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1585: $ Loop over quadrature points (q):
1586: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1587: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1588: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1589: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1590: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1591: Level: intermediate
1593: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1594: @*/
1595: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1596: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1597: {
1598: PetscFE fe;
1603: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1604: if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1605: return(0);
1606: }
1608: /*@C
1609: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1611: Not collective
1613: Input Parameters:
1614: + prob - The PetscDS specifying the discretizations and continuum functions
1615: . jtype - The type of matrix pointwise functions that should be used
1616: . fieldI - The test field being integrated
1617: . fieldJ - The basis field being integrated
1618: . Ne - The number of elements in the chunk
1619: . fgeom - The face geometry for each cell in the chunk
1620: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1621: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1622: . probAux - The PetscDS specifying the auxiliary discretizations
1623: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1624: . t - The time
1625: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1627: Output Parameter
1628: . elemMat - the element matrices for the Jacobian from each element
1630: Note:
1631: $ Loop over batch of elements (e):
1632: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1633: $ Loop over quadrature points (q):
1634: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1635: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1636: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1637: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1638: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1639: Level: developer
1641: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1642: @*/
1643: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1644: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1645: {
1646: PetscFE fe;
1651: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1652: if (fe->ops->integratehybridjacobian) {(*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1653: return(0);
1654: }
1656: /*@
1657: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1659: Input Parameters:
1660: + fe - The finite element space
1661: - height - The height of the Plex point
1663: Output Parameter:
1664: . subfe - The subspace of this FE space
1666: Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1668: Level: advanced
1670: .seealso: PetscFECreateDefault()
1671: @*/
1672: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1673: {
1674: PetscSpace P, subP;
1675: PetscDualSpace Q, subQ;
1676: PetscQuadrature subq;
1677: PetscFEType fetype;
1678: PetscInt dim, Nc;
1679: PetscErrorCode ierr;
1684: if (height == 0) {
1685: *subfe = fe;
1686: return(0);
1687: }
1688: PetscFEGetBasisSpace(fe, &P);
1689: PetscFEGetDualSpace(fe, &Q);
1690: PetscFEGetNumComponents(fe, &Nc);
1691: PetscFEGetFaceQuadrature(fe, &subq);
1692: PetscDualSpaceGetDimension(Q, &dim);
1693: if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
1694: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1695: if (height <= dim) {
1696: if (!fe->subspaces[height-1]) {
1697: PetscFE sub = NULL;
1698: const char *name;
1700: PetscSpaceGetHeightSubspace(P, height, &subP);
1701: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1702: if (subQ) {
1703: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1704: PetscObjectGetName((PetscObject) fe, &name);
1705: PetscObjectSetName((PetscObject) sub, name);
1706: PetscFEGetType(fe, &fetype);
1707: PetscFESetType(sub, fetype);
1708: PetscFESetBasisSpace(sub, subP);
1709: PetscFESetDualSpace(sub, subQ);
1710: PetscFESetNumComponents(sub, Nc);
1711: PetscFESetUp(sub);
1712: PetscFESetQuadrature(sub, subq);
1713: }
1714: fe->subspaces[height-1] = sub;
1715: }
1716: *subfe = fe->subspaces[height-1];
1717: } else {
1718: *subfe = NULL;
1719: }
1720: return(0);
1721: }
1723: /*@
1724: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1725: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1726: sparsity). It is also used to create an interpolation between regularly refined meshes.
1728: Collective on fem
1730: Input Parameter:
1731: . fe - The initial PetscFE
1733: Output Parameter:
1734: . feRef - The refined PetscFE
1736: Level: advanced
1738: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1739: @*/
1740: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1741: {
1742: PetscSpace P, Pref;
1743: PetscDualSpace Q, Qref;
1744: DM K, Kref;
1745: PetscQuadrature q, qref;
1746: const PetscReal *v0, *jac;
1747: PetscInt numComp, numSubelements;
1748: PetscInt cStart, cEnd, c;
1749: PetscDualSpace *cellSpaces;
1750: PetscErrorCode ierr;
1753: PetscFEGetBasisSpace(fe, &P);
1754: PetscFEGetDualSpace(fe, &Q);
1755: PetscFEGetQuadrature(fe, &q);
1756: PetscDualSpaceGetDM(Q, &K);
1757: /* Create space */
1758: PetscObjectReference((PetscObject) P);
1759: Pref = P;
1760: /* Create dual space */
1761: PetscDualSpaceDuplicate(Q, &Qref);
1762: PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1763: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1764: PetscDualSpaceSetDM(Qref, Kref);
1765: DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1766: PetscMalloc1(cEnd - cStart, &cellSpaces);
1767: /* TODO: fix for non-uniform refinement */
1768: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1769: PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1770: PetscFree(cellSpaces);
1771: DMDestroy(&Kref);
1772: PetscDualSpaceSetUp(Qref);
1773: /* Create element */
1774: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1775: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1776: PetscFESetBasisSpace(*feRef, Pref);
1777: PetscFESetDualSpace(*feRef, Qref);
1778: PetscFEGetNumComponents(fe, &numComp);
1779: PetscFESetNumComponents(*feRef, numComp);
1780: PetscFESetUp(*feRef);
1781: PetscSpaceDestroy(&Pref);
1782: PetscDualSpaceDestroy(&Qref);
1783: /* Create quadrature */
1784: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1785: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1786: PetscFESetQuadrature(*feRef, qref);
1787: PetscQuadratureDestroy(&qref);
1788: return(0);
1789: }
1791: /*@C
1792: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1794: Collective
1796: Input Parameters:
1797: + comm - The MPI comm
1798: . dim - The spatial dimension
1799: . Nc - The number of components
1800: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1801: . prefix - The options prefix, or NULL
1802: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1804: Output Parameter:
1805: . fem - The PetscFE object
1807: Note:
1808: Each object is SetFromOption() during creation, so that the object may be customized from the command line.
1810: Level: beginner
1812: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1813: @*/
1814: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1815: {
1816: PetscQuadrature q, fq;
1817: DM K;
1818: PetscSpace P;
1819: PetscDualSpace Q;
1820: PetscInt order, quadPointsPerEdge;
1821: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1822: PetscErrorCode ierr;
1825: /* Create space */
1826: PetscSpaceCreate(comm, &P);
1827: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1828: PetscSpacePolynomialSetTensor(P, tensor);
1829: PetscSpaceSetNumComponents(P, Nc);
1830: PetscSpaceSetNumVariables(P, dim);
1831: PetscSpaceSetFromOptions(P);
1832: PetscSpaceSetUp(P);
1833: PetscSpaceGetDegree(P, &order, NULL);
1834: PetscSpacePolynomialGetTensor(P, &tensor);
1835: /* Create dual space */
1836: PetscDualSpaceCreate(comm, &Q);
1837: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1838: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1839: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1840: PetscDualSpaceSetDM(Q, K);
1841: DMDestroy(&K);
1842: PetscDualSpaceSetNumComponents(Q, Nc);
1843: PetscDualSpaceSetOrder(Q, order);
1844: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1845: PetscDualSpaceSetFromOptions(Q);
1846: PetscDualSpaceSetUp(Q);
1847: /* Create element */
1848: PetscFECreate(comm, fem);
1849: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1850: PetscFESetBasisSpace(*fem, P);
1851: PetscFESetDualSpace(*fem, Q);
1852: PetscFESetNumComponents(*fem, Nc);
1853: PetscFESetFromOptions(*fem);
1854: PetscFESetUp(*fem);
1855: PetscSpaceDestroy(&P);
1856: PetscDualSpaceDestroy(&Q);
1857: /* Create quadrature (with specified order if given) */
1858: qorder = qorder >= 0 ? qorder : order;
1859: PetscObjectOptionsBegin((PetscObject)*fem);
1860: PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1861: PetscOptionsEnd();
1862: quadPointsPerEdge = PetscMax(qorder + 1,1);
1863: if (isSimplex) {
1864: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1865: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1866: } else {
1867: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1868: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1869: }
1870: PetscFESetQuadrature(*fem, q);
1871: PetscFESetFaceQuadrature(*fem, fq);
1872: PetscQuadratureDestroy(&q);
1873: PetscQuadratureDestroy(&fq);
1874: return(0);
1875: }
1877: /*@
1878: PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k
1880: Collective
1882: Input Parameters:
1883: + comm - The MPI comm
1884: . dim - The spatial dimension
1885: . Nc - The number of components
1886: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1887: . k - The degree k of the space
1888: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1890: Output Parameter:
1891: . fem - The PetscFE object
1893: Level: beginner
1895: Notes:
1896: For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k.
1898: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1899: @*/
1900: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1901: {
1902: PetscQuadrature q, fq;
1903: DM K;
1904: PetscSpace P;
1905: PetscDualSpace Q;
1906: PetscInt quadPointsPerEdge;
1907: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1908: char name[64];
1909: PetscErrorCode ierr;
1912: /* Create space */
1913: PetscSpaceCreate(comm, &P);
1914: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1915: PetscSpacePolynomialSetTensor(P, tensor);
1916: PetscSpaceSetNumComponents(P, Nc);
1917: PetscSpaceSetNumVariables(P, dim);
1918: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1919: PetscSpaceSetUp(P);
1920: /* Create dual space */
1921: PetscDualSpaceCreate(comm, &Q);
1922: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1923: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1924: PetscDualSpaceSetDM(Q, K);
1925: DMDestroy(&K);
1926: PetscDualSpaceSetNumComponents(Q, Nc);
1927: PetscDualSpaceSetOrder(Q, k);
1928: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1929: PetscDualSpaceSetUp(Q);
1930: /* Create finite element */
1931: PetscFECreate(comm, fem);
1932: PetscFESetType(*fem, PETSCFEBASIC);
1933: PetscFESetBasisSpace(*fem, P);
1934: PetscFESetDualSpace(*fem, Q);
1935: PetscFESetNumComponents(*fem, Nc);
1936: PetscFESetUp(*fem);
1937: PetscSpaceDestroy(&P);
1938: PetscDualSpaceDestroy(&Q);
1939: /* Create quadrature (with specified order if given) */
1940: qorder = qorder >= 0 ? qorder : k;
1941: quadPointsPerEdge = PetscMax(qorder + 1,1);
1942: if (isSimplex) {
1943: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1944: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1945: } else {
1946: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1947: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1948: }
1949: PetscFESetQuadrature(*fem, q);
1950: PetscFESetFaceQuadrature(*fem, fq);
1951: PetscQuadratureDestroy(&q);
1952: PetscQuadratureDestroy(&fq);
1953: /* Set finite element name */
1954: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1955: PetscFESetName(*fem, name);
1956: return(0);
1957: }
1959: /*@C
1960: PetscFESetName - Names the FE and its subobjects
1962: Not collective
1964: Input Parameters:
1965: + fe - The PetscFE
1966: - name - The name
1968: Level: intermediate
1970: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1971: @*/
1972: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1973: {
1974: PetscSpace P;
1975: PetscDualSpace Q;
1979: PetscFEGetBasisSpace(fe, &P);
1980: PetscFEGetDualSpace(fe, &Q);
1981: PetscObjectSetName((PetscObject) fe, name);
1982: PetscObjectSetName((PetscObject) P, name);
1983: PetscObjectSetName((PetscObject) Q, name);
1984: return(0);
1985: }
1987: 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[])
1988: {
1989: PetscInt dOffset = 0, fOffset = 0, f;
1992: for (f = 0; f < Nf; ++f) {
1993: PetscFE fe;
1994: const PetscInt cdim = T[f]->cdim;
1995: const PetscInt Nq = T[f]->Np;
1996: const PetscInt Nbf = T[f]->Nb;
1997: const PetscInt Ncf = T[f]->Nc;
1998: const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
1999: const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2000: PetscInt b, c, d;
2002: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2003: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2004: for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2005: for (b = 0; b < Nbf; ++b) {
2006: for (c = 0; c < Ncf; ++c) {
2007: const PetscInt cidx = b*Ncf+c;
2009: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2010: for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2011: }
2012: }
2013: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2014: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2015: if (u_t) {
2016: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2017: for (b = 0; b < Nbf; ++b) {
2018: for (c = 0; c < Ncf; ++c) {
2019: const PetscInt cidx = b*Ncf+c;
2021: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2022: }
2023: }
2024: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2025: }
2026: fOffset += Ncf;
2027: dOffset += Nbf;
2028: }
2029: return 0;
2030: }
2032: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2033: {
2034: PetscInt dOffset = 0, fOffset = 0, g;
2037: for (g = 0; g < 2*Nf-1; ++g) {
2038: if (!T[g/2]) continue;
2039: {
2040: PetscFE fe;
2041: const PetscInt f = g/2;
2042: const PetscInt cdim = T[f]->cdim;
2043: const PetscInt Nq = T[f]->Np;
2044: const PetscInt Nbf = T[f]->Nb;
2045: const PetscInt Ncf = T[f]->Nc;
2046: const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2047: const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2048: PetscInt b, c, d;
2050: fe = (PetscFE) ds->disc[f];
2051: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2052: for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2053: for (b = 0; b < Nbf; ++b) {
2054: for (c = 0; c < Ncf; ++c) {
2055: const PetscInt cidx = b*Ncf+c;
2057: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2058: for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2059: }
2060: }
2061: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2062: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2063: if (u_t) {
2064: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2065: for (b = 0; b < Nbf; ++b) {
2066: for (c = 0; c < Ncf; ++c) {
2067: const PetscInt cidx = b*Ncf+c;
2069: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2070: }
2071: }
2072: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2073: }
2074: fOffset += Ncf;
2075: dOffset += Nbf;
2076: }
2077: }
2078: return 0;
2079: }
2081: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2082: {
2083: PetscFE fe;
2084: PetscTabulation Tc;
2085: PetscInt b, c;
2086: PetscErrorCode ierr;
2088: if (!prob) return 0;
2089: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2090: PetscFEGetFaceCentroidTabulation(fe, &Tc);
2091: {
2092: const PetscReal *faceBasis = Tc->T[0];
2093: const PetscInt Nb = Tc->Nb;
2094: const PetscInt Nc = Tc->Nc;
2096: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2097: for (b = 0; b < Nb; ++b) {
2098: for (c = 0; c < Nc; ++c) {
2099: u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2100: }
2101: }
2102: }
2103: return 0;
2104: }
2106: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2107: {
2108: const PetscInt dE = T->cdim; /* fegeom->dimEmbed */
2109: const PetscInt Nq = T->Np;
2110: const PetscInt Nb = T->Nb;
2111: const PetscInt Nc = T->Nc;
2112: const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2113: const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2114: PetscInt q, b, c, d;
2115: PetscErrorCode ierr;
2117: for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2118: for (q = 0; q < Nq; ++q) {
2119: for (b = 0; b < Nb; ++b) {
2120: for (c = 0; c < Nc; ++c) {
2121: const PetscInt bcidx = b*Nc+c;
2123: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2124: for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2125: }
2126: }
2127: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2128: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2129: for (b = 0; b < Nb; ++b) {
2130: for (c = 0; c < Nc; ++c) {
2131: const PetscInt bcidx = b*Nc+c;
2132: const PetscInt qcidx = q*Nc+c;
2134: elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2135: for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2136: }
2137: }
2138: }
2139: return(0);
2140: }
2142: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2143: {
2144: const PetscInt dE = T->cdim;
2145: const PetscInt Nq = T->Np;
2146: const PetscInt Nb = T->Nb;
2147: const PetscInt Nc = T->Nc;
2148: const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2149: const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2150: PetscInt q, b, c, d, s;
2151: PetscErrorCode ierr;
2153: for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0;
2154: for (q = 0; q < Nq; ++q) {
2155: for (b = 0; b < Nb; ++b) {
2156: for (c = 0; c < Nc; ++c) {
2157: const PetscInt bcidx = b*Nc+c;
2159: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2160: for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2161: }
2162: }
2163: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2164: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2165: for (s = 0; s < 2; ++s) {
2166: for (b = 0; b < Nb; ++b) {
2167: for (c = 0; c < Nc; ++c) {
2168: const PetscInt bcidx = b*Nc+c;
2169: const PetscInt qcidx = (q*2+s)*Nc+c;
2171: elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2172: for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2173: }
2174: }
2175: }
2176: }
2177: return(0);
2178: }
2180: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2181: {
2182: const PetscInt dE = TI->cdim;
2183: const PetscInt NqI = TI->Np;
2184: const PetscInt NbI = TI->Nb;
2185: const PetscInt NcI = TI->Nc;
2186: const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2187: const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2188: const PetscInt NqJ = TJ->Np;
2189: const PetscInt NbJ = TJ->Nb;
2190: const PetscInt NcJ = TJ->Nc;
2191: const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2192: const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2193: PetscInt f, fc, g, gc, df, dg;
2194: PetscErrorCode ierr;
2196: for (f = 0; f < NbI; ++f) {
2197: for (fc = 0; fc < NcI; ++fc) {
2198: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2200: tmpBasisI[fidx] = basisI[fidx];
2201: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2202: }
2203: }
2204: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2205: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2206: for (g = 0; g < NbJ; ++g) {
2207: for (gc = 0; gc < NcJ; ++gc) {
2208: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2210: tmpBasisJ[gidx] = basisJ[gidx];
2211: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2212: }
2213: }
2214: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2215: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2216: for (f = 0; f < NbI; ++f) {
2217: for (fc = 0; fc < NcI; ++fc) {
2218: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2219: const PetscInt i = offsetI+f; /* Element matrix row */
2220: for (g = 0; g < NbJ; ++g) {
2221: for (gc = 0; gc < NcJ; ++gc) {
2222: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2223: const PetscInt j = offsetJ+g; /* Element matrix column */
2224: const PetscInt fOff = eOffset+i*totDim+j;
2226: elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2227: for (df = 0; df < dE; ++df) {
2228: elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2229: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2230: for (dg = 0; dg < dE; ++dg) {
2231: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2232: }
2233: }
2234: }
2235: }
2236: }
2237: }
2238: return(0);
2239: }
2241: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2242: {
2243: const PetscInt dE = TI->cdim;
2244: const PetscInt NqI = TI->Np;
2245: const PetscInt NbI = TI->Nb;
2246: const PetscInt NcI = TI->Nc;
2247: const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2248: const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2249: const PetscInt NqJ = TJ->Np;
2250: const PetscInt NbJ = TJ->Nb;
2251: const PetscInt NcJ = TJ->Nc;
2252: const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2253: const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2254: const PetscInt Ns = isHybridI ? 1 : 2;
2255: const PetscInt Nt = isHybridJ ? 1 : 2;
2256: PetscInt f, fc, g, gc, df, dg, s, t;
2257: PetscErrorCode ierr;
2259: for (f = 0; f < NbI; ++f) {
2260: for (fc = 0; fc < NcI; ++fc) {
2261: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2263: tmpBasisI[fidx] = basisI[fidx];
2264: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2265: }
2266: }
2267: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2268: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2269: for (g = 0; g < NbJ; ++g) {
2270: for (gc = 0; gc < NcJ; ++gc) {
2271: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2273: tmpBasisJ[gidx] = basisJ[gidx];
2274: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2275: }
2276: }
2277: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2278: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2279: for (s = 0; s < Ns; ++s) {
2280: for (f = 0; f < NbI; ++f) {
2281: for (fc = 0; fc < NcI; ++fc) {
2282: const PetscInt sc = NcI*s+fc; /* components from each side of the surface */
2283: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2284: const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */
2285: for (t = 0; t < Nt; ++t) {
2286: for (g = 0; g < NbJ; ++g) {
2287: for (gc = 0; gc < NcJ; ++gc) {
2288: const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */
2289: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2290: const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */
2291: const PetscInt fOff = eOffset+i*totDim+j;
2293: elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2294: for (df = 0; df < dE; ++df) {
2295: elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2296: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2297: for (dg = 0; dg < dE; ++dg) {
2298: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2299: }
2300: }
2301: }
2302: }
2303: }
2304: }
2305: }
2306: }
2307: return(0);
2308: }
2310: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2311: {
2312: PetscDualSpace dsp;
2313: DM dm;
2314: PetscQuadrature quadDef;
2315: PetscInt dim, cdim, Nq;
2316: PetscErrorCode ierr;
2319: PetscFEGetDualSpace(fe, &dsp);
2320: PetscDualSpaceGetDM(dsp, &dm);
2321: DMGetDimension(dm, &dim);
2322: DMGetCoordinateDim(dm, &cdim);
2323: PetscFEGetQuadrature(fe, &quadDef);
2324: quad = quad ? quad : quadDef;
2325: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2326: PetscMalloc1(Nq*cdim, &cgeom->v);
2327: PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2328: PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2329: PetscMalloc1(Nq, &cgeom->detJ);
2330: cgeom->dim = dim;
2331: cgeom->dimEmbed = cdim;
2332: cgeom->numCells = 1;
2333: cgeom->numPoints = Nq;
2334: DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2335: return(0);
2336: }
2338: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2339: {
2343: PetscFree(cgeom->v);
2344: PetscFree(cgeom->J);
2345: PetscFree(cgeom->invJ);
2346: PetscFree(cgeom->detJ);
2347: return(0);
2348: }