Actual source code: dualspace.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: PetscClassId PETSCDUALSPACE_CLASSID = 0;
6: PetscLogEvent PETSCDUALSPACE_SetUp;
8: PetscFunctionList PetscDualSpaceList = NULL;
9: PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE;
11: /*
12: PetscDualSpaceLatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
13: Ordering is lexicographic with lowest index as least significant in ordering.
14: e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {0,2}.
16: Input Parameters:
17: + len - The length of the tuple
18: . max - The maximum sum
19: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
21: Output Parameter:
22: . tup - A tuple of len integers whos sum is at most 'max'
24: Level: developer
26: .seealso: `PetscDualSpaceTensorPointLexicographic_Internal()`
27: */
28: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
29: {
30: while (len--) {
31: max -= tup[len];
32: if (!max) {
33: tup[len] = 0;
34: break;
35: }
36: }
37: tup[++len]++;
38: return 0;
39: }
41: /*
42: PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
43: Ordering is lexicographic with lowest index as least significant in ordering.
44: e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.
46: Input Parameters:
47: + len - The length of the tuple
48: . max - The maximum value
49: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
51: Output Parameter:
52: . tup - A tuple of len integers whos sum is at most 'max'
54: Level: developer
56: .seealso: `PetscDualSpaceLatticePointLexicographic_Internal()`
57: */
58: PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
59: {
60: PetscInt i;
62: for (i = 0; i < len; i++) {
63: if (tup[i] < max) {
64: break;
65: } else {
66: tup[i] = 0;
67: }
68: }
69: tup[i]++;
70: return 0;
71: }
73: /*@C
74: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
76: Not Collective
78: Input Parameters:
79: + name - The name of a new user-defined creation routine
80: - create_func - The creation routine itself
82: Notes:
83: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
85: Sample usage:
86: .vb
87: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
88: .ve
90: Then, your PetscDualSpace type can be chosen with the procedural interface via
91: .vb
92: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
93: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
94: .ve
95: or at runtime via the option
96: .vb
97: -petscdualspace_type my_dual_space
98: .ve
100: Level: advanced
102: .seealso: `PetscDualSpaceRegisterAll()`, `PetscDualSpaceRegisterDestroy()`
104: @*/
105: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
106: {
107: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
108: return 0;
109: }
111: /*@C
112: PetscDualSpaceSetType - Builds a particular PetscDualSpace
114: Collective on sp
116: Input Parameters:
117: + sp - The PetscDualSpace object
118: - name - The kind of space
120: Options Database Key:
121: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
123: Level: intermediate
125: .seealso: `PetscDualSpaceGetType()`, `PetscDualSpaceCreate()`
126: @*/
127: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
128: {
129: PetscErrorCode (*r)(PetscDualSpace);
130: PetscBool match;
133: PetscObjectTypeCompare((PetscObject)sp, name, &match);
134: if (match) return 0;
136: if (!PetscDualSpaceRegisterAllCalled) PetscDualSpaceRegisterAll();
137: PetscFunctionListFind(PetscDualSpaceList, name, &r);
140: PetscTryTypeMethod(sp, destroy);
141: sp->ops->destroy = NULL;
143: (*r)(sp);
144: PetscObjectChangeTypeName((PetscObject)sp, name);
145: return 0;
146: }
148: /*@C
149: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
151: Not Collective
153: Input Parameter:
154: . sp - The PetscDualSpace
156: Output Parameter:
157: . name - The PetscDualSpace type name
159: Level: intermediate
161: .seealso: `PetscDualSpaceSetType()`, `PetscDualSpaceCreate()`
162: @*/
163: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
164: {
167: if (!PetscDualSpaceRegisterAllCalled) PetscDualSpaceRegisterAll();
168: *name = ((PetscObject)sp)->type_name;
169: return 0;
170: }
172: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
173: {
174: PetscViewerFormat format;
175: PetscInt pdim, f;
177: PetscDualSpaceGetDimension(sp, &pdim);
178: PetscObjectPrintClassNamePrefixType((PetscObject)sp, v);
179: PetscViewerASCIIPushTab(v);
180: if (sp->k) {
181: PetscViewerASCIIPrintf(v, "Dual space for %" PetscInt_FMT "-forms %swith %" PetscInt_FMT " components, size %" PetscInt_FMT "\n", PetscAbsInt(sp->k), sp->k < 0 ? "(stored in dual form) " : "", sp->Nc, pdim);
182: } else {
183: PetscViewerASCIIPrintf(v, "Dual space with %" PetscInt_FMT " components, size %" PetscInt_FMT "\n", sp->Nc, pdim);
184: }
185: PetscTryTypeMethod(sp, view, v);
186: PetscViewerGetFormat(v, &format);
187: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
188: PetscViewerASCIIPushTab(v);
189: for (f = 0; f < pdim; ++f) {
190: PetscViewerASCIIPrintf(v, "Dual basis vector %" PetscInt_FMT "\n", f);
191: PetscViewerASCIIPushTab(v);
192: PetscQuadratureView(sp->functional[f], v);
193: PetscViewerASCIIPopTab(v);
194: }
195: PetscViewerASCIIPopTab(v);
196: }
197: PetscViewerASCIIPopTab(v);
198: return 0;
199: }
201: /*@C
202: PetscDualSpaceViewFromOptions - View from Options
204: Collective on PetscDualSpace
206: Input Parameters:
207: + A - the PetscDualSpace object
208: . obj - Optional object, proivides prefix
209: - name - command line option
211: Level: intermediate
212: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscObjectViewFromOptions()`, `PetscDualSpaceCreate()`
213: @*/
214: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A, PetscObject obj, const char name[])
215: {
217: PetscObjectViewFromOptions((PetscObject)A, obj, name);
218: return 0;
219: }
221: /*@
222: PetscDualSpaceView - Views a PetscDualSpace
224: Collective on sp
226: Input Parameters:
227: + sp - the PetscDualSpace object to view
228: - v - the viewer
230: Level: beginner
232: .seealso `PetscDualSpaceDestroy()`, `PetscDualSpace`
233: @*/
234: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
235: {
236: PetscBool iascii;
240: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &v);
241: PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii);
242: if (iascii) PetscDualSpaceView_ASCII(sp, v);
243: return 0;
244: }
246: /*@
247: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
249: Collective on sp
251: Input Parameter:
252: . sp - the PetscDualSpace object to set options for
254: Options Database:
255: + -petscdualspace_order <order> - the approximation order of the space
256: . -petscdualspace_form_degree <deg> - the form degree, say 0 for point evaluations, or 2 for area integrals
257: . -petscdualspace_components <c> - the number of components, say d for a vector field
258: - -petscdualspace_refcell <celltype> - Reference cell type name
260: Level: intermediate
262: .seealso `PetscDualSpaceView()`, `PetscDualSpace`, `PetscObjectSetFromOptions()`
263: @*/
264: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
265: {
266: DMPolytopeType refCell = DM_POLYTOPE_TRIANGLE;
267: const char *defaultType;
268: char name[256];
269: PetscBool flg;
272: if (!((PetscObject)sp)->type_name) {
273: defaultType = PETSCDUALSPACELAGRANGE;
274: } else {
275: defaultType = ((PetscObject)sp)->type_name;
276: }
277: if (!PetscSpaceRegisterAllCalled) PetscSpaceRegisterAll();
279: PetscObjectOptionsBegin((PetscObject)sp);
280: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
281: if (flg) {
282: PetscDualSpaceSetType(sp, name);
283: } else if (!((PetscObject)sp)->type_name) {
284: PetscDualSpaceSetType(sp, defaultType);
285: }
286: PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL, 0);
287: PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);
288: PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL, 1);
289: PetscTryTypeMethod(sp, setfromoptions, PetscOptionsObject);
290: PetscOptionsEnum("-petscdualspace_refcell", "Reference cell shape", "PetscDualSpaceSetReferenceCell", DMPolytopeTypes, (PetscEnum)refCell, (PetscEnum *)&refCell, &flg);
291: if (flg) {
292: DM K;
294: DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K);
295: PetscDualSpaceSetDM(sp, K);
296: DMDestroy(&K);
297: }
299: /* process any options handlers added with PetscObjectAddOptionsHandler() */
300: PetscObjectProcessOptionsHandlers((PetscObject)sp, PetscOptionsObject);
301: PetscOptionsEnd();
302: sp->setfromoptionscalled = PETSC_TRUE;
303: return 0;
304: }
306: /*@
307: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
309: Collective on sp
311: Input Parameter:
312: . sp - the PetscDualSpace object to setup
314: Level: intermediate
316: .seealso `PetscDualSpaceView()`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
317: @*/
318: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
319: {
321: if (sp->setupcalled) return 0;
322: PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
323: sp->setupcalled = PETSC_TRUE;
324: PetscTryTypeMethod(sp, setup);
325: PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
326: if (sp->setfromoptionscalled) PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
327: return 0;
328: }
330: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
331: {
332: PetscInt pStart = -1, pEnd = -1, depth = -1;
334: if (!dm) return 0;
335: DMPlexGetChart(dm, &pStart, &pEnd);
336: DMPlexGetDepth(dm, &depth);
338: if (sp->pointSpaces) {
339: PetscInt i;
341: for (i = 0; i < pEnd - pStart; i++) PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
342: }
343: PetscFree(sp->pointSpaces);
345: if (sp->heightSpaces) {
346: PetscInt i;
348: for (i = 0; i <= depth; i++) PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
349: }
350: PetscFree(sp->heightSpaces);
352: PetscSectionDestroy(&(sp->pointSection));
353: PetscQuadratureDestroy(&(sp->intNodes));
354: VecDestroy(&(sp->intDofValues));
355: VecDestroy(&(sp->intNodeValues));
356: MatDestroy(&(sp->intMat));
357: PetscQuadratureDestroy(&(sp->allNodes));
358: VecDestroy(&(sp->allDofValues));
359: VecDestroy(&(sp->allNodeValues));
360: MatDestroy(&(sp->allMat));
361: PetscFree(sp->numDof);
362: return 0;
363: }
365: /*@
366: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
368: Collective on sp
370: Input Parameter:
371: . sp - the PetscDualSpace object to destroy
373: Level: beginner
375: .seealso `PetscDualSpaceView()`, `PetscDualSpace()`, `PetscDualSpaceCreate()`
376: @*/
377: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
378: {
379: PetscInt dim, f;
380: DM dm;
382: if (!*sp) return 0;
385: if (--((PetscObject)(*sp))->refct > 0) {
386: *sp = NULL;
387: return 0;
388: }
389: ((PetscObject)(*sp))->refct = 0;
391: PetscDualSpaceGetDimension(*sp, &dim);
392: dm = (*sp)->dm;
394: PetscTryTypeMethod((*sp), destroy);
395: PetscDualSpaceClearDMData_Internal(*sp, dm);
397: for (f = 0; f < dim; ++f) PetscQuadratureDestroy(&(*sp)->functional[f]);
398: PetscFree((*sp)->functional);
399: DMDestroy(&(*sp)->dm);
400: PetscHeaderDestroy(sp);
401: return 0;
402: }
404: /*@
405: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
407: Collective
409: Input Parameter:
410: . comm - The communicator for the PetscDualSpace object
412: Output Parameter:
413: . sp - The PetscDualSpace object
415: Level: beginner
417: .seealso: `PetscDualSpaceSetType()`, `PETSCDUALSPACELAGRANGE`
418: @*/
419: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
420: {
421: PetscDualSpace s;
424: PetscCitationsRegister(FECitation, &FEcite);
425: *sp = NULL;
426: PetscFEInitializePackage();
428: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
430: s->order = 0;
431: s->Nc = 1;
432: s->k = 0;
433: s->spdim = -1;
434: s->spintdim = -1;
435: s->uniform = PETSC_TRUE;
436: s->setupcalled = PETSC_FALSE;
438: *sp = s;
439: return 0;
440: }
442: /*@
443: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
445: Collective on sp
447: Input Parameter:
448: . sp - The original PetscDualSpace
450: Output Parameter:
451: . spNew - The duplicate PetscDualSpace
453: Level: beginner
455: .seealso: `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
456: @*/
457: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
458: {
459: DM dm;
460: PetscDualSpaceType type;
461: const char *name;
465: PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
466: PetscObjectGetName((PetscObject)sp, &name);
467: PetscObjectSetName((PetscObject)*spNew, name);
468: PetscDualSpaceGetType(sp, &type);
469: PetscDualSpaceSetType(*spNew, type);
470: PetscDualSpaceGetDM(sp, &dm);
471: PetscDualSpaceSetDM(*spNew, dm);
473: (*spNew)->order = sp->order;
474: (*spNew)->k = sp->k;
475: (*spNew)->Nc = sp->Nc;
476: (*spNew)->uniform = sp->uniform;
477: PetscTryTypeMethod(sp, duplicate, *spNew);
478: return 0;
479: }
481: /*@
482: PetscDualSpaceGetDM - Get the DM representing the reference cell
484: Not collective
486: Input Parameter:
487: . sp - The PetscDualSpace
489: Output Parameter:
490: . dm - The reference cell
492: Level: intermediate
494: .seealso: `PetscDualSpaceSetDM()`, `PetscDualSpaceCreate()`
495: @*/
496: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
497: {
500: *dm = sp->dm;
501: return 0;
502: }
504: /*@
505: PetscDualSpaceSetDM - Get the DM representing the reference cell
507: Not collective
509: Input Parameters:
510: + sp - The PetscDualSpace
511: - dm - The reference cell
513: Level: intermediate
515: .seealso: `PetscDualSpaceGetDM()`, `PetscDualSpaceCreate()`
516: @*/
517: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
518: {
522: PetscObjectReference((PetscObject)dm);
523: if (sp->dm && sp->dm != dm) PetscDualSpaceClearDMData_Internal(sp, sp->dm);
524: DMDestroy(&sp->dm);
525: sp->dm = dm;
526: return 0;
527: }
529: /*@
530: PetscDualSpaceGetOrder - Get the order of the dual space
532: Not collective
534: Input Parameter:
535: . sp - The PetscDualSpace
537: Output Parameter:
538: . order - The order
540: Level: intermediate
542: .seealso: `PetscDualSpaceSetOrder()`, `PetscDualSpaceCreate()`
543: @*/
544: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
545: {
548: *order = sp->order;
549: return 0;
550: }
552: /*@
553: PetscDualSpaceSetOrder - Set the order of the dual space
555: Not collective
557: Input Parameters:
558: + sp - The PetscDualSpace
559: - order - The order
561: Level: intermediate
563: .seealso: `PetscDualSpaceGetOrder()`, `PetscDualSpaceCreate()`
564: @*/
565: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
566: {
569: sp->order = order;
570: return 0;
571: }
573: /*@
574: PetscDualSpaceGetNumComponents - Return the number of components for this space
576: Input Parameter:
577: . sp - The PetscDualSpace
579: Output Parameter:
580: . Nc - The number of components
582: Note: A vector space, for example, will have d components, where d is the spatial dimension
584: Level: intermediate
586: .seealso: `PetscDualSpaceSetNumComponents()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
587: @*/
588: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
589: {
592: *Nc = sp->Nc;
593: return 0;
594: }
596: /*@
597: PetscDualSpaceSetNumComponents - Set the number of components for this space
599: Input Parameters:
600: + sp - The PetscDualSpace
601: - order - The number of components
603: Level: intermediate
605: .seealso: `PetscDualSpaceGetNumComponents()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
606: @*/
607: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
608: {
611: sp->Nc = Nc;
612: return 0;
613: }
615: /*@
616: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
618: Not collective
620: Input Parameters:
621: + sp - The PetscDualSpace
622: - i - The basis number
624: Output Parameter:
625: . functional - The basis functional
627: Level: intermediate
629: .seealso: `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`
630: @*/
631: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
632: {
633: PetscInt dim;
637: PetscDualSpaceGetDimension(sp, &dim);
639: *functional = sp->functional[i];
640: return 0;
641: }
643: /*@
644: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
646: Not collective
648: Input Parameter:
649: . sp - The PetscDualSpace
651: Output Parameter:
652: . dim - The dimension
654: Level: intermediate
656: .seealso: `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
657: @*/
658: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
659: {
662: if (sp->spdim < 0) {
663: PetscSection section;
665: PetscDualSpaceGetSection(sp, §ion);
666: if (section) {
667: PetscSectionGetStorageSize(section, &(sp->spdim));
668: } else sp->spdim = 0;
669: }
670: *dim = sp->spdim;
671: return 0;
672: }
674: /*@
675: PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain
677: Not collective
679: Input Parameter:
680: . sp - The PetscDualSpace
682: Output Parameter:
683: . dim - The dimension
685: Level: intermediate
687: .seealso: `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
688: @*/
689: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
690: {
693: if (sp->spintdim < 0) {
694: PetscSection section;
696: PetscDualSpaceGetSection(sp, §ion);
697: if (section) {
698: PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
699: } else sp->spintdim = 0;
700: }
701: *intdim = sp->spintdim;
702: return 0;
703: }
705: /*@
706: PetscDualSpaceGetUniform - Whether this dual space is uniform
708: Not collective
710: Input Parameters:
711: . sp - A dual space
713: Output Parameters:
714: . uniform - PETSC_TRUE if (a) the dual space is the same for each point in a stratum of the reference DMPlex, and
715: (b) every symmetry of each point in the reference DMPlex is also a symmetry of the point's dual space.
717: Level: advanced
719: Note: all of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells
720: with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform.
722: .seealso: `PetscDualSpaceGetPointSubspace()`, `PetscDualSpaceGetSymmetries()`
723: @*/
724: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
725: {
728: *uniform = sp->uniform;
729: return 0;
730: }
732: /*@C
733: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
735: Not collective
737: Input Parameter:
738: . sp - The PetscDualSpace
740: Output Parameter:
741: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
743: Level: intermediate
745: .seealso: `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
746: @*/
747: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
748: {
752: if (!sp->numDof) {
753: DM dm;
754: PetscInt depth, d;
755: PetscSection section;
757: PetscDualSpaceGetDM(sp, &dm);
758: DMPlexGetDepth(dm, &depth);
759: PetscCalloc1(depth + 1, &(sp->numDof));
760: PetscDualSpaceGetSection(sp, §ion);
761: for (d = 0; d <= depth; d++) {
762: PetscInt dStart, dEnd;
764: DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
765: if (dEnd <= dStart) continue;
766: PetscSectionGetDof(section, dStart, &(sp->numDof[d]));
767: }
768: }
769: *numDof = sp->numDof;
771: return 0;
772: }
774: /* create the section of the right size and set a permutation for topological ordering */
775: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
776: {
777: DM dm;
778: PetscInt pStart, pEnd, cStart, cEnd, c, depth, count, i;
779: PetscInt *seen, *perm;
780: PetscSection section;
782: dm = sp->dm;
783: PetscSectionCreate(PETSC_COMM_SELF, §ion);
784: DMPlexGetChart(dm, &pStart, &pEnd);
785: PetscSectionSetChart(section, pStart, pEnd);
786: PetscCalloc1(pEnd - pStart, &seen);
787: PetscMalloc1(pEnd - pStart, &perm);
788: DMPlexGetDepth(dm, &depth);
789: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
790: for (c = cStart, count = 0; c < cEnd; c++) {
791: PetscInt closureSize = -1, e;
792: PetscInt *closure = NULL;
794: perm[count++] = c;
795: seen[c - pStart] = 1;
796: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
797: for (e = 0; e < closureSize; e++) {
798: PetscInt point = closure[2 * e];
800: if (seen[point - pStart]) continue;
801: perm[count++] = point;
802: seen[point - pStart] = 1;
803: }
804: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
805: }
807: for (i = 0; i < pEnd - pStart; i++)
808: if (perm[i] != i) break;
809: if (i < pEnd - pStart) {
810: IS permIS;
812: ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
813: ISSetPermutation(permIS);
814: PetscSectionSetPermutation(section, permIS);
815: ISDestroy(&permIS);
816: } else {
817: PetscFree(perm);
818: }
819: PetscFree(seen);
820: *topSection = section;
821: return 0;
822: }
824: /* mark boundary points and set up */
825: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
826: {
827: DM dm;
828: DMLabel boundary;
829: PetscInt pStart, pEnd, p;
831: dm = sp->dm;
832: DMLabelCreate(PETSC_COMM_SELF, "boundary", &boundary);
833: PetscDualSpaceGetDM(sp, &dm);
834: DMPlexMarkBoundaryFaces(dm, 1, boundary);
835: DMPlexLabelComplete(dm, boundary);
836: DMPlexGetChart(dm, &pStart, &pEnd);
837: for (p = pStart; p < pEnd; p++) {
838: PetscInt bval;
840: DMLabelGetValue(boundary, p, &bval);
841: if (bval == 1) {
842: PetscInt dof;
844: PetscSectionGetDof(section, p, &dof);
845: PetscSectionSetConstraintDof(section, p, dof);
846: }
847: }
848: DMLabelDestroy(&boundary);
849: PetscSectionSetUp(section);
850: return 0;
851: }
853: /*@
854: PetscDualSpaceGetSection - Create a PetscSection over the reference cell with the layout from this space
856: Collective on sp
858: Input Parameters:
859: . sp - The PetscDualSpace
861: Output Parameter:
862: . section - The section
864: Level: advanced
866: .seealso: `PetscDualSpaceCreate()`, `DMPLEX`
867: @*/
868: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
869: {
870: PetscInt pStart, pEnd, p;
872: if (!sp->dm) {
873: *section = NULL;
874: return 0;
875: }
876: if (!sp->pointSection) {
877: /* mark the boundary */
878: PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
879: DMPlexGetChart(sp->dm, &pStart, &pEnd);
880: for (p = pStart; p < pEnd; p++) {
881: PetscDualSpace psp;
883: PetscDualSpaceGetPointSubspace(sp, p, &psp);
884: if (psp) {
885: PetscInt dof;
887: PetscDualSpaceGetInteriorDimension(psp, &dof);
888: PetscSectionSetDof(sp->pointSection, p, dof);
889: }
890: }
891: PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection);
892: }
893: *section = sp->pointSection;
894: return 0;
895: }
897: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
898: * have one cell */
899: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
900: {
901: PetscReal *sv0, *v0, *J;
902: PetscSection section;
903: PetscInt dim, s, k;
904: DM dm;
906: PetscDualSpaceGetDM(sp, &dm);
907: DMGetDimension(dm, &dim);
908: PetscDualSpaceGetSection(sp, §ion);
909: PetscMalloc3(dim, &v0, dim, &sv0, dim * dim, &J);
910: PetscDualSpaceGetFormDegree(sp, &k);
911: for (s = sStart; s < sEnd; s++) {
912: PetscReal detJ, hdetJ;
913: PetscDualSpace ssp;
914: PetscInt dof, off, f, sdim;
915: PetscInt i, j;
916: DM sdm;
918: PetscDualSpaceGetPointSubspace(sp, s, &ssp);
919: if (!ssp) continue;
920: PetscSectionGetDof(section, s, &dof);
921: PetscSectionGetOffset(section, s, &off);
922: /* get the first vertex of the reference cell */
923: PetscDualSpaceGetDM(ssp, &sdm);
924: DMGetDimension(sdm, &sdim);
925: DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);
926: DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);
927: /* compactify Jacobian */
928: for (i = 0; i < dim; i++)
929: for (j = 0; j < sdim; j++) J[i * sdim + j] = J[i * dim + j];
930: for (f = 0; f < dof; f++) {
931: PetscQuadrature fn;
933: PetscDualSpaceGetFunctional(ssp, f, &fn);
934: PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off + f]));
935: }
936: }
937: PetscFree3(v0, sv0, J);
938: return 0;
939: }
941: /*@C
942: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
944: Input Parameters:
945: + sp - The PetscDualSpace object
946: . f - The basis functional index
947: . time - The time
948: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian) (or evaluated at the coordinates of the functional)
949: . numComp - The number of components for the function
950: . func - The input function
951: - ctx - A context for the function
953: Output Parameter:
954: . value - numComp output values
956: Note: The calling sequence for the callback func is given by:
958: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
959: $ PetscInt numComponents, PetscScalar values[], void *ctx)
961: Level: beginner
963: .seealso: `PetscDualSpaceCreate()`
964: @*/
965: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
966: {
970: PetscUseTypeMethod(sp, apply, f, time, cgeom, numComp, func, ctx, value);
971: return 0;
972: }
974: /*@C
975: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
977: Input Parameters:
978: + sp - The PetscDualSpace object
979: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
981: Output Parameter:
982: . spValue - The values of all dual space functionals
984: Level: beginner
986: .seealso: `PetscDualSpaceCreate()`
987: @*/
988: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
989: {
991: PetscUseTypeMethod(sp, applyall, pointEval, spValue);
992: return 0;
993: }
995: /*@C
996: PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
998: Input Parameters:
999: + sp - The PetscDualSpace object
1000: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1002: Output Parameter:
1003: . spValue - The values of interior dual space functionals
1005: Level: beginner
1007: .seealso: `PetscDualSpaceCreate()`
1008: @*/
1009: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1010: {
1012: PetscUseTypeMethod(sp, applyint, pointEval, spValue);
1013: return 0;
1014: }
1016: /*@C
1017: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
1019: Input Parameters:
1020: + sp - The PetscDualSpace object
1021: . f - The basis functional index
1022: . time - The time
1023: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1024: . Nc - The number of components for the function
1025: . func - The input function
1026: - ctx - A context for the function
1028: Output Parameter:
1029: . value - The output value
1031: Note: The calling sequence for the callback func is given by:
1033: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1034: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1036: and the idea is to evaluate the functional as an integral
1038: $ n(f) = int dx n(x) . f(x)
1040: where both n and f have Nc components.
1042: Level: beginner
1044: .seealso: `PetscDualSpaceCreate()`
1045: @*/
1046: PetscErrorCode PetscDualSpaceApplyDefault(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1047: {
1048: DM dm;
1049: PetscQuadrature n;
1050: const PetscReal *points, *weights;
1051: PetscReal x[3];
1052: PetscScalar *val;
1053: PetscInt dim, dE, qNc, c, Nq, q;
1054: PetscBool isAffine;
1058: PetscDualSpaceGetDM(sp, &dm);
1059: PetscDualSpaceGetFunctional(sp, f, &n);
1060: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
1063: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1064: *value = 0.0;
1065: isAffine = cgeom->isAffine;
1066: dE = cgeom->dimEmbed;
1067: for (q = 0; q < Nq; ++q) {
1068: if (isAffine) {
1069: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q * dim], x);
1070: (*func)(dE, time, x, Nc, val, ctx);
1071: } else {
1072: (*func)(dE, time, &cgeom->v[dE * q], Nc, val, ctx);
1073: }
1074: for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1075: }
1076: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1077: return 0;
1078: }
1080: /*@C
1081: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
1083: Input Parameters:
1084: + sp - The PetscDualSpace object
1085: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
1087: Output Parameter:
1088: . spValue - The values of all dual space functionals
1090: Level: beginner
1092: .seealso: `PetscDualSpaceCreate()`
1093: @*/
1094: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1095: {
1096: Vec pointValues, dofValues;
1097: Mat allMat;
1102: PetscDualSpaceGetAllData(sp, NULL, &allMat);
1103: if (!(sp->allNodeValues)) MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1104: pointValues = sp->allNodeValues;
1105: if (!(sp->allDofValues)) MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1106: dofValues = sp->allDofValues;
1107: VecPlaceArray(pointValues, pointEval);
1108: VecPlaceArray(dofValues, spValue);
1109: MatMult(allMat, pointValues, dofValues);
1110: VecResetArray(dofValues);
1111: VecResetArray(pointValues);
1112: return 0;
1113: }
1115: /*@C
1116: PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
1118: Input Parameters:
1119: + sp - The PetscDualSpace object
1120: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1122: Output Parameter:
1123: . spValue - The values of interior dual space functionals
1125: Level: beginner
1127: .seealso: `PetscDualSpaceCreate()`
1128: @*/
1129: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1130: {
1131: Vec pointValues, dofValues;
1132: Mat intMat;
1137: PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1138: if (!(sp->intNodeValues)) MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1139: pointValues = sp->intNodeValues;
1140: if (!(sp->intDofValues)) MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1141: dofValues = sp->intDofValues;
1142: VecPlaceArray(pointValues, pointEval);
1143: VecPlaceArray(dofValues, spValue);
1144: MatMult(intMat, pointValues, dofValues);
1145: VecResetArray(dofValues);
1146: VecResetArray(pointValues);
1147: return 0;
1148: }
1150: /*@
1151: PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values
1153: Input Parameter:
1154: . sp - The dualspace
1156: Output Parameters:
1157: + allNodes - A PetscQuadrature object containing all evaluation nodes
1158: - allMat - A Mat for the node-to-dof transformation
1160: Level: advanced
1162: .seealso: `PetscDualSpaceCreate()`
1163: @*/
1164: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1165: {
1169: if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1170: PetscQuadrature qpoints;
1171: Mat amat;
1173: PetscUseTypeMethod(sp, createalldata, &qpoints, &amat);
1174: PetscQuadratureDestroy(&(sp->allNodes));
1175: MatDestroy(&(sp->allMat));
1176: sp->allNodes = qpoints;
1177: sp->allMat = amat;
1178: }
1179: if (allNodes) *allNodes = sp->allNodes;
1180: if (allMat) *allMat = sp->allMat;
1181: return 0;
1182: }
1184: /*@
1185: PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals
1187: Input Parameter:
1188: . sp - The dualspace
1190: Output Parameters:
1191: + allNodes - A PetscQuadrature object containing all evaluation nodes
1192: - allMat - A Mat for the node-to-dof transformation
1194: Level: advanced
1196: .seealso: `PetscDualSpaceCreate()`
1197: @*/
1198: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1199: {
1200: PetscInt spdim;
1201: PetscInt numPoints, offset;
1202: PetscReal *points;
1203: PetscInt f, dim;
1204: PetscInt Nc, nrows, ncols;
1205: PetscInt maxNumPoints;
1206: PetscQuadrature q;
1207: Mat A;
1209: PetscDualSpaceGetNumComponents(sp, &Nc);
1210: PetscDualSpaceGetDimension(sp, &spdim);
1211: if (!spdim) {
1212: PetscQuadratureCreate(PETSC_COMM_SELF, allNodes);
1213: PetscQuadratureSetData(*allNodes, 0, 0, 0, NULL, NULL);
1214: }
1215: nrows = spdim;
1216: PetscDualSpaceGetFunctional(sp, 0, &q);
1217: PetscQuadratureGetData(q, &dim, NULL, &numPoints, NULL, NULL);
1218: maxNumPoints = numPoints;
1219: for (f = 1; f < spdim; f++) {
1220: PetscInt Np;
1222: PetscDualSpaceGetFunctional(sp, f, &q);
1223: PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL);
1224: numPoints += Np;
1225: maxNumPoints = PetscMax(maxNumPoints, Np);
1226: }
1227: ncols = numPoints * Nc;
1228: PetscMalloc1(dim * numPoints, &points);
1229: MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1230: for (f = 0, offset = 0; f < spdim; f++) {
1231: const PetscReal *p, *w;
1232: PetscInt Np, i;
1233: PetscInt fnc;
1235: PetscDualSpaceGetFunctional(sp, f, &q);
1236: PetscQuadratureGetData(q, NULL, &fnc, &Np, &p, &w);
1238: for (i = 0; i < Np * dim; i++) points[offset * dim + i] = p[i];
1239: for (i = 0; i < Np * Nc; i++) MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1240: offset += Np;
1241: }
1242: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1243: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1244: PetscQuadratureCreate(PETSC_COMM_SELF, allNodes);
1245: PetscQuadratureSetData(*allNodes, dim, 0, numPoints, points, NULL);
1246: *allMat = A;
1247: return 0;
1248: }
1250: /*@
1251: PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1252: this space, as well as the matrix that computes the degrees of freedom from the quadrature values. Degrees of
1253: freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the
1254: reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by
1255: PetscDualSpaceGetSection()).
1257: Input Parameter:
1258: . sp - The dualspace
1260: Output Parameters:
1261: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1262: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1263: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1264: npoints is the number of points in intNodes and nc is PetscDualSpaceGetNumComponents().
1266: Level: advanced
1268: .seealso: `PetscDualSpaceCreate()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceGetNumComponents()`, `PetscQuadratureGetData()`
1269: @*/
1270: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1271: {
1275: if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1276: PetscQuadrature qpoints;
1277: Mat imat;
1279: PetscUseTypeMethod(sp, createintdata, &qpoints, &imat);
1280: PetscQuadratureDestroy(&(sp->intNodes));
1281: MatDestroy(&(sp->intMat));
1282: sp->intNodes = qpoints;
1283: sp->intMat = imat;
1284: }
1285: if (intNodes) *intNodes = sp->intNodes;
1286: if (intMat) *intMat = sp->intMat;
1287: return 0;
1288: }
1290: /*@
1291: PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values
1293: Input Parameter:
1294: . sp - The dualspace
1296: Output Parameters:
1297: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1298: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1299: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1300: npoints is the number of points in allNodes and nc is PetscDualSpaceGetNumComponents().
1302: Level: advanced
1304: .seealso: `PetscDualSpaceCreate()`, `PetscDualSpaceGetInteriorData()`
1305: @*/
1306: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1307: {
1308: DM dm;
1309: PetscInt spdim0;
1310: PetscInt Nc;
1311: PetscInt pStart, pEnd, p, f;
1312: PetscSection section;
1313: PetscInt numPoints, offset, matoffset;
1314: PetscReal *points;
1315: PetscInt dim;
1316: PetscInt *nnz;
1317: PetscQuadrature q;
1318: Mat imat;
1321: PetscDualSpaceGetSection(sp, §ion);
1322: PetscSectionGetConstrainedStorageSize(section, &spdim0);
1323: if (!spdim0) {
1324: *intNodes = NULL;
1325: *intMat = NULL;
1326: return 0;
1327: }
1328: PetscDualSpaceGetNumComponents(sp, &Nc);
1329: PetscSectionGetChart(section, &pStart, &pEnd);
1330: PetscDualSpaceGetDM(sp, &dm);
1331: DMGetDimension(dm, &dim);
1332: PetscMalloc1(spdim0, &nnz);
1333: for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1334: PetscInt dof, cdof, off, d;
1336: PetscSectionGetDof(section, p, &dof);
1337: PetscSectionGetConstraintDof(section, p, &cdof);
1338: if (!(dof - cdof)) continue;
1339: PetscSectionGetOffset(section, p, &off);
1340: for (d = 0; d < dof; d++, off++, f++) {
1341: PetscInt Np;
1343: PetscDualSpaceGetFunctional(sp, off, &q);
1344: PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL);
1345: nnz[f] = Np * Nc;
1346: numPoints += Np;
1347: }
1348: }
1349: MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1350: PetscFree(nnz);
1351: PetscMalloc1(dim * numPoints, &points);
1352: for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1353: PetscInt dof, cdof, off, d;
1355: PetscSectionGetDof(section, p, &dof);
1356: PetscSectionGetConstraintDof(section, p, &cdof);
1357: if (!(dof - cdof)) continue;
1358: PetscSectionGetOffset(section, p, &off);
1359: for (d = 0; d < dof; d++, off++, f++) {
1360: const PetscReal *p;
1361: const PetscReal *w;
1362: PetscInt Np, i;
1364: PetscDualSpaceGetFunctional(sp, off, &q);
1365: PetscQuadratureGetData(q, NULL, NULL, &Np, &p, &w);
1366: for (i = 0; i < Np * dim; i++) points[offset + i] = p[i];
1367: for (i = 0; i < Np * Nc; i++) MatSetValue(imat, f, matoffset + i, w[i], INSERT_VALUES);
1368: offset += Np * dim;
1369: matoffset += Np * Nc;
1370: }
1371: }
1372: PetscQuadratureCreate(PETSC_COMM_SELF, intNodes);
1373: PetscQuadratureSetData(*intNodes, dim, 0, numPoints, points, NULL);
1374: MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1375: MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1376: *intMat = imat;
1377: return 0;
1378: }
1380: /*@
1381: PetscDualSpaceEqual - Determine if a dual space is equivalent
1383: Input Parameters:
1384: + A - A PetscDualSpace object
1385: - B - Another PetscDualSpace object
1387: Output Parameter:
1388: . equal - PETSC_TRUE if the dual spaces are equivalent
1390: Level: advanced
1392: .seealso: `PetscDualSpaceCreate()`
1393: @*/
1394: PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal)
1395: {
1396: PetscInt sizeA, sizeB, dimA, dimB;
1397: const PetscInt *dofA, *dofB;
1398: PetscQuadrature quadA, quadB;
1399: Mat matA, matB;
1404: *equal = PETSC_FALSE;
1405: PetscDualSpaceGetDimension(A, &sizeA);
1406: PetscDualSpaceGetDimension(B, &sizeB);
1407: if (sizeB != sizeA) return 0;
1408: DMGetDimension(A->dm, &dimA);
1409: DMGetDimension(B->dm, &dimB);
1410: if (dimA != dimB) return 0;
1412: PetscDualSpaceGetNumDof(A, &dofA);
1413: PetscDualSpaceGetNumDof(B, &dofB);
1414: for (PetscInt d = 0; d < dimA; d++) {
1415: if (dofA[d] != dofB[d]) return 0;
1416: }
1418: PetscDualSpaceGetInteriorData(A, &quadA, &matA);
1419: PetscDualSpaceGetInteriorData(B, &quadB, &matB);
1420: if (!quadA && !quadB) {
1421: *equal = PETSC_TRUE;
1422: } else if (quadA && quadB) {
1423: PetscQuadratureEqual(quadA, quadB, equal);
1424: if (*equal == PETSC_FALSE) return 0;
1425: if (!matA && !matB) return 0;
1426: if (matA && matB) MatEqual(matA, matB, equal);
1427: else *equal = PETSC_FALSE;
1428: }
1429: return 0;
1430: }
1432: /*@C
1433: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
1435: Input Parameters:
1436: + sp - The PetscDualSpace object
1437: . f - The basis functional index
1438: . time - The time
1439: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1440: . Nc - The number of components for the function
1441: . func - The input function
1442: - ctx - A context for the function
1444: Output Parameter:
1445: . value - The output value (scalar)
1447: Note: The calling sequence for the callback func is given by:
1449: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1450: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1452: and the idea is to evaluate the functional as an integral
1454: $ n(f) = int dx n(x) . f(x)
1456: where both n and f have Nc components.
1458: Level: beginner
1460: .seealso: `PetscDualSpaceCreate()`
1461: @*/
1462: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1463: {
1464: DM dm;
1465: PetscQuadrature n;
1466: const PetscReal *points, *weights;
1467: PetscScalar *val;
1468: PetscInt dimEmbed, qNc, c, Nq, q;
1472: PetscDualSpaceGetDM(sp, &dm);
1473: DMGetCoordinateDim(dm, &dimEmbed);
1474: PetscDualSpaceGetFunctional(sp, f, &n);
1475: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1477: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1478: *value = 0.;
1479: for (q = 0; q < Nq; ++q) {
1480: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1481: for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1482: }
1483: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1484: return 0;
1485: }
1487: /*@
1488: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
1489: given height. This assumes that the reference cell is symmetric over points of this height.
1491: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1492: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1493: support extracting subspaces, then NULL is returned.
1495: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1497: Not collective
1499: Input Parameters:
1500: + sp - the PetscDualSpace object
1501: - height - the height of the mesh point for which the subspace is desired
1503: Output Parameter:
1504: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1505: point, which will be of lesser dimension if height > 0.
1507: Level: advanced
1509: .seealso: `PetscSpaceGetHeightSubspace()`, `PetscDualSpace`
1510: @*/
1511: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1512: {
1513: PetscInt depth = -1, cStart, cEnd;
1514: DM dm;
1519: *subsp = NULL;
1520: dm = sp->dm;
1521: DMPlexGetDepth(dm, &depth);
1523: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1524: if (height == 0 && cEnd == cStart + 1) {
1525: *subsp = sp;
1526: return 0;
1527: }
1528: if (!sp->heightSpaces) {
1529: PetscInt h;
1530: PetscCalloc1(depth + 1, &(sp->heightSpaces));
1532: for (h = 0; h <= depth; h++) {
1533: if (h == 0 && cEnd == cStart + 1) continue;
1534: if (sp->ops->createheightsubspace) (*sp->ops->createheightsubspace)(sp, height, &(sp->heightSpaces[h]));
1535: else if (sp->pointSpaces) {
1536: PetscInt hStart, hEnd;
1538: DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
1539: if (hEnd > hStart) {
1540: const char *name;
1542: PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1543: if (sp->pointSpaces[hStart]) {
1544: PetscObjectGetName((PetscObject)sp, &name);
1545: PetscObjectSetName((PetscObject)sp->pointSpaces[hStart], name);
1546: }
1547: sp->heightSpaces[h] = sp->pointSpaces[hStart];
1548: }
1549: }
1550: }
1551: }
1552: *subsp = sp->heightSpaces[height];
1553: return 0;
1554: }
1556: /*@
1557: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
1559: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
1560: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
1561: subspaces, then NULL is returned.
1563: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1565: Not collective
1567: Input Parameters:
1568: + sp - the PetscDualSpace object
1569: - point - the point (in the dual space's DM) for which the subspace is desired
1571: Output Parameters:
1572: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1573: point, which will be of lesser dimension if height > 0.
1575: Level: advanced
1577: .seealso: `PetscDualSpace`
1578: @*/
1579: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1580: {
1581: PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1582: DM dm;
1586: *bdsp = NULL;
1587: dm = sp->dm;
1588: DMPlexGetChart(dm, &pStart, &pEnd);
1590: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1591: if (point == cStart && cEnd == cStart + 1) { /* the dual space is only equivalent to the dual space on a cell if the reference mesh has just one cell */
1592: *bdsp = sp;
1593: return 0;
1594: }
1595: if (!sp->pointSpaces) {
1596: PetscInt p;
1597: PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));
1599: for (p = 0; p < pEnd - pStart; p++) {
1600: if (p + pStart == cStart && cEnd == cStart + 1) continue;
1601: if (sp->ops->createpointsubspace) (*sp->ops->createpointsubspace)(sp, p + pStart, &(sp->pointSpaces[p]));
1602: else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1603: PetscInt dim, depth, height;
1604: DMLabel label;
1606: DMPlexGetDepth(dm, &dim);
1607: DMPlexGetDepthLabel(dm, &label);
1608: DMLabelGetValue(label, p + pStart, &depth);
1609: height = dim - depth;
1610: PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1611: PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1612: }
1613: }
1614: }
1615: *bdsp = sp->pointSpaces[point - pStart];
1616: return 0;
1617: }
1619: /*@C
1620: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
1622: Not collective
1624: Input Parameter:
1625: . sp - the PetscDualSpace object
1627: Output Parameters:
1628: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1629: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation
1631: Note: The permutation and flip arrays are organized in the following way
1632: $ perms[p][ornt][dof # on point] = new local dof #
1633: $ flips[p][ornt][dof # on point] = reversal or not
1635: Level: developer
1637: @*/
1638: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1639: {
1641: if (perms) {
1643: *perms = NULL;
1644: }
1645: if (flips) {
1647: *flips = NULL;
1648: }
1649: if (sp->ops->getsymmetries) (sp->ops->getsymmetries)(sp, perms, flips);
1650: return 0;
1651: }
1653: /*@
1654: PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1655: dual space's functionals.
1657: Input Parameter:
1658: . dsp - The PetscDualSpace
1660: Output Parameter:
1661: . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1662: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1663: the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz).
1664: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1665: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1666: but are stored as 1-forms.
1668: Level: developer
1670: .seealso: `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1671: @*/
1672: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1673: {
1677: *k = dsp->k;
1678: return 0;
1679: }
1681: /*@
1682: PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1683: dual space's functionals.
1685: Input Parameters:
1686: + dsp - The PetscDualSpace
1687: - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1688: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1689: the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz).
1690: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1691: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1692: but are stored as 1-forms.
1694: Level: developer
1696: .seealso: `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1697: @*/
1698: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1699: {
1700: PetscInt dim;
1705: dim = dsp->dm->dim;
1707: dsp->k = k;
1708: return 0;
1709: }
1711: /*@
1712: PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space
1714: Input Parameter:
1715: . dsp - The PetscDualSpace
1717: Output Parameter:
1718: . k - The simplex dimension
1720: Level: developer
1722: Note: Currently supported values are
1723: $ 0: These are H_1 methods that only transform coordinates
1724: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1725: $ 2: These are the same as 1
1726: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1728: .seealso: `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1729: @*/
1730: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1731: {
1732: PetscInt dim;
1737: dim = dsp->dm->dim;
1738: if (!dsp->k) *k = IDENTITY_TRANSFORM;
1739: else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1740: else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1741: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1742: return 0;
1743: }
1745: /*@C
1746: PetscDualSpaceTransform - Transform the function values
1748: Input Parameters:
1749: + dsp - The PetscDualSpace
1750: . trans - The type of transform
1751: . isInverse - Flag to invert the transform
1752: . fegeom - The cell geometry
1753: . Nv - The number of function samples
1754: . Nc - The number of function components
1755: - vals - The function values
1757: Output Parameter:
1758: . vals - The transformed function values
1760: Level: intermediate
1762: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1764: .seealso: `PetscDualSpaceTransformGradient()`, `PetscDualSpaceTransformHessian()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1765: @*/
1766: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1767: {
1768: PetscReal Jstar[9] = {0};
1769: PetscInt dim, v, c, Nk;
1775: /* TODO: not handling dimEmbed != dim right now */
1776: dim = dsp->dm->dim;
1777: /* No change needed for 0-forms */
1778: if (!dsp->k) return 0;
1779: PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1780: /* TODO: use fegeom->isAffine */
1781: PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1782: for (v = 0; v < Nv; ++v) {
1783: switch (Nk) {
1784: case 1:
1785: for (c = 0; c < Nc; c++) vals[v * Nc + c] *= Jstar[0];
1786: break;
1787: case 2:
1788: for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1789: break;
1790: case 3:
1791: for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1792: break;
1793: default:
1794: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %" PetscInt_FMT " for transformation", Nk);
1795: }
1796: }
1797: return 0;
1798: }
1800: /*@C
1801: PetscDualSpaceTransformGradient - Transform the function gradient values
1803: Input Parameters:
1804: + dsp - The PetscDualSpace
1805: . trans - The type of transform
1806: . isInverse - Flag to invert the transform
1807: . fegeom - The cell geometry
1808: . Nv - The number of function gradient samples
1809: . Nc - The number of function components
1810: - vals - The function gradient values
1812: Output Parameter:
1813: . vals - The transformed function gradient values
1815: Level: intermediate
1817: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1819: .seealso: `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1820: @*/
1821: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1822: {
1823: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1824: PetscInt v, c, d;
1830: #ifdef PETSC_USE_DEBUG
1832: #endif
1833: /* Transform gradient */
1834: if (dim == dE) {
1835: for (v = 0; v < Nv; ++v) {
1836: for (c = 0; c < Nc; ++c) {
1837: switch (dim) {
1838: case 1:
1839: vals[(v * Nc + c) * dim] *= fegeom->invJ[0];
1840: break;
1841: case 2:
1842: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1843: break;
1844: case 3:
1845: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1846: break;
1847: default:
1848: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1849: }
1850: }
1851: }
1852: } else {
1853: for (v = 0; v < Nv; ++v) {
1854: for (c = 0; c < Nc; ++c) DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v * Nc + c) * dE], &vals[(v * Nc + c) * dE]);
1855: }
1856: }
1857: /* Assume its a vector, otherwise assume its a bunch of scalars */
1858: if (Nc == 1 || Nc != dim) return 0;
1859: switch (trans) {
1860: case IDENTITY_TRANSFORM:
1861: break;
1862: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1863: if (isInverse) {
1864: for (v = 0; v < Nv; ++v) {
1865: for (d = 0; d < dim; ++d) {
1866: switch (dim) {
1867: case 2:
1868: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1869: break;
1870: case 3:
1871: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1872: break;
1873: default:
1874: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1875: }
1876: }
1877: }
1878: } else {
1879: for (v = 0; v < Nv; ++v) {
1880: for (d = 0; d < dim; ++d) {
1881: switch (dim) {
1882: case 2:
1883: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1884: break;
1885: case 3:
1886: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1887: break;
1888: default:
1889: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1890: }
1891: }
1892: }
1893: }
1894: break;
1895: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1896: if (isInverse) {
1897: for (v = 0; v < Nv; ++v) {
1898: for (d = 0; d < dim; ++d) {
1899: switch (dim) {
1900: case 2:
1901: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1902: break;
1903: case 3:
1904: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1905: break;
1906: default:
1907: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1908: }
1909: for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] *= fegeom->detJ[0];
1910: }
1911: }
1912: } else {
1913: for (v = 0; v < Nv; ++v) {
1914: for (d = 0; d < dim; ++d) {
1915: switch (dim) {
1916: case 2:
1917: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1918: break;
1919: case 3:
1920: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1921: break;
1922: default:
1923: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1924: }
1925: for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] /= fegeom->detJ[0];
1926: }
1927: }
1928: }
1929: break;
1930: }
1931: return 0;
1932: }
1934: /*@C
1935: PetscDualSpaceTransformHessian - Transform the function Hessian values
1937: Input Parameters:
1938: + dsp - The PetscDualSpace
1939: . trans - The type of transform
1940: . isInverse - Flag to invert the transform
1941: . fegeom - The cell geometry
1942: . Nv - The number of function Hessian samples
1943: . Nc - The number of function components
1944: - vals - The function gradient values
1946: Output Parameter:
1947: . vals - The transformed function Hessian values
1949: Level: intermediate
1951: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1953: .seealso: `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1954: @*/
1955: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1956: {
1957: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1958: PetscInt v, c;
1964: #ifdef PETSC_USE_DEBUG
1966: #endif
1967: /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
1968: if (dim == dE) {
1969: for (v = 0; v < Nv; ++v) {
1970: for (c = 0; c < Nc; ++c) {
1971: switch (dim) {
1972: case 1:
1973: vals[(v * Nc + c) * dim * dim] *= PetscSqr(fegeom->invJ[0]);
1974: break;
1975: case 2:
1976: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
1977: break;
1978: case 3:
1979: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
1980: break;
1981: default:
1982: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1983: }
1984: }
1985: }
1986: } else {
1987: for (v = 0; v < Nv; ++v) {
1988: for (c = 0; c < Nc; ++c) DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v * Nc + c) * dE * dE], &vals[(v * Nc + c) * dE * dE]);
1989: }
1990: }
1991: /* Assume its a vector, otherwise assume its a bunch of scalars */
1992: if (Nc == 1 || Nc != dim) return 0;
1993: switch (trans) {
1994: case IDENTITY_TRANSFORM:
1995: break;
1996: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1997: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
1998: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1999: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2000: }
2001: return 0;
2002: }
2004: /*@C
2005: PetscDualSpacePullback - Transform the given functional so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.
2007: Input Parameters:
2008: + dsp - The PetscDualSpace
2009: . fegeom - The geometry for this cell
2010: . Nq - The number of function samples
2011: . Nc - The number of function components
2012: - pointEval - The function values
2014: Output Parameter:
2015: . pointEval - The transformed function values
2017: Level: advanced
2019: Note: Functions transform in a complementary way (pushforward) to functionals, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.
2021: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2023: .seealso: `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2024: @*/
2025: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2026: {
2027: PetscDualSpaceTransformType trans;
2028: PetscInt k;
2034: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2035: This determines their transformation properties. */
2036: PetscDualSpaceGetDeRahm(dsp, &k);
2037: switch (k) {
2038: case 0: /* H^1 point evaluations */
2039: trans = IDENTITY_TRANSFORM;
2040: break;
2041: case 1: /* Hcurl preserves tangential edge traces */
2042: trans = COVARIANT_PIOLA_TRANSFORM;
2043: break;
2044: case 2:
2045: case 3: /* Hdiv preserve normal traces */
2046: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2047: break;
2048: default:
2049: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2050: }
2051: PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2052: return 0;
2053: }
2055: /*@C
2056: PetscDualSpacePushforward - Transform the given function so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.
2058: Input Parameters:
2059: + dsp - The PetscDualSpace
2060: . fegeom - The geometry for this cell
2061: . Nq - The number of function samples
2062: . Nc - The number of function components
2063: - pointEval - The function values
2065: Output Parameter:
2066: . pointEval - The transformed function values
2068: Level: advanced
2070: Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.
2072: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2074: .seealso: `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2075: @*/
2076: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2077: {
2078: PetscDualSpaceTransformType trans;
2079: PetscInt k;
2085: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2086: This determines their transformation properties. */
2087: PetscDualSpaceGetDeRahm(dsp, &k);
2088: switch (k) {
2089: case 0: /* H^1 point evaluations */
2090: trans = IDENTITY_TRANSFORM;
2091: break;
2092: case 1: /* Hcurl preserves tangential edge traces */
2093: trans = COVARIANT_PIOLA_TRANSFORM;
2094: break;
2095: case 2:
2096: case 3: /* Hdiv preserve normal traces */
2097: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2098: break;
2099: default:
2100: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2101: }
2102: PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2103: return 0;
2104: }
2106: /*@C
2107: PetscDualSpacePushforwardGradient - Transform the given function gradient so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.
2109: Input Parameters:
2110: + dsp - The PetscDualSpace
2111: . fegeom - The geometry for this cell
2112: . Nq - The number of function gradient samples
2113: . Nc - The number of function components
2114: - pointEval - The function gradient values
2116: Output Parameter:
2117: . pointEval - The transformed function gradient values
2119: Level: advanced
2121: Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.
2123: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2125: .seealso: `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2126: @*/
2127: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2128: {
2129: PetscDualSpaceTransformType trans;
2130: PetscInt k;
2136: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2137: This determines their transformation properties. */
2138: PetscDualSpaceGetDeRahm(dsp, &k);
2139: switch (k) {
2140: case 0: /* H^1 point evaluations */
2141: trans = IDENTITY_TRANSFORM;
2142: break;
2143: case 1: /* Hcurl preserves tangential edge traces */
2144: trans = COVARIANT_PIOLA_TRANSFORM;
2145: break;
2146: case 2:
2147: case 3: /* Hdiv preserve normal traces */
2148: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2149: break;
2150: default:
2151: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2152: }
2153: PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2154: return 0;
2155: }
2157: /*@C
2158: PetscDualSpacePushforwardHessian - Transform the given function Hessian so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.
2160: Input Parameters:
2161: + dsp - The PetscDualSpace
2162: . fegeom - The geometry for this cell
2163: . Nq - The number of function Hessian samples
2164: . Nc - The number of function components
2165: - pointEval - The function gradient values
2167: Output Parameter:
2168: . pointEval - The transformed function Hessian values
2170: Level: advanced
2172: Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.
2174: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2176: .seealso: `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2177: @*/
2178: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2179: {
2180: PetscDualSpaceTransformType trans;
2181: PetscInt k;
2187: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2188: This determines their transformation properties. */
2189: PetscDualSpaceGetDeRahm(dsp, &k);
2190: switch (k) {
2191: case 0: /* H^1 point evaluations */
2192: trans = IDENTITY_TRANSFORM;
2193: break;
2194: case 1: /* Hcurl preserves tangential edge traces */
2195: trans = COVARIANT_PIOLA_TRANSFORM;
2196: break;
2197: case 2:
2198: case 3: /* Hdiv preserve normal traces */
2199: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2200: break;
2201: default:
2202: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2203: }
2204: PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2205: return 0;
2206: }