Actual source code: fe.c
petsc-3.11.4 2019-09-28
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: .keywords: PetscFE, register
86: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
88: @*/
89: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
90: {
94: PetscFunctionListAdd(&PetscFEList, sname, function);
95: return(0);
96: }
98: /*@C
99: PetscFESetType - Builds a particular PetscFE
101: Collective on PetscFE
103: Input Parameters:
104: + fem - The PetscFE object
105: - name - The kind of FEM space
107: Options Database Key:
108: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
110: Level: intermediate
112: .keywords: PetscFE, set, type
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: .keywords: PetscFE, get, type, name
153: .seealso: PetscFESetType(), PetscFECreate()
154: @*/
155: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
156: {
162: if (!PetscFERegisterAllCalled) {
163: PetscFERegisterAll();
164: }
165: *name = ((PetscObject) fem)->type_name;
166: return(0);
167: }
169: /*@C
170: PetscFEView - Views a PetscFE
172: Collective on PetscFE
174: Input Parameter:
175: + fem - the PetscFE object to view
176: - viewer - the viewer
178: Level: developer
180: .seealso PetscFEDestroy()
181: @*/
182: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
183: {
184: PetscBool iascii;
190: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
191: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
192: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
193: if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
194: return(0);
195: }
197: /*@
198: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
200: Collective on PetscFE
202: Input Parameter:
203: . fem - the PetscFE object to set options for
205: Options Database:
206: . -petscfe_num_blocks the number of cell blocks to integrate concurrently
207: . -petscfe_num_batches the number of cell batches to integrate serially
209: Level: developer
211: .seealso PetscFEView()
212: @*/
213: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
214: {
215: const char *defaultType;
216: char name[256];
217: PetscBool flg;
222: if (!((PetscObject) fem)->type_name) {
223: defaultType = PETSCFEBASIC;
224: } else {
225: defaultType = ((PetscObject) fem)->type_name;
226: }
227: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
229: PetscObjectOptionsBegin((PetscObject) fem);
230: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
231: if (flg) {
232: PetscFESetType(fem, name);
233: } else if (!((PetscObject) fem)->type_name) {
234: PetscFESetType(fem, defaultType);
235: }
236: PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
237: PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
238: if (fem->ops->setfromoptions) {
239: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
240: }
241: /* process any options handlers added with PetscObjectAddOptionsHandler() */
242: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
243: PetscOptionsEnd();
244: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
245: return(0);
246: }
248: /*@C
249: PetscFESetUp - Construct data structures for the PetscFE
251: Collective on PetscFE
253: Input Parameter:
254: . fem - the PetscFE object to setup
256: Level: developer
258: .seealso PetscFEView(), PetscFEDestroy()
259: @*/
260: PetscErrorCode PetscFESetUp(PetscFE fem)
261: {
266: if (fem->setupcalled) return(0);
267: fem->setupcalled = PETSC_TRUE;
268: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
269: return(0);
270: }
272: /*@
273: PetscFEDestroy - Destroys a PetscFE object
275: Collective on PetscFE
277: Input Parameter:
278: . fem - the PetscFE object to destroy
280: Level: developer
282: .seealso PetscFEView()
283: @*/
284: PetscErrorCode PetscFEDestroy(PetscFE *fem)
285: {
289: if (!*fem) return(0);
292: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
293: ((PetscObject) (*fem))->refct = 0;
295: if ((*fem)->subspaces) {
296: PetscInt dim, d;
298: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
299: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
300: }
301: PetscFree((*fem)->subspaces);
302: PetscFree((*fem)->invV);
303: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
304: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
305: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
306: PetscSpaceDestroy(&(*fem)->basisSpace);
307: PetscDualSpaceDestroy(&(*fem)->dualSpace);
308: PetscQuadratureDestroy(&(*fem)->quadrature);
309: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
311: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
312: PetscHeaderDestroy(fem);
313: return(0);
314: }
316: /*@
317: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
319: Collective on MPI_Comm
321: Input Parameter:
322: . comm - The communicator for the PetscFE object
324: Output Parameter:
325: . fem - The PetscFE object
327: Level: beginner
329: .seealso: PetscFESetType(), PETSCFEGALERKIN
330: @*/
331: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
332: {
333: PetscFE f;
338: PetscCitationsRegister(FECitation,&FEcite);
339: *fem = NULL;
340: PetscFEInitializePackage();
342: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
344: f->basisSpace = NULL;
345: f->dualSpace = NULL;
346: f->numComponents = 1;
347: f->subspaces = NULL;
348: f->invV = NULL;
349: f->B = NULL;
350: f->D = NULL;
351: f->H = NULL;
352: f->Bf = NULL;
353: f->Df = NULL;
354: f->Hf = NULL;
355: PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
356: PetscMemzero(&f->faceQuadrature, sizeof(PetscQuadrature));
357: f->blockSize = 0;
358: f->numBlocks = 1;
359: f->batchSize = 0;
360: f->numBatches = 1;
362: *fem = f;
363: return(0);
364: }
366: /*@
367: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
369: Not collective
371: Input Parameter:
372: . fem - The PetscFE object
374: Output Parameter:
375: . dim - The spatial dimension
377: Level: intermediate
379: .seealso: PetscFECreate()
380: @*/
381: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
382: {
383: DM dm;
389: PetscDualSpaceGetDM(fem->dualSpace, &dm);
390: DMGetDimension(dm, dim);
391: return(0);
392: }
394: /*@
395: PetscFESetNumComponents - Sets the number of components in the element
397: Not collective
399: Input Parameters:
400: + fem - The PetscFE object
401: - comp - The number of field components
403: Level: intermediate
405: .seealso: PetscFECreate()
406: @*/
407: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
408: {
411: fem->numComponents = comp;
412: return(0);
413: }
415: /*@
416: PetscFEGetNumComponents - Returns the number of components in the element
418: Not collective
420: Input Parameter:
421: . fem - The PetscFE object
423: Output Parameter:
424: . comp - The number of field components
426: Level: intermediate
428: .seealso: PetscFECreate()
429: @*/
430: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
431: {
435: *comp = fem->numComponents;
436: return(0);
437: }
439: /*@
440: PetscFESetTileSizes - Sets the tile sizes for evaluation
442: Not collective
444: Input Parameters:
445: + fem - The PetscFE object
446: . blockSize - The number of elements in a block
447: . numBlocks - The number of blocks in a batch
448: . batchSize - The number of elements in a batch
449: - numBatches - The number of batches in a chunk
451: Level: intermediate
453: .seealso: PetscFECreate()
454: @*/
455: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
456: {
459: fem->blockSize = blockSize;
460: fem->numBlocks = numBlocks;
461: fem->batchSize = batchSize;
462: fem->numBatches = numBatches;
463: return(0);
464: }
466: /*@
467: PetscFEGetTileSizes - Returns the tile sizes for evaluation
469: Not collective
471: Input Parameter:
472: . fem - The PetscFE object
474: Output Parameters:
475: + blockSize - The number of elements in a block
476: . numBlocks - The number of blocks in a batch
477: . batchSize - The number of elements in a batch
478: - numBatches - The number of batches in a chunk
480: Level: intermediate
482: .seealso: PetscFECreate()
483: @*/
484: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
485: {
492: if (blockSize) *blockSize = fem->blockSize;
493: if (numBlocks) *numBlocks = fem->numBlocks;
494: if (batchSize) *batchSize = fem->batchSize;
495: if (numBatches) *numBatches = fem->numBatches;
496: return(0);
497: }
499: /*@
500: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
502: Not collective
504: Input Parameter:
505: . fem - The PetscFE object
507: Output Parameter:
508: . sp - The PetscSpace object
510: Level: intermediate
512: .seealso: PetscFECreate()
513: @*/
514: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
515: {
519: *sp = fem->basisSpace;
520: return(0);
521: }
523: /*@
524: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
526: Not collective
528: Input Parameters:
529: + fem - The PetscFE object
530: - sp - The PetscSpace object
532: Level: intermediate
534: .seealso: PetscFECreate()
535: @*/
536: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
537: {
543: PetscSpaceDestroy(&fem->basisSpace);
544: fem->basisSpace = sp;
545: PetscObjectReference((PetscObject) fem->basisSpace);
546: return(0);
547: }
549: /*@
550: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
552: Not collective
554: Input Parameter:
555: . fem - The PetscFE object
557: Output Parameter:
558: . sp - The PetscDualSpace object
560: Level: intermediate
562: .seealso: PetscFECreate()
563: @*/
564: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
565: {
569: *sp = fem->dualSpace;
570: return(0);
571: }
573: /*@
574: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
576: Not collective
578: Input Parameters:
579: + fem - The PetscFE object
580: - sp - The PetscDualSpace object
582: Level: intermediate
584: .seealso: PetscFECreate()
585: @*/
586: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
587: {
593: PetscDualSpaceDestroy(&fem->dualSpace);
594: fem->dualSpace = sp;
595: PetscObjectReference((PetscObject) fem->dualSpace);
596: return(0);
597: }
599: /*@
600: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
602: Not collective
604: Input Parameter:
605: . fem - The PetscFE object
607: Output Parameter:
608: . q - The PetscQuadrature object
610: Level: intermediate
612: .seealso: PetscFECreate()
613: @*/
614: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
615: {
619: *q = fem->quadrature;
620: return(0);
621: }
623: /*@
624: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
626: Not collective
628: Input Parameters:
629: + fem - The PetscFE object
630: - q - The PetscQuadrature object
632: Level: intermediate
634: .seealso: PetscFECreate()
635: @*/
636: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
637: {
638: PetscInt Nc, qNc;
643: PetscFEGetNumComponents(fem, &Nc);
644: PetscQuadratureGetNumComponents(q, &qNc);
645: 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);
646: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
647: PetscQuadratureDestroy(&fem->quadrature);
648: fem->quadrature = q;
649: PetscObjectReference((PetscObject) q);
650: return(0);
651: }
653: /*@
654: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
656: Not collective
658: Input Parameter:
659: . fem - The PetscFE object
661: Output Parameter:
662: . q - The PetscQuadrature object
664: Level: intermediate
666: .seealso: PetscFECreate()
667: @*/
668: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
669: {
673: *q = fem->faceQuadrature;
674: return(0);
675: }
677: /*@
678: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
680: Not collective
682: Input Parameters:
683: + fem - The PetscFE object
684: - q - The PetscQuadrature object
686: Level: intermediate
688: .seealso: PetscFECreate()
689: @*/
690: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
691: {
696: PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
697: PetscQuadratureDestroy(&fem->faceQuadrature);
698: fem->faceQuadrature = q;
699: PetscObjectReference((PetscObject) q);
700: return(0);
701: }
703: /*@C
704: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
706: Not collective
708: Input Parameter:
709: . fem - The PetscFE object
711: Output Parameter:
712: . numDof - Array with the number of dofs per dimension
714: Level: intermediate
716: .seealso: PetscFECreate()
717: @*/
718: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
719: {
725: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
726: return(0);
727: }
729: /*@C
730: PetscFEGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points
732: Not collective
734: Input Parameter:
735: . fem - The PetscFE object
737: Output Parameters:
738: + B - The basis function values at quadrature points
739: . D - The basis function derivatives at quadrature points
740: - H - The basis function second derivatives at quadrature points
742: Note:
743: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
744: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
745: $ 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
747: Level: intermediate
749: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
750: @*/
751: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
752: {
753: PetscInt npoints;
754: const PetscReal *points;
755: PetscErrorCode ierr;
762: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
763: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
764: if (B) *B = fem->B;
765: if (D) *D = fem->D;
766: if (H) *H = fem->H;
767: return(0);
768: }
770: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
771: {
772: PetscErrorCode ierr;
779: if (!fem->Bf) {
780: const PetscReal xi0[3] = {-1., -1., -1.};
781: PetscReal v0[3], J[9], detJ;
782: PetscQuadrature fq;
783: PetscDualSpace sp;
784: DM dm;
785: const PetscInt *faces;
786: PetscInt dim, numFaces, f, npoints, q;
787: const PetscReal *points;
788: PetscReal *facePoints;
790: PetscFEGetDualSpace(fem, &sp);
791: PetscDualSpaceGetDM(sp, &dm);
792: DMGetDimension(dm, &dim);
793: DMPlexGetConeSize(dm, 0, &numFaces);
794: DMPlexGetCone(dm, 0, &faces);
795: PetscFEGetFaceQuadrature(fem, &fq);
796: if (fq) {
797: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
798: PetscMalloc1(numFaces*npoints*dim, &facePoints);
799: for (f = 0; f < numFaces; ++f) {
800: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
801: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
802: }
803: PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
804: PetscFree(facePoints);
805: }
806: }
807: if (Bf) *Bf = fem->Bf;
808: if (Df) *Df = fem->Df;
809: if (Hf) *Hf = fem->Hf;
810: return(0);
811: }
813: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
814: {
815: PetscErrorCode ierr;
820: if (!fem->F) {
821: PetscDualSpace sp;
822: DM dm;
823: const PetscInt *cone;
824: PetscReal *centroids;
825: PetscInt dim, numFaces, f;
827: PetscFEGetDualSpace(fem, &sp);
828: PetscDualSpaceGetDM(sp, &dm);
829: DMGetDimension(dm, &dim);
830: DMPlexGetConeSize(dm, 0, &numFaces);
831: DMPlexGetCone(dm, 0, &cone);
832: PetscMalloc1(numFaces*dim, ¢roids);
833: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
834: PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
835: PetscFree(centroids);
836: }
837: *F = fem->F;
838: return(0);
839: }
841: /*@C
842: PetscFEGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
844: Not collective
846: Input Parameters:
847: + fem - The PetscFE object
848: . npoints - The number of tabulation points
849: - points - The tabulation point coordinates
851: Output Parameters:
852: + B - The basis function values at tabulation points
853: . D - The basis function derivatives at tabulation points
854: - H - The basis function second derivatives at tabulation points
856: Note:
857: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
858: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
859: $ 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
861: Level: intermediate
863: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
864: @*/
865: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
866: {
867: DM dm;
868: PetscInt pdim; /* Dimension of FE space P */
869: PetscInt dim; /* Spatial dimension */
870: PetscInt comp; /* Field components */
871: PetscErrorCode ierr;
874: if (!npoints) {
875: if (B) *B = NULL;
876: if (D) *D = NULL;
877: if (H) *H = NULL;
878: return(0);
879: }
885: PetscDualSpaceGetDM(fem->dualSpace, &dm);
886: DMGetDimension(dm, &dim);
887: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
888: PetscFEGetNumComponents(fem, &comp);
889: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
890: if (!dim) {
891: if (D) *D = NULL;
892: if (H) *H = NULL;
893: } else {
894: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
895: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
896: }
897: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
898: return(0);
899: }
901: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
902: {
903: DM dm;
908: PetscDualSpaceGetDM(fem->dualSpace, &dm);
909: if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
910: if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
911: if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
912: return(0);
913: }
915: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
916: {
917: PetscSpace bsp, bsubsp;
918: PetscDualSpace dsp, dsubsp;
919: PetscInt dim, depth, numComp, i, j, coneSize, order;
920: PetscFEType type;
921: DM dm;
922: DMLabel label;
923: PetscReal *xi, *v, *J, detJ;
924: const char *name;
925: PetscQuadrature origin, fullQuad, subQuad;
931: PetscFEGetBasisSpace(fe,&bsp);
932: PetscFEGetDualSpace(fe,&dsp);
933: PetscDualSpaceGetDM(dsp,&dm);
934: DMGetDimension(dm,&dim);
935: DMPlexGetDepthLabel(dm,&label);
936: DMLabelGetValue(label,refPoint,&depth);
937: PetscCalloc1(depth,&xi);
938: PetscMalloc1(dim,&v);
939: PetscMalloc1(dim*dim,&J);
940: for (i = 0; i < depth; i++) xi[i] = 0.;
941: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
942: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
943: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
944: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
945: for (i = 1; i < dim; i++) {
946: for (j = 0; j < depth; j++) {
947: J[i * depth + j] = J[i * dim + j];
948: }
949: }
950: PetscQuadratureDestroy(&origin);
951: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
952: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
953: PetscSpaceSetUp(bsubsp);
954: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
955: PetscFEGetType(fe,&type);
956: PetscFESetType(*trFE,type);
957: PetscFEGetNumComponents(fe,&numComp);
958: PetscFESetNumComponents(*trFE,numComp);
959: PetscFESetBasisSpace(*trFE,bsubsp);
960: PetscFESetDualSpace(*trFE,dsubsp);
961: PetscObjectGetName((PetscObject) fe, &name);
962: if (name) {PetscFESetName(*trFE, name);}
963: PetscFEGetQuadrature(fe,&fullQuad);
964: PetscQuadratureGetOrder(fullQuad,&order);
965: DMPlexGetConeSize(dm,refPoint,&coneSize);
966: if (coneSize == 2 * depth) {
967: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
968: } else {
969: PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
970: }
971: PetscFESetQuadrature(*trFE,subQuad);
972: PetscFESetUp(*trFE);
973: PetscQuadratureDestroy(&subQuad);
974: PetscSpaceDestroy(&bsubsp);
975: return(0);
976: }
978: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
979: {
980: PetscInt hStart, hEnd;
981: PetscDualSpace dsp;
982: DM dm;
988: *trFE = NULL;
989: PetscFEGetDualSpace(fe,&dsp);
990: PetscDualSpaceGetDM(dsp,&dm);
991: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
992: if (hEnd <= hStart) return(0);
993: PetscFECreatePointTrace(fe,hStart,trFE);
994: return(0);
995: }
998: /*@
999: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1001: Not collective
1003: Input Parameter:
1004: . fe - The PetscFE
1006: Output Parameter:
1007: . dim - The dimension
1009: Level: intermediate
1011: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1012: @*/
1013: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1014: {
1020: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1021: return(0);
1022: }
1024: /*
1025: Purpose: Compute element vector for chunk of elements
1027: Input:
1028: Sizes:
1029: Ne: number of elements
1030: Nf: number of fields
1031: PetscFE
1032: dim: spatial dimension
1033: Nb: number of basis functions
1034: Nc: number of field components
1035: PetscQuadrature
1036: Nq: number of quadrature points
1038: Geometry:
1039: PetscFEGeom[Ne] possibly *Nq
1040: PetscReal v0s[dim]
1041: PetscReal n[dim]
1042: PetscReal jacobians[dim*dim]
1043: PetscReal jacobianInverses[dim*dim]
1044: PetscReal jacobianDeterminants
1045: FEM:
1046: PetscFE
1047: PetscQuadrature
1048: PetscReal quadPoints[Nq*dim]
1049: PetscReal quadWeights[Nq]
1050: PetscReal basis[Nq*Nb*Nc]
1051: PetscReal basisDer[Nq*Nb*Nc*dim]
1052: PetscScalar coefficients[Ne*Nb*Nc]
1053: PetscScalar elemVec[Ne*Nb*Nc]
1055: Problem:
1056: PetscInt f: the active field
1057: f0, f1
1059: Work Space:
1060: PetscFE
1061: PetscScalar f0[Nq*dim];
1062: PetscScalar f1[Nq*dim*dim];
1063: PetscScalar u[Nc];
1064: PetscScalar gradU[Nc*dim];
1065: PetscReal x[dim];
1066: PetscScalar realSpaceDer[dim];
1068: Purpose: Compute element vector for N_cb batches of elements
1070: Input:
1071: Sizes:
1072: N_cb: Number of serial cell batches
1074: Geometry:
1075: PetscReal v0s[Ne*dim]
1076: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1077: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1078: PetscReal jacobianDeterminants[Ne] possibly *Nq
1079: FEM:
1080: static PetscReal quadPoints[Nq*dim]
1081: static PetscReal quadWeights[Nq]
1082: static PetscReal basis[Nq*Nb*Nc]
1083: static PetscReal basisDer[Nq*Nb*Nc*dim]
1084: PetscScalar coefficients[Ne*Nb*Nc]
1085: PetscScalar elemVec[Ne*Nb*Nc]
1087: ex62.c:
1088: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1089: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1090: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1091: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1093: ex52.c:
1094: 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)
1095: 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)
1097: ex52_integrateElement.cu
1098: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1100: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1101: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1102: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1104: ex52_integrateElementOpenCL.c:
1105: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1106: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1107: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1109: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1110: */
1112: /*@C
1113: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1115: Not collective
1117: Input Parameters:
1118: + fem - The PetscFE object for the field being integrated
1119: . prob - The PetscDS specifying the discretizations and continuum functions
1120: . field - The field being integrated
1121: . Ne - The number of elements in the chunk
1122: . cgeom - The cell geometry for each cell in the chunk
1123: . coefficients - The array of FEM basis coefficients for the elements
1124: . probAux - The PetscDS specifying the auxiliary discretizations
1125: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1127: Output Parameter
1128: . integral - the integral for this field
1130: Level: developer
1132: .seealso: PetscFEIntegrateResidual()
1133: @*/
1134: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1135: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1136: {
1142: if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1143: return(0);
1144: }
1146: /*@C
1147: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1149: Not collective
1151: Input Parameters:
1152: + fem - The PetscFE object for the field being integrated
1153: . prob - The PetscDS specifying the discretizations and continuum functions
1154: . field - The field being integrated
1155: . obj_func - The function to be integrated
1156: . Ne - The number of elements in the chunk
1157: . fgeom - The face geometry for each face in the chunk
1158: . coefficients - The array of FEM basis coefficients for the elements
1159: . probAux - The PetscDS specifying the auxiliary discretizations
1160: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1162: Output Parameter
1163: . integral - the integral for this field
1165: Level: developer
1167: .seealso: PetscFEIntegrateResidual()
1168: @*/
1169: PetscErrorCode PetscFEIntegrateBd(PetscFE fem, PetscDS prob, PetscInt field,
1170: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1171: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1172: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1173: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1174: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1175: {
1181: if (fem->ops->integratebd) {(*fem->ops->integratebd)(fem, prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1182: return(0);
1183: }
1185: /*@C
1186: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1188: Not collective
1190: Input Parameters:
1191: + fem - The PetscFE object for the field being integrated
1192: . prob - The PetscDS specifying the discretizations and continuum functions
1193: . field - The field being integrated
1194: . Ne - The number of elements in the chunk
1195: . cgeom - The cell geometry for each cell in the chunk
1196: . coefficients - The array of FEM basis coefficients for the elements
1197: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1198: . probAux - The PetscDS specifying the auxiliary discretizations
1199: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1200: - t - The time
1202: Output Parameter
1203: . elemVec - the element residual vectors from each element
1205: Note:
1206: $ Loop over batch of elements (e):
1207: $ Loop over quadrature points (q):
1208: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1209: $ Call f_0 and f_1
1210: $ Loop over element vector entries (f,fc --> i):
1211: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1213: Level: developer
1215: .seealso: PetscFEIntegrateResidual()
1216: @*/
1217: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1218: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1219: {
1225: if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1226: return(0);
1227: }
1229: /*@C
1230: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1232: Not collective
1234: Input Parameters:
1235: + fem - The PetscFE object for the field being integrated
1236: . prob - The PetscDS specifying the discretizations and continuum functions
1237: . field - The field being integrated
1238: . Ne - The number of elements in the chunk
1239: . fgeom - The face geometry for each cell in the chunk
1240: . coefficients - The array of FEM basis coefficients for the elements
1241: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1242: . probAux - The PetscDS specifying the auxiliary discretizations
1243: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1244: - t - The time
1246: Output Parameter
1247: . elemVec - the element residual vectors from each element
1249: Level: developer
1251: .seealso: PetscFEIntegrateResidual()
1252: @*/
1253: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1254: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1255: {
1260: if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1261: return(0);
1262: }
1264: /*@C
1265: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1267: Not collective
1269: Input Parameters:
1270: + fem - The PetscFE object for the field being integrated
1271: . prob - The PetscDS specifying the discretizations and continuum functions
1272: . jtype - The type of matrix pointwise functions that should be used
1273: . fieldI - The test field being integrated
1274: . fieldJ - The basis field being integrated
1275: . Ne - The number of elements in the chunk
1276: . cgeom - The cell geometry for each cell in the chunk
1277: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1278: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1279: . probAux - The PetscDS specifying the auxiliary discretizations
1280: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1281: . t - The time
1282: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1284: Output Parameter
1285: . elemMat - the element matrices for the Jacobian from each element
1287: Note:
1288: $ Loop over batch of elements (e):
1289: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1290: $ Loop over quadrature points (q):
1291: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1292: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1293: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1294: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1295: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1296: Level: developer
1298: .seealso: PetscFEIntegrateResidual()
1299: @*/
1300: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1301: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1302: {
1307: if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1308: return(0);
1309: }
1311: /*@C
1312: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1314: Not collective
1316: Input Parameters:
1317: + fem = The PetscFE object for the field being integrated
1318: . prob - The PetscDS specifying the discretizations and continuum functions
1319: . fieldI - The test field being integrated
1320: . fieldJ - The basis field being integrated
1321: . Ne - The number of elements in the chunk
1322: . fgeom - The face geometry for each cell in the chunk
1323: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1324: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1325: . probAux - The PetscDS specifying the auxiliary discretizations
1326: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1327: . t - The time
1328: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1330: Output Parameter
1331: . elemMat - the element matrices for the Jacobian from each element
1333: Note:
1334: $ Loop over batch of elements (e):
1335: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1336: $ Loop over quadrature points (q):
1337: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1338: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1339: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1340: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1341: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1342: Level: developer
1344: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1345: @*/
1346: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1347: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1348: {
1353: if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1354: return(0);
1355: }
1357: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1358: {
1359: PetscSpace P, subP;
1360: PetscDualSpace Q, subQ;
1361: PetscQuadrature subq;
1362: PetscFEType fetype;
1363: PetscInt dim, Nc;
1364: PetscErrorCode ierr;
1369: if (height == 0) {
1370: *subfe = fe;
1371: return(0);
1372: }
1373: PetscFEGetBasisSpace(fe, &P);
1374: PetscFEGetDualSpace(fe, &Q);
1375: PetscFEGetNumComponents(fe, &Nc);
1376: PetscFEGetFaceQuadrature(fe, &subq);
1377: PetscDualSpaceGetDimension(Q, &dim);
1378: 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);}
1379: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1380: if (height <= dim) {
1381: if (!fe->subspaces[height-1]) {
1382: PetscFE sub;
1383: const char *name;
1385: PetscSpaceGetHeightSubspace(P, height, &subP);
1386: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1387: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1388: PetscObjectGetName((PetscObject) fe, &name);
1389: PetscObjectSetName((PetscObject) sub, name);
1390: PetscFEGetType(fe, &fetype);
1391: PetscFESetType(sub, fetype);
1392: PetscFESetBasisSpace(sub, subP);
1393: PetscFESetDualSpace(sub, subQ);
1394: PetscFESetNumComponents(sub, Nc);
1395: PetscFESetUp(sub);
1396: PetscFESetQuadrature(sub, subq);
1397: fe->subspaces[height-1] = sub;
1398: }
1399: *subfe = fe->subspaces[height-1];
1400: } else {
1401: *subfe = NULL;
1402: }
1403: return(0);
1404: }
1406: /*@
1407: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1408: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1409: sparsity). It is also used to create an interpolation between regularly refined meshes.
1411: Collective on PetscFE
1413: Input Parameter:
1414: . fe - The initial PetscFE
1416: Output Parameter:
1417: . feRef - The refined PetscFE
1419: Level: developer
1421: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1422: @*/
1423: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1424: {
1425: PetscSpace P, Pref;
1426: PetscDualSpace Q, Qref;
1427: DM K, Kref;
1428: PetscQuadrature q, qref;
1429: const PetscReal *v0, *jac;
1430: PetscInt numComp, numSubelements;
1431: PetscErrorCode ierr;
1434: PetscFEGetBasisSpace(fe, &P);
1435: PetscFEGetDualSpace(fe, &Q);
1436: PetscFEGetQuadrature(fe, &q);
1437: PetscDualSpaceGetDM(Q, &K);
1438: /* Create space */
1439: PetscObjectReference((PetscObject) P);
1440: Pref = P;
1441: /* Create dual space */
1442: PetscDualSpaceDuplicate(Q, &Qref);
1443: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1444: PetscDualSpaceSetDM(Qref, Kref);
1445: DMDestroy(&Kref);
1446: PetscDualSpaceSetUp(Qref);
1447: /* Create element */
1448: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1449: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1450: PetscFESetBasisSpace(*feRef, Pref);
1451: PetscFESetDualSpace(*feRef, Qref);
1452: PetscFEGetNumComponents(fe, &numComp);
1453: PetscFESetNumComponents(*feRef, numComp);
1454: PetscFESetUp(*feRef);
1455: PetscSpaceDestroy(&Pref);
1456: PetscDualSpaceDestroy(&Qref);
1457: /* Create quadrature */
1458: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1459: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1460: PetscFESetQuadrature(*feRef, qref);
1461: PetscQuadratureDestroy(&qref);
1462: return(0);
1463: }
1465: /*@C
1466: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1468: Collective on DM
1470: Input Parameters:
1471: + comm - The MPI comm
1472: . dim - The spatial dimension
1473: . Nc - The number of components
1474: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1475: . prefix - The options prefix, or NULL
1476: - qorder - The quadrature order
1478: Output Parameter:
1479: . fem - The PetscFE object
1481: Level: beginner
1483: .keywords: PetscFE, finite element
1484: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1485: @*/
1486: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1487: {
1488: PetscQuadrature q, fq;
1489: DM K;
1490: PetscSpace P;
1491: PetscDualSpace Q;
1492: PetscInt order, quadPointsPerEdge;
1493: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1494: PetscErrorCode ierr;
1497: /* Create space */
1498: PetscSpaceCreate(comm, &P);
1499: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1500: PetscSpacePolynomialSetTensor(P, tensor);
1501: PetscSpaceSetNumComponents(P, Nc);
1502: PetscSpaceSetNumVariables(P, dim);
1503: PetscSpaceSetFromOptions(P);
1504: PetscSpaceSetUp(P);
1505: PetscSpaceGetDegree(P, &order, NULL);
1506: PetscSpacePolynomialGetTensor(P, &tensor);
1507: /* Create dual space */
1508: PetscDualSpaceCreate(comm, &Q);
1509: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1510: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1511: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1512: PetscDualSpaceSetDM(Q, K);
1513: DMDestroy(&K);
1514: PetscDualSpaceSetNumComponents(Q, Nc);
1515: PetscDualSpaceSetOrder(Q, order);
1516: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1517: PetscDualSpaceSetFromOptions(Q);
1518: PetscDualSpaceSetUp(Q);
1519: /* Create element */
1520: PetscFECreate(comm, fem);
1521: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1522: PetscFESetFromOptions(*fem);
1523: PetscFESetBasisSpace(*fem, P);
1524: PetscFESetDualSpace(*fem, Q);
1525: PetscFESetNumComponents(*fem, Nc);
1526: PetscFESetUp(*fem);
1527: PetscSpaceDestroy(&P);
1528: PetscDualSpaceDestroy(&Q);
1529: /* Create quadrature (with specified order if given) */
1530: qorder = qorder >= 0 ? qorder : order;
1531: PetscObjectOptionsBegin((PetscObject)*fem);
1532: PetscOptionsInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadture points per edge","PetscFECreateDefault",qorder,&qorder,NULL);
1533: PetscOptionsEnd();
1534: quadPointsPerEdge = PetscMax(qorder + 1,1);
1535: if (isSimplex) {
1536: PetscDTGaussJacobiQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1537: PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1538: } else {
1539: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1540: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1541: }
1542: PetscFESetQuadrature(*fem, q);
1543: PetscFESetFaceQuadrature(*fem, fq);
1544: PetscQuadratureDestroy(&q);
1545: PetscQuadratureDestroy(&fq);
1546: return(0);
1547: }
1549: /*@C
1550: PetscFESetName - Names the FE and its subobjects
1552: Not collective
1554: Input Parameters:
1555: + fe - The PetscFE
1556: - name - The name
1558: Level: beginner
1560: .keywords: PetscFE, finite element
1561: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1562: @*/
1563: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1564: {
1565: PetscSpace P;
1566: PetscDualSpace Q;
1570: PetscFEGetBasisSpace(fe, &P);
1571: PetscFEGetDualSpace(fe, &Q);
1572: PetscObjectSetName((PetscObject) fe, name);
1573: PetscObjectSetName((PetscObject) P, name);
1574: PetscObjectSetName((PetscObject) Q, name);
1575: return(0);
1576: }