Actual source code: fe.c
petsc-3.13.6 2020-09-29
1: /* Basis Jet Tabulation
3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6: as a prime basis.
8: \psi_i = \sum_k \alpha_{ki} \phi_k
10: Our nodal basis is defined in terms of the dual basis $n_j$
12: n_j \cdot \psi_i = \delta_{ji}
14: and we may act on the first equation to obtain
16: n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17: \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18: I = V \alpha
20: so the coefficients of the nodal basis in the prime basis are
22: \alpha = V^{-1}
24: We will define the dual basis vectors $n_j$ using a quadrature rule.
26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28: be implemented exactly as in FIAT using functionals $L_j$.
30: I will have to count the degrees correctly for the Legendre product when we are on simplices.
32: We will have three objects:
33: - Space, P: this just need point evaluation I think
34: - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35: - FEM: This keeps {P, P', Q}
36: */
37: #include <petsc/private/petscfeimpl.h>
38: #include <petscdmplex.h>
40: PetscBool FEcite = PETSC_FALSE;
41: const char FECitation[] = "@article{kirby2004,\n"
42: " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
43: " journal = {ACM Transactions on Mathematical Software},\n"
44: " author = {Robert C. Kirby},\n"
45: " volume = {30},\n"
46: " number = {4},\n"
47: " pages = {502--516},\n"
48: " doi = {10.1145/1039813.1039820},\n"
49: " year = {2004}\n}\n";
51: PetscClassId PETSCFE_CLASSID = 0;
53: PetscFunctionList PetscFEList = NULL;
54: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
56: /*@C
57: PetscFERegister - Adds a new PetscFE implementation
59: Not Collective
61: Input Parameters:
62: + name - The name of a new user-defined creation routine
63: - create_func - The creation routine itself
65: Notes:
66: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
68: Sample usage:
69: .vb
70: PetscFERegister("my_fe", MyPetscFECreate);
71: .ve
73: Then, your PetscFE type can be chosen with the procedural interface via
74: .vb
75: PetscFECreate(MPI_Comm, PetscFE *);
76: PetscFESetType(PetscFE, "my_fe");
77: .ve
78: or at runtime via the option
79: .vb
80: -petscfe_type my_fe
81: .ve
83: Level: advanced
85: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
87: @*/
88: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
89: {
93: PetscFunctionListAdd(&PetscFEList, sname, function);
94: return(0);
95: }
97: /*@C
98: PetscFESetType - Builds a particular PetscFE
100: Collective on fem
102: Input Parameters:
103: + fem - The PetscFE object
104: - name - The kind of FEM space
106: Options Database Key:
107: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
109: Level: intermediate
111: .seealso: PetscFEGetType(), PetscFECreate()
112: @*/
113: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
114: {
115: PetscErrorCode (*r)(PetscFE);
116: PetscBool match;
121: PetscObjectTypeCompare((PetscObject) fem, name, &match);
122: if (match) return(0);
124: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
125: PetscFunctionListFind(PetscFEList, name, &r);
126: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
128: if (fem->ops->destroy) {
129: (*fem->ops->destroy)(fem);
130: fem->ops->destroy = NULL;
131: }
132: (*r)(fem);
133: PetscObjectChangeTypeName((PetscObject) fem, name);
134: return(0);
135: }
137: /*@C
138: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
140: Not Collective
142: Input Parameter:
143: . fem - The PetscFE
145: Output Parameter:
146: . name - The PetscFE type name
148: Level: intermediate
150: .seealso: PetscFESetType(), PetscFECreate()
151: @*/
152: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
153: {
159: if (!PetscFERegisterAllCalled) {
160: PetscFERegisterAll();
161: }
162: *name = ((PetscObject) fem)->type_name;
163: return(0);
164: }
166: /*@C
167: PetscFEViewFromOptions - View from Options
169: Collective on PetscFE
171: Input Parameters:
172: + A - the PetscFE object
173: . obj - Optional object
174: - name - command line option
176: Level: intermediate
177: .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
178: @*/
179: PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
180: {
185: PetscObjectViewFromOptions((PetscObject)A,obj,name);
186: return(0);
187: }
189: /*@C
190: PetscFEView - Views a PetscFE
192: Collective on fem
194: Input Parameter:
195: + fem - the PetscFE object to view
196: - viewer - the viewer
198: Level: beginner
200: .seealso PetscFEDestroy()
201: @*/
202: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
203: {
204: PetscBool iascii;
210: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
211: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
212: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
213: if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
214: return(0);
215: }
217: /*@
218: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
220: Collective on fem
222: Input Parameter:
223: . fem - the PetscFE object to set options for
225: Options Database:
226: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
227: - -petscfe_num_batches - the number of cell batches to integrate serially
229: Level: intermediate
231: .seealso PetscFEView()
232: @*/
233: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
234: {
235: const char *defaultType;
236: char name[256];
237: PetscBool flg;
242: if (!((PetscObject) fem)->type_name) {
243: defaultType = PETSCFEBASIC;
244: } else {
245: defaultType = ((PetscObject) fem)->type_name;
246: }
247: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
249: PetscObjectOptionsBegin((PetscObject) fem);
250: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
251: if (flg) {
252: PetscFESetType(fem, name);
253: } else if (!((PetscObject) fem)->type_name) {
254: PetscFESetType(fem, defaultType);
255: }
256: PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
257: PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
258: if (fem->ops->setfromoptions) {
259: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
260: }
261: /* process any options handlers added with PetscObjectAddOptionsHandler() */
262: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
263: PetscOptionsEnd();
264: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
265: return(0);
266: }
268: /*@C
269: PetscFESetUp - Construct data structures for the PetscFE
271: Collective on fem
273: Input Parameter:
274: . fem - the PetscFE object to setup
276: Level: intermediate
278: .seealso PetscFEView(), PetscFEDestroy()
279: @*/
280: PetscErrorCode PetscFESetUp(PetscFE fem)
281: {
286: if (fem->setupcalled) return(0);
287: fem->setupcalled = PETSC_TRUE;
288: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
289: return(0);
290: }
292: /*@
293: PetscFEDestroy - Destroys a PetscFE object
295: Collective on fem
297: Input Parameter:
298: . fem - the PetscFE object to destroy
300: Level: beginner
302: .seealso PetscFEView()
303: @*/
304: PetscErrorCode PetscFEDestroy(PetscFE *fem)
305: {
309: if (!*fem) return(0);
312: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
313: ((PetscObject) (*fem))->refct = 0;
315: if ((*fem)->subspaces) {
316: PetscInt dim, d;
318: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
319: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
320: }
321: PetscFree((*fem)->subspaces);
322: PetscFree((*fem)->invV);
323: PetscTabulationDestroy(&(*fem)->T);
324: PetscTabulationDestroy(&(*fem)->Tf);
325: PetscTabulationDestroy(&(*fem)->Tc);
326: PetscSpaceDestroy(&(*fem)->basisSpace);
327: PetscDualSpaceDestroy(&(*fem)->dualSpace);
328: PetscQuadratureDestroy(&(*fem)->quadrature);
329: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
331: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
332: PetscHeaderDestroy(fem);
333: return(0);
334: }
336: /*@
337: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
339: Collective
341: Input Parameter:
342: . comm - The communicator for the PetscFE object
344: Output Parameter:
345: . fem - The PetscFE object
347: Level: beginner
349: .seealso: PetscFESetType(), PETSCFEGALERKIN
350: @*/
351: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
352: {
353: PetscFE f;
358: PetscCitationsRegister(FECitation,&FEcite);
359: *fem = NULL;
360: PetscFEInitializePackage();
362: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
364: f->basisSpace = NULL;
365: f->dualSpace = NULL;
366: f->numComponents = 1;
367: f->subspaces = NULL;
368: f->invV = NULL;
369: f->T = NULL;
370: f->Tf = NULL;
371: f->Tc = NULL;
372: PetscArrayzero(&f->quadrature, 1);
373: PetscArrayzero(&f->faceQuadrature, 1);
374: f->blockSize = 0;
375: f->numBlocks = 1;
376: f->batchSize = 0;
377: f->numBatches = 1;
379: *fem = f;
380: return(0);
381: }
383: /*@
384: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
386: Not collective
388: Input Parameter:
389: . fem - The PetscFE object
391: Output Parameter:
392: . dim - The spatial dimension
394: Level: intermediate
396: .seealso: PetscFECreate()
397: @*/
398: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
399: {
400: DM dm;
406: PetscDualSpaceGetDM(fem->dualSpace, &dm);
407: DMGetDimension(dm, dim);
408: return(0);
409: }
411: /*@
412: PetscFESetNumComponents - Sets the number of components in the element
414: Not collective
416: Input Parameters:
417: + fem - The PetscFE object
418: - comp - The number of field components
420: Level: intermediate
422: .seealso: PetscFECreate()
423: @*/
424: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
425: {
428: fem->numComponents = comp;
429: return(0);
430: }
432: /*@
433: PetscFEGetNumComponents - Returns the number of components in the element
435: Not collective
437: Input Parameter:
438: . fem - The PetscFE object
440: Output Parameter:
441: . comp - The number of field components
443: Level: intermediate
445: .seealso: PetscFECreate()
446: @*/
447: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
448: {
452: *comp = fem->numComponents;
453: return(0);
454: }
456: /*@
457: PetscFESetTileSizes - Sets the tile sizes for evaluation
459: Not collective
461: Input Parameters:
462: + fem - The PetscFE object
463: . blockSize - The number of elements in a block
464: . numBlocks - The number of blocks in a batch
465: . batchSize - The number of elements in a batch
466: - numBatches - The number of batches in a chunk
468: Level: intermediate
470: .seealso: PetscFECreate()
471: @*/
472: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
473: {
476: fem->blockSize = blockSize;
477: fem->numBlocks = numBlocks;
478: fem->batchSize = batchSize;
479: fem->numBatches = numBatches;
480: return(0);
481: }
483: /*@
484: PetscFEGetTileSizes - Returns the tile sizes for evaluation
486: Not collective
488: Input Parameter:
489: . fem - The PetscFE object
491: Output Parameters:
492: + blockSize - The number of elements in a block
493: . numBlocks - The number of blocks in a batch
494: . batchSize - The number of elements in a batch
495: - numBatches - The number of batches in a chunk
497: Level: intermediate
499: .seealso: PetscFECreate()
500: @*/
501: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
502: {
509: if (blockSize) *blockSize = fem->blockSize;
510: if (numBlocks) *numBlocks = fem->numBlocks;
511: if (batchSize) *batchSize = fem->batchSize;
512: if (numBatches) *numBatches = fem->numBatches;
513: return(0);
514: }
516: /*@
517: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
519: Not collective
521: Input Parameter:
522: . fem - The PetscFE object
524: Output Parameter:
525: . sp - The PetscSpace object
527: Level: intermediate
529: .seealso: PetscFECreate()
530: @*/
531: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
532: {
536: *sp = fem->basisSpace;
537: return(0);
538: }
540: /*@
541: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
543: Not collective
545: Input Parameters:
546: + fem - The PetscFE object
547: - sp - The PetscSpace object
549: Level: intermediate
551: .seealso: PetscFECreate()
552: @*/
553: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
554: {
560: PetscSpaceDestroy(&fem->basisSpace);
561: fem->basisSpace = sp;
562: PetscObjectReference((PetscObject) fem->basisSpace);
563: return(0);
564: }
566: /*@
567: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
569: Not collective
571: Input Parameter:
572: . fem - The PetscFE object
574: Output Parameter:
575: . sp - The PetscDualSpace object
577: Level: intermediate
579: .seealso: PetscFECreate()
580: @*/
581: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
582: {
586: *sp = fem->dualSpace;
587: return(0);
588: }
590: /*@
591: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
593: Not collective
595: Input Parameters:
596: + fem - The PetscFE object
597: - sp - The PetscDualSpace object
599: Level: intermediate
601: .seealso: PetscFECreate()
602: @*/
603: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
604: {
610: PetscDualSpaceDestroy(&fem->dualSpace);
611: fem->dualSpace = sp;
612: PetscObjectReference((PetscObject) fem->dualSpace);
613: return(0);
614: }
616: /*@
617: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
619: Not collective
621: Input Parameter:
622: . fem - The PetscFE object
624: Output Parameter:
625: . q - The PetscQuadrature object
627: Level: intermediate
629: .seealso: PetscFECreate()
630: @*/
631: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
632: {
636: *q = fem->quadrature;
637: return(0);
638: }
640: /*@
641: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
643: Not collective
645: Input Parameters:
646: + fem - The PetscFE object
647: - q - The PetscQuadrature object
649: Level: intermediate
651: .seealso: PetscFECreate()
652: @*/
653: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
654: {
655: PetscInt Nc, qNc;
660: if (q == fem->quadrature) return(0);
661: PetscFEGetNumComponents(fem, &Nc);
662: PetscQuadratureGetNumComponents(q, &qNc);
663: 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);
664: PetscTabulationDestroy(&fem->T);
665: PetscTabulationDestroy(&fem->Tc);
666: PetscObjectReference((PetscObject) q);
667: PetscQuadratureDestroy(&fem->quadrature);
668: fem->quadrature = q;
669: return(0);
670: }
672: /*@
673: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
675: Not collective
677: Input Parameter:
678: . fem - The PetscFE object
680: Output Parameter:
681: . q - The PetscQuadrature object
683: Level: intermediate
685: .seealso: PetscFECreate()
686: @*/
687: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
688: {
692: *q = fem->faceQuadrature;
693: return(0);
694: }
696: /*@
697: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
699: Not collective
701: Input Parameters:
702: + fem - The PetscFE object
703: - q - The PetscQuadrature object
705: Level: intermediate
707: .seealso: PetscFECreate()
708: @*/
709: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
710: {
711: PetscInt Nc, qNc;
716: PetscFEGetNumComponents(fem, &Nc);
717: PetscQuadratureGetNumComponents(q, &qNc);
718: 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);
719: PetscTabulationDestroy(&fem->Tf);
720: PetscQuadratureDestroy(&fem->faceQuadrature);
721: fem->faceQuadrature = q;
722: PetscObjectReference((PetscObject) q);
723: return(0);
724: }
726: /*@
727: PetscFECopyQuadrature - Copy both volumetric and surface quadrature
729: Not collective
731: Input Parameters:
732: + sfe - The PetscFE source for the quadratures
733: - tfe - The PetscFE target for the quadratures
735: Level: intermediate
737: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
738: @*/
739: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
740: {
741: PetscQuadrature q;
742: PetscErrorCode ierr;
747: PetscFEGetQuadrature(sfe, &q);
748: PetscFESetQuadrature(tfe, q);
749: PetscFEGetFaceQuadrature(sfe, &q);
750: PetscFESetFaceQuadrature(tfe, q);
751: return(0);
752: }
754: /*@C
755: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
757: Not collective
759: Input Parameter:
760: . fem - The PetscFE object
762: Output Parameter:
763: . numDof - Array with the number of dofs per dimension
765: Level: intermediate
767: .seealso: PetscFECreate()
768: @*/
769: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
770: {
776: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
777: return(0);
778: }
780: /*@C
781: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
783: Not collective
785: Input Parameter:
786: . fem - The PetscFE object
788: Output Parameter:
789: . T - The basis function values and derivatives at quadrature points
791: Note:
792: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
793: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
794: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
796: Level: intermediate
798: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
799: @*/
800: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T)
801: {
802: PetscInt npoints;
803: const PetscReal *points;
804: PetscErrorCode ierr;
809: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
810: if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);}
811: *T = fem->T;
812: return(0);
813: }
815: /*@C
816: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
818: Not collective
820: Input Parameter:
821: . fem - The PetscFE object
823: Output Parameters:
824: . Tf - The basis function values and derviatives at face quadrature points
826: Note:
827: $ T->T[0] = Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
828: $ T->T[1] = Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d
829: $ T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e
831: Level: intermediate
833: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
834: @*/
835: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf)
836: {
837: PetscErrorCode ierr;
842: if (!fem->Tf) {
843: const PetscReal xi0[3] = {-1., -1., -1.};
844: PetscReal v0[3], J[9], detJ;
845: PetscQuadrature fq;
846: PetscDualSpace sp;
847: DM dm;
848: const PetscInt *faces;
849: PetscInt dim, numFaces, f, npoints, q;
850: const PetscReal *points;
851: PetscReal *facePoints;
853: PetscFEGetDualSpace(fem, &sp);
854: PetscDualSpaceGetDM(sp, &dm);
855: DMGetDimension(dm, &dim);
856: DMPlexGetConeSize(dm, 0, &numFaces);
857: DMPlexGetCone(dm, 0, &faces);
858: PetscFEGetFaceQuadrature(fem, &fq);
859: if (fq) {
860: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
861: PetscMalloc1(numFaces*npoints*dim, &facePoints);
862: for (f = 0; f < numFaces; ++f) {
863: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
864: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
865: }
866: PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);
867: PetscFree(facePoints);
868: }
869: }
870: *Tf = fem->Tf;
871: return(0);
872: }
874: /*@C
875: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
877: Not collective
879: Input Parameter:
880: . fem - The PetscFE object
882: Output Parameters:
883: . Tc - The basis function values at face centroid points
885: Note:
886: $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
888: Level: intermediate
890: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
891: @*/
892: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
893: {
894: PetscErrorCode ierr;
899: if (!fem->Tc) {
900: PetscDualSpace sp;
901: DM dm;
902: const PetscInt *cone;
903: PetscReal *centroids;
904: PetscInt dim, numFaces, f;
906: PetscFEGetDualSpace(fem, &sp);
907: PetscDualSpaceGetDM(sp, &dm);
908: DMGetDimension(dm, &dim);
909: DMPlexGetConeSize(dm, 0, &numFaces);
910: DMPlexGetCone(dm, 0, &cone);
911: PetscMalloc1(numFaces*dim, ¢roids);
912: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
913: PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
914: PetscFree(centroids);
915: }
916: *Tc = fem->Tc;
917: return(0);
918: }
920: /*@C
921: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
923: Not collective
925: Input Parameters:
926: + fem - The PetscFE object
927: . nrepl - The number of replicas
928: . npoints - The number of tabulation points in a replica
929: . points - The tabulation point coordinates
930: - K - The number of derivatives calculated
932: Output Parameter:
933: . T - The basis function values and derivatives at tabulation points
935: Note:
936: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
937: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
938: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
940: Level: intermediate
942: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
943: @*/
944: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
945: {
946: DM dm;
947: PetscDualSpace Q;
948: PetscInt Nb; /* Dimension of FE space P */
949: PetscInt Nc; /* Field components */
950: PetscInt cdim; /* Reference coordinate dimension */
951: PetscInt k;
952: PetscErrorCode ierr;
955: if (!npoints || !fem->dualSpace || K < 0) {
956: *T = NULL;
957: return(0);
958: }
962: PetscFEGetDualSpace(fem, &Q);
963: PetscDualSpaceGetDM(Q, &dm);
964: DMGetDimension(dm, &cdim);
965: PetscDualSpaceGetDimension(Q, &Nb);
966: PetscFEGetNumComponents(fem, &Nc);
967: PetscMalloc1(1, T);
968: (*T)->K = !cdim ? 0 : K;
969: (*T)->Nr = nrepl;
970: (*T)->Np = npoints;
971: (*T)->Nb = Nb;
972: (*T)->Nc = Nc;
973: (*T)->cdim = cdim;
974: PetscMalloc1((*T)->K+1, &(*T)->T);
975: for (k = 0; k <= (*T)->K; ++k) {
976: PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
977: }
978: (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
979: return(0);
980: }
982: /*@C
983: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
985: Not collective
987: Input Parameters:
988: + fem - The PetscFE object
989: . npoints - The number of tabulation points
990: . points - The tabulation point coordinates
991: . K - The number of derivatives calculated
992: - T - An existing tabulation object with enough allocated space
994: Output Parameter:
995: . T - The basis function values and derivatives at tabulation points
997: Note:
998: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
999: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1000: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1002: Level: intermediate
1004: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1005: @*/
1006: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1007: {
1011: if (!npoints || !fem->dualSpace || K < 0) return(0);
1015: #ifdef PETSC_USE_DEBUG
1016: {
1017: DM dm;
1018: PetscDualSpace Q;
1019: PetscInt Nb; /* Dimension of FE space P */
1020: PetscInt Nc; /* Field components */
1021: PetscInt cdim; /* Reference coordinate dimension */
1023: PetscFEGetDualSpace(fem, &Q);
1024: PetscDualSpaceGetDM(Q, &dm);
1025: DMGetDimension(dm, &cdim);
1026: PetscDualSpaceGetDimension(Q, &Nb);
1027: PetscFEGetNumComponents(fem, &Nc);
1028: if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K);
1029: if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1030: if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1031: if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1032: }
1033: #endif
1034: T->Nr = 1;
1035: T->Np = npoints;
1036: (*fem->ops->createtabulation)(fem, npoints, points, K, T);
1037: return(0);
1038: }
1040: /*@C
1041: PetscTabulationDestroy - Frees memory from the associated tabulation.
1043: Not collective
1045: Input Parameter:
1046: . T - The tabulation
1048: Level: intermediate
1050: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1051: @*/
1052: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1053: {
1054: PetscInt k;
1059: if (!T || !(*T)) return(0);
1060: for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);}
1061: PetscFree((*T)->T);
1062: PetscFree(*T);
1063: *T = NULL;
1064: return(0);
1065: }
1067: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1068: {
1069: PetscSpace bsp, bsubsp;
1070: PetscDualSpace dsp, dsubsp;
1071: PetscInt dim, depth, numComp, i, j, coneSize, order;
1072: PetscFEType type;
1073: DM dm;
1074: DMLabel label;
1075: PetscReal *xi, *v, *J, detJ;
1076: const char *name;
1077: PetscQuadrature origin, fullQuad, subQuad;
1083: PetscFEGetBasisSpace(fe,&bsp);
1084: PetscFEGetDualSpace(fe,&dsp);
1085: PetscDualSpaceGetDM(dsp,&dm);
1086: DMGetDimension(dm,&dim);
1087: DMPlexGetDepthLabel(dm,&label);
1088: DMLabelGetValue(label,refPoint,&depth);
1089: PetscCalloc1(depth,&xi);
1090: PetscMalloc1(dim,&v);
1091: PetscMalloc1(dim*dim,&J);
1092: for (i = 0; i < depth; i++) xi[i] = 0.;
1093: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1094: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1095: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1096: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1097: for (i = 1; i < dim; i++) {
1098: for (j = 0; j < depth; j++) {
1099: J[i * depth + j] = J[i * dim + j];
1100: }
1101: }
1102: PetscQuadratureDestroy(&origin);
1103: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1104: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1105: PetscSpaceSetUp(bsubsp);
1106: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1107: PetscFEGetType(fe,&type);
1108: PetscFESetType(*trFE,type);
1109: PetscFEGetNumComponents(fe,&numComp);
1110: PetscFESetNumComponents(*trFE,numComp);
1111: PetscFESetBasisSpace(*trFE,bsubsp);
1112: PetscFESetDualSpace(*trFE,dsubsp);
1113: PetscObjectGetName((PetscObject) fe, &name);
1114: if (name) {PetscFESetName(*trFE, name);}
1115: PetscFEGetQuadrature(fe,&fullQuad);
1116: PetscQuadratureGetOrder(fullQuad,&order);
1117: DMPlexGetConeSize(dm,refPoint,&coneSize);
1118: if (coneSize == 2 * depth) {
1119: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1120: } else {
1121: PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1122: }
1123: PetscFESetQuadrature(*trFE,subQuad);
1124: PetscFESetUp(*trFE);
1125: PetscQuadratureDestroy(&subQuad);
1126: PetscSpaceDestroy(&bsubsp);
1127: return(0);
1128: }
1130: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1131: {
1132: PetscInt hStart, hEnd;
1133: PetscDualSpace dsp;
1134: DM dm;
1140: *trFE = NULL;
1141: PetscFEGetDualSpace(fe,&dsp);
1142: PetscDualSpaceGetDM(dsp,&dm);
1143: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1144: if (hEnd <= hStart) return(0);
1145: PetscFECreatePointTrace(fe,hStart,trFE);
1146: return(0);
1147: }
1150: /*@
1151: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1153: Not collective
1155: Input Parameter:
1156: . fe - The PetscFE
1158: Output Parameter:
1159: . dim - The dimension
1161: Level: intermediate
1163: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1164: @*/
1165: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1166: {
1172: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1173: return(0);
1174: }
1176: /*@C
1177: PetscFEPushforward - Map the reference element function to real space
1179: Input Parameters:
1180: + fe - The PetscFE
1181: . fegeom - The cell geometry
1182: . Nv - The number of function values
1183: - vals - The function values
1185: Output Parameter:
1186: . vals - The transformed function values
1188: Level: advanced
1190: Note: This just forwards the call onto PetscDualSpacePushforward().
1192: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1194: .seealso: PetscDualSpacePushforward()
1195: @*/
1196: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1197: {
1201: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1202: return(0);
1203: }
1205: /*@C
1206: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1208: Input Parameters:
1209: + fe - The PetscFE
1210: . fegeom - The cell geometry
1211: . Nv - The number of function gradient values
1212: - vals - The function gradient values
1214: Output Parameter:
1215: . vals - The transformed function gradient values
1217: Level: advanced
1219: Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1221: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1223: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1224: @*/
1225: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1226: {
1230: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1231: return(0);
1232: }
1234: /*
1235: Purpose: Compute element vector for chunk of elements
1237: Input:
1238: Sizes:
1239: Ne: number of elements
1240: Nf: number of fields
1241: PetscFE
1242: dim: spatial dimension
1243: Nb: number of basis functions
1244: Nc: number of field components
1245: PetscQuadrature
1246: Nq: number of quadrature points
1248: Geometry:
1249: PetscFEGeom[Ne] possibly *Nq
1250: PetscReal v0s[dim]
1251: PetscReal n[dim]
1252: PetscReal jacobians[dim*dim]
1253: PetscReal jacobianInverses[dim*dim]
1254: PetscReal jacobianDeterminants
1255: FEM:
1256: PetscFE
1257: PetscQuadrature
1258: PetscReal quadPoints[Nq*dim]
1259: PetscReal quadWeights[Nq]
1260: PetscReal basis[Nq*Nb*Nc]
1261: PetscReal basisDer[Nq*Nb*Nc*dim]
1262: PetscScalar coefficients[Ne*Nb*Nc]
1263: PetscScalar elemVec[Ne*Nb*Nc]
1265: Problem:
1266: PetscInt f: the active field
1267: f0, f1
1269: Work Space:
1270: PetscFE
1271: PetscScalar f0[Nq*dim];
1272: PetscScalar f1[Nq*dim*dim];
1273: PetscScalar u[Nc];
1274: PetscScalar gradU[Nc*dim];
1275: PetscReal x[dim];
1276: PetscScalar realSpaceDer[dim];
1278: Purpose: Compute element vector for N_cb batches of elements
1280: Input:
1281: Sizes:
1282: N_cb: Number of serial cell batches
1284: Geometry:
1285: PetscReal v0s[Ne*dim]
1286: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1287: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1288: PetscReal jacobianDeterminants[Ne] possibly *Nq
1289: FEM:
1290: static PetscReal quadPoints[Nq*dim]
1291: static PetscReal quadWeights[Nq]
1292: static PetscReal basis[Nq*Nb*Nc]
1293: static PetscReal basisDer[Nq*Nb*Nc*dim]
1294: PetscScalar coefficients[Ne*Nb*Nc]
1295: PetscScalar elemVec[Ne*Nb*Nc]
1297: ex62.c:
1298: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1299: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1300: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1301: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1303: ex52.c:
1304: 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)
1305: 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)
1307: ex52_integrateElement.cu
1308: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1310: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1311: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1312: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1314: ex52_integrateElementOpenCL.c:
1315: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1316: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1317: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1319: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1320: */
1322: /*@C
1323: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1325: Not collective
1327: Input Parameters:
1328: + fem - The PetscFE object for the field being integrated
1329: . prob - The PetscDS specifying the discretizations and continuum functions
1330: . field - The field being integrated
1331: . Ne - The number of elements in the chunk
1332: . cgeom - The cell geometry for each cell in the chunk
1333: . coefficients - The array of FEM basis coefficients for the elements
1334: . probAux - The PetscDS specifying the auxiliary discretizations
1335: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1337: Output Parameter:
1338: . integral - the integral for this field
1340: Level: intermediate
1342: .seealso: PetscFEIntegrateResidual()
1343: @*/
1344: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1345: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1346: {
1347: PetscFE fe;
1352: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1353: if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1354: return(0);
1355: }
1357: /*@C
1358: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1360: Not collective
1362: Input Parameters:
1363: + fem - The PetscFE object for the field being integrated
1364: . prob - The PetscDS specifying the discretizations and continuum functions
1365: . field - The field being integrated
1366: . obj_func - The function to be integrated
1367: . Ne - The number of elements in the chunk
1368: . fgeom - The face geometry for each face in the chunk
1369: . coefficients - The array of FEM basis coefficients for the elements
1370: . probAux - The PetscDS specifying the auxiliary discretizations
1371: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1373: Output Parameter:
1374: . integral - the integral for this field
1376: Level: intermediate
1378: .seealso: PetscFEIntegrateResidual()
1379: @*/
1380: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1381: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1382: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1383: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1384: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1385: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1386: {
1387: PetscFE fe;
1392: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1393: if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1394: return(0);
1395: }
1397: /*@C
1398: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1400: Not collective
1402: Input Parameters:
1403: + fem - The PetscFE object for the field being integrated
1404: . prob - The PetscDS specifying the discretizations and continuum functions
1405: . field - The field being integrated
1406: . Ne - The number of elements in the chunk
1407: . cgeom - The cell geometry for each cell in the chunk
1408: . coefficients - The array of FEM basis coefficients for the elements
1409: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1410: . probAux - The PetscDS specifying the auxiliary discretizations
1411: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1412: - t - The time
1414: Output Parameter:
1415: . elemVec - the element residual vectors from each element
1417: Note:
1418: $ Loop over batch of elements (e):
1419: $ Loop over quadrature points (q):
1420: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1421: $ Call f_0 and f_1
1422: $ Loop over element vector entries (f,fc --> i):
1423: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1425: Level: intermediate
1427: .seealso: PetscFEIntegrateResidual()
1428: @*/
1429: PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1430: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1431: {
1432: PetscFE fe;
1437: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1438: if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1439: return(0);
1440: }
1442: /*@C
1443: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1445: Not collective
1447: Input Parameters:
1448: + fem - The PetscFE object for the field being integrated
1449: . prob - The PetscDS specifying the discretizations and continuum functions
1450: . field - The field being integrated
1451: . Ne - The number of elements in the chunk
1452: . fgeom - The face geometry for each cell in the chunk
1453: . coefficients - The array of FEM basis coefficients for the elements
1454: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1455: . probAux - The PetscDS specifying the auxiliary discretizations
1456: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1457: - t - The time
1459: Output Parameter:
1460: . elemVec - the element residual vectors from each element
1462: Level: intermediate
1464: .seealso: PetscFEIntegrateResidual()
1465: @*/
1466: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1467: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1468: {
1469: PetscFE fe;
1474: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1475: if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1476: return(0);
1477: }
1479: /*@C
1480: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1482: Not collective
1484: Input Parameters:
1485: + fem - The PetscFE object for the field being integrated
1486: . prob - The PetscDS specifying the discretizations and continuum functions
1487: . jtype - The type of matrix pointwise functions that should be used
1488: . fieldI - The test field being integrated
1489: . fieldJ - The basis field being integrated
1490: . Ne - The number of elements in the chunk
1491: . cgeom - The cell geometry for each cell in the chunk
1492: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1493: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1494: . probAux - The PetscDS specifying the auxiliary discretizations
1495: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1496: . t - The time
1497: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1499: Output Parameter:
1500: . elemMat - the element matrices for the Jacobian from each element
1502: Note:
1503: $ Loop over batch of elements (e):
1504: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1505: $ Loop over quadrature points (q):
1506: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1507: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1508: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1509: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1510: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1511: Level: intermediate
1513: .seealso: PetscFEIntegrateResidual()
1514: @*/
1515: PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1516: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1517: {
1518: PetscFE fe;
1523: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1524: if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1525: return(0);
1526: }
1528: /*@C
1529: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1531: Not collective
1533: Input Parameters:
1534: + prob - The PetscDS specifying the discretizations and continuum functions
1535: . fieldI - The test field being integrated
1536: . fieldJ - The basis field being integrated
1537: . Ne - The number of elements in the chunk
1538: . fgeom - The face geometry for each cell in the chunk
1539: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1540: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1541: . probAux - The PetscDS specifying the auxiliary discretizations
1542: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1543: . t - The time
1544: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1546: Output Parameter:
1547: . elemMat - the element matrices for the Jacobian from each element
1549: Note:
1550: $ Loop over batch of elements (e):
1551: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1552: $ Loop over quadrature points (q):
1553: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1554: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1555: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1556: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1557: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1558: Level: intermediate
1560: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1561: @*/
1562: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1563: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1564: {
1565: PetscFE fe;
1570: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1571: if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1572: return(0);
1573: }
1575: /*@
1576: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1578: Input Parameters:
1579: + fe - The finite element space
1580: - height - The height of the Plex point
1582: Output Parameter:
1583: . subfe - The subspace of this FE space
1585: Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1587: Level: advanced
1589: .seealso: PetscFECreateDefault()
1590: @*/
1591: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1592: {
1593: PetscSpace P, subP;
1594: PetscDualSpace Q, subQ;
1595: PetscQuadrature subq;
1596: PetscFEType fetype;
1597: PetscInt dim, Nc;
1598: PetscErrorCode ierr;
1603: if (height == 0) {
1604: *subfe = fe;
1605: return(0);
1606: }
1607: PetscFEGetBasisSpace(fe, &P);
1608: PetscFEGetDualSpace(fe, &Q);
1609: PetscFEGetNumComponents(fe, &Nc);
1610: PetscFEGetFaceQuadrature(fe, &subq);
1611: PetscDualSpaceGetDimension(Q, &dim);
1612: 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);}
1613: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1614: if (height <= dim) {
1615: if (!fe->subspaces[height-1]) {
1616: PetscFE sub;
1617: const char *name;
1619: PetscSpaceGetHeightSubspace(P, height, &subP);
1620: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1621: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1622: PetscObjectGetName((PetscObject) fe, &name);
1623: PetscObjectSetName((PetscObject) sub, name);
1624: PetscFEGetType(fe, &fetype);
1625: PetscFESetType(sub, fetype);
1626: PetscFESetBasisSpace(sub, subP);
1627: PetscFESetDualSpace(sub, subQ);
1628: PetscFESetNumComponents(sub, Nc);
1629: PetscFESetUp(sub);
1630: PetscFESetQuadrature(sub, subq);
1631: fe->subspaces[height-1] = sub;
1632: }
1633: *subfe = fe->subspaces[height-1];
1634: } else {
1635: *subfe = NULL;
1636: }
1637: return(0);
1638: }
1640: /*@
1641: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1642: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1643: sparsity). It is also used to create an interpolation between regularly refined meshes.
1645: Collective on fem
1647: Input Parameter:
1648: . fe - The initial PetscFE
1650: Output Parameter:
1651: . feRef - The refined PetscFE
1653: Level: advanced
1655: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1656: @*/
1657: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1658: {
1659: PetscSpace P, Pref;
1660: PetscDualSpace Q, Qref;
1661: DM K, Kref;
1662: PetscQuadrature q, qref;
1663: const PetscReal *v0, *jac;
1664: PetscInt numComp, numSubelements;
1665: PetscInt cStart, cEnd, c;
1666: PetscDualSpace *cellSpaces;
1667: PetscErrorCode ierr;
1670: PetscFEGetBasisSpace(fe, &P);
1671: PetscFEGetDualSpace(fe, &Q);
1672: PetscFEGetQuadrature(fe, &q);
1673: PetscDualSpaceGetDM(Q, &K);
1674: /* Create space */
1675: PetscObjectReference((PetscObject) P);
1676: Pref = P;
1677: /* Create dual space */
1678: PetscDualSpaceDuplicate(Q, &Qref);
1679: PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1680: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1681: PetscDualSpaceSetDM(Qref, Kref);
1682: DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1683: PetscMalloc1(cEnd - cStart, &cellSpaces);
1684: /* TODO: fix for non-uniform refinement */
1685: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1686: PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1687: PetscFree(cellSpaces);
1688: DMDestroy(&Kref);
1689: PetscDualSpaceSetUp(Qref);
1690: /* Create element */
1691: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1692: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1693: PetscFESetBasisSpace(*feRef, Pref);
1694: PetscFESetDualSpace(*feRef, Qref);
1695: PetscFEGetNumComponents(fe, &numComp);
1696: PetscFESetNumComponents(*feRef, numComp);
1697: PetscFESetUp(*feRef);
1698: PetscSpaceDestroy(&Pref);
1699: PetscDualSpaceDestroy(&Qref);
1700: /* Create quadrature */
1701: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1702: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1703: PetscFESetQuadrature(*feRef, qref);
1704: PetscQuadratureDestroy(&qref);
1705: return(0);
1706: }
1708: /*@C
1709: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1711: Collective
1713: Input Parameters:
1714: + comm - The MPI comm
1715: . dim - The spatial dimension
1716: . Nc - The number of components
1717: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1718: . prefix - The options prefix, or NULL
1719: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1721: Output Parameter:
1722: . fem - The PetscFE object
1724: Note:
1725: Each object is SetFromOption() during creation, so that the object may be customized from the command line.
1727: Level: beginner
1729: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1730: @*/
1731: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1732: {
1733: PetscQuadrature q, fq;
1734: DM K;
1735: PetscSpace P;
1736: PetscDualSpace Q;
1737: PetscInt order, quadPointsPerEdge;
1738: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1739: PetscErrorCode ierr;
1742: /* Create space */
1743: PetscSpaceCreate(comm, &P);
1744: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1745: PetscSpacePolynomialSetTensor(P, tensor);
1746: PetscSpaceSetNumComponents(P, Nc);
1747: PetscSpaceSetNumVariables(P, dim);
1748: PetscSpaceSetFromOptions(P);
1749: PetscSpaceSetUp(P);
1750: PetscSpaceGetDegree(P, &order, NULL);
1751: PetscSpacePolynomialGetTensor(P, &tensor);
1752: /* Create dual space */
1753: PetscDualSpaceCreate(comm, &Q);
1754: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1755: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1756: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1757: PetscDualSpaceSetDM(Q, K);
1758: DMDestroy(&K);
1759: PetscDualSpaceSetNumComponents(Q, Nc);
1760: PetscDualSpaceSetOrder(Q, order);
1761: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1762: PetscDualSpaceSetFromOptions(Q);
1763: PetscDualSpaceSetUp(Q);
1764: /* Create element */
1765: PetscFECreate(comm, fem);
1766: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1767: PetscFESetBasisSpace(*fem, P);
1768: PetscFESetDualSpace(*fem, Q);
1769: PetscFESetNumComponents(*fem, Nc);
1770: PetscFESetFromOptions(*fem);
1771: PetscFESetUp(*fem);
1772: PetscSpaceDestroy(&P);
1773: PetscDualSpaceDestroy(&Q);
1774: /* Create quadrature (with specified order if given) */
1775: qorder = qorder >= 0 ? qorder : order;
1776: PetscObjectOptionsBegin((PetscObject)*fem);
1777: PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1778: PetscOptionsEnd();
1779: quadPointsPerEdge = PetscMax(qorder + 1,1);
1780: if (isSimplex) {
1781: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1782: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1783: } else {
1784: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1785: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1786: }
1787: PetscFESetQuadrature(*fem, q);
1788: PetscFESetFaceQuadrature(*fem, fq);
1789: PetscQuadratureDestroy(&q);
1790: PetscQuadratureDestroy(&fq);
1791: return(0);
1792: }
1794: /*@
1795: PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k
1797: Collective
1799: Input Parameters:
1800: + comm - The MPI comm
1801: . dim - The spatial dimension
1802: . Nc - The number of components
1803: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1804: . k - The degree k of the space
1805: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1807: Output Parameter:
1808: . fem - The PetscFE object
1810: Level: beginner
1812: Notes:
1813: For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k.
1815: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1816: @*/
1817: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1818: {
1819: PetscQuadrature q, fq;
1820: DM K;
1821: PetscSpace P;
1822: PetscDualSpace Q;
1823: PetscInt quadPointsPerEdge;
1824: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1825: char name[64];
1826: PetscErrorCode ierr;
1829: /* Create space */
1830: PetscSpaceCreate(comm, &P);
1831: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1832: PetscSpacePolynomialSetTensor(P, tensor);
1833: PetscSpaceSetNumComponents(P, Nc);
1834: PetscSpaceSetNumVariables(P, dim);
1835: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1836: PetscSpaceSetUp(P);
1837: /* Create dual space */
1838: PetscDualSpaceCreate(comm, &Q);
1839: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1840: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1841: PetscDualSpaceSetDM(Q, K);
1842: DMDestroy(&K);
1843: PetscDualSpaceSetNumComponents(Q, Nc);
1844: PetscDualSpaceSetOrder(Q, k);
1845: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1846: PetscDualSpaceSetUp(Q);
1847: /* Create element */
1848: PetscFECreate(comm, fem);
1849: PetscSNPrintf(name, 64, "P%d", (int) k);
1850: PetscObjectSetName((PetscObject) *fem, name);
1851: PetscFESetType(*fem, PETSCFEBASIC);
1852: PetscFESetBasisSpace(*fem, P);
1853: PetscFESetDualSpace(*fem, Q);
1854: PetscFESetNumComponents(*fem, Nc);
1855: PetscFESetUp(*fem);
1856: PetscSpaceDestroy(&P);
1857: PetscDualSpaceDestroy(&Q);
1858: /* Create quadrature (with specified order if given) */
1859: qorder = qorder >= 0 ? qorder : k;
1860: quadPointsPerEdge = PetscMax(qorder + 1,1);
1861: if (isSimplex) {
1862: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1863: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1864: } else {
1865: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1866: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1867: }
1868: PetscFESetQuadrature(*fem, q);
1869: PetscFESetFaceQuadrature(*fem, fq);
1870: PetscQuadratureDestroy(&q);
1871: PetscQuadratureDestroy(&fq);
1872: return(0);
1873: }
1875: /*@C
1876: PetscFESetName - Names the FE and its subobjects
1878: Not collective
1880: Input Parameters:
1881: + fe - The PetscFE
1882: - name - The name
1884: Level: intermediate
1886: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1887: @*/
1888: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1889: {
1890: PetscSpace P;
1891: PetscDualSpace Q;
1895: PetscFEGetBasisSpace(fe, &P);
1896: PetscFEGetDualSpace(fe, &Q);
1897: PetscObjectSetName((PetscObject) fe, name);
1898: PetscObjectSetName((PetscObject) P, name);
1899: PetscObjectSetName((PetscObject) Q, name);
1900: return(0);
1901: }
1903: PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1904: {
1905: PetscInt dOffset = 0, fOffset = 0, f;
1908: for (f = 0; f < Nf; ++f) {
1909: PetscFE fe;
1910: const PetscInt cdim = T[f]->cdim;
1911: const PetscInt Nq = T[f]->Np;
1912: const PetscInt Nbf = T[f]->Nb;
1913: const PetscInt Ncf = T[f]->Nc;
1914: const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
1915: const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
1916: PetscInt b, c, d;
1918: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
1919: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
1920: for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
1921: for (b = 0; b < Nbf; ++b) {
1922: for (c = 0; c < Ncf; ++c) {
1923: const PetscInt cidx = b*Ncf+c;
1925: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
1926: for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
1927: }
1928: }
1929: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
1930: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
1931: if (u_t) {
1932: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
1933: for (b = 0; b < Nbf; ++b) {
1934: for (c = 0; c < Ncf; ++c) {
1935: const PetscInt cidx = b*Ncf+c;
1937: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
1938: }
1939: }
1940: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
1941: }
1942: fOffset += Ncf;
1943: dOffset += Nbf;
1944: }
1945: return 0;
1946: }
1948: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
1949: {
1950: PetscFE fe;
1951: PetscTabulation Tc;
1952: PetscInt b, c;
1953: PetscErrorCode ierr;
1955: if (!prob) return 0;
1956: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1957: PetscFEGetFaceCentroidTabulation(fe, &Tc);
1958: {
1959: const PetscReal *faceBasis = Tc->T[0];
1960: const PetscInt Nb = Tc->Nb;
1961: const PetscInt Nc = Tc->Nc;
1963: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
1964: for (b = 0; b < Nb; ++b) {
1965: for (c = 0; c < Nc; ++c) {
1966: const PetscInt cidx = b*Nc+c;
1968: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
1969: }
1970: }
1971: }
1972: return 0;
1973: }
1975: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
1976: {
1977: const PetscInt dim = T->cdim;
1978: const PetscInt Nq = T->Np;
1979: const PetscInt Nb = T->Nb;
1980: const PetscInt Nc = T->Nc;
1981: const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
1982: const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dim];
1983: PetscInt q, b, c, d;
1984: PetscErrorCode ierr;
1986: for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
1987: for (q = 0; q < Nq; ++q) {
1988: for (b = 0; b < Nb; ++b) {
1989: for (c = 0; c < Nc; ++c) {
1990: const PetscInt bcidx = b*Nc+c;
1992: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
1993: for (d = 0; d < dim; ++d) tmpBasisDer[bcidx*dim+d] = basisDer[q*Nb*Nc*dim+bcidx*dim+d];
1994: }
1995: }
1996: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
1997: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
1998: for (b = 0; b < Nb; ++b) {
1999: for (c = 0; c < Nc; ++c) {
2000: const PetscInt bcidx = b*Nc+c;
2001: const PetscInt qcidx = q*Nc+c;
2003: elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2004: for (d = 0; d < dim; ++d) elemVec[b] += tmpBasisDer[bcidx*dim+d]*f1[qcidx*dim+d];
2005: }
2006: }
2007: }
2008: return(0);
2009: }
2011: PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2012: {
2013: const PetscInt dim = TI->cdim;
2014: const PetscInt NqI = TI->Np;
2015: const PetscInt NbI = TI->Nb;
2016: const PetscInt NcI = TI->Nc;
2017: const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2018: const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dim];
2019: const PetscInt NqJ = TJ->Np;
2020: const PetscInt NbJ = TJ->Nb;
2021: const PetscInt NcJ = TJ->Nc;
2022: const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2023: const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dim];
2024: PetscInt f, fc, g, gc, df, dg;
2025: PetscErrorCode ierr;
2027: for (f = 0; f < NbI; ++f) {
2028: for (fc = 0; fc < NcI; ++fc) {
2029: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2031: tmpBasisI[fidx] = basisI[fidx];
2032: for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dim+df] = basisDerI[fidx*dim+df];
2033: }
2034: }
2035: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2036: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2037: for (g = 0; g < NbJ; ++g) {
2038: for (gc = 0; gc < NcJ; ++gc) {
2039: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2041: tmpBasisJ[gidx] = basisJ[gidx];
2042: for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dim+dg] = basisDerJ[gidx*dim+dg];
2043: }
2044: }
2045: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2046: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2047: for (f = 0; f < NbI; ++f) {
2048: for (fc = 0; fc < NcI; ++fc) {
2049: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2050: const PetscInt i = offsetI+f; /* Element matrix row */
2051: for (g = 0; g < NbJ; ++g) {
2052: for (gc = 0; gc < NcJ; ++gc) {
2053: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2054: const PetscInt j = offsetJ+g; /* Element matrix column */
2055: const PetscInt fOff = eOffset+i*totDim+j;
2057: elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2058: for (df = 0; df < dim; ++df) {
2059: elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dim+df]*tmpBasisDerJ[gidx*dim+df];
2060: elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g2[(fc*NcJ+gc)*dim+df]*tmpBasisJ[gidx];
2061: for (dg = 0; dg < dim; ++dg) {
2062: elemMat[fOff] += tmpBasisDerI[fidx*dim+df]*g3[((fc*NcJ+gc)*dim+df)*dim+dg]*tmpBasisDerJ[gidx*dim+dg];
2063: }
2064: }
2065: }
2066: }
2067: }
2068: }
2069: return(0);
2070: }
2072: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2073: {
2074: PetscDualSpace dsp;
2075: DM dm;
2076: PetscQuadrature quadDef;
2077: PetscInt dim, cdim, Nq;
2078: PetscErrorCode ierr;
2081: PetscFEGetDualSpace(fe, &dsp);
2082: PetscDualSpaceGetDM(dsp, &dm);
2083: DMGetDimension(dm, &dim);
2084: DMGetCoordinateDim(dm, &cdim);
2085: PetscFEGetQuadrature(fe, &quadDef);
2086: quad = quad ? quad : quadDef;
2087: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2088: PetscMalloc1(Nq*cdim, &cgeom->v);
2089: PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2090: PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2091: PetscMalloc1(Nq, &cgeom->detJ);
2092: cgeom->dim = dim;
2093: cgeom->dimEmbed = cdim;
2094: cgeom->numCells = 1;
2095: cgeom->numPoints = Nq;
2096: DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2097: return(0);
2098: }
2100: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2101: {
2105: PetscFree(cgeom->v);
2106: PetscFree(cgeom->J);
2107: PetscFree(cgeom->invJ);
2108: PetscFree(cgeom->detJ);
2109: return(0);
2110: }