Actual source code: fe.c
petsc-3.10.5 2019-03-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: - v - the viewer
178: Level: developer
180: .seealso PetscFEDestroy()
181: @*/
182: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
183: {
188: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);}
189: if (fem->ops->view) {(*fem->ops->view)(fem, v);}
190: return(0);
191: }
193: /*@
194: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
196: Collective on PetscFE
198: Input Parameter:
199: . fem - the PetscFE object to set options for
201: Options Database:
202: . -petscfe_num_blocks the number of cell blocks to integrate concurrently
203: . -petscfe_num_batches the number of cell batches to integrate serially
205: Level: developer
207: .seealso PetscFEView()
208: @*/
209: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
210: {
211: const char *defaultType;
212: char name[256];
213: PetscBool flg;
218: if (!((PetscObject) fem)->type_name) {
219: defaultType = PETSCFEBASIC;
220: } else {
221: defaultType = ((PetscObject) fem)->type_name;
222: }
223: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
225: PetscObjectOptionsBegin((PetscObject) fem);
226: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
227: if (flg) {
228: PetscFESetType(fem, name);
229: } else if (!((PetscObject) fem)->type_name) {
230: PetscFESetType(fem, defaultType);
231: }
232: PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
233: PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
234: if (fem->ops->setfromoptions) {
235: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
236: }
237: /* process any options handlers added with PetscObjectAddOptionsHandler() */
238: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
239: PetscOptionsEnd();
240: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
241: return(0);
242: }
244: /*@C
245: PetscFESetUp - Construct data structures for the PetscFE
247: Collective on PetscFE
249: Input Parameter:
250: . fem - the PetscFE object to setup
252: Level: developer
254: .seealso PetscFEView(), PetscFEDestroy()
255: @*/
256: PetscErrorCode PetscFESetUp(PetscFE fem)
257: {
262: if (fem->setupcalled) return(0);
263: fem->setupcalled = PETSC_TRUE;
264: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
265: return(0);
266: }
268: /*@
269: PetscFEDestroy - Destroys a PetscFE object
271: Collective on PetscFE
273: Input Parameter:
274: . fem - the PetscFE object to destroy
276: Level: developer
278: .seealso PetscFEView()
279: @*/
280: PetscErrorCode PetscFEDestroy(PetscFE *fem)
281: {
285: if (!*fem) return(0);
288: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
289: ((PetscObject) (*fem))->refct = 0;
291: if ((*fem)->subspaces) {
292: PetscInt dim, d;
294: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
295: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
296: }
297: PetscFree((*fem)->subspaces);
298: PetscFree((*fem)->invV);
299: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
300: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->Bf, &(*fem)->Df, NULL /*&(*fem)->Hf*/);
301: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
302: PetscSpaceDestroy(&(*fem)->basisSpace);
303: PetscDualSpaceDestroy(&(*fem)->dualSpace);
304: PetscQuadratureDestroy(&(*fem)->quadrature);
305: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
307: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
308: PetscHeaderDestroy(fem);
309: return(0);
310: }
312: /*@
313: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
315: Collective on MPI_Comm
317: Input Parameter:
318: . comm - The communicator for the PetscFE object
320: Output Parameter:
321: . fem - The PetscFE object
323: Level: beginner
325: .seealso: PetscFESetType(), PETSCFEGALERKIN
326: @*/
327: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
328: {
329: PetscFE f;
334: PetscCitationsRegister(FECitation,&FEcite);
335: *fem = NULL;
336: PetscFEInitializePackage();
338: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
340: f->basisSpace = NULL;
341: f->dualSpace = NULL;
342: f->numComponents = 1;
343: f->subspaces = NULL;
344: f->invV = NULL;
345: f->B = NULL;
346: f->D = NULL;
347: f->H = NULL;
348: f->Bf = NULL;
349: f->Df = NULL;
350: f->Hf = NULL;
351: PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
352: PetscMemzero(&f->faceQuadrature, sizeof(PetscQuadrature));
353: f->blockSize = 0;
354: f->numBlocks = 1;
355: f->batchSize = 0;
356: f->numBatches = 1;
358: *fem = f;
359: return(0);
360: }
362: /*@
363: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
365: Not collective
367: Input Parameter:
368: . fem - The PetscFE object
370: Output Parameter:
371: . dim - The spatial dimension
373: Level: intermediate
375: .seealso: PetscFECreate()
376: @*/
377: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
378: {
379: DM dm;
385: PetscDualSpaceGetDM(fem->dualSpace, &dm);
386: DMGetDimension(dm, dim);
387: return(0);
388: }
390: /*@
391: PetscFESetNumComponents - Sets the number of components in the element
393: Not collective
395: Input Parameters:
396: + fem - The PetscFE object
397: - comp - The number of field components
399: Level: intermediate
401: .seealso: PetscFECreate()
402: @*/
403: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
404: {
407: fem->numComponents = comp;
408: return(0);
409: }
411: /*@
412: PetscFEGetNumComponents - Returns the number of components in the element
414: Not collective
416: Input Parameter:
417: . fem - The PetscFE object
419: Output Parameter:
420: . comp - The number of field components
422: Level: intermediate
424: .seealso: PetscFECreate()
425: @*/
426: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
427: {
431: *comp = fem->numComponents;
432: return(0);
433: }
435: /*@
436: PetscFESetTileSizes - Sets the tile sizes for evaluation
438: Not collective
440: Input Parameters:
441: + fem - The PetscFE object
442: . blockSize - The number of elements in a block
443: . numBlocks - The number of blocks in a batch
444: . batchSize - The number of elements in a batch
445: - numBatches - The number of batches in a chunk
447: Level: intermediate
449: .seealso: PetscFECreate()
450: @*/
451: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
452: {
455: fem->blockSize = blockSize;
456: fem->numBlocks = numBlocks;
457: fem->batchSize = batchSize;
458: fem->numBatches = numBatches;
459: return(0);
460: }
462: /*@
463: PetscFEGetTileSizes - Returns the tile sizes for evaluation
465: Not collective
467: Input Parameter:
468: . fem - The PetscFE object
470: Output Parameters:
471: + blockSize - The number of elements in a block
472: . numBlocks - The number of blocks in a batch
473: . batchSize - The number of elements in a batch
474: - numBatches - The number of batches in a chunk
476: Level: intermediate
478: .seealso: PetscFECreate()
479: @*/
480: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
481: {
488: if (blockSize) *blockSize = fem->blockSize;
489: if (numBlocks) *numBlocks = fem->numBlocks;
490: if (batchSize) *batchSize = fem->batchSize;
491: if (numBatches) *numBatches = fem->numBatches;
492: return(0);
493: }
495: /*@
496: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
498: Not collective
500: Input Parameter:
501: . fem - The PetscFE object
503: Output Parameter:
504: . sp - The PetscSpace object
506: Level: intermediate
508: .seealso: PetscFECreate()
509: @*/
510: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
511: {
515: *sp = fem->basisSpace;
516: return(0);
517: }
519: /*@
520: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
522: Not collective
524: Input Parameters:
525: + fem - The PetscFE object
526: - sp - The PetscSpace object
528: Level: intermediate
530: .seealso: PetscFECreate()
531: @*/
532: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
533: {
539: PetscSpaceDestroy(&fem->basisSpace);
540: fem->basisSpace = sp;
541: PetscObjectReference((PetscObject) fem->basisSpace);
542: return(0);
543: }
545: /*@
546: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
548: Not collective
550: Input Parameter:
551: . fem - The PetscFE object
553: Output Parameter:
554: . sp - The PetscDualSpace object
556: Level: intermediate
558: .seealso: PetscFECreate()
559: @*/
560: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
561: {
565: *sp = fem->dualSpace;
566: return(0);
567: }
569: /*@
570: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
572: Not collective
574: Input Parameters:
575: + fem - The PetscFE object
576: - sp - The PetscDualSpace object
578: Level: intermediate
580: .seealso: PetscFECreate()
581: @*/
582: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
583: {
589: PetscDualSpaceDestroy(&fem->dualSpace);
590: fem->dualSpace = sp;
591: PetscObjectReference((PetscObject) fem->dualSpace);
592: return(0);
593: }
595: /*@
596: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
598: Not collective
600: Input Parameter:
601: . fem - The PetscFE object
603: Output Parameter:
604: . q - The PetscQuadrature object
606: Level: intermediate
608: .seealso: PetscFECreate()
609: @*/
610: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
611: {
615: *q = fem->quadrature;
616: return(0);
617: }
619: /*@
620: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
622: Not collective
624: Input Parameters:
625: + fem - The PetscFE object
626: - q - The PetscQuadrature object
628: Level: intermediate
630: .seealso: PetscFECreate()
631: @*/
632: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
633: {
634: PetscInt Nc, qNc;
639: PetscFEGetNumComponents(fem, &Nc);
640: PetscQuadratureGetNumComponents(q, &qNc);
641: 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);
642: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
643: PetscQuadratureDestroy(&fem->quadrature);
644: fem->quadrature = q;
645: PetscObjectReference((PetscObject) q);
646: return(0);
647: }
649: /*@
650: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
652: Not collective
654: Input Parameter:
655: . fem - The PetscFE object
657: Output Parameter:
658: . q - The PetscQuadrature object
660: Level: intermediate
662: .seealso: PetscFECreate()
663: @*/
664: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
665: {
669: *q = fem->faceQuadrature;
670: return(0);
671: }
673: /*@
674: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
676: Not collective
678: Input Parameters:
679: + fem - The PetscFE object
680: - q - The PetscQuadrature object
682: Level: intermediate
684: .seealso: PetscFECreate()
685: @*/
686: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
687: {
692: PetscFERestoreTabulation(fem, 0, NULL, &fem->Bf, &fem->Df, NULL /*&(*fem)->Hf*/);
693: PetscQuadratureDestroy(&fem->faceQuadrature);
694: fem->faceQuadrature = q;
695: PetscObjectReference((PetscObject) q);
696: return(0);
697: }
699: /*@C
700: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
702: Not collective
704: Input Parameter:
705: . fem - The PetscFE object
707: Output Parameter:
708: . numDof - Array with the number of dofs per dimension
710: Level: intermediate
712: .seealso: PetscFECreate()
713: @*/
714: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
715: {
721: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
722: return(0);
723: }
725: /*@C
726: PetscFEGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points
728: Not collective
730: Input Parameter:
731: . fem - The PetscFE object
733: Output Parameters:
734: + B - The basis function values at quadrature points
735: . D - The basis function derivatives at quadrature points
736: - H - The basis function second derivatives at quadrature points
738: Note:
739: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
740: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
741: $ 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
743: Level: intermediate
745: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation()
746: @*/
747: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
748: {
749: PetscInt npoints;
750: const PetscReal *points;
751: PetscErrorCode ierr;
758: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
759: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
760: if (B) *B = fem->B;
761: if (D) *D = fem->D;
762: if (H) *H = fem->H;
763: return(0);
764: }
766: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **Bf, PetscReal **Df, PetscReal **Hf)
767: {
768: PetscErrorCode ierr;
775: if (!fem->Bf) {
776: const PetscReal xi0[3] = {-1., -1., -1.};
777: PetscReal v0[3], J[9], detJ;
778: PetscQuadrature fq;
779: PetscDualSpace sp;
780: DM dm;
781: const PetscInt *faces;
782: PetscInt dim, numFaces, f, npoints, q;
783: const PetscReal *points;
784: PetscReal *facePoints;
786: PetscFEGetDualSpace(fem, &sp);
787: PetscDualSpaceGetDM(sp, &dm);
788: DMGetDimension(dm, &dim);
789: DMPlexGetConeSize(dm, 0, &numFaces);
790: DMPlexGetCone(dm, 0, &faces);
791: PetscFEGetFaceQuadrature(fem, &fq);
792: if (fq) {
793: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
794: PetscMalloc1(numFaces*npoints*dim, &facePoints);
795: for (f = 0; f < numFaces; ++f) {
796: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
797: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
798: }
799: PetscFEGetTabulation(fem, numFaces*npoints, facePoints, &fem->Bf, &fem->Df, NULL/*&fem->Hf*/);
800: PetscFree(facePoints);
801: }
802: }
803: if (Bf) *Bf = fem->Bf;
804: if (Df) *Df = fem->Df;
805: if (Hf) *Hf = fem->Hf;
806: return(0);
807: }
809: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscReal **F)
810: {
811: PetscErrorCode ierr;
816: if (!fem->F) {
817: PetscDualSpace sp;
818: DM dm;
819: const PetscInt *cone;
820: PetscReal *centroids;
821: PetscInt dim, numFaces, f;
823: PetscFEGetDualSpace(fem, &sp);
824: PetscDualSpaceGetDM(sp, &dm);
825: DMGetDimension(dm, &dim);
826: DMPlexGetConeSize(dm, 0, &numFaces);
827: DMPlexGetCone(dm, 0, &cone);
828: PetscMalloc1(numFaces*dim, ¢roids);
829: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
830: PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
831: PetscFree(centroids);
832: }
833: *F = fem->F;
834: return(0);
835: }
837: /*@C
838: PetscFEGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
840: Not collective
842: Input Parameters:
843: + fem - The PetscFE object
844: . npoints - The number of tabulation points
845: - points - The tabulation point coordinates
847: Output Parameters:
848: + B - The basis function values at tabulation points
849: . D - The basis function derivatives at tabulation points
850: - H - The basis function second derivatives at tabulation points
852: Note:
853: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
854: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
855: $ 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
857: Level: intermediate
859: .seealso: PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
860: @*/
861: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
862: {
863: DM dm;
864: PetscInt pdim; /* Dimension of FE space P */
865: PetscInt dim; /* Spatial dimension */
866: PetscInt comp; /* Field components */
867: PetscErrorCode ierr;
870: if (!npoints) {
871: if (B) *B = NULL;
872: if (D) *D = NULL;
873: if (H) *H = NULL;
874: return(0);
875: }
881: PetscDualSpaceGetDM(fem->dualSpace, &dm);
882: DMGetDimension(dm, &dim);
883: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
884: PetscFEGetNumComponents(fem, &comp);
885: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, MPIU_REAL, B);}
886: if (!dim) {
887: if (D) *D = NULL;
888: if (H) *H = NULL;
889: } else {
890: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, MPIU_REAL, D);}
891: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, MPIU_REAL, H);}
892: }
893: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
894: return(0);
895: }
897: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
898: {
899: DM dm;
904: PetscDualSpaceGetDM(fem->dualSpace, &dm);
905: if (B && *B) {DMRestoreWorkArray(dm, 0, MPIU_REAL, B);}
906: if (D && *D) {DMRestoreWorkArray(dm, 0, MPIU_REAL, D);}
907: if (H && *H) {DMRestoreWorkArray(dm, 0, MPIU_REAL, H);}
908: return(0);
909: }
911: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
912: {
913: PetscSpace bsp, bsubsp;
914: PetscDualSpace dsp, dsubsp;
915: PetscInt dim, depth, numComp, i, j, coneSize, order;
916: PetscFEType type;
917: DM dm;
918: DMLabel label;
919: PetscReal *xi, *v, *J, detJ;
920: PetscQuadrature origin, fullQuad, subQuad;
926: PetscFEGetBasisSpace(fe,&bsp);
927: PetscFEGetDualSpace(fe,&dsp);
928: PetscDualSpaceGetDM(dsp,&dm);
929: DMGetDimension(dm,&dim);
930: DMPlexGetDepthLabel(dm,&label);
931: DMLabelGetValue(label,refPoint,&depth);
932: PetscCalloc1(depth,&xi);
933: PetscMalloc1(dim,&v);
934: PetscMalloc1(dim*dim,&J);
935: for (i = 0; i < depth; i++) xi[i] = 0.;
936: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
937: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
938: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
939: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
940: for (i = 1; i < dim; i++) {
941: for (j = 0; j < depth; j++) {
942: J[i * depth + j] = J[i * dim + j];
943: }
944: }
945: PetscQuadratureDestroy(&origin);
946: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
947: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
948: PetscSpaceSetUp(bsubsp);
949: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
950: PetscFEGetType(fe,&type);
951: PetscFESetType(*trFE,type);
952: PetscFEGetNumComponents(fe,&numComp);
953: PetscFESetNumComponents(*trFE,numComp);
954: PetscFESetBasisSpace(*trFE,bsubsp);
955: PetscFESetDualSpace(*trFE,dsubsp);
956: PetscFEGetQuadrature(fe,&fullQuad);
957: PetscQuadratureGetOrder(fullQuad,&order);
958: DMPlexGetConeSize(dm,refPoint,&coneSize);
959: if (coneSize == 2 * depth) {
960: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
961: } else {
962: PetscDTGaussJacobiQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
963: }
964: PetscFESetQuadrature(*trFE,subQuad);
965: PetscFESetUp(*trFE);
966: PetscQuadratureDestroy(&subQuad);
967: PetscSpaceDestroy(&bsubsp);
968: return(0);
969: }
971: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
972: {
973: PetscInt hStart, hEnd;
974: PetscDualSpace dsp;
975: DM dm;
981: *trFE = NULL;
982: PetscFEGetDualSpace(fe,&dsp);
983: PetscDualSpaceGetDM(dsp,&dm);
984: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
985: if (hEnd <= hStart) return(0);
986: PetscFECreatePointTrace(fe,hStart,trFE);
987: return(0);
988: }
991: /*@
992: PetscFEGetDimension - Get the dimension of the finite element space on a cell
994: Not collective
996: Input Parameter:
997: . fe - The PetscFE
999: Output Parameter:
1000: . dim - The dimension
1002: Level: intermediate
1004: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1005: @*/
1006: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1007: {
1013: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1014: return(0);
1015: }
1017: /*
1018: Purpose: Compute element vector for chunk of elements
1020: Input:
1021: Sizes:
1022: Ne: number of elements
1023: Nf: number of fields
1024: PetscFE
1025: dim: spatial dimension
1026: Nb: number of basis functions
1027: Nc: number of field components
1028: PetscQuadrature
1029: Nq: number of quadrature points
1031: Geometry:
1032: PetscFEGeom[Ne] possibly *Nq
1033: PetscReal v0s[dim]
1034: PetscReal n[dim]
1035: PetscReal jacobians[dim*dim]
1036: PetscReal jacobianInverses[dim*dim]
1037: PetscReal jacobianDeterminants
1038: FEM:
1039: PetscFE
1040: PetscQuadrature
1041: PetscReal quadPoints[Nq*dim]
1042: PetscReal quadWeights[Nq]
1043: PetscReal basis[Nq*Nb*Nc]
1044: PetscReal basisDer[Nq*Nb*Nc*dim]
1045: PetscScalar coefficients[Ne*Nb*Nc]
1046: PetscScalar elemVec[Ne*Nb*Nc]
1048: Problem:
1049: PetscInt f: the active field
1050: f0, f1
1052: Work Space:
1053: PetscFE
1054: PetscScalar f0[Nq*dim];
1055: PetscScalar f1[Nq*dim*dim];
1056: PetscScalar u[Nc];
1057: PetscScalar gradU[Nc*dim];
1058: PetscReal x[dim];
1059: PetscScalar realSpaceDer[dim];
1061: Purpose: Compute element vector for N_cb batches of elements
1063: Input:
1064: Sizes:
1065: N_cb: Number of serial cell batches
1067: Geometry:
1068: PetscReal v0s[Ne*dim]
1069: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1070: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1071: PetscReal jacobianDeterminants[Ne] possibly *Nq
1072: FEM:
1073: static PetscReal quadPoints[Nq*dim]
1074: static PetscReal quadWeights[Nq]
1075: static PetscReal basis[Nq*Nb*Nc]
1076: static PetscReal basisDer[Nq*Nb*Nc*dim]
1077: PetscScalar coefficients[Ne*Nb*Nc]
1078: PetscScalar elemVec[Ne*Nb*Nc]
1080: ex62.c:
1081: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1082: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1083: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1084: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1086: ex52.c:
1087: 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)
1088: 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)
1090: ex52_integrateElement.cu
1091: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1093: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1094: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1095: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1097: ex52_integrateElementOpenCL.c:
1098: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1099: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1100: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1102: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1103: */
1105: /*@C
1106: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1108: Not collective
1110: Input Parameters:
1111: + fem - The PetscFE object for the field being integrated
1112: . prob - The PetscDS specifying the discretizations and continuum functions
1113: . field - The field being integrated
1114: . Ne - The number of elements in the chunk
1115: . cgeom - The cell geometry for each cell in the chunk
1116: . coefficients - The array of FEM basis coefficients for the elements
1117: . probAux - The PetscDS specifying the auxiliary discretizations
1118: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1120: Output Parameter
1121: . integral - the integral for this field
1123: Level: developer
1125: .seealso: PetscFEIntegrateResidual()
1126: @*/
1127: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1128: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1129: {
1135: if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1136: return(0);
1137: }
1139: /*@C
1140: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1142: Not collective
1144: Input Parameters:
1145: + fem - The PetscFE object for the field being integrated
1146: . prob - The PetscDS specifying the discretizations and continuum functions
1147: . field - The field being integrated
1148: . obj_func - The function to be integrated
1149: . Ne - The number of elements in the chunk
1150: . fgeom - The face geometry for each face in the chunk
1151: . coefficients - The array of FEM basis coefficients for the elements
1152: . probAux - The PetscDS specifying the auxiliary discretizations
1153: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1155: Output Parameter
1156: . integral - the integral for this field
1158: Level: developer
1160: .seealso: PetscFEIntegrateResidual()
1161: @*/
1162: PetscErrorCode PetscFEIntegrateBd(PetscFE fem, PetscDS prob, PetscInt field,
1163: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1164: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1165: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1166: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1167: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1168: {
1174: if (fem->ops->integratebd) {(*fem->ops->integratebd)(fem, prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1175: return(0);
1176: }
1178: /*@C
1179: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1181: Not collective
1183: Input Parameters:
1184: + fem - The PetscFE object for the field being integrated
1185: . prob - The PetscDS specifying the discretizations and continuum functions
1186: . field - The field being integrated
1187: . Ne - The number of elements in the chunk
1188: . cgeom - The cell geometry for each cell in the chunk
1189: . coefficients - The array of FEM basis coefficients for the elements
1190: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1191: . probAux - The PetscDS specifying the auxiliary discretizations
1192: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1193: - t - The time
1195: Output Parameter
1196: . elemVec - the element residual vectors from each element
1198: Note:
1199: $ Loop over batch of elements (e):
1200: $ Loop over quadrature points (q):
1201: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1202: $ Call f_0 and f_1
1203: $ Loop over element vector entries (f,fc --> i):
1204: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1206: Level: developer
1208: .seealso: PetscFEIntegrateResidual()
1209: @*/
1210: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1211: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1212: {
1218: if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1219: return(0);
1220: }
1222: /*@C
1223: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1225: Not collective
1227: Input Parameters:
1228: + fem - The PetscFE object for the field being integrated
1229: . prob - The PetscDS specifying the discretizations and continuum functions
1230: . field - The field being integrated
1231: . Ne - The number of elements in the chunk
1232: . fgeom - The face geometry for each cell in the chunk
1233: . coefficients - The array of FEM basis coefficients for the elements
1234: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1235: . probAux - The PetscDS specifying the auxiliary discretizations
1236: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1237: - t - The time
1239: Output Parameter
1240: . elemVec - the element residual vectors from each element
1242: Level: developer
1244: .seealso: PetscFEIntegrateResidual()
1245: @*/
1246: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1247: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1248: {
1253: if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1254: return(0);
1255: }
1257: /*@C
1258: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1260: Not collective
1262: Input Parameters:
1263: + fem - The PetscFE object for the field being integrated
1264: . prob - The PetscDS specifying the discretizations and continuum functions
1265: . jtype - The type of matrix pointwise functions that should be used
1266: . fieldI - The test field being integrated
1267: . fieldJ - The basis field being integrated
1268: . Ne - The number of elements in the chunk
1269: . cgeom - The cell geometry for each cell in the chunk
1270: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1271: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1272: . probAux - The PetscDS specifying the auxiliary discretizations
1273: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1274: . t - The time
1275: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1277: Output Parameter
1278: . elemMat - the element matrices for the Jacobian from each element
1280: Note:
1281: $ Loop over batch of elements (e):
1282: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1283: $ Loop over quadrature points (q):
1284: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1285: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1286: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1287: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1288: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1289: Level: developer
1291: .seealso: PetscFEIntegrateResidual()
1292: @*/
1293: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1294: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1295: {
1300: if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1301: return(0);
1302: }
1304: /*@C
1305: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1307: Not collective
1309: Input Parameters:
1310: + fem = The PetscFE object for the field being integrated
1311: . prob - The PetscDS specifying the discretizations and continuum functions
1312: . fieldI - The test field being integrated
1313: . fieldJ - The basis field being integrated
1314: . Ne - The number of elements in the chunk
1315: . fgeom - The face geometry for each cell in the chunk
1316: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1317: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1318: . probAux - The PetscDS specifying the auxiliary discretizations
1319: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1320: . t - The time
1321: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1323: Output Parameter
1324: . elemMat - the element matrices for the Jacobian from each element
1326: Note:
1327: $ Loop over batch of elements (e):
1328: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1329: $ Loop over quadrature points (q):
1330: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1331: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1332: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1333: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1334: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1335: Level: developer
1337: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1338: @*/
1339: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1340: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1341: {
1346: if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1347: return(0);
1348: }
1350: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1351: {
1352: PetscSpace P, subP;
1353: PetscDualSpace Q, subQ;
1354: PetscQuadrature subq;
1355: PetscFEType fetype;
1356: PetscInt dim, Nc;
1357: PetscErrorCode ierr;
1362: if (height == 0) {
1363: *subfe = fe;
1364: return(0);
1365: }
1366: PetscFEGetBasisSpace(fe, &P);
1367: PetscFEGetDualSpace(fe, &Q);
1368: PetscFEGetNumComponents(fe, &Nc);
1369: PetscFEGetFaceQuadrature(fe, &subq);
1370: PetscDualSpaceGetDimension(Q, &dim);
1371: 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);}
1372: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1373: if (height <= dim) {
1374: if (!fe->subspaces[height-1]) {
1375: PetscFE sub;
1377: PetscSpaceGetHeightSubspace(P, height, &subP);
1378: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1379: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1380: PetscFEGetType(fe, &fetype);
1381: PetscFESetType(sub, fetype);
1382: PetscFESetBasisSpace(sub, subP);
1383: PetscFESetDualSpace(sub, subQ);
1384: PetscFESetNumComponents(sub, Nc);
1385: PetscFESetUp(sub);
1386: PetscFESetQuadrature(sub, subq);
1387: fe->subspaces[height-1] = sub;
1388: }
1389: *subfe = fe->subspaces[height-1];
1390: } else {
1391: *subfe = NULL;
1392: }
1393: return(0);
1394: }
1396: /*@
1397: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1398: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1399: sparsity). It is also used to create an interpolation between regularly refined meshes.
1401: Collective on PetscFE
1403: Input Parameter:
1404: . fe - The initial PetscFE
1406: Output Parameter:
1407: . feRef - The refined PetscFE
1409: Level: developer
1411: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1412: @*/
1413: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1414: {
1415: PetscSpace P, Pref;
1416: PetscDualSpace Q, Qref;
1417: DM K, Kref;
1418: PetscQuadrature q, qref;
1419: const PetscReal *v0, *jac;
1420: PetscInt numComp, numSubelements;
1421: PetscErrorCode ierr;
1424: PetscFEGetBasisSpace(fe, &P);
1425: PetscFEGetDualSpace(fe, &Q);
1426: PetscFEGetQuadrature(fe, &q);
1427: PetscDualSpaceGetDM(Q, &K);
1428: /* Create space */
1429: PetscObjectReference((PetscObject) P);
1430: Pref = P;
1431: /* Create dual space */
1432: PetscDualSpaceDuplicate(Q, &Qref);
1433: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1434: PetscDualSpaceSetDM(Qref, Kref);
1435: DMDestroy(&Kref);
1436: PetscDualSpaceSetUp(Qref);
1437: /* Create element */
1438: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1439: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1440: PetscFESetBasisSpace(*feRef, Pref);
1441: PetscFESetDualSpace(*feRef, Qref);
1442: PetscFEGetNumComponents(fe, &numComp);
1443: PetscFESetNumComponents(*feRef, numComp);
1444: PetscFESetUp(*feRef);
1445: PetscSpaceDestroy(&Pref);
1446: PetscDualSpaceDestroy(&Qref);
1447: /* Create quadrature */
1448: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1449: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1450: PetscFESetQuadrature(*feRef, qref);
1451: PetscQuadratureDestroy(&qref);
1452: return(0);
1453: }
1455: /*@C
1456: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1458: Collective on DM
1460: Input Parameters:
1461: + comm - The MPI comm
1462: . dim - The spatial dimension
1463: . Nc - The number of components
1464: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1465: . prefix - The options prefix, or NULL
1466: - qorder - The quadrature order
1468: Output Parameter:
1469: . fem - The PetscFE object
1471: Level: beginner
1473: .keywords: PetscFE, finite element
1474: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1475: @*/
1476: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1477: {
1478: PetscQuadrature q, fq;
1479: DM K;
1480: PetscSpace P;
1481: PetscDualSpace Q;
1482: PetscInt order, quadPointsPerEdge;
1483: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1484: PetscErrorCode ierr;
1487: /* Create space */
1488: PetscSpaceCreate(comm, &P);
1489: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1490: PetscSpacePolynomialSetTensor(P, tensor);
1491: PetscSpaceSetNumComponents(P, Nc);
1492: PetscSpaceSetNumVariables(P, dim);
1493: PetscSpaceSetFromOptions(P);
1494: PetscSpaceSetUp(P);
1495: PetscSpaceGetDegree(P, &order, NULL);
1496: PetscSpacePolynomialGetTensor(P, &tensor);
1497: /* Create dual space */
1498: PetscDualSpaceCreate(comm, &Q);
1499: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1500: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1501: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1502: PetscDualSpaceSetDM(Q, K);
1503: DMDestroy(&K);
1504: PetscDualSpaceSetNumComponents(Q, Nc);
1505: PetscDualSpaceSetOrder(Q, order);
1506: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1507: PetscDualSpaceSetFromOptions(Q);
1508: PetscDualSpaceSetUp(Q);
1509: /* Create element */
1510: PetscFECreate(comm, fem);
1511: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1512: PetscFESetFromOptions(*fem);
1513: PetscFESetBasisSpace(*fem, P);
1514: PetscFESetDualSpace(*fem, Q);
1515: PetscFESetNumComponents(*fem, Nc);
1516: PetscFESetUp(*fem);
1517: PetscSpaceDestroy(&P);
1518: PetscDualSpaceDestroy(&Q);
1519: /* Create quadrature (with specified order if given) */
1520: qorder = qorder >= 0 ? qorder : order;
1521: PetscObjectOptionsBegin((PetscObject)*fem);
1522: PetscOptionsInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadture points per edge","PetscFECreateDefault",qorder,&qorder,NULL);
1523: PetscOptionsEnd();
1524: quadPointsPerEdge = PetscMax(qorder + 1,1);
1525: if (isSimplex) {
1526: PetscDTGaussJacobiQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1527: PetscDTGaussJacobiQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1528: }
1529: else {
1530: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1531: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1532: }
1533: PetscFESetQuadrature(*fem, q);
1534: PetscFESetFaceQuadrature(*fem, fq);
1535: PetscQuadratureDestroy(&q);
1536: PetscQuadratureDestroy(&fq);
1537: return(0);
1538: }