Actual source code: fe.c
1: /* Basis Jet Tabulation
3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6: as a prime basis.
8: \psi_i = \sum_k \alpha_{ki} \phi_k
10: Our nodal basis is defined in terms of the dual basis $n_j$
12: n_j \cdot \psi_i = \delta_{ji}
14: and we may act on the first equation to obtain
16: n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17: \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18: I = V \alpha
20: so the coefficients of the nodal basis in the prime basis are
22: \alpha = V^{-1}
24: We will define the dual basis vectors $n_j$ using a quadrature rule.
26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28: be implemented exactly as in FIAT using functionals $L_j$.
30: I will have to count the degrees correctly for the Legendre product when we are on simplices.
32: We will have three objects:
33: - Space, P: this just need point evaluation I think
34: - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35: - FEM: This keeps {P, P', Q}
36: */
37: #include <petsc/private/petscfeimpl.h>
38: #include <petscdmplex.h>
40: PetscBool FEcite = PETSC_FALSE;
41: const char FECitation[] = "@article{kirby2004,\n"
42: " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
43: " journal = {ACM Transactions on Mathematical Software},\n"
44: " author = {Robert C. Kirby},\n"
45: " volume = {30},\n"
46: " number = {4},\n"
47: " pages = {502--516},\n"
48: " doi = {10.1145/1039813.1039820},\n"
49: " year = {2004}\n}\n";
51: PetscClassId PETSCFE_CLASSID = 0;
53: PetscLogEvent PETSCFE_SetUp;
55: PetscFunctionList PetscFEList = NULL;
56: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
58: /*@C
59: PetscFERegister - Adds a new PetscFE implementation
61: Not Collective
63: Input Parameters:
64: + name - The name of a new user-defined creation routine
65: - create_func - The creation routine itself
67: Notes:
68: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
70: Sample usage:
71: .vb
72: PetscFERegister("my_fe", MyPetscFECreate);
73: .ve
75: Then, your PetscFE type can be chosen with the procedural interface via
76: .vb
77: PetscFECreate(MPI_Comm, PetscFE *);
78: PetscFESetType(PetscFE, "my_fe");
79: .ve
80: or at runtime via the option
81: .vb
82: -petscfe_type my_fe
83: .ve
85: Level: advanced
87: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
89: @*/
90: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
91: {
95: PetscFunctionListAdd(&PetscFEList, sname, function);
96: return(0);
97: }
99: /*@C
100: PetscFESetType - Builds a particular PetscFE
102: Collective on fem
104: Input Parameters:
105: + fem - The PetscFE object
106: - name - The kind of FEM space
108: Options Database Key:
109: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
111: Level: intermediate
113: .seealso: PetscFEGetType(), PetscFECreate()
114: @*/
115: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
116: {
117: PetscErrorCode (*r)(PetscFE);
118: PetscBool match;
123: PetscObjectTypeCompare((PetscObject) fem, name, &match);
124: if (match) return(0);
126: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
127: PetscFunctionListFind(PetscFEList, name, &r);
128: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
130: if (fem->ops->destroy) {
131: (*fem->ops->destroy)(fem);
132: fem->ops->destroy = NULL;
133: }
134: (*r)(fem);
135: PetscObjectChangeTypeName((PetscObject) fem, name);
136: return(0);
137: }
139: /*@C
140: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
142: Not Collective
144: Input Parameter:
145: . fem - The PetscFE
147: Output Parameter:
148: . name - The PetscFE type name
150: Level: intermediate
152: .seealso: PetscFESetType(), PetscFECreate()
153: @*/
154: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
155: {
161: if (!PetscFERegisterAllCalled) {
162: PetscFERegisterAll();
163: }
164: *name = ((PetscObject) fem)->type_name;
165: return(0);
166: }
168: /*@C
169: PetscFEViewFromOptions - View from Options
171: Collective on PetscFE
173: Input Parameters:
174: + A - the PetscFE object
175: . obj - Optional object
176: - name - command line option
178: Level: intermediate
179: .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
180: @*/
181: PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
182: {
187: PetscObjectViewFromOptions((PetscObject)A,obj,name);
188: return(0);
189: }
191: /*@C
192: PetscFEView - Views a PetscFE
194: Collective on fem
196: Input Parameters:
197: + fem - the PetscFE object to view
198: - viewer - the viewer
200: Level: beginner
202: .seealso PetscFEDestroy()
203: @*/
204: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
205: {
206: PetscBool iascii;
212: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);}
213: PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);
214: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
215: if (fem->ops->view) {(*fem->ops->view)(fem, viewer);}
216: return(0);
217: }
219: /*@
220: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
222: Collective on fem
224: Input Parameter:
225: . fem - the PetscFE object to set options for
227: Options Database:
228: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
229: - -petscfe_num_batches - the number of cell batches to integrate serially
231: Level: intermediate
233: .seealso PetscFEView()
234: @*/
235: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
236: {
237: const char *defaultType;
238: char name[256];
239: PetscBool flg;
244: if (!((PetscObject) fem)->type_name) {
245: defaultType = PETSCFEBASIC;
246: } else {
247: defaultType = ((PetscObject) fem)->type_name;
248: }
249: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
251: PetscObjectOptionsBegin((PetscObject) fem);
252: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
253: if (flg) {
254: PetscFESetType(fem, name);
255: } else if (!((PetscObject) fem)->type_name) {
256: PetscFESetType(fem, defaultType);
257: }
258: PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);
259: PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);
260: if (fem->ops->setfromoptions) {
261: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
262: }
263: /* process any options handlers added with PetscObjectAddOptionsHandler() */
264: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
265: PetscOptionsEnd();
266: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
267: return(0);
268: }
270: /*@C
271: PetscFESetUp - Construct data structures for the PetscFE
273: Collective on fem
275: Input Parameter:
276: . fem - the PetscFE object to setup
278: Level: intermediate
280: .seealso PetscFEView(), PetscFEDestroy()
281: @*/
282: PetscErrorCode PetscFESetUp(PetscFE fem)
283: {
288: if (fem->setupcalled) return(0);
289: PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);
290: fem->setupcalled = PETSC_TRUE;
291: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
292: PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);
293: return(0);
294: }
296: /*@
297: PetscFEDestroy - Destroys a PetscFE object
299: Collective on fem
301: Input Parameter:
302: . fem - the PetscFE object to destroy
304: Level: beginner
306: .seealso PetscFEView()
307: @*/
308: PetscErrorCode PetscFEDestroy(PetscFE *fem)
309: {
313: if (!*fem) return(0);
316: if (--((PetscObject)(*fem))->refct > 0) {*fem = NULL; return(0);}
317: ((PetscObject) (*fem))->refct = 0;
319: if ((*fem)->subspaces) {
320: PetscInt dim, d;
322: PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);
323: for (d = 0; d < dim; ++d) {PetscFEDestroy(&(*fem)->subspaces[d]);}
324: }
325: PetscFree((*fem)->subspaces);
326: PetscFree((*fem)->invV);
327: PetscTabulationDestroy(&(*fem)->T);
328: PetscTabulationDestroy(&(*fem)->Tf);
329: PetscTabulationDestroy(&(*fem)->Tc);
330: PetscSpaceDestroy(&(*fem)->basisSpace);
331: PetscDualSpaceDestroy(&(*fem)->dualSpace);
332: PetscQuadratureDestroy(&(*fem)->quadrature);
333: PetscQuadratureDestroy(&(*fem)->faceQuadrature);
334: #ifdef PETSC_HAVE_LIBCEED
335: CeedBasisDestroy(&(*fem)->ceedBasis);
336: CeedDestroy(&(*fem)->ceed);
337: #endif
339: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
340: PetscHeaderDestroy(fem);
341: return(0);
342: }
344: /*@
345: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
347: Collective
349: Input Parameter:
350: . comm - The communicator for the PetscFE object
352: Output Parameter:
353: . fem - The PetscFE object
355: Level: beginner
357: .seealso: PetscFESetType(), PETSCFEGALERKIN
358: @*/
359: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
360: {
361: PetscFE f;
366: PetscCitationsRegister(FECitation,&FEcite);
367: *fem = NULL;
368: PetscFEInitializePackage();
370: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
372: f->basisSpace = NULL;
373: f->dualSpace = NULL;
374: f->numComponents = 1;
375: f->subspaces = NULL;
376: f->invV = NULL;
377: f->T = NULL;
378: f->Tf = NULL;
379: f->Tc = NULL;
380: PetscArrayzero(&f->quadrature, 1);
381: PetscArrayzero(&f->faceQuadrature, 1);
382: f->blockSize = 0;
383: f->numBlocks = 1;
384: f->batchSize = 0;
385: f->numBatches = 1;
387: *fem = f;
388: return(0);
389: }
391: /*@
392: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
394: Not collective
396: Input Parameter:
397: . fem - The PetscFE object
399: Output Parameter:
400: . dim - The spatial dimension
402: Level: intermediate
404: .seealso: PetscFECreate()
405: @*/
406: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
407: {
408: DM dm;
414: PetscDualSpaceGetDM(fem->dualSpace, &dm);
415: DMGetDimension(dm, dim);
416: return(0);
417: }
419: /*@
420: PetscFESetNumComponents - Sets the number of components in the element
422: Not collective
424: Input Parameters:
425: + fem - The PetscFE object
426: - comp - The number of field components
428: Level: intermediate
430: .seealso: PetscFECreate()
431: @*/
432: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
433: {
436: fem->numComponents = comp;
437: return(0);
438: }
440: /*@
441: PetscFEGetNumComponents - Returns the number of components in the element
443: Not collective
445: Input Parameter:
446: . fem - The PetscFE object
448: Output Parameter:
449: . comp - The number of field components
451: Level: intermediate
453: .seealso: PetscFECreate()
454: @*/
455: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
456: {
460: *comp = fem->numComponents;
461: return(0);
462: }
464: /*@
465: PetscFESetTileSizes - Sets the tile sizes for evaluation
467: Not collective
469: Input Parameters:
470: + fem - The PetscFE object
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 PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
481: {
484: fem->blockSize = blockSize;
485: fem->numBlocks = numBlocks;
486: fem->batchSize = batchSize;
487: fem->numBatches = numBatches;
488: return(0);
489: }
491: /*@
492: PetscFEGetTileSizes - Returns the tile sizes for evaluation
494: Not collective
496: Input Parameter:
497: . fem - The PetscFE object
499: Output Parameters:
500: + blockSize - The number of elements in a block
501: . numBlocks - The number of blocks in a batch
502: . batchSize - The number of elements in a batch
503: - numBatches - The number of batches in a chunk
505: Level: intermediate
507: .seealso: PetscFECreate()
508: @*/
509: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
510: {
517: if (blockSize) *blockSize = fem->blockSize;
518: if (numBlocks) *numBlocks = fem->numBlocks;
519: if (batchSize) *batchSize = fem->batchSize;
520: if (numBatches) *numBatches = fem->numBatches;
521: return(0);
522: }
524: /*@
525: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
527: Not collective
529: Input Parameter:
530: . fem - The PetscFE object
532: Output Parameter:
533: . sp - The PetscSpace object
535: Level: intermediate
537: .seealso: PetscFECreate()
538: @*/
539: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
540: {
544: *sp = fem->basisSpace;
545: return(0);
546: }
548: /*@
549: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
551: Not collective
553: Input Parameters:
554: + fem - The PetscFE object
555: - sp - The PetscSpace object
557: Level: intermediate
559: .seealso: PetscFECreate()
560: @*/
561: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
562: {
568: PetscSpaceDestroy(&fem->basisSpace);
569: fem->basisSpace = sp;
570: PetscObjectReference((PetscObject) fem->basisSpace);
571: return(0);
572: }
574: /*@
575: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
577: Not collective
579: Input Parameter:
580: . fem - The PetscFE object
582: Output Parameter:
583: . sp - The PetscDualSpace object
585: Level: intermediate
587: .seealso: PetscFECreate()
588: @*/
589: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
590: {
594: *sp = fem->dualSpace;
595: return(0);
596: }
598: /*@
599: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
601: Not collective
603: Input Parameters:
604: + fem - The PetscFE object
605: - sp - The PetscDualSpace object
607: Level: intermediate
609: .seealso: PetscFECreate()
610: @*/
611: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
612: {
618: PetscDualSpaceDestroy(&fem->dualSpace);
619: fem->dualSpace = sp;
620: PetscObjectReference((PetscObject) fem->dualSpace);
621: return(0);
622: }
624: /*@
625: PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
627: Not collective
629: Input Parameter:
630: . fem - The PetscFE object
632: Output Parameter:
633: . q - The PetscQuadrature object
635: Level: intermediate
637: .seealso: PetscFECreate()
638: @*/
639: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
640: {
644: *q = fem->quadrature;
645: return(0);
646: }
648: /*@
649: PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
651: Not collective
653: Input Parameters:
654: + fem - The PetscFE object
655: - q - The PetscQuadrature object
657: Level: intermediate
659: .seealso: PetscFECreate()
660: @*/
661: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
662: {
663: PetscInt Nc, qNc;
668: if (q == fem->quadrature) return(0);
669: PetscFEGetNumComponents(fem, &Nc);
670: PetscQuadratureGetNumComponents(q, &qNc);
671: 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);
672: PetscTabulationDestroy(&fem->T);
673: PetscTabulationDestroy(&fem->Tc);
674: PetscObjectReference((PetscObject) q);
675: PetscQuadratureDestroy(&fem->quadrature);
676: fem->quadrature = q;
677: return(0);
678: }
680: /*@
681: PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
683: Not collective
685: Input Parameter:
686: . fem - The PetscFE object
688: Output Parameter:
689: . q - The PetscQuadrature object
691: Level: intermediate
693: .seealso: PetscFECreate()
694: @*/
695: PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
696: {
700: *q = fem->faceQuadrature;
701: return(0);
702: }
704: /*@
705: PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
707: Not collective
709: Input Parameters:
710: + fem - The PetscFE object
711: - q - The PetscQuadrature object
713: Level: intermediate
715: .seealso: PetscFECreate()
716: @*/
717: PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
718: {
719: PetscInt Nc, qNc;
724: PetscFEGetNumComponents(fem, &Nc);
725: PetscQuadratureGetNumComponents(q, &qNc);
726: 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);
727: PetscTabulationDestroy(&fem->Tf);
728: PetscQuadratureDestroy(&fem->faceQuadrature);
729: fem->faceQuadrature = q;
730: PetscObjectReference((PetscObject) q);
731: return(0);
732: }
734: /*@
735: PetscFECopyQuadrature - Copy both volumetric and surface quadrature
737: Not collective
739: Input Parameters:
740: + sfe - The PetscFE source for the quadratures
741: - tfe - The PetscFE target for the quadratures
743: Level: intermediate
745: .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
746: @*/
747: PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
748: {
749: PetscQuadrature q;
750: PetscErrorCode ierr;
755: PetscFEGetQuadrature(sfe, &q);
756: PetscFESetQuadrature(tfe, q);
757: PetscFEGetFaceQuadrature(sfe, &q);
758: PetscFESetFaceQuadrature(tfe, q);
759: return(0);
760: }
762: /*@C
763: PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
765: Not collective
767: Input Parameter:
768: . fem - The PetscFE object
770: Output Parameter:
771: . numDof - Array with the number of dofs per dimension
773: Level: intermediate
775: .seealso: PetscFECreate()
776: @*/
777: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
778: {
784: PetscDualSpaceGetNumDof(fem->dualSpace, numDof);
785: return(0);
786: }
788: /*@C
789: PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
791: Not collective
793: Input Parameters:
794: + fem - The PetscFE object
795: - k - The highest derivative we need to tabulate, very often 1
797: Output Parameter:
798: . T - The basis function values and derivatives at quadrature points
800: Note:
801: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
802: $ 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
803: $ 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
805: Level: intermediate
807: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
808: @*/
809: PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T)
810: {
811: PetscInt npoints;
812: const PetscReal *points;
813: PetscErrorCode ierr;
818: PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);
819: if (!fem->T) {PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);}
820: if (fem->T && k > fem->T->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->T->K);
821: *T = fem->T;
822: return(0);
823: }
825: /*@C
826: PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
828: Not collective
830: Input Parameters:
831: + fem - The PetscFE object
832: - k - The highest derivative we need to tabulate, very often 1
834: Output Parameters:
835: . Tf - The basis function values and derivatives at face quadrature points
837: Note:
838: $ 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
839: $ 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
840: $ 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
842: Level: intermediate
844: .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
845: @*/
846: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf)
847: {
848: PetscErrorCode ierr;
853: if (!fem->Tf) {
854: const PetscReal xi0[3] = {-1., -1., -1.};
855: PetscReal v0[3], J[9], detJ;
856: PetscQuadrature fq;
857: PetscDualSpace sp;
858: DM dm;
859: const PetscInt *faces;
860: PetscInt dim, numFaces, f, npoints, q;
861: const PetscReal *points;
862: PetscReal *facePoints;
864: PetscFEGetDualSpace(fem, &sp);
865: PetscDualSpaceGetDM(sp, &dm);
866: DMGetDimension(dm, &dim);
867: DMPlexGetConeSize(dm, 0, &numFaces);
868: DMPlexGetCone(dm, 0, &faces);
869: PetscFEGetFaceQuadrature(fem, &fq);
870: if (fq) {
871: PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);
872: PetscMalloc1(numFaces*npoints*dim, &facePoints);
873: for (f = 0; f < numFaces; ++f) {
874: DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);
875: for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
876: }
877: PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);
878: PetscFree(facePoints);
879: }
880: }
881: if (fem->Tf && k > fem->Tf->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->Tf->K);
882: *Tf = fem->Tf;
883: return(0);
884: }
886: /*@C
887: PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
889: Not collective
891: Input Parameter:
892: . fem - The PetscFE object
894: Output Parameters:
895: . Tc - The basis function values at face centroid points
897: Note:
898: $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
900: Level: intermediate
902: .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
903: @*/
904: PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
905: {
906: PetscErrorCode ierr;
911: if (!fem->Tc) {
912: PetscDualSpace sp;
913: DM dm;
914: const PetscInt *cone;
915: PetscReal *centroids;
916: PetscInt dim, numFaces, f;
918: PetscFEGetDualSpace(fem, &sp);
919: PetscDualSpaceGetDM(sp, &dm);
920: DMGetDimension(dm, &dim);
921: DMPlexGetConeSize(dm, 0, &numFaces);
922: DMPlexGetCone(dm, 0, &cone);
923: PetscMalloc1(numFaces*dim, ¢roids);
924: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
925: PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);
926: PetscFree(centroids);
927: }
928: *Tc = fem->Tc;
929: return(0);
930: }
932: /*@C
933: PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
935: Not collective
937: Input Parameters:
938: + fem - The PetscFE object
939: . nrepl - The number of replicas
940: . npoints - The number of tabulation points in a replica
941: . points - The tabulation point coordinates
942: - K - The number of derivatives calculated
944: Output Parameter:
945: . T - The basis function values and derivatives at tabulation points
947: Note:
948: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
949: $ 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
950: $ 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
952: Level: intermediate
954: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
955: @*/
956: PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
957: {
958: DM dm;
959: PetscDualSpace Q;
960: PetscInt Nb; /* Dimension of FE space P */
961: PetscInt Nc; /* Field components */
962: PetscInt cdim; /* Reference coordinate dimension */
963: PetscInt k;
964: PetscErrorCode ierr;
967: if (!npoints || !fem->dualSpace || K < 0) {
968: *T = NULL;
969: return(0);
970: }
974: PetscFEGetDualSpace(fem, &Q);
975: PetscDualSpaceGetDM(Q, &dm);
976: DMGetDimension(dm, &cdim);
977: PetscDualSpaceGetDimension(Q, &Nb);
978: PetscFEGetNumComponents(fem, &Nc);
979: PetscMalloc1(1, T);
980: (*T)->K = !cdim ? 0 : K;
981: (*T)->Nr = nrepl;
982: (*T)->Np = npoints;
983: (*T)->Nb = Nb;
984: (*T)->Nc = Nc;
985: (*T)->cdim = cdim;
986: PetscMalloc1((*T)->K+1, &(*T)->T);
987: for (k = 0; k <= (*T)->K; ++k) {
988: PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
989: }
990: (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);
991: return(0);
992: }
994: /*@C
995: PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
997: Not collective
999: Input Parameters:
1000: + fem - The PetscFE object
1001: . npoints - The number of tabulation points
1002: . points - The tabulation point coordinates
1003: . K - The number of derivatives calculated
1004: - T - An existing tabulation object with enough allocated space
1006: Output Parameter:
1007: . T - The basis function values and derivatives at tabulation points
1009: Note:
1010: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1011: $ 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
1012: $ 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
1014: Level: intermediate
1016: .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1017: @*/
1018: PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1019: {
1023: if (!npoints || !fem->dualSpace || K < 0) return(0);
1027: if (PetscDefined(USE_DEBUG)) {
1028: DM dm;
1029: PetscDualSpace Q;
1030: PetscInt Nb; /* Dimension of FE space P */
1031: PetscInt Nc; /* Field components */
1032: PetscInt cdim; /* Reference coordinate dimension */
1034: PetscFEGetDualSpace(fem, &Q);
1035: PetscDualSpaceGetDM(Q, &dm);
1036: DMGetDimension(dm, &cdim);
1037: PetscDualSpaceGetDimension(Q, &Nb);
1038: PetscFEGetNumComponents(fem, &Nc);
1039: 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);
1040: if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1041: if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1042: if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1043: }
1044: T->Nr = 1;
1045: T->Np = npoints;
1046: (*fem->ops->createtabulation)(fem, npoints, points, K, T);
1047: return(0);
1048: }
1050: /*@C
1051: PetscTabulationDestroy - Frees memory from the associated tabulation.
1053: Not collective
1055: Input Parameter:
1056: . T - The tabulation
1058: Level: intermediate
1060: .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1061: @*/
1062: PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1063: {
1064: PetscInt k;
1069: if (!T || !(*T)) return(0);
1070: for (k = 0; k <= (*T)->K; ++k) {PetscFree((*T)->T[k]);}
1071: PetscFree((*T)->T);
1072: PetscFree(*T);
1073: *T = NULL;
1074: return(0);
1075: }
1077: PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1078: {
1079: PetscSpace bsp, bsubsp;
1080: PetscDualSpace dsp, dsubsp;
1081: PetscInt dim, depth, numComp, i, j, coneSize, order;
1082: PetscFEType type;
1083: DM dm;
1084: DMLabel label;
1085: PetscReal *xi, *v, *J, detJ;
1086: const char *name;
1087: PetscQuadrature origin, fullQuad, subQuad;
1093: PetscFEGetBasisSpace(fe,&bsp);
1094: PetscFEGetDualSpace(fe,&dsp);
1095: PetscDualSpaceGetDM(dsp,&dm);
1096: DMGetDimension(dm,&dim);
1097: DMPlexGetDepthLabel(dm,&label);
1098: DMLabelGetValue(label,refPoint,&depth);
1099: PetscCalloc1(depth,&xi);
1100: PetscMalloc1(dim,&v);
1101: PetscMalloc1(dim*dim,&J);
1102: for (i = 0; i < depth; i++) xi[i] = 0.;
1103: PetscQuadratureCreate(PETSC_COMM_SELF,&origin);
1104: PetscQuadratureSetData(origin,depth,0,1,xi,NULL);
1105: DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);
1106: /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1107: for (i = 1; i < dim; i++) {
1108: for (j = 0; j < depth; j++) {
1109: J[i * depth + j] = J[i * dim + j];
1110: }
1111: }
1112: PetscQuadratureDestroy(&origin);
1113: PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);
1114: PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);
1115: PetscSpaceSetUp(bsubsp);
1116: PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);
1117: PetscFEGetType(fe,&type);
1118: PetscFESetType(*trFE,type);
1119: PetscFEGetNumComponents(fe,&numComp);
1120: PetscFESetNumComponents(*trFE,numComp);
1121: PetscFESetBasisSpace(*trFE,bsubsp);
1122: PetscFESetDualSpace(*trFE,dsubsp);
1123: PetscObjectGetName((PetscObject) fe, &name);
1124: if (name) {PetscFESetName(*trFE, name);}
1125: PetscFEGetQuadrature(fe,&fullQuad);
1126: PetscQuadratureGetOrder(fullQuad,&order);
1127: DMPlexGetConeSize(dm,refPoint,&coneSize);
1128: if (coneSize == 2 * depth) {
1129: PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1130: } else {
1131: PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);
1132: }
1133: PetscFESetQuadrature(*trFE,subQuad);
1134: PetscFESetUp(*trFE);
1135: PetscQuadratureDestroy(&subQuad);
1136: PetscSpaceDestroy(&bsubsp);
1137: return(0);
1138: }
1140: PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1141: {
1142: PetscInt hStart, hEnd;
1143: PetscDualSpace dsp;
1144: DM dm;
1150: *trFE = NULL;
1151: PetscFEGetDualSpace(fe,&dsp);
1152: PetscDualSpaceGetDM(dsp,&dm);
1153: DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);
1154: if (hEnd <= hStart) return(0);
1155: PetscFECreatePointTrace(fe,hStart,trFE);
1156: return(0);
1157: }
1159: /*@
1160: PetscFEGetDimension - Get the dimension of the finite element space on a cell
1162: Not collective
1164: Input Parameter:
1165: . fe - The PetscFE
1167: Output Parameter:
1168: . dim - The dimension
1170: Level: intermediate
1172: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1173: @*/
1174: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1175: {
1181: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
1182: return(0);
1183: }
1185: /*@C
1186: PetscFEPushforward - Map the reference element function to real space
1188: Input Parameters:
1189: + fe - The PetscFE
1190: . fegeom - The cell geometry
1191: . Nv - The number of function values
1192: - vals - The function values
1194: Output Parameter:
1195: . vals - The transformed function values
1197: Level: advanced
1199: Note: This just forwards the call onto PetscDualSpacePushforward().
1201: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1203: .seealso: PetscDualSpacePushforward()
1204: @*/
1205: PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1206: {
1210: PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1211: return(0);
1212: }
1214: /*@C
1215: PetscFEPushforwardGradient - Map the reference element function gradient to real space
1217: Input Parameters:
1218: + fe - The PetscFE
1219: . fegeom - The cell geometry
1220: . Nv - The number of function gradient values
1221: - vals - The function gradient values
1223: Output Parameter:
1224: . vals - The transformed function gradient values
1226: Level: advanced
1228: Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1230: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1232: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1233: @*/
1234: PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1235: {
1239: PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1240: return(0);
1241: }
1243: /*@C
1244: PetscFEPushforwardHessian - Map the reference element function Hessian to real space
1246: Input Parameters:
1247: + fe - The PetscFE
1248: . fegeom - The cell geometry
1249: . Nv - The number of function Hessian values
1250: - vals - The function Hessian values
1252: Output Parameter:
1253: . vals - The transformed function Hessian values
1255: Level: advanced
1257: Note: This just forwards the call onto PetscDualSpacePushforwardHessian().
1259: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1261: .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward()
1262: @*/
1263: PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1264: {
1268: PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);
1269: return(0);
1270: }
1272: /*
1273: Purpose: Compute element vector for chunk of elements
1275: Input:
1276: Sizes:
1277: Ne: number of elements
1278: Nf: number of fields
1279: PetscFE
1280: dim: spatial dimension
1281: Nb: number of basis functions
1282: Nc: number of field components
1283: PetscQuadrature
1284: Nq: number of quadrature points
1286: Geometry:
1287: PetscFEGeom[Ne] possibly *Nq
1288: PetscReal v0s[dim]
1289: PetscReal n[dim]
1290: PetscReal jacobians[dim*dim]
1291: PetscReal jacobianInverses[dim*dim]
1292: PetscReal jacobianDeterminants
1293: FEM:
1294: PetscFE
1295: PetscQuadrature
1296: PetscReal quadPoints[Nq*dim]
1297: PetscReal quadWeights[Nq]
1298: PetscReal basis[Nq*Nb*Nc]
1299: PetscReal basisDer[Nq*Nb*Nc*dim]
1300: PetscScalar coefficients[Ne*Nb*Nc]
1301: PetscScalar elemVec[Ne*Nb*Nc]
1303: Problem:
1304: PetscInt f: the active field
1305: f0, f1
1307: Work Space:
1308: PetscFE
1309: PetscScalar f0[Nq*dim];
1310: PetscScalar f1[Nq*dim*dim];
1311: PetscScalar u[Nc];
1312: PetscScalar gradU[Nc*dim];
1313: PetscReal x[dim];
1314: PetscScalar realSpaceDer[dim];
1316: Purpose: Compute element vector for N_cb batches of elements
1318: Input:
1319: Sizes:
1320: N_cb: Number of serial cell batches
1322: Geometry:
1323: PetscReal v0s[Ne*dim]
1324: PetscReal jacobians[Ne*dim*dim] possibly *Nq
1325: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1326: PetscReal jacobianDeterminants[Ne] possibly *Nq
1327: FEM:
1328: static PetscReal quadPoints[Nq*dim]
1329: static PetscReal quadWeights[Nq]
1330: static PetscReal basis[Nq*Nb*Nc]
1331: static PetscReal basisDer[Nq*Nb*Nc*dim]
1332: PetscScalar coefficients[Ne*Nb*Nc]
1333: PetscScalar elemVec[Ne*Nb*Nc]
1335: ex62.c:
1336: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1337: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1338: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1339: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1341: ex52.c:
1342: 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)
1343: 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)
1345: ex52_integrateElement.cu
1346: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1348: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1349: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1350: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1352: ex52_integrateElementOpenCL.c:
1353: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1354: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1355: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1357: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1358: */
1360: /*@C
1361: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1363: Not collective
1365: Input Parameters:
1366: + prob - The PetscDS specifying the discretizations and continuum functions
1367: . field - The field being integrated
1368: . Ne - The number of elements in the chunk
1369: . cgeom - The cell geometry for each cell in the chunk
1370: . coefficients - The array of FEM basis coefficients for the elements
1371: . probAux - The PetscDS specifying the auxiliary discretizations
1372: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1374: Output Parameter:
1375: . integral - the integral for this field
1377: Level: intermediate
1379: .seealso: PetscFEIntegrateResidual()
1380: @*/
1381: PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1382: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1383: {
1384: PetscFE fe;
1389: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1390: if (fe->ops->integrate) {(*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);}
1391: return(0);
1392: }
1394: /*@C
1395: PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1397: Not collective
1399: Input Parameters:
1400: + prob - The PetscDS specifying the discretizations and continuum functions
1401: . field - The field being integrated
1402: . obj_func - The function to be integrated
1403: . Ne - The number of elements in the chunk
1404: . fgeom - The face geometry for each face in the chunk
1405: . coefficients - The array of FEM basis coefficients for the elements
1406: . probAux - The PetscDS specifying the auxiliary discretizations
1407: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1409: Output Parameter:
1410: . integral - the integral for this field
1412: Level: intermediate
1414: .seealso: PetscFEIntegrateResidual()
1415: @*/
1416: PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1417: void (*obj_func)(PetscInt, PetscInt, PetscInt,
1418: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1419: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1420: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1421: PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1422: {
1423: PetscFE fe;
1428: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1429: if (fe->ops->integratebd) {(*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
1430: return(0);
1431: }
1433: /*@C
1434: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1436: Not collective
1438: Input Parameters:
1439: + ds - The PetscDS specifying the discretizations and continuum functions
1440: . key - The (label+value, field) being integrated
1441: . Ne - The number of elements in the chunk
1442: . cgeom - The cell geometry for each cell in the chunk
1443: . coefficients - The array of FEM basis coefficients for the elements
1444: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1445: . probAux - The PetscDS specifying the auxiliary discretizations
1446: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1447: - t - The time
1449: Output Parameter:
1450: . elemVec - the element residual vectors from each element
1452: Note:
1453: $ Loop over batch of elements (e):
1454: $ Loop over quadrature points (q):
1455: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1456: $ Call f_0 and f_1
1457: $ Loop over element vector entries (f,fc --> i):
1458: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1460: Level: intermediate
1462: .seealso: PetscFEIntegrateResidual()
1463: @*/
1464: PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1465: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1466: {
1467: PetscFE fe;
1472: PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);
1473: if (fe->ops->integrateresidual) {(*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1474: return(0);
1475: }
1477: /*@C
1478: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1480: Not collective
1482: Input Parameters:
1483: + ds - The PetscDS specifying the discretizations and continuum functions
1484: . wf - The PetscWeakForm object holding the pointwise functions
1485: . key - The (label+value, field) being integrated
1486: . Ne - The number of elements in the chunk
1487: . fgeom - The face geometry for each cell in the chunk
1488: . coefficients - The array of FEM basis coefficients for the elements
1489: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1490: . probAux - The PetscDS specifying the auxiliary discretizations
1491: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1492: - t - The time
1494: Output Parameter:
1495: . elemVec - the element residual vectors from each element
1497: Level: intermediate
1499: .seealso: PetscFEIntegrateResidual()
1500: @*/
1501: PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1502: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1503: {
1504: PetscFE fe;
1509: PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);
1510: if (fe->ops->integratebdresidual) {(*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1511: return(0);
1512: }
1514: /*@C
1515: PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1517: Not collective
1519: Input Parameters:
1520: + prob - The PetscDS specifying the discretizations and continuum functions
1521: . key - The (label+value, field) being integrated
1522: . Ne - The number of elements in the chunk
1523: . fgeom - The face geometry for each cell in the chunk
1524: . coefficients - The array of FEM basis coefficients for the elements
1525: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1526: . probAux - The PetscDS specifying the auxiliary discretizations
1527: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1528: - t - The time
1530: Output Parameter
1531: . elemVec - the element residual vectors from each element
1533: Level: developer
1535: .seealso: PetscFEIntegrateResidual()
1536: @*/
1537: PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1538: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1539: {
1540: PetscFE fe;
1545: PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe);
1546: if (fe->ops->integratehybridresidual) {(*fe->ops->integratehybridresidual)(prob, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
1547: return(0);
1548: }
1550: /*@C
1551: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1553: Not collective
1555: Input Parameters:
1556: + ds - The PetscDS specifying the discretizations and continuum functions
1557: . jtype - The type of matrix pointwise functions that should be used
1558: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1559: . Ne - The number of elements in the chunk
1560: . cgeom - The cell geometry for each cell in the chunk
1561: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1562: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1563: . probAux - The PetscDS specifying the auxiliary discretizations
1564: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1565: . t - The time
1566: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1568: Output Parameter:
1569: . elemMat - the element matrices for the Jacobian from each element
1571: Note:
1572: $ Loop over batch of elements (e):
1573: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1574: $ Loop over quadrature points (q):
1575: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1576: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1577: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1578: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1579: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1580: Level: intermediate
1582: .seealso: PetscFEIntegrateResidual()
1583: @*/
1584: PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
1585: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1586: {
1587: PetscFE fe;
1588: PetscInt Nf;
1593: PetscDSGetNumFields(ds, &Nf);
1594: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1595: if (fe->ops->integratejacobian) {(*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1596: return(0);
1597: }
1599: /*@C
1600: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1602: Not collective
1604: Input Parameters:
1605: + ds - The PetscDS specifying the discretizations and continuum functions
1606: . wf - The PetscWeakForm holding the pointwise functions
1607: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1608: . Ne - The number of elements in the chunk
1609: . fgeom - The face geometry for each cell in the chunk
1610: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1611: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1612: . probAux - The PetscDS specifying the auxiliary discretizations
1613: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1614: . t - The time
1615: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1617: Output Parameter:
1618: . elemMat - the element matrices for the Jacobian from each element
1620: Note:
1621: $ Loop over batch of elements (e):
1622: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1623: $ Loop over quadrature points (q):
1624: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1625: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1626: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1627: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1628: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1629: Level: intermediate
1631: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1632: @*/
1633: PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1634: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1635: {
1636: PetscFE fe;
1637: PetscInt Nf;
1642: PetscDSGetNumFields(ds, &Nf);
1643: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1644: if (fe->ops->integratebdjacobian) {(*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1645: return(0);
1646: }
1648: /*@C
1649: PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1651: Not collective
1653: Input Parameters:
1654: + ds - The PetscDS specifying the discretizations and continuum functions
1655: . jtype - The type of matrix pointwise functions that should be used
1656: . key - The (label+value, fieldI*Nf + fieldJ) being integrated
1657: . Ne - The number of elements in the chunk
1658: . fgeom - The face geometry for each cell in the chunk
1659: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1660: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1661: . probAux - The PetscDS specifying the auxiliary discretizations
1662: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1663: . t - The time
1664: - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1666: Output Parameter
1667: . elemMat - the element matrices for the Jacobian from each element
1669: Note:
1670: $ Loop over batch of elements (e):
1671: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1672: $ Loop over quadrature points (q):
1673: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1674: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1675: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1676: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1677: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1678: Level: developer
1680: .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1681: @*/
1682: PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
1683: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1684: {
1685: PetscFE fe;
1686: PetscInt Nf;
1691: PetscDSGetNumFields(ds, &Nf);
1692: PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);
1693: if (fe->ops->integratehybridjacobian) {(*fe->ops->integratehybridjacobian)(ds, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
1694: return(0);
1695: }
1697: /*@
1698: PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1700: Input Parameters:
1701: + fe - The finite element space
1702: - height - The height of the Plex point
1704: Output Parameter:
1705: . subfe - The subspace of this FE space
1707: Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1709: Level: advanced
1711: .seealso: PetscFECreateDefault()
1712: @*/
1713: PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1714: {
1715: PetscSpace P, subP;
1716: PetscDualSpace Q, subQ;
1717: PetscQuadrature subq;
1718: PetscFEType fetype;
1719: PetscInt dim, Nc;
1720: PetscErrorCode ierr;
1725: if (height == 0) {
1726: *subfe = fe;
1727: return(0);
1728: }
1729: PetscFEGetBasisSpace(fe, &P);
1730: PetscFEGetDualSpace(fe, &Q);
1731: PetscFEGetNumComponents(fe, &Nc);
1732: PetscFEGetFaceQuadrature(fe, &subq);
1733: PetscDualSpaceGetDimension(Q, &dim);
1734: 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);
1735: if (!fe->subspaces) {PetscCalloc1(dim, &fe->subspaces);}
1736: if (height <= dim) {
1737: if (!fe->subspaces[height-1]) {
1738: PetscFE sub = NULL;
1739: const char *name;
1741: PetscSpaceGetHeightSubspace(P, height, &subP);
1742: PetscDualSpaceGetHeightSubspace(Q, height, &subQ);
1743: if (subQ) {
1744: PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);
1745: PetscObjectGetName((PetscObject) fe, &name);
1746: PetscObjectSetName((PetscObject) sub, name);
1747: PetscFEGetType(fe, &fetype);
1748: PetscFESetType(sub, fetype);
1749: PetscFESetBasisSpace(sub, subP);
1750: PetscFESetDualSpace(sub, subQ);
1751: PetscFESetNumComponents(sub, Nc);
1752: PetscFESetUp(sub);
1753: PetscFESetQuadrature(sub, subq);
1754: }
1755: fe->subspaces[height-1] = sub;
1756: }
1757: *subfe = fe->subspaces[height-1];
1758: } else {
1759: *subfe = NULL;
1760: }
1761: return(0);
1762: }
1764: /*@
1765: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1766: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1767: sparsity). It is also used to create an interpolation between regularly refined meshes.
1769: Collective on fem
1771: Input Parameter:
1772: . fe - The initial PetscFE
1774: Output Parameter:
1775: . feRef - The refined PetscFE
1777: Level: advanced
1779: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1780: @*/
1781: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1782: {
1783: PetscSpace P, Pref;
1784: PetscDualSpace Q, Qref;
1785: DM K, Kref;
1786: PetscQuadrature q, qref;
1787: const PetscReal *v0, *jac;
1788: PetscInt numComp, numSubelements;
1789: PetscInt cStart, cEnd, c;
1790: PetscDualSpace *cellSpaces;
1791: PetscErrorCode ierr;
1794: PetscFEGetBasisSpace(fe, &P);
1795: PetscFEGetDualSpace(fe, &Q);
1796: PetscFEGetQuadrature(fe, &q);
1797: PetscDualSpaceGetDM(Q, &K);
1798: /* Create space */
1799: PetscObjectReference((PetscObject) P);
1800: Pref = P;
1801: /* Create dual space */
1802: PetscDualSpaceDuplicate(Q, &Qref);
1803: PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);
1804: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
1805: PetscDualSpaceSetDM(Qref, Kref);
1806: DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);
1807: PetscMalloc1(cEnd - cStart, &cellSpaces);
1808: /* TODO: fix for non-uniform refinement */
1809: for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1810: PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);
1811: PetscFree(cellSpaces);
1812: DMDestroy(&Kref);
1813: PetscDualSpaceSetUp(Qref);
1814: /* Create element */
1815: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
1816: PetscFESetType(*feRef, PETSCFECOMPOSITE);
1817: PetscFESetBasisSpace(*feRef, Pref);
1818: PetscFESetDualSpace(*feRef, Qref);
1819: PetscFEGetNumComponents(fe, &numComp);
1820: PetscFESetNumComponents(*feRef, numComp);
1821: PetscFESetUp(*feRef);
1822: PetscSpaceDestroy(&Pref);
1823: PetscDualSpaceDestroy(&Qref);
1824: /* Create quadrature */
1825: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
1826: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1827: PetscFESetQuadrature(*feRef, qref);
1828: PetscQuadratureDestroy(&qref);
1829: return(0);
1830: }
1832: /*@C
1833: PetscFECreateDefault - Create a PetscFE for basic FEM computation
1835: Collective
1837: Input Parameters:
1838: + comm - The MPI comm
1839: . dim - The spatial dimension
1840: . Nc - The number of components
1841: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1842: . prefix - The options prefix, or NULL
1843: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1845: Output Parameter:
1846: . fem - The PetscFE object
1848: Note:
1849: Each subobject is SetFromOption() during creation, so that the object may be customized from the command line, using the prefix specified above. See the links below for the particular options available.
1851: Level: beginner
1853: .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1854: @*/
1855: PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1856: {
1857: PetscQuadrature q, fq;
1858: DM K;
1859: PetscSpace P;
1860: PetscDualSpace Q;
1861: PetscInt order, quadPointsPerEdge;
1862: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1863: PetscErrorCode ierr;
1866: /* Create space */
1867: PetscSpaceCreate(comm, &P);
1868: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1869: PetscSpacePolynomialSetTensor(P, tensor);
1870: PetscSpaceSetNumComponents(P, Nc);
1871: PetscSpaceSetNumVariables(P, dim);
1872: PetscSpaceSetFromOptions(P);
1873: PetscSpaceSetUp(P);
1874: PetscSpaceGetDegree(P, &order, NULL);
1875: PetscSpacePolynomialGetTensor(P, &tensor);
1876: /* Create dual space */
1877: PetscDualSpaceCreate(comm, &Q);
1878: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1879: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1880: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1881: PetscDualSpaceSetDM(Q, K);
1882: DMDestroy(&K);
1883: PetscDualSpaceSetNumComponents(Q, Nc);
1884: PetscDualSpaceSetOrder(Q, order);
1885: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1886: PetscDualSpaceSetFromOptions(Q);
1887: PetscDualSpaceSetUp(Q);
1888: /* Create element */
1889: PetscFECreate(comm, fem);
1890: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1891: PetscFESetBasisSpace(*fem, P);
1892: PetscFESetDualSpace(*fem, Q);
1893: PetscFESetNumComponents(*fem, Nc);
1894: PetscFESetFromOptions(*fem);
1895: PetscFESetUp(*fem);
1896: PetscSpaceDestroy(&P);
1897: PetscDualSpaceDestroy(&Q);
1898: /* Create quadrature (with specified order if given) */
1899: qorder = qorder >= 0 ? qorder : order;
1900: PetscObjectOptionsBegin((PetscObject)*fem);
1901: PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);
1902: PetscOptionsEnd();
1903: quadPointsPerEdge = PetscMax(qorder + 1,1);
1904: if (isSimplex) {
1905: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1906: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1907: } else {
1908: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1909: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1910: }
1911: PetscFESetQuadrature(*fem, q);
1912: PetscFESetFaceQuadrature(*fem, fq);
1913: PetscQuadratureDestroy(&q);
1914: PetscQuadratureDestroy(&fq);
1915: return(0);
1916: }
1918: /*@
1919: PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k
1921: Collective
1923: Input Parameters:
1924: + comm - The MPI comm
1925: . dim - The spatial dimension
1926: . Nc - The number of components
1927: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1928: . k - The degree k of the space
1929: - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1931: Output Parameter:
1932: . fem - The PetscFE object
1934: Level: beginner
1936: Notes:
1937: 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.
1939: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1940: @*/
1941: PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1942: {
1943: PetscQuadrature q, fq;
1944: DM K;
1945: PetscSpace P;
1946: PetscDualSpace Q;
1947: PetscInt quadPointsPerEdge;
1948: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1949: char name[64];
1950: PetscErrorCode ierr;
1953: /* Create space */
1954: PetscSpaceCreate(comm, &P);
1955: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1956: PetscSpacePolynomialSetTensor(P, tensor);
1957: PetscSpaceSetNumComponents(P, Nc);
1958: PetscSpaceSetNumVariables(P, dim);
1959: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1960: PetscSpaceSetUp(P);
1961: /* Create dual space */
1962: PetscDualSpaceCreate(comm, &Q);
1963: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1964: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1965: PetscDualSpaceSetDM(Q, K);
1966: DMDestroy(&K);
1967: PetscDualSpaceSetNumComponents(Q, Nc);
1968: PetscDualSpaceSetOrder(Q, k);
1969: PetscDualSpaceLagrangeSetTensor(Q, tensor);
1970: PetscDualSpaceSetUp(Q);
1971: /* Create finite element */
1972: PetscFECreate(comm, fem);
1973: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1974: PetscObjectSetName((PetscObject) *fem, name);
1975: PetscFESetType(*fem, PETSCFEBASIC);
1976: PetscFESetBasisSpace(*fem, P);
1977: PetscFESetDualSpace(*fem, Q);
1978: PetscFESetNumComponents(*fem, Nc);
1979: PetscFESetUp(*fem);
1980: PetscSpaceDestroy(&P);
1981: PetscDualSpaceDestroy(&Q);
1982: /* Create quadrature (with specified order if given) */
1983: qorder = qorder >= 0 ? qorder : k;
1984: quadPointsPerEdge = PetscMax(qorder + 1,1);
1985: if (isSimplex) {
1986: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1987: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1988: } else {
1989: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
1990: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
1991: }
1992: PetscFESetQuadrature(*fem, q);
1993: PetscFESetFaceQuadrature(*fem, fq);
1994: PetscQuadratureDestroy(&q);
1995: PetscQuadratureDestroy(&fq);
1996: /* Set finite element name */
1997: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1998: PetscFESetName(*fem, name);
1999: return(0);
2000: }
2002: /*@C
2003: PetscFESetName - Names the FE and its subobjects
2005: Not collective
2007: Input Parameters:
2008: + fe - The PetscFE
2009: - name - The name
2011: Level: intermediate
2013: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
2014: @*/
2015: PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
2016: {
2017: PetscSpace P;
2018: PetscDualSpace Q;
2022: PetscFEGetBasisSpace(fe, &P);
2023: PetscFEGetDualSpace(fe, &Q);
2024: PetscObjectSetName((PetscObject) fe, name);
2025: PetscObjectSetName((PetscObject) P, name);
2026: PetscObjectSetName((PetscObject) Q, name);
2027: return(0);
2028: }
2030: 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[])
2031: {
2032: PetscInt dOffset = 0, fOffset = 0, f, g;
2035: for (f = 0; f < Nf; ++f) {
2036: PetscFE fe;
2037: const PetscInt k = ds->jetDegree[f];
2038: const PetscInt cdim = T[f]->cdim;
2039: const PetscInt Nq = T[f]->Np;
2040: const PetscInt Nbf = T[f]->Nb;
2041: const PetscInt Ncf = T[f]->Nc;
2042: const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2043: const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2044: const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL;
2045: PetscInt hOffset = 0, b, c, d;
2047: PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);
2048: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2049: for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2050: for (b = 0; b < Nbf; ++b) {
2051: for (c = 0; c < Ncf; ++c) {
2052: const PetscInt cidx = b*Ncf+c;
2054: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2055: for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2056: }
2057: }
2058: if (k > 1) {
2059: for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim;
2060: for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0;
2061: for (b = 0; b < Nbf; ++b) {
2062: for (c = 0; c < Ncf; ++c) {
2063: const PetscInt cidx = b*Ncf+c;
2065: for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b];
2066: }
2067: }
2068: PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);
2069: }
2070: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2071: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2072: if (u_t) {
2073: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2074: for (b = 0; b < Nbf; ++b) {
2075: for (c = 0; c < Ncf; ++c) {
2076: const PetscInt cidx = b*Ncf+c;
2078: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2079: }
2080: }
2081: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2082: }
2083: fOffset += Ncf;
2084: dOffset += Nbf;
2085: }
2086: return 0;
2087: }
2089: PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2090: {
2091: PetscInt dOffset = 0, fOffset = 0, g;
2094: for (g = 0; g < 2*Nf-1; ++g) {
2095: if (!T[g/2]) continue;
2096: {
2097: PetscFE fe;
2098: const PetscInt f = g/2;
2099: const PetscInt cdim = T[f]->cdim;
2100: const PetscInt Nq = T[f]->Np;
2101: const PetscInt Nbf = T[f]->Nb;
2102: const PetscInt Ncf = T[f]->Nc;
2103: const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2104: const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2105: PetscInt b, c, d;
2107: fe = (PetscFE) ds->disc[f];
2108: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2109: for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2110: for (b = 0; b < Nbf; ++b) {
2111: for (c = 0; c < Ncf; ++c) {
2112: const PetscInt cidx = b*Ncf+c;
2114: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2115: for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2116: }
2117: }
2118: PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);
2119: PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);
2120: if (u_t) {
2121: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2122: for (b = 0; b < Nbf; ++b) {
2123: for (c = 0; c < Ncf; ++c) {
2124: const PetscInt cidx = b*Ncf+c;
2126: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2127: }
2128: }
2129: PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);
2130: }
2131: fOffset += Ncf;
2132: dOffset += Nbf;
2133: }
2134: }
2135: return 0;
2136: }
2138: PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2139: {
2140: PetscFE fe;
2141: PetscTabulation Tc;
2142: PetscInt b, c;
2143: PetscErrorCode ierr;
2145: if (!prob) return 0;
2146: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2147: PetscFEGetFaceCentroidTabulation(fe, &Tc);
2148: {
2149: const PetscReal *faceBasis = Tc->T[0];
2150: const PetscInt Nb = Tc->Nb;
2151: const PetscInt Nc = Tc->Nc;
2153: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2154: for (b = 0; b < Nb; ++b) {
2155: for (c = 0; c < Nc; ++c) {
2156: u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2157: }
2158: }
2159: }
2160: return 0;
2161: }
2163: PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2164: {
2165: PetscFEGeom pgeom;
2166: const PetscInt dEt = T->cdim;
2167: const PetscInt dE = fegeom->dimEmbed;
2168: const PetscInt Nq = T->Np;
2169: const PetscInt Nb = T->Nb;
2170: const PetscInt Nc = T->Nc;
2171: const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2172: const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dEt];
2173: PetscInt q, b, c, d;
2174: PetscErrorCode ierr;
2176: for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2177: for (q = 0; q < Nq; ++q) {
2178: for (b = 0; b < Nb; ++b) {
2179: for (c = 0; c < Nc; ++c) {
2180: const PetscInt bcidx = b*Nc+c;
2182: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2183: for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dEt+bcidx*dEt+d];
2184: }
2185: }
2186: PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);
2187: PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);
2188: PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer);
2189: for (b = 0; b < Nb; ++b) {
2190: for (c = 0; c < Nc; ++c) {
2191: const PetscInt bcidx = b*Nc+c;
2192: const PetscInt qcidx = q*Nc+c;
2194: elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2195: for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2196: }
2197: }
2198: }
2199: return(0);
2200: }
2202: PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2203: {
2204: const PetscInt dE = T->cdim;
2205: const PetscInt Nq = T->Np;
2206: const PetscInt Nb = T->Nb;
2207: const PetscInt Nc = T->Nc;
2208: const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2209: const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2210: PetscInt q, b, c, d, s;
2211: PetscErrorCode ierr;
2213: for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0;
2214: for (q = 0; q < Nq; ++q) {
2215: for (b = 0; b < Nb; ++b) {
2216: for (c = 0; c < Nc; ++c) {
2217: const PetscInt bcidx = b*Nc+c;
2219: tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2220: for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2221: }
2222: }
2223: PetscFEPushforward(fe, fegeom, Nb, tmpBasis);
2224: PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);
2225: for (s = 0; s < 2; ++s) {
2226: for (b = 0; b < Nb; ++b) {
2227: for (c = 0; c < Nc; ++c) {
2228: const PetscInt bcidx = b*Nc+c;
2229: const PetscInt qcidx = (q*2+s)*Nc+c;
2231: elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2232: for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2233: }
2234: }
2235: }
2236: }
2237: return(0);
2238: }
2240: 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[])
2241: {
2242: const PetscInt dE = TI->cdim;
2243: const PetscInt NqI = TI->Np;
2244: const PetscInt NbI = TI->Nb;
2245: const PetscInt NcI = TI->Nc;
2246: const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2247: const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2248: const PetscInt NqJ = TJ->Np;
2249: const PetscInt NbJ = TJ->Nb;
2250: const PetscInt NcJ = TJ->Nc;
2251: const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2252: const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2253: PetscInt f, fc, g, gc, df, dg;
2254: PetscErrorCode ierr;
2256: for (f = 0; f < NbI; ++f) {
2257: for (fc = 0; fc < NcI; ++fc) {
2258: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2260: tmpBasisI[fidx] = basisI[fidx];
2261: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2262: }
2263: }
2264: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2265: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2266: for (g = 0; g < NbJ; ++g) {
2267: for (gc = 0; gc < NcJ; ++gc) {
2268: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2270: tmpBasisJ[gidx] = basisJ[gidx];
2271: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2272: }
2273: }
2274: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2275: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2276: for (f = 0; f < NbI; ++f) {
2277: for (fc = 0; fc < NcI; ++fc) {
2278: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2279: const PetscInt i = offsetI+f; /* Element matrix row */
2280: for (g = 0; g < NbJ; ++g) {
2281: for (gc = 0; gc < NcJ; ++gc) {
2282: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2283: const PetscInt j = offsetJ+g; /* Element matrix column */
2284: const PetscInt fOff = eOffset+i*totDim+j;
2286: elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2287: for (df = 0; df < dE; ++df) {
2288: elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2289: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2290: for (dg = 0; dg < dE; ++dg) {
2291: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2292: }
2293: }
2294: }
2295: }
2296: }
2297: }
2298: return(0);
2299: }
2301: PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2302: {
2303: const PetscInt dE = TI->cdim;
2304: const PetscInt NqI = TI->Np;
2305: const PetscInt NbI = TI->Nb;
2306: const PetscInt NcI = TI->Nc;
2307: const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2308: const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2309: const PetscInt NqJ = TJ->Np;
2310: const PetscInt NbJ = TJ->Nb;
2311: const PetscInt NcJ = TJ->Nc;
2312: const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2313: const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2314: const PetscInt Ns = isHybridI ? 1 : 2;
2315: const PetscInt Nt = isHybridJ ? 1 : 2;
2316: PetscInt f, fc, g, gc, df, dg, s, t;
2317: PetscErrorCode ierr;
2319: for (f = 0; f < NbI; ++f) {
2320: for (fc = 0; fc < NcI; ++fc) {
2321: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2323: tmpBasisI[fidx] = basisI[fidx];
2324: for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2325: }
2326: }
2327: PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);
2328: PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);
2329: for (g = 0; g < NbJ; ++g) {
2330: for (gc = 0; gc < NcJ; ++gc) {
2331: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2333: tmpBasisJ[gidx] = basisJ[gidx];
2334: for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2335: }
2336: }
2337: PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);
2338: PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);
2339: for (s = 0; s < Ns; ++s) {
2340: for (f = 0; f < NbI; ++f) {
2341: for (fc = 0; fc < NcI; ++fc) {
2342: const PetscInt sc = NcI*s+fc; /* components from each side of the surface */
2343: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2344: const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */
2345: for (t = 0; t < Nt; ++t) {
2346: for (g = 0; g < NbJ; ++g) {
2347: for (gc = 0; gc < NcJ; ++gc) {
2348: const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */
2349: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2350: const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */
2351: const PetscInt fOff = eOffset+i*totDim+j;
2353: elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2354: for (df = 0; df < dE; ++df) {
2355: elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2356: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2357: for (dg = 0; dg < dE; ++dg) {
2358: elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2359: }
2360: }
2361: }
2362: }
2363: }
2364: }
2365: }
2366: }
2367: return(0);
2368: }
2370: PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2371: {
2372: PetscDualSpace dsp;
2373: DM dm;
2374: PetscQuadrature quadDef;
2375: PetscInt dim, cdim, Nq;
2376: PetscErrorCode ierr;
2379: PetscFEGetDualSpace(fe, &dsp);
2380: PetscDualSpaceGetDM(dsp, &dm);
2381: DMGetDimension(dm, &dim);
2382: DMGetCoordinateDim(dm, &cdim);
2383: PetscFEGetQuadrature(fe, &quadDef);
2384: quad = quad ? quad : quadDef;
2385: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
2386: PetscMalloc1(Nq*cdim, &cgeom->v);
2387: PetscMalloc1(Nq*cdim*cdim, &cgeom->J);
2388: PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);
2389: PetscMalloc1(Nq, &cgeom->detJ);
2390: cgeom->dim = dim;
2391: cgeom->dimEmbed = cdim;
2392: cgeom->numCells = 1;
2393: cgeom->numPoints = Nq;
2394: DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);
2395: return(0);
2396: }
2398: PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2399: {
2403: PetscFree(cgeom->v);
2404: PetscFree(cgeom->J);
2405: PetscFree(cgeom->invJ);
2406: PetscFree(cgeom->detJ);
2407: return(0);
2408: }