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 whose sum is at most `max`
24: Level: developer
26: .seealso: `PetscDualSpaceType`, `PetscDualSpaceTensorPointLexicographic_Internal()`
27: */
28: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
29: {
30: PetscFunctionBegin;
31: while (len--) {
32: max -= tup[len];
33: if (!max) {
34: tup[len] = 0;
35: break;
36: }
37: }
38: tup[++len]++;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: /*
43: PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
44: Ordering is lexicographic with lowest index as least significant in ordering.
45: 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}.
47: Input Parameters:
48: + len - The length of the tuple
49: . max - The maximum value
50: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
52: Output Parameter:
53: . tup - A tuple of `len` integers whose entries are at most `max`
55: Level: developer
57: .seealso: `PetscDualSpaceType`, `PetscDualSpaceLatticePointLexicographic_Internal()`
58: */
59: PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
60: {
61: PetscInt i;
63: PetscFunctionBegin;
64: for (i = 0; i < len; i++) {
65: if (tup[i] < max) {
66: break;
67: } else {
68: tup[i] = 0;
69: }
70: }
71: tup[i]++;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: /*@C
76: PetscDualSpaceRegister - Adds a new `PetscDualSpaceType`
78: Not Collective, No Fortran Support
80: Input Parameters:
81: + sname - The name of a new user-defined creation routine
82: - function - The creation routine
84: Example Usage:
85: .vb
86: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
87: .ve
89: Then, your PetscDualSpace type can be chosen with the procedural interface via
90: .vb
91: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
92: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
93: .ve
94: or at runtime via the option
95: .vb
96: -petscdualspace_type my_dual_space
97: .ve
99: Level: advanced
101: Note:
102: `PetscDualSpaceRegister()` may be called multiple times to add several user-defined `PetscDualSpace`
104: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceRegisterAll()`, `PetscDualSpaceRegisterDestroy()`
105: @*/
106: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
107: {
108: PetscFunctionBegin;
109: PetscCall(PetscFunctionListAdd(&PetscDualSpaceList, sname, function));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /*@C
114: PetscDualSpaceSetType - Builds a particular `PetscDualSpace` based on its `PetscDualSpaceType`
116: Collective
118: Input Parameters:
119: + sp - The `PetscDualSpace` object
120: - name - The kind of space
122: Options Database Key:
123: . -petscdualspace_type type - Sets the `PetscDualSpace` type; see `PetscDualSpaceType` for the choices
125: Level: intermediate
127: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceGetType()`, `PetscDualSpaceCreate()`
128: @*/
129: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
130: {
131: PetscErrorCode (*r)(PetscDualSpace);
132: PetscBool match;
134: PetscFunctionBegin;
136: PetscCall(PetscObjectTypeCompare((PetscObject)sp, name, &match));
137: if (match) PetscFunctionReturn(PETSC_SUCCESS);
139: if (!PetscDualSpaceRegisterAllCalled) PetscCall(PetscDualSpaceRegisterAll());
140: PetscCall(PetscFunctionListFind(PetscDualSpaceList, name, &r));
141: PetscCheck(r, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
143: PetscTryTypeMethod(sp, destroy);
144: sp->ops->destroy = NULL;
146: PetscCall((*r)(sp));
147: PetscCall(PetscObjectChangeTypeName((PetscObject)sp, name));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@C
152: PetscDualSpaceGetType - Gets the `PetscDualSpaceType` name (as a string) from the object.
154: Not Collective
156: Input Parameter:
157: . sp - The `PetscDualSpace`
159: Output Parameter:
160: . name - The `PetscDualSpaceType` name
162: Level: intermediate
164: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceSetType()`, `PetscDualSpaceCreate()`
165: @*/
166: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
167: {
168: PetscFunctionBegin;
170: PetscAssertPointer(name, 2);
171: if (!PetscDualSpaceRegisterAllCalled) PetscCall(PetscDualSpaceRegisterAll());
172: *name = ((PetscObject)sp)->type_name;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
177: {
178: PetscViewerFormat format;
179: PetscInt pdim;
181: PetscFunctionBegin;
182: PetscCall(PetscDualSpaceGetDimension(sp, &pdim));
183: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sp, v));
184: PetscCall(PetscViewerASCIIPushTab(v));
185: if (sp->k != 0 && sp->k != PETSC_FORM_DEGREE_UNDEFINED) {
186: PetscCall(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));
187: } else {
188: PetscCall(PetscViewerASCIIPrintf(v, "Dual space with %" PetscInt_FMT " components, size %" PetscInt_FMT "\n", sp->Nc, pdim));
189: }
190: PetscTryTypeMethod(sp, view, v);
191: PetscCall(PetscViewerGetFormat(v, &format));
192: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
193: PetscCall(PetscViewerASCIIPushTab(v));
194: for (PetscInt f = 0; f < pdim; ++f) {
195: PetscCall(PetscViewerASCIIPrintf(v, "Dual basis vector %" PetscInt_FMT "\n", f));
196: PetscCall(PetscViewerASCIIPushTab(v));
197: PetscCall(PetscQuadratureView(sp->functional[f], v));
198: PetscCall(PetscViewerASCIIPopTab(v));
199: }
200: PetscCall(PetscViewerASCIIPopTab(v));
201: }
202: PetscCall(PetscViewerASCIIPopTab(v));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*@C
207: PetscDualSpaceViewFromOptions - View a `PetscDualSpace` based on values in the options database
209: Collective
211: Input Parameters:
212: + A - the `PetscDualSpace` object
213: . obj - Optional object, provides the options prefix
214: - name - command line option name
216: Options Database Key:
217: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
219: Level: intermediate
221: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscObjectViewFromOptions()`, `PetscDualSpaceCreate()`
222: @*/
223: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A, PeOp PetscObject obj, const char name[])
224: {
225: PetscFunctionBegin;
227: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@C
232: PetscDualSpaceView - Views a `PetscDualSpace`
234: Collective
236: Input Parameters:
237: + sp - the `PetscDualSpace` object to view
238: - v - the viewer
240: Level: beginner
242: .seealso: `PetscViewer`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
243: @*/
244: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
245: {
246: PetscBool isascii;
248: PetscFunctionBegin;
251: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &v));
252: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
253: if (isascii) PetscCall(PetscDualSpaceView_ASCII(sp, v));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@C
258: PetscDualSpaceSetFromOptions - sets parameters in a `PetscDualSpace` from the options database
260: Collective
262: Input Parameter:
263: . sp - the `PetscDualSpace` object to set options for
265: Options Database Keys:
266: + -petscdualspace_order order - the approximation order of the space
267: . -petscdualspace_form_degree deg - the form degree, say 0 for point evaluations, or 2 for area integrals
268: . -petscdualspace_components c - the number of components, say d for a vector field
269: . -petscdualspace_refcell celltype - Reference cell type name
270: . -petscdualspace_lagrange_continuity (true|false) - Flag for continuous element
271: . -petscdualspace_lagrange_tensor (true|false) - Flag for tensor dual space
272: . -petscdualspace_lagrange_trimmed (true|false) - Flag for trimmed dual space
273: . -petscdualspace_lagrange_node_type nodetype - Lagrange node location type
274: . -petscdualspace_lagrange_node_endpoints (true|false) - Flag for nodes that include endpoints
275: . -petscdualspace_lagrange_node_exponent exponent - Gauss-Jacobi weight function exponent
276: . -petscdualspace_lagrange_use_moments (true|false) - Use moments (where appropriate) for functionals
277: - -petscdualspace_lagrange_moment_order order - Quadrature order for moment functionals
279: Level: intermediate
281: .seealso: `PetscDualSpaceView()`, `PetscDualSpace`, `PetscObjectSetFromOptions()`
282: @*/
283: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
284: {
285: DMPolytopeType refCell = DM_POLYTOPE_TRIANGLE;
286: const char *defaultType;
287: char name[256];
288: PetscBool flg;
290: PetscFunctionBegin;
292: if (!((PetscObject)sp)->type_name) defaultType = PETSCDUALSPACELAGRANGE;
293: else defaultType = ((PetscObject)sp)->type_name;
294: if (!PetscSpaceRegisterAllCalled) PetscCall(PetscSpaceRegisterAll());
296: PetscObjectOptionsBegin((PetscObject)sp);
297: PetscCall(PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg));
298: if (flg) PetscCall(PetscDualSpaceSetType(sp, name));
299: else if (!((PetscObject)sp)->type_name) PetscCall(PetscDualSpaceSetType(sp, defaultType));
300: PetscCall(PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL, 0));
301: PetscCall(PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL));
302: PetscCall(PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL, 1));
303: PetscTryTypeMethod(sp, setfromoptions, PetscOptionsObject);
304: PetscCall(PetscOptionsEnum("-petscdualspace_refcell", "Reference cell shape", "PetscDualSpaceSetReferenceCell", DMPolytopeTypes, (PetscEnum)refCell, (PetscEnum *)&refCell, &flg));
305: if (flg) {
306: DM K;
308: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K));
309: PetscCall(PetscDualSpaceSetDM(sp, K));
310: PetscCall(DMDestroy(&K));
311: }
313: /* process any options handlers added with PetscObjectAddOptionsHandler() */
314: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)sp, PetscOptionsObject));
315: PetscOptionsEnd();
316: sp->setfromoptionscalled = PETSC_TRUE;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@C
321: PetscDualSpaceSetUp - Construct a basis for a `PetscDualSpace`
323: Collective
325: Input Parameter:
326: . sp - the `PetscDualSpace` object to setup
328: Level: intermediate
330: .seealso: `PetscDualSpaceView()`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
331: @*/
332: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
333: {
334: PetscFunctionBegin;
336: if (sp->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
337: PetscCall(PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0));
338: sp->setupcalled = PETSC_TRUE;
339: PetscTryTypeMethod(sp, setup);
340: PetscCall(PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0));
341: if (sp->setfromoptionscalled) PetscCall(PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view"));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
346: {
347: PetscInt pStart = -1, pEnd = -1, depth = -1;
349: PetscFunctionBegin;
350: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
351: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
352: PetscCall(DMPlexGetDepth(dm, &depth));
354: if (sp->pointSpaces) {
355: for (PetscInt i = 0; i < pEnd - pStart; i++) PetscCall(PetscDualSpaceDestroy(&sp->pointSpaces[i]));
356: }
357: PetscCall(PetscFree(sp->pointSpaces));
359: if (sp->heightSpaces) {
360: for (PetscInt i = 0; i <= depth; i++) PetscCall(PetscDualSpaceDestroy(&sp->heightSpaces[i]));
361: }
362: PetscCall(PetscFree(sp->heightSpaces));
364: PetscCall(PetscSectionDestroy(&sp->pointSection));
365: PetscCall(PetscSectionDestroy(&sp->intPointSection));
366: PetscCall(PetscQuadratureDestroy(&sp->intNodes));
367: PetscCall(VecDestroy(&sp->intDofValues));
368: PetscCall(VecDestroy(&sp->intNodeValues));
369: PetscCall(MatDestroy(&sp->intMat));
370: PetscCall(PetscQuadratureDestroy(&sp->allNodes));
371: PetscCall(VecDestroy(&sp->allDofValues));
372: PetscCall(VecDestroy(&sp->allNodeValues));
373: PetscCall(MatDestroy(&sp->allMat));
374: PetscCall(PetscFree(sp->numDof));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@C
379: PetscDualSpaceDestroy - Destroys a `PetscDualSpace` object
381: Collective
383: Input Parameter:
384: . sp - the `PetscDualSpace` object to destroy
386: Level: beginner
388: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscDualSpace()`, `PetscDualSpaceCreate()`
389: @*/
390: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
391: {
392: PetscInt dim;
393: DM dm;
395: PetscFunctionBegin;
396: if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);
399: if (--((PetscObject)*sp)->refct > 0) {
400: *sp = NULL;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
403: ((PetscObject)*sp)->refct = 0;
405: PetscCall(PetscDualSpaceGetDimension(*sp, &dim));
406: dm = (*sp)->dm;
408: PetscTryTypeMethod(*sp, destroy);
409: PetscCall(PetscDualSpaceClearDMData_Internal(*sp, dm));
411: for (PetscInt f = 0; f < dim; ++f) PetscCall(PetscQuadratureDestroy(&(*sp)->functional[f]));
412: PetscCall(PetscFree((*sp)->functional));
413: PetscCall(DMDestroy(&(*sp)->dm));
414: PetscCall(PetscHeaderDestroy(sp));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /*@C
419: PetscDualSpaceCreate - Creates an empty `PetscDualSpace` object. The type can then be set with `PetscDualSpaceSetType()`.
421: Collective
423: Input Parameter:
424: . comm - The communicator for the `PetscDualSpace` object
426: Output Parameter:
427: . sp - The `PetscDualSpace` object
429: Level: beginner
431: .seealso: `PetscDualSpace`, `PetscDualSpaceSetType()`, `PETSCDUALSPACELAGRANGE`
432: @*/
433: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
434: {
435: PetscDualSpace s;
437: PetscFunctionBegin;
438: PetscAssertPointer(sp, 2);
439: PetscCall(PetscCitationsRegister(FECitation, &FEcite));
440: PetscCall(PetscFEInitializePackage());
442: PetscCall(PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView));
443: s->order = 0;
444: s->Nc = 1;
445: s->k = 0;
446: s->spdim = -1;
447: s->spintdim = -1;
448: s->uniform = PETSC_TRUE;
449: s->setupcalled = PETSC_FALSE;
450: *sp = s;
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*@C
455: PetscDualSpaceDuplicate - Creates a duplicate `PetscDualSpace` object that is not setup.
457: Collective
459: Input Parameter:
460: . sp - The original `PetscDualSpace`
462: Output Parameter:
463: . spNew - The duplicate `PetscDualSpace`
465: Level: beginner
467: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
468: @*/
469: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
470: {
471: DM dm;
472: PetscDualSpaceType type;
473: const char *name;
475: PetscFunctionBegin;
477: PetscAssertPointer(spNew, 2);
478: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew));
479: name = ((PetscObject)sp)->name;
480: if (name) PetscCall(PetscObjectSetName((PetscObject)*spNew, name));
481: PetscCall(PetscDualSpaceGetType(sp, &type));
482: PetscCall(PetscDualSpaceSetType(*spNew, type));
483: PetscCall(PetscDualSpaceGetDM(sp, &dm));
484: PetscCall(PetscDualSpaceSetDM(*spNew, dm));
486: (*spNew)->order = sp->order;
487: (*spNew)->k = sp->k;
488: (*spNew)->Nc = sp->Nc;
489: (*spNew)->uniform = sp->uniform;
490: PetscTryTypeMethod(sp, duplicate, *spNew);
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@C
495: PetscDualSpaceGetDM - Get the `DM` representing the reference cell of a `PetscDualSpace`
497: Not Collective
499: Input Parameter:
500: . sp - The `PetscDualSpace`
502: Output Parameter:
503: . dm - The reference cell, that is a `DM` that consists of a single cell
505: Level: intermediate
507: .seealso: `PetscDualSpace`, `PetscDualSpaceSetDM()`, `PetscDualSpaceCreate()`
508: @*/
509: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
510: {
511: PetscFunctionBegin;
513: PetscAssertPointer(dm, 2);
514: *dm = sp->dm;
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@C
519: PetscDualSpaceSetDM - Get the `DM` representing the reference cell
521: Not Collective
523: Input Parameters:
524: + sp - The `PetscDual`Space
525: - dm - The reference cell
527: Level: intermediate
529: .seealso: `PetscDualSpace`, `DM`, `PetscDualSpaceGetDM()`, `PetscDualSpaceCreate()`
530: @*/
531: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
532: {
533: PetscFunctionBegin;
536: PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up");
537: PetscCall(PetscObjectReference((PetscObject)dm));
538: if (sp->dm && sp->dm != dm) PetscCall(PetscDualSpaceClearDMData_Internal(sp, sp->dm));
539: PetscCall(DMDestroy(&sp->dm));
540: sp->dm = dm;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@C
545: PetscDualSpaceGetOrder - Get the order of the dual space
547: Not Collective
549: Input Parameter:
550: . sp - The `PetscDualSpace`
552: Output Parameter:
553: . order - The order
555: Level: intermediate
557: .seealso: `PetscDualSpace`, `PetscDualSpaceSetOrder()`, `PetscDualSpaceCreate()`
558: @*/
559: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
560: {
561: PetscFunctionBegin;
563: PetscAssertPointer(order, 2);
564: *order = sp->order;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@C
569: PetscDualSpaceSetOrder - Set the order of the dual space
571: Not Collective
573: Input Parameters:
574: + sp - The `PetscDualSpace`
575: - order - The order
577: Level: intermediate
579: .seealso: `PetscDualSpace`, `PetscDualSpaceGetOrder()`, `PetscDualSpaceCreate()`
580: @*/
581: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
582: {
583: PetscFunctionBegin;
585: PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up");
586: sp->order = order;
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /*@C
591: PetscDualSpaceGetNumComponents - Return the number of components for this space
593: Input Parameter:
594: . sp - The `PetscDualSpace`
596: Output Parameter:
597: . Nc - The number of components
599: Level: intermediate
601: Note:
602: A vector space, for example, will have d components, where d is the spatial dimension
604: .seealso: `PetscDualSpaceSetNumComponents()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
605: @*/
606: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
607: {
608: PetscFunctionBegin;
610: PetscAssertPointer(Nc, 2);
611: *Nc = sp->Nc;
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@C
616: PetscDualSpaceSetNumComponents - Set the number of components for this space
618: Input Parameters:
619: + sp - The `PetscDualSpace`
620: - Nc - The number of components
622: Level: intermediate
624: .seealso: `PetscDualSpaceGetNumComponents()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
625: @*/
626: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
627: {
628: PetscFunctionBegin;
630: PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
631: sp->Nc = Nc;
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*@C
636: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
638: Not Collective
640: Input Parameters:
641: + sp - The `PetscDualSpace`
642: - i - The basis number
644: Output Parameter:
645: . functional - The basis functional
647: Level: intermediate
649: .seealso: `PetscDualSpace`, `PetscQuadrature`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`
650: @*/
651: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
652: {
653: PetscInt dim;
655: PetscFunctionBegin;
657: PetscAssertPointer(functional, 3);
658: PetscCall(PetscDualSpaceGetDimension(sp, &dim));
659: PetscCheck(!(i < 0) && !(i >= dim), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", i, dim);
660: *functional = sp->functional[i];
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*@C
665: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
667: Not Collective
669: Input Parameter:
670: . sp - The `PetscDualSpace`
672: Output Parameter:
673: . dim - The dimension
675: Level: intermediate
677: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
678: @*/
679: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
680: {
681: PetscFunctionBegin;
683: PetscAssertPointer(dim, 2);
684: if (sp->spdim < 0) {
685: PetscSection section;
687: PetscCall(PetscDualSpaceGetSection(sp, §ion));
688: if (section) PetscCall(PetscSectionGetStorageSize(section, &sp->spdim));
689: else sp->spdim = 0;
690: }
691: *dim = sp->spdim;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@C
696: PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain
698: Not Collective
700: Input Parameter:
701: . sp - The `PetscDualSpace`
703: Output Parameter:
704: . intdim - The dimension
706: Level: intermediate
708: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
709: @*/
710: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
711: {
712: PetscFunctionBegin;
714: PetscAssertPointer(intdim, 2);
715: if (sp->spintdim < 0) {
716: PetscSection section;
718: PetscCall(PetscDualSpaceGetSection(sp, §ion));
719: if (section) PetscCall(PetscSectionGetConstrainedStorageSize(section, &sp->spintdim));
720: else sp->spintdim = 0;
721: }
722: *intdim = sp->spintdim;
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@C
727: PetscDualSpaceGetUniform - Whether this dual space is uniform
729: Not Collective
731: Input Parameter:
732: . sp - A dual space
734: Output Parameter:
735: . uniform - `PETSC_TRUE` if (a) the dual space is the same for each point in a stratum of the reference `DMPLEX`, and
736: (b) every symmetry of each point in the reference `DMPLEX` is also a symmetry of the point's dual space.
738: Level: advanced
740: Note:
741: All of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells
742: with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform.
744: .seealso: `PetscDualSpace`, `PetscDualSpaceGetPointSubspace()`, `PetscDualSpaceGetSymmetries()`
745: @*/
746: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
747: {
748: PetscFunctionBegin;
750: PetscAssertPointer(uniform, 2);
751: *uniform = sp->uniform;
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /*@CC
756: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
758: Not Collective
760: Input Parameter:
761: . sp - The `PetscDualSpace`
763: Output Parameter:
764: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
766: Level: intermediate
768: Note:
769: Do not free `numDof`
771: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
772: @*/
773: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt *numDof[])
774: {
775: PetscFunctionBegin;
777: PetscAssertPointer(numDof, 2);
778: PetscCheck(sp->uniform, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform space does not have a fixed number of dofs for each height");
779: if (!sp->numDof) {
780: DM dm;
781: PetscInt depth;
782: PetscSection section;
784: PetscCall(PetscDualSpaceGetDM(sp, &dm));
785: PetscCall(DMPlexGetDepth(dm, &depth));
786: PetscCall(PetscCalloc1(depth + 1, &sp->numDof));
787: PetscCall(PetscDualSpaceGetSection(sp, §ion));
788: for (PetscInt d = 0; d <= depth; d++) {
789: PetscInt dStart, dEnd;
791: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
792: if (dEnd <= dStart) continue;
793: PetscCall(PetscSectionGetDof(section, dStart, &sp->numDof[d]));
794: }
795: }
796: *numDof = sp->numDof;
797: PetscCheck(*numDof, PetscObjectComm((PetscObject)sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: /* create the section of the right size and set a permutation for topological ordering */
802: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
803: {
804: DM dm;
805: PetscInt pStart, pEnd, cStart, cEnd, c, depth, count, i;
806: PetscInt *seen, *perm;
807: PetscSection section;
809: PetscFunctionBegin;
810: dm = sp->dm;
811: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
812: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
813: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
814: PetscCall(PetscCalloc1(pEnd - pStart, &seen));
815: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
816: PetscCall(DMPlexGetDepth(dm, &depth));
817: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
818: for (c = cStart, count = 0; c < cEnd; c++) {
819: PetscInt closureSize = -1, e;
820: PetscInt *closure = NULL;
822: perm[count++] = c;
823: seen[c - pStart] = 1;
824: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
825: for (e = 0; e < closureSize; e++) {
826: PetscInt point = closure[2 * e];
828: if (seen[point - pStart]) continue;
829: perm[count++] = point;
830: seen[point - pStart] = 1;
831: }
832: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
833: }
834: PetscCheck(count == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering");
835: for (i = 0; i < pEnd - pStart; i++)
836: if (perm[i] != i) break;
837: if (i < pEnd - pStart) {
838: IS permIS;
840: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
841: PetscCall(ISSetPermutation(permIS));
842: PetscCall(PetscSectionSetPermutation(section, permIS));
843: PetscCall(ISDestroy(&permIS));
844: } else {
845: PetscCall(PetscFree(perm));
846: }
847: PetscCall(PetscFree(seen));
848: *topSection = section;
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /* mark boundary points and set up */
853: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
854: {
855: DM dm;
856: DMLabel boundary;
857: PetscInt pStart, pEnd, p;
859: PetscFunctionBegin;
860: dm = sp->dm;
861: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "boundary", &boundary));
862: PetscCall(PetscDualSpaceGetDM(sp, &dm));
863: PetscCall(DMPlexMarkBoundaryFaces(dm, 1, boundary));
864: PetscCall(DMPlexLabelComplete(dm, boundary));
865: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
866: for (p = pStart; p < pEnd; p++) {
867: PetscInt bval;
869: PetscCall(DMLabelGetValue(boundary, p, &bval));
870: if (bval == 1) {
871: PetscInt dof;
873: PetscCall(PetscSectionGetDof(section, p, &dof));
874: PetscCall(PetscSectionSetConstraintDof(section, p, dof));
875: }
876: }
877: PetscCall(DMLabelDestroy(&boundary));
878: PetscCall(PetscSectionSetUp(section));
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: /*@C
883: PetscDualSpaceGetSection - Create a `PetscSection` over the reference cell with the layout from this space
885: Collective
887: Input Parameter:
888: . sp - The `PetscDualSpace`
890: Output Parameter:
891: . section - The section
893: Level: advanced
895: .seealso: `PetscDualSpace`, `PetscSection`, `PetscDualSpaceCreate()`, `DMPLEX`
896: @*/
897: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
898: {
899: PetscInt pStart, pEnd, p;
901: PetscFunctionBegin;
902: if (!sp->dm) {
903: *section = NULL;
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
906: if (!sp->pointSection) {
907: /* mark the boundary */
908: PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->pointSection));
909: PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd));
910: for (p = pStart; p < pEnd; p++) {
911: PetscDualSpace psp;
913: PetscCall(PetscDualSpaceGetPointSubspace(sp, p, &psp));
914: if (psp) {
915: PetscInt dof;
917: PetscCall(PetscDualSpaceGetInteriorDimension(psp, &dof));
918: PetscCall(PetscSectionSetDof(sp->pointSection, p, dof));
919: }
920: }
921: PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));
922: }
923: *section = sp->pointSection;
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: /*@C
928: PetscDualSpaceGetInteriorSection - Create a `PetscSection` over the reference cell with the layout from this space
929: for interior degrees of freedom
931: Collective
933: Input Parameter:
934: . sp - The `PetscDualSpace`
936: Output Parameter:
937: . section - The interior section
939: Level: advanced
941: Note:
942: Most reference domains have one cell, in which case the only cell will have
943: all of the interior degrees of freedom in the interior section. But
944: for `PETSCDUALSPACEREFINED` there may be other mesh points in the interior,
945: and this section describes their layout.
947: .seealso: `PetscDualSpace`, `PetscSection`, `PetscDualSpaceCreate()`, `DMPLEX`
948: @*/
949: PetscErrorCode PetscDualSpaceGetInteriorSection(PetscDualSpace sp, PetscSection *section)
950: {
951: PetscInt pStart, pEnd, p;
953: PetscFunctionBegin;
954: if (!sp->dm) {
955: *section = NULL;
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
958: if (!sp->intPointSection) {
959: PetscSection full_section;
961: PetscCall(PetscDualSpaceGetSection(sp, &full_section));
962: PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->intPointSection));
963: PetscCall(PetscSectionGetChart(full_section, &pStart, &pEnd));
964: for (p = pStart; p < pEnd; p++) {
965: PetscInt dof, cdof;
967: PetscCall(PetscSectionGetDof(full_section, p, &dof));
968: PetscCall(PetscSectionGetConstraintDof(full_section, p, &cdof));
969: PetscCall(PetscSectionSetDof(sp->intPointSection, p, dof - cdof));
970: }
971: PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->intPointSection));
972: }
973: *section = sp->intPointSection;
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
978: * have one cell */
979: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
980: {
981: PetscReal *sv0, *v0, *J;
982: PetscSection section;
983: PetscInt dim, s, k;
984: DM dm;
986: PetscFunctionBegin;
987: PetscCall(PetscDualSpaceGetDM(sp, &dm));
988: PetscCall(DMGetDimension(dm, &dim));
989: PetscCall(PetscDualSpaceGetSection(sp, §ion));
990: PetscCall(PetscMalloc3(dim, &v0, dim, &sv0, dim * dim, &J));
991: PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
992: for (s = sStart; s < sEnd; s++) {
993: PetscReal detJ, hdetJ;
994: PetscDualSpace ssp;
995: PetscInt dof, off, f, sdim;
996: PetscInt i;
997: DM sdm;
999: PetscCall(PetscDualSpaceGetPointSubspace(sp, s, &ssp));
1000: if (!ssp) continue;
1001: PetscCall(PetscSectionGetDof(section, s, &dof));
1002: PetscCall(PetscSectionGetOffset(section, s, &off));
1003: /* get the first vertex of the reference cell */
1004: PetscCall(PetscDualSpaceGetDM(ssp, &sdm));
1005: PetscCall(DMGetDimension(sdm, &sdim));
1006: PetscCall(DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ));
1007: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ));
1008: /* compactify Jacobian */
1009: for (i = 0; i < dim; i++)
1010: for (PetscInt j = 0; j < sdim; j++) J[i * sdim + j] = J[i * dim + j];
1011: for (f = 0; f < dof; f++) {
1012: PetscQuadrature fn;
1014: PetscCall(PetscDualSpaceGetFunctional(ssp, f, &fn));
1015: PetscCall(PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &sp->functional[off + f]));
1016: }
1017: }
1018: PetscCall(PetscFree3(v0, sv0, J));
1019: PetscFunctionReturn(PETSC_SUCCESS);
1020: }
1022: /*@C
1023: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1025: Input Parameters:
1026: + sp - The `PetscDualSpace` object
1027: . f - The basis functional index
1028: . time - The time
1029: . 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)
1030: . numComp - The number of components for the function
1031: . func - The input function
1032: - ctx - A context for the function
1034: Output Parameter:
1035: . value - numComp output values
1037: Calling sequence:
1038: .vb
1039: PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], PetscCtx ctx)
1040: .ve
1042: Level: beginner
1044: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1045: @*/
1046: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscCtx ctx, PetscScalar *value)
1047: {
1048: PetscFunctionBegin;
1050: PetscAssertPointer(cgeom, 4);
1051: PetscAssertPointer(value, 8);
1052: PetscUseTypeMethod(sp, apply, f, time, cgeom, numComp, func, ctx, value);
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1056: /*@C
1057: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by `PetscDualSpaceGetAllData()`
1059: Input Parameters:
1060: + sp - The `PetscDualSpace` object
1061: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`
1063: Output Parameter:
1064: . spValue - The values of all dual space functionals
1066: Level: advanced
1068: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1069: @*/
1070: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1071: {
1072: PetscFunctionBegin;
1074: PetscUseTypeMethod(sp, applyall, pointEval, spValue);
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@C
1079: PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by `PetscDualSpaceGetInteriorData()`
1081: Input Parameters:
1082: + sp - The `PetscDualSpace` object
1083: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`
1085: Output Parameter:
1086: . spValue - The values of interior dual space functionals
1088: Level: advanced
1090: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1091: @*/
1092: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1093: {
1094: PetscFunctionBegin;
1096: PetscUseTypeMethod(sp, applyint, pointEval, spValue);
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /*@C
1101: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
1103: Input Parameters:
1104: + sp - The `PetscDualSpace` object
1105: . f - The basis functional index
1106: . time - The time
1107: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1108: . Nc - The number of components for the function
1109: . func - The input function
1110: - ctx - A context for the function
1112: Output Parameter:
1113: . value - The output value
1115: Calling sequence:
1116: .vb
1117: PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[],PetscInt numComponents, PetscScalar values[], PetscCtx ctx)
1118: .ve
1120: Level: advanced
1122: Note:
1123: The idea is to evaluate the functional as an integral $ n(f) = \int dx n(x) . f(x) $ where both n and f have Nc components.
1125: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1126: @*/
1127: PetscErrorCode PetscDualSpaceApplyDefault(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscCtx ctx, PetscScalar *value)
1128: {
1129: DM dm;
1130: PetscQuadrature n;
1131: const PetscReal *points, *weights;
1132: PetscReal x[3];
1133: PetscScalar *val;
1134: PetscInt dim, dE, qNc, c, Nq, q;
1135: PetscBool isAffine;
1137: PetscFunctionBegin;
1139: PetscAssertPointer(value, 8);
1140: PetscCall(PetscDualSpaceGetDM(sp, &dm));
1141: PetscCall(PetscDualSpaceGetFunctional(sp, f, &n));
1142: PetscCall(PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights));
1143: PetscCheck(dim == cgeom->dim, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %" PetscInt_FMT " != cell geometry dimension %" PetscInt_FMT, dim, cgeom->dim);
1144: PetscCheck(qNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature components %" PetscInt_FMT " != function components %" PetscInt_FMT, qNc, Nc);
1145: PetscCall(DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val));
1146: *value = 0.0;
1147: isAffine = cgeom->isAffine;
1148: dE = cgeom->dimEmbed;
1149: for (q = 0; q < Nq; ++q) {
1150: if (isAffine) {
1151: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q * dim], x);
1152: PetscCall((*func)(dE, time, x, Nc, val, ctx));
1153: } else {
1154: PetscCall((*func)(dE, time, &cgeom->v[dE * q], Nc, val, ctx));
1155: }
1156: for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1157: }
1158: PetscCall(DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val));
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: /*@C
1163: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by `PetscDualSpaceGetAllData()`
1165: Input Parameters:
1166: + sp - The `PetscDualSpace` object
1167: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`
1169: Output Parameter:
1170: . spValue - The values of all dual space functionals
1172: Level: advanced
1174: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1175: @*/
1176: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1177: {
1178: Vec pointValues, dofValues;
1179: Mat allMat;
1181: PetscFunctionBegin;
1183: PetscAssertPointer(pointEval, 2);
1184: PetscAssertPointer(spValue, 3);
1185: PetscCall(PetscDualSpaceGetAllData(sp, NULL, &allMat));
1186: if (!sp->allNodeValues) PetscCall(MatCreateVecs(allMat, &sp->allNodeValues, NULL));
1187: pointValues = sp->allNodeValues;
1188: if (!sp->allDofValues) PetscCall(MatCreateVecs(allMat, NULL, &sp->allDofValues));
1189: dofValues = sp->allDofValues;
1190: PetscCall(VecPlaceArray(pointValues, pointEval));
1191: PetscCall(VecPlaceArray(dofValues, spValue));
1192: PetscCall(MatMult(allMat, pointValues, dofValues));
1193: PetscCall(VecResetArray(dofValues));
1194: PetscCall(VecResetArray(pointValues));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*@C
1199: PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by `PetscDualSpaceGetInteriorData()`
1201: Input Parameters:
1202: + sp - The `PetscDualSpace` object
1203: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`
1205: Output Parameter:
1206: . spValue - The values of interior dual space functionals
1208: Level: advanced
1210: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1211: @*/
1212: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1213: {
1214: Vec pointValues, dofValues;
1215: Mat intMat;
1217: PetscFunctionBegin;
1219: PetscAssertPointer(pointEval, 2);
1220: PetscAssertPointer(spValue, 3);
1221: PetscCall(PetscDualSpaceGetInteriorData(sp, NULL, &intMat));
1222: if (!sp->intNodeValues) PetscCall(MatCreateVecs(intMat, &sp->intNodeValues, NULL));
1223: pointValues = sp->intNodeValues;
1224: if (!sp->intDofValues) PetscCall(MatCreateVecs(intMat, NULL, &sp->intDofValues));
1225: dofValues = sp->intDofValues;
1226: PetscCall(VecPlaceArray(pointValues, pointEval));
1227: PetscCall(VecPlaceArray(dofValues, spValue));
1228: PetscCall(MatMult(intMat, pointValues, dofValues));
1229: PetscCall(VecResetArray(dofValues));
1230: PetscCall(VecResetArray(pointValues));
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: /*@C
1235: PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values
1237: Input Parameter:
1238: . sp - The dualspace
1240: Output Parameters:
1241: + allNodes - A `PetscQuadrature` object containing all evaluation nodes, pass `NULL` if not needed
1242: - allMat - A `Mat` for the node-to-dof transformation, pass `NULL` if not needed
1244: Level: advanced
1246: .seealso: `PetscQuadrature`, `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`
1247: @*/
1248: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PeOp PetscQuadrature *allNodes, PeOp Mat *allMat)
1249: {
1250: PetscFunctionBegin;
1252: if (allNodes) PetscAssertPointer(allNodes, 2);
1253: if (allMat) PetscAssertPointer(allMat, 3);
1254: if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1255: PetscQuadrature qpoints;
1256: Mat amat;
1258: PetscUseTypeMethod(sp, createalldata, &qpoints, &amat);
1259: PetscCall(PetscQuadratureDestroy(&sp->allNodes));
1260: PetscCall(MatDestroy(&sp->allMat));
1261: sp->allNodes = qpoints;
1262: sp->allMat = amat;
1263: }
1264: if (allNodes) *allNodes = sp->allNodes;
1265: if (allMat) *allMat = sp->allMat;
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: /*@C
1270: PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals
1272: Input Parameter:
1273: . sp - The dualspace
1275: Output Parameters:
1276: + allNodes - A `PetscQuadrature` object containing all evaluation nodes
1277: - allMat - A `Mat` for the node-to-dof transformation
1279: Level: advanced
1281: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`, `PetscQuadrature`
1282: @*/
1283: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1284: {
1285: PetscInt spdim;
1286: PetscInt numPoints, offset;
1287: PetscReal *points;
1288: PetscInt f, dim;
1289: PetscInt Nc, nrows, ncols;
1290: PetscInt maxNumPoints;
1291: PetscQuadrature q;
1292: Mat A;
1294: PetscFunctionBegin;
1295: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1296: PetscCall(PetscDualSpaceGetDimension(sp, &spdim));
1297: if (!spdim) {
1298: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allNodes));
1299: PetscCall(PetscQuadratureSetData(*allNodes, 0, 0, 0, NULL, NULL));
1300: }
1301: nrows = spdim;
1302: PetscCall(PetscDualSpaceGetFunctional(sp, 0, &q));
1303: PetscCall(PetscQuadratureGetData(q, &dim, NULL, &numPoints, NULL, NULL));
1304: maxNumPoints = numPoints;
1305: for (f = 1; f < spdim; f++) {
1306: PetscInt Np;
1308: PetscCall(PetscDualSpaceGetFunctional(sp, f, &q));
1309: PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL));
1310: numPoints += Np;
1311: maxNumPoints = PetscMax(maxNumPoints, Np);
1312: }
1313: ncols = numPoints * Nc;
1314: PetscCall(PetscMalloc1(dim * numPoints, &points));
1315: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A));
1316: for (f = 0, offset = 0; f < spdim; f++) {
1317: const PetscReal *p, *w;
1318: PetscInt Np;
1319: PetscInt fnc;
1321: PetscCall(PetscDualSpaceGetFunctional(sp, f, &q));
1322: PetscCall(PetscQuadratureGetData(q, NULL, &fnc, &Np, &p, &w));
1323: PetscCheck(fnc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch");
1324: for (PetscInt i = 0; i < Np * dim; i++) points[offset * dim + i] = p[i];
1325: for (PetscInt i = 0; i < Np * Nc; i++) PetscCall(MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES));
1326: offset += Np;
1327: }
1328: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1329: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1330: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allNodes));
1331: PetscCall(PetscQuadratureSetData(*allNodes, dim, 0, numPoints, points, NULL));
1332: *allMat = A;
1333: PetscFunctionReturn(PETSC_SUCCESS);
1334: }
1336: /*@C
1337: PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1338: this space, as well as the matrix that computes the degrees of freedom from the quadrature
1339: values.
1341: Input Parameter:
1342: . sp - The dualspace
1344: Output Parameters:
1345: + intNodes - A `PetscQuadrature` object containing all evaluation points needed to evaluate interior degrees of freedom,
1346: pass `NULL` if not needed
1347: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1348: the size of the constrained layout (`PetscSectionGetConstrainStorageSize()`) of the dual space section,
1349: npoints is the number of points in intNodes and nc is `PetscDualSpaceGetNumComponents()`.
1350: Pass `NULL` if not needed
1352: Level: advanced
1354: Notes:
1355: Degrees of freedom are interior degrees of freedom if they belong (by
1356: `PetscDualSpaceGetSection()`) to interior points in the references, complementary boundary
1357: degrees of freedom are marked as constrained in the section returned by
1358: `PetscDualSpaceGetSection()`).
1360: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceGetNumComponents()`, `PetscQuadratureGetData()`
1361: @*/
1362: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PeOp PetscQuadrature *intNodes, PeOp Mat *intMat)
1363: {
1364: PetscFunctionBegin;
1366: if (intNodes) PetscAssertPointer(intNodes, 2);
1367: if (intMat) PetscAssertPointer(intMat, 3);
1368: if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1369: PetscQuadrature qpoints;
1370: Mat imat;
1372: PetscUseTypeMethod(sp, createintdata, &qpoints, &imat);
1373: PetscCall(PetscQuadratureDestroy(&sp->intNodes));
1374: PetscCall(MatDestroy(&sp->intMat));
1375: sp->intNodes = qpoints;
1376: sp->intMat = imat;
1377: }
1378: if (intNodes) *intNodes = sp->intNodes;
1379: if (intMat) *intMat = sp->intMat;
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: /*@C
1384: PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values
1386: Input Parameter:
1387: . sp - The dualspace
1389: Output Parameters:
1390: + intNodes - A `PetscQuadrature` object containing all evaluation points needed to evaluate interior degrees of freedom
1391: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1392: the size of the constrained layout (`PetscSectionGetConstrainStorageSize()`) of the dual space section,
1393: npoints is the number of points in allNodes and nc is `PetscDualSpaceGetNumComponents()`.
1395: Level: advanced
1397: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetInteriorData()`
1398: @*/
1399: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1400: {
1401: DM dm;
1402: PetscInt spdim0;
1403: PetscInt Nc;
1404: PetscInt pStart, pEnd, p, f;
1405: PetscSection section;
1406: PetscInt numPoints, offset, matoffset;
1407: PetscReal *points;
1408: PetscInt dim;
1409: PetscInt *nnz;
1410: PetscQuadrature q;
1411: Mat imat;
1413: PetscFunctionBegin;
1415: PetscCall(PetscDualSpaceGetSection(sp, §ion));
1416: PetscCall(PetscSectionGetConstrainedStorageSize(section, &spdim0));
1417: if (!spdim0) {
1418: *intNodes = NULL;
1419: *intMat = NULL;
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1422: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1423: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
1424: PetscCall(PetscDualSpaceGetDM(sp, &dm));
1425: PetscCall(DMGetDimension(dm, &dim));
1426: PetscCall(PetscMalloc1(spdim0, &nnz));
1427: for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1428: PetscInt dof, cdof, off, d;
1430: PetscCall(PetscSectionGetDof(section, p, &dof));
1431: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
1432: if (!(dof - cdof)) continue;
1433: PetscCall(PetscSectionGetOffset(section, p, &off));
1434: for (d = 0; d < dof; d++, off++, f++) {
1435: PetscInt Np;
1437: PetscCall(PetscDualSpaceGetFunctional(sp, off, &q));
1438: PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL));
1439: nnz[f] = Np * Nc;
1440: numPoints += Np;
1441: }
1442: }
1443: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat));
1444: PetscCall(PetscFree(nnz));
1445: PetscCall(PetscMalloc1(dim * numPoints, &points));
1446: for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1447: PetscInt dof, cdof, off, d;
1449: PetscCall(PetscSectionGetDof(section, p, &dof));
1450: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
1451: if (!(dof - cdof)) continue;
1452: PetscCall(PetscSectionGetOffset(section, p, &off));
1453: for (d = 0; d < dof; d++, off++, f++) {
1454: const PetscReal *p;
1455: const PetscReal *w;
1456: PetscInt Np;
1458: PetscCall(PetscDualSpaceGetFunctional(sp, off, &q));
1459: PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, &p, &w));
1460: for (PetscInt i = 0; i < Np * dim; i++) points[offset + i] = p[i];
1461: for (PetscInt i = 0; i < Np * Nc; i++) PetscCall(MatSetValue(imat, f, matoffset + i, w[i], INSERT_VALUES));
1462: offset += Np * dim;
1463: matoffset += Np * Nc;
1464: }
1465: }
1466: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, intNodes));
1467: PetscCall(PetscQuadratureSetData(*intNodes, dim, 0, numPoints, points, NULL));
1468: PetscCall(MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY));
1469: PetscCall(MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY));
1470: *intMat = imat;
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: /*@C
1475: PetscDualSpaceEqual - Determine if two dual spaces are equivalent
1477: Input Parameters:
1478: + A - A `PetscDualSpace` object
1479: - B - Another `PetscDualSpace` object
1481: Output Parameter:
1482: . equal - `PETSC_TRUE` if the dual spaces are equivalent
1484: Level: advanced
1486: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1487: @*/
1488: PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal)
1489: {
1490: PetscInt sizeA, sizeB, dimA, dimB;
1491: const PetscInt *dofA, *dofB;
1492: PetscQuadrature quadA, quadB;
1493: Mat matA, matB;
1495: PetscFunctionBegin;
1498: PetscAssertPointer(equal, 3);
1499: *equal = PETSC_FALSE;
1500: PetscCall(PetscDualSpaceGetDimension(A, &sizeA));
1501: PetscCall(PetscDualSpaceGetDimension(B, &sizeB));
1502: if (sizeB != sizeA) PetscFunctionReturn(PETSC_SUCCESS);
1503: PetscCall(DMGetDimension(A->dm, &dimA));
1504: PetscCall(DMGetDimension(B->dm, &dimB));
1505: if (dimA != dimB) PetscFunctionReturn(PETSC_SUCCESS);
1507: PetscCall(PetscDualSpaceGetNumDof(A, &dofA));
1508: PetscCall(PetscDualSpaceGetNumDof(B, &dofB));
1509: for (PetscInt d = 0; d < dimA; d++) {
1510: if (dofA[d] != dofB[d]) PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: PetscCall(PetscDualSpaceGetInteriorData(A, &quadA, &matA));
1514: PetscCall(PetscDualSpaceGetInteriorData(B, &quadB, &matB));
1515: if (!quadA && !quadB) {
1516: *equal = PETSC_TRUE;
1517: } else if (quadA && quadB) {
1518: PetscCall(PetscQuadratureEqual(quadA, quadB, equal));
1519: if (*equal == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1520: if (!matA && !matB) PetscFunctionReturn(PETSC_SUCCESS);
1521: if (matA && matB) PetscCall(MatEqual(matA, matB, equal));
1522: else *equal = PETSC_FALSE;
1523: }
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: /*@C
1528: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
1530: Input Parameters:
1531: + sp - The `PetscDualSpace` object
1532: . f - The basis functional index
1533: . time - The time
1534: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1535: . Nc - The number of components for the function
1536: . func - The input function
1537: - ctx - A context for the function
1539: Output Parameter:
1540: . value - The output value (scalar)
1542: Calling sequence:
1543: .vb
1544: PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], PetscCtx ctx)
1545: .ve
1547: Level: advanced
1549: Note:
1550: The idea is to evaluate the functional as an integral $ n(f) = \int dx n(x) . f(x)$ where both n and f have Nc components.
1552: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1553: @*/
1554: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscCtx ctx, PetscScalar *value)
1555: {
1556: DM dm;
1557: PetscQuadrature n;
1558: const PetscReal *points, *weights;
1559: PetscScalar *val;
1560: PetscInt dimEmbed, qNc, c, Nq, q;
1562: PetscFunctionBegin;
1564: PetscAssertPointer(value, 8);
1565: PetscCall(PetscDualSpaceGetDM(sp, &dm));
1566: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1567: PetscCall(PetscDualSpaceGetFunctional(sp, f, &n));
1568: PetscCall(PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights));
1569: PetscCheck(qNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature components %" PetscInt_FMT " != function components %" PetscInt_FMT, qNc, Nc);
1570: PetscCall(DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val));
1571: *value = 0.;
1572: for (q = 0; q < Nq; ++q) {
1573: PetscCall((*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx));
1574: for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1575: }
1576: PetscCall(DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val));
1577: PetscFunctionReturn(PETSC_SUCCESS);
1578: }
1580: /*@C
1581: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
1582: given height. This assumes that the reference cell is symmetric over points of this height.
1584: Not Collective
1586: Input Parameters:
1587: + sp - the `PetscDualSpace` object
1588: - height - the height of the mesh point for which the subspace is desired
1590: Output Parameter:
1591: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1592: point, which will be of lesser dimension if height > 0.
1594: Level: advanced
1596: Notes:
1597: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1598: pointwise values are not defined on the element boundaries), or if the implementation of `PetscDualSpace` does not
1599: support extracting subspaces, then `NULL` is returned.
1601: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1603: .seealso: `PetscDualSpace`, `PetscSpaceGetHeightSubspace()`, `PetscDualSpaceGetPointSubspace()`
1604: @*/
1605: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1606: {
1607: PetscInt depth = -1, cStart, cEnd;
1608: DM dm;
1610: PetscFunctionBegin;
1612: PetscAssertPointer(subsp, 3);
1613: PetscCheck(sp->uniform, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform dual space does not have a single dual space at each height");
1614: *subsp = NULL;
1615: dm = sp->dm;
1616: PetscCall(DMPlexGetDepth(dm, &depth));
1617: PetscCheck(height >= 0 && height <= depth, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height");
1618: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1619: if (height == 0 && cEnd == cStart + 1) {
1620: *subsp = sp;
1621: PetscFunctionReturn(PETSC_SUCCESS);
1622: }
1623: if (!sp->heightSpaces) {
1624: PetscCall(PetscCalloc1(depth + 1, &sp->heightSpaces));
1626: for (PetscInt h = 0; h <= depth; h++) {
1627: if (h == 0 && cEnd == cStart + 1) continue;
1628: if (sp->ops->createheightsubspace) PetscUseTypeMethod(sp, createheightsubspace, height, &sp->heightSpaces[h]);
1629: else if (sp->pointSpaces) {
1630: PetscInt hStart, hEnd;
1632: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
1633: if (hEnd > hStart) {
1634: const char *name;
1636: PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[hStart]));
1637: if (sp->pointSpaces[hStart]) {
1638: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
1639: PetscCall(PetscObjectSetName((PetscObject)sp->pointSpaces[hStart], name));
1640: }
1641: sp->heightSpaces[h] = sp->pointSpaces[hStart];
1642: }
1643: }
1644: }
1645: }
1646: *subsp = sp->heightSpaces[height];
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: /*@C
1651: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
1653: Not Collective
1655: Input Parameters:
1656: + sp - the `PetscDualSpace` object
1657: - point - the point (in the dual space's DM) for which the subspace is desired
1659: Output Parameter:
1660: . bdsp - the subspace.
1662: Level: advanced
1664: Notes:
1665: The functionals in the subspace are with respect to the intrinsic geometry of the point,
1666: which will be of lesser dimension if height > 0.
1668: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
1669: defined on the element boundaries), or if the implementation of `PetscDualSpace` does not support extracting
1670: subspaces, then `NULL` is returned.
1672: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1674: .seealso: `PetscDualSpace`, `PetscDualSpaceGetHeightSubspace()`
1675: @*/
1676: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1677: {
1678: PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1679: DM dm;
1681: PetscFunctionBegin;
1683: PetscAssertPointer(bdsp, 3);
1684: *bdsp = NULL;
1685: dm = sp->dm;
1686: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1687: PetscCheck(point >= pStart && point <= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point");
1688: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1689: 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 */
1690: *bdsp = sp;
1691: PetscFunctionReturn(PETSC_SUCCESS);
1692: }
1693: if (!sp->pointSpaces) {
1694: PetscCall(PetscCalloc1(pEnd - pStart, &sp->pointSpaces));
1696: for (PetscInt p = 0; p < pEnd - pStart; p++) {
1697: if (p + pStart == cStart && cEnd == cStart + 1) continue;
1698: if (sp->ops->createpointsubspace) PetscUseTypeMethod(sp, createpointsubspace, p + pStart, &sp->pointSpaces[p]);
1699: else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1700: PetscInt dim, depth, height;
1701: DMLabel label;
1703: PetscCall(DMPlexGetDepth(dm, &dim));
1704: PetscCall(DMPlexGetDepthLabel(dm, &label));
1705: PetscCall(DMLabelGetValue(label, p + pStart, &depth));
1706: height = dim - depth;
1707: PetscCall(PetscDualSpaceGetHeightSubspace(sp, height, &sp->pointSpaces[p]));
1708: PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[p]));
1709: }
1710: }
1711: }
1712: *bdsp = sp->pointSpaces[point - pStart];
1713: PetscFunctionReturn(PETSC_SUCCESS);
1714: }
1716: /*@C
1717: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
1719: Not Collective
1721: Input Parameter:
1722: . sp - the `PetscDualSpace` object
1724: Output Parameters:
1725: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1726: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation
1728: Level: developer
1730: Note:
1731: The permutation and flip arrays are organized in the following way
1732: .vb
1733: perms[p][ornt][dof # on point] = new local dof #
1734: flips[p][ornt][dof # on point] = reversal or not
1735: .ve
1737: .seealso: `PetscDualSpace`
1738: @*/
1739: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1740: {
1741: PetscFunctionBegin;
1743: if (perms) {
1744: PetscAssertPointer(perms, 2);
1745: *perms = NULL;
1746: }
1747: if (flips) {
1748: PetscAssertPointer(flips, 3);
1749: *flips = NULL;
1750: }
1751: PetscTryTypeMethod(sp, getsymmetries, perms, flips);
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: /*@C
1756: PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1757: dual space's functionals.
1759: Input Parameter:
1760: . dsp - The `PetscDualSpace`
1762: Output Parameter:
1763: . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1764: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1765: 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).
1766: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1767: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1768: but are stored as 1-forms.
1770: Level: developer
1772: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1773: @*/
1774: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1775: {
1776: PetscFunctionBeginHot;
1778: PetscAssertPointer(k, 2);
1779: *k = dsp->k;
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: /*@C
1784: PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1785: dual space's functionals.
1787: Input Parameters:
1788: + dsp - The `PetscDualSpace`
1789: - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1790: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1791: 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).
1792: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1793: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1794: but are stored as 1-forms.
1796: Level: developer
1798: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1799: @*/
1800: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1801: {
1802: PetscInt dim;
1804: PetscFunctionBeginHot;
1806: PetscCheck(!dsp->setupcalled, PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
1807: dim = dsp->dm->dim;
1808: PetscCheck((k >= -dim && k <= dim) || k == PETSC_FORM_DEGREE_UNDEFINED, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %" PetscInt_FMT "-form on %" PetscInt_FMT "-dimensional reference cell", PetscAbsInt(k), dim);
1809: dsp->k = k;
1810: PetscFunctionReturn(PETSC_SUCCESS);
1811: }
1813: /*@C
1814: PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space
1816: Input Parameter:
1817: . dsp - The `PetscDualSpace`
1819: Output Parameter:
1820: . k - The simplex dimension
1822: Level: developer
1824: Note:
1825: Currently supported values are
1826: .vb
1827: 0: These are H_1 methods that only transform coordinates
1828: 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1829: 2: These are the same as 1
1830: 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1831: .ve
1833: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1834: @*/
1835: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1836: {
1837: PetscInt dim;
1839: PetscFunctionBeginHot;
1841: PetscAssertPointer(k, 2);
1842: dim = dsp->dm->dim;
1843: if (!dsp->k) *k = IDENTITY_TRANSFORM;
1844: else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1845: else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1846: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1847: PetscFunctionReturn(PETSC_SUCCESS);
1848: }
1850: /*@C
1851: PetscDualSpaceTransform - Transform the function values
1853: Input Parameters:
1854: + dsp - The `PetscDualSpace`
1855: . trans - The type of transform
1856: . isInverse - Flag to invert the transform
1857: . fegeom - The cell geometry
1858: . Nv - The number of function samples
1859: . Nc - The number of function components
1860: - vals - The function values
1862: Output Parameter:
1863: . vals - The transformed function values
1865: Level: intermediate
1867: Note:
1868: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1870: .seealso: `PetscDualSpace`, `PetscDualSpaceTransformGradient()`, `PetscDualSpaceTransformHessian()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1871: @*/
1872: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1873: {
1874: PetscReal Jstar[9] = {0};
1875: PetscInt dim, v, c, Nk;
1877: PetscFunctionBeginHot;
1879: PetscAssertPointer(fegeom, 4);
1880: PetscAssertPointer(vals, 7);
1881: /* TODO: not handling dimEmbed != dim right now */
1882: dim = dsp->dm->dim;
1883: /* No change needed for 0-forms */
1884: if (!dsp->k) PetscFunctionReturn(PETSC_SUCCESS);
1885: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk));
1886: /* TODO: use fegeom->isAffine */
1887: PetscCall(PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar));
1888: for (v = 0; v < Nv; ++v) {
1889: switch (Nk) {
1890: case 1:
1891: for (c = 0; c < Nc; c++) vals[v * Nc + c] *= Jstar[0];
1892: break;
1893: case 2:
1894: for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1895: break;
1896: case 3:
1897: for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1898: break;
1899: default:
1900: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %" PetscInt_FMT " for transformation", Nk);
1901: }
1902: }
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: /*@C
1907: PetscDualSpaceTransformGradient - Transform the function gradient values
1909: Input Parameters:
1910: + dsp - The `PetscDualSpace`
1911: . trans - The type of transform
1912: . isInverse - Flag to invert the transform
1913: . fegeom - The cell geometry
1914: . Nv - The number of function gradient samples
1915: . Nc - The number of function components
1916: - vals - The function gradient values
1918: Output Parameter:
1919: . vals - The transformed function gradient values
1921: Level: intermediate
1923: Note:
1924: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1926: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1927: @*/
1928: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1929: {
1930: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1931: PetscInt v, c, d;
1933: PetscFunctionBeginHot;
1935: PetscAssertPointer(fegeom, 4);
1936: PetscAssertPointer(vals, 7);
1937: PetscAssert(dE > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %" PetscInt_FMT, dE);
1938: /* Transform gradient */
1939: if (dim == dE) {
1940: for (v = 0; v < Nv; ++v) {
1941: for (c = 0; c < Nc; ++c) {
1942: switch (dim) {
1943: case 1:
1944: vals[(v * Nc + c) * dim] *= fegeom->invJ[0];
1945: break;
1946: case 2:
1947: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1948: break;
1949: case 3:
1950: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1951: break;
1952: default:
1953: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1954: }
1955: }
1956: }
1957: } else {
1958: for (v = 0; v < Nv; ++v) {
1959: for (c = 0; c < Nc; ++c) DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v * Nc + c) * dE], &vals[(v * Nc + c) * dE]);
1960: }
1961: }
1962: /* Assume its a vector, otherwise assume its a bunch of scalars */
1963: if (Nc == 1 || Nc != dim) PetscFunctionReturn(PETSC_SUCCESS);
1964: switch (trans) {
1965: case IDENTITY_TRANSFORM:
1966: break;
1967: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1968: if (isInverse) {
1969: for (v = 0; v < Nv; ++v) {
1970: for (d = 0; d < dim; ++d) {
1971: switch (dim) {
1972: case 2:
1973: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1974: break;
1975: case 3:
1976: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1977: break;
1978: default:
1979: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1980: }
1981: }
1982: }
1983: } else {
1984: for (v = 0; v < Nv; ++v) {
1985: for (d = 0; d < dim; ++d) {
1986: switch (dim) {
1987: case 2:
1988: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1989: break;
1990: case 3:
1991: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1992: break;
1993: default:
1994: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1995: }
1996: }
1997: }
1998: }
1999: break;
2000: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2001: if (isInverse) {
2002: for (v = 0; v < Nv; ++v) {
2003: for (d = 0; d < dim; ++d) {
2004: switch (dim) {
2005: case 2:
2006: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2007: break;
2008: case 3:
2009: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2010: break;
2011: default:
2012: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2013: }
2014: for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] *= fegeom->detJ[0];
2015: }
2016: }
2017: } else {
2018: for (v = 0; v < Nv; ++v) {
2019: for (d = 0; d < dim; ++d) {
2020: switch (dim) {
2021: case 2:
2022: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2023: break;
2024: case 3:
2025: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2026: break;
2027: default:
2028: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2029: }
2030: for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] /= fegeom->detJ[0];
2031: }
2032: }
2033: }
2034: break;
2035: }
2036: PetscFunctionReturn(PETSC_SUCCESS);
2037: }
2039: /*@C
2040: PetscDualSpaceTransformHessian - Transform the function Hessian values
2042: Input Parameters:
2043: + dsp - The `PetscDualSpace`
2044: . trans - The type of transform
2045: . isInverse - Flag to invert the transform
2046: . fegeom - The cell geometry
2047: . Nv - The number of function Hessian samples
2048: . Nc - The number of function components
2049: - vals - The function gradient values
2051: Output Parameter:
2052: . vals - The transformed function Hessian values
2054: Level: intermediate
2056: Note:
2057: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2059: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
2060: @*/
2061: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
2062: {
2063: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
2065: PetscFunctionBeginHot;
2067: PetscAssertPointer(fegeom, 4);
2068: PetscAssertPointer(vals, 7);
2069: PetscAssert(dE > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %" PetscInt_FMT, dE);
2070: /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
2071: if (dim == dE) {
2072: for (PetscInt v = 0; v < Nv; ++v) {
2073: for (PetscInt c = 0; c < Nc; ++c) {
2074: switch (dim) {
2075: case 1:
2076: vals[(v * Nc + c) * dim * dim] *= PetscSqr(fegeom->invJ[0]);
2077: break;
2078: case 2:
2079: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
2080: break;
2081: case 3:
2082: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
2083: break;
2084: default:
2085: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2086: }
2087: }
2088: }
2089: } else {
2090: for (PetscInt v = 0; v < Nv; ++v) {
2091: for (PetscInt c = 0; c < Nc; ++c) DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v * Nc + c) * dE * dE], &vals[(v * Nc + c) * dE * dE]);
2092: }
2093: }
2094: /* Assume its a vector, otherwise assume its a bunch of scalars */
2095: if (Nc == 1 || Nc != dim) PetscFunctionReturn(PETSC_SUCCESS);
2096: switch (trans) {
2097: case IDENTITY_TRANSFORM:
2098: break;
2099: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
2100: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2101: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2102: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2103: }
2104: PetscFunctionReturn(PETSC_SUCCESS);
2105: }
2107: /*@C
2108: 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.
2110: Input Parameters:
2111: + dsp - The `PetscDualSpace`
2112: . fegeom - The geometry for this cell
2113: . Nq - The number of function samples
2114: . Nc - The number of function components
2115: - pointEval - The function values
2117: Output Parameter:
2118: . pointEval - The transformed function values
2120: Level: advanced
2122: Notes:
2123: 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.
2125: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2127: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2128: @*/
2129: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2130: {
2131: PetscDualSpaceTransformType trans;
2132: PetscInt k;
2134: PetscFunctionBeginHot;
2136: PetscAssertPointer(fegeom, 2);
2137: PetscAssertPointer(pointEval, 5);
2138: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2139: This determines their transformation properties. */
2140: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2141: switch (k) {
2142: case 0: /* H^1 point evaluations */
2143: trans = IDENTITY_TRANSFORM;
2144: break;
2145: case 1: /* Hcurl preserves tangential edge traces */
2146: trans = COVARIANT_PIOLA_TRANSFORM;
2147: break;
2148: case 2:
2149: case 3: /* Hdiv preserve normal traces */
2150: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2151: break;
2152: default:
2153: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2154: }
2155: PetscCall(PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval));
2156: PetscFunctionReturn(PETSC_SUCCESS);
2157: }
2159: /*@C
2160: 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.
2162: Input Parameters:
2163: + dsp - The `PetscDualSpace`
2164: . fegeom - The geometry for this cell
2165: . Nq - The number of function samples
2166: . Nc - The number of function components
2167: - pointEval - The function values
2169: Output Parameter:
2170: . pointEval - The transformed function values
2172: Level: advanced
2174: Notes:
2175: 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.
2177: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2179: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2180: @*/
2181: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2182: {
2183: PetscDualSpaceTransformType trans;
2184: PetscInt k;
2186: PetscFunctionBeginHot;
2188: PetscAssertPointer(fegeom, 2);
2189: PetscAssertPointer(pointEval, 5);
2190: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2191: This determines their transformation properties. */
2192: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2193: switch (k) {
2194: case 0: /* H^1 point evaluations */
2195: trans = IDENTITY_TRANSFORM;
2196: break;
2197: case 1: /* Hcurl preserves tangential edge traces */
2198: trans = COVARIANT_PIOLA_TRANSFORM;
2199: break;
2200: case 2:
2201: case 3: /* Hdiv preserve normal traces */
2202: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2203: break;
2204: default:
2205: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2206: }
2207: PetscCall(PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2208: PetscFunctionReturn(PETSC_SUCCESS);
2209: }
2211: /*@C
2212: 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.
2214: Input Parameters:
2215: + dsp - The `PetscDualSpace`
2216: . fegeom - The geometry for this cell
2217: . Nq - The number of function gradient samples
2218: . Nc - The number of function components
2219: - pointEval - The function gradient values
2221: Output Parameter:
2222: . pointEval - The transformed function gradient values
2224: Level: advanced
2226: Notes:
2227: 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.
2229: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2231: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2232: @*/
2233: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2234: {
2235: PetscDualSpaceTransformType trans;
2236: PetscInt k;
2238: PetscFunctionBeginHot;
2240: PetscAssertPointer(fegeom, 2);
2241: PetscAssertPointer(pointEval, 5);
2242: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2243: This determines their transformation properties. */
2244: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2245: switch (k) {
2246: case 0: /* H^1 point evaluations */
2247: trans = IDENTITY_TRANSFORM;
2248: break;
2249: case 1: /* Hcurl preserves tangential edge traces */
2250: trans = COVARIANT_PIOLA_TRANSFORM;
2251: break;
2252: case 2:
2253: case 3: /* Hdiv preserve normal traces */
2254: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2255: break;
2256: default:
2257: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2258: }
2259: PetscCall(PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2260: PetscFunctionReturn(PETSC_SUCCESS);
2261: }
2263: /*@C
2264: 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.
2266: Input Parameters:
2267: + dsp - The `PetscDualSpace`
2268: . fegeom - The geometry for this cell
2269: . Nq - The number of function Hessian samples
2270: . Nc - The number of function components
2271: - pointEval - The function gradient values
2273: Output Parameter:
2274: . pointEval - The transformed function Hessian values
2276: Level: advanced
2278: Notes:
2279: 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.
2281: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2283: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2284: @*/
2285: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2286: {
2287: PetscDualSpaceTransformType trans;
2288: PetscInt k;
2290: PetscFunctionBeginHot;
2292: PetscAssertPointer(fegeom, 2);
2293: PetscAssertPointer(pointEval, 5);
2294: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2295: This determines their transformation properties. */
2296: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2297: switch (k) {
2298: case 0: /* H^1 point evaluations */
2299: trans = IDENTITY_TRANSFORM;
2300: break;
2301: case 1: /* Hcurl preserves tangential edge traces */
2302: trans = COVARIANT_PIOLA_TRANSFORM;
2303: break;
2304: case 2:
2305: case 3: /* Hdiv preserve normal traces */
2306: trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2307: break;
2308: default:
2309: SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2310: }
2311: PetscCall(PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2312: PetscFunctionReturn(PETSC_SUCCESS);
2313: }