Actual source code: dualspace.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: PetscClassId PETSCDUALSPACE_CLASSID = 0;
6: PetscFunctionList PetscDualSpaceList = NULL;
7: PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE;
9: const char *const PetscDualSpaceReferenceCells[] = {"SIMPLEX", "TENSOR", "PetscDualSpaceReferenceCell", "PETSCDUALSPACE_REFCELL_",0};
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}, {2,0}.
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: {
31: while (len--) {
32: max -= tup[len];
33: if (!max) {
34: tup[len] = 0;
35: break;
36: }
37: }
38: tup[++len]++;
39: return(0);
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 whos sum is at most 'max'
55: Level: developer
57: .seealso: PetscDualSpaceLatticePointLexicographic_Internal()
58: */
59: PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
60: {
61: PetscInt i;
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: return(0);
73: }
75: /*@C
76: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
78: Not Collective
80: Input Parameters:
81: + name - The name of a new user-defined creation routine
82: - create_func - The creation routine itself
84: Notes:
85: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
87: Sample usage:
88: .vb
89: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
90: .ve
92: Then, your PetscDualSpace type can be chosen with the procedural interface via
93: .vb
94: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
95: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
96: .ve
97: or at runtime via the option
98: .vb
99: -petscdualspace_type my_dual_space
100: .ve
102: Level: advanced
104: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()
106: @*/
107: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
108: {
112: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
113: return(0);
114: }
116: /*@C
117: PetscDualSpaceSetType - Builds a particular PetscDualSpace
119: Collective on sp
121: Input Parameters:
122: + sp - The PetscDualSpace object
123: - name - The kind of space
125: Options Database Key:
126: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
128: Level: intermediate
130: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
131: @*/
132: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
133: {
134: PetscErrorCode (*r)(PetscDualSpace);
135: PetscBool match;
140: PetscObjectTypeCompare((PetscObject) sp, name, &match);
141: if (match) return(0);
143: if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
144: PetscFunctionListFind(PetscDualSpaceList, name, &r);
145: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
147: if (sp->ops->destroy) {
148: (*sp->ops->destroy)(sp);
149: sp->ops->destroy = NULL;
150: }
151: (*r)(sp);
152: PetscObjectChangeTypeName((PetscObject) sp, name);
153: return(0);
154: }
156: /*@C
157: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
159: Not Collective
161: Input Parameter:
162: . sp - The PetscDualSpace
164: Output Parameter:
165: . name - The PetscDualSpace type name
167: Level: intermediate
169: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
170: @*/
171: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
172: {
178: if (!PetscDualSpaceRegisterAllCalled) {
179: PetscDualSpaceRegisterAll();
180: }
181: *name = ((PetscObject) sp)->type_name;
182: return(0);
183: }
185: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
186: {
187: PetscViewerFormat format;
188: PetscInt pdim, f;
189: PetscErrorCode ierr;
192: PetscDualSpaceGetDimension(sp, &pdim);
193: PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);
194: PetscViewerASCIIPushTab(v);
195: PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
196: if (sp->ops->view) {(*sp->ops->view)(sp, v);}
197: PetscViewerGetFormat(v, &format);
198: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
199: PetscViewerASCIIPushTab(v);
200: for (f = 0; f < pdim; ++f) {
201: PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
202: PetscViewerASCIIPushTab(v);
203: PetscQuadratureView(sp->functional[f], v);
204: PetscViewerASCIIPopTab(v);
205: }
206: PetscViewerASCIIPopTab(v);
207: }
208: PetscViewerASCIIPopTab(v);
209: return(0);
210: }
212: /*@
213: PetscDualSpaceView - Views a PetscDualSpace
215: Collective on sp
217: Input Parameter:
218: + sp - the PetscDualSpace object to view
219: - v - the viewer
221: Level: beginner
223: .seealso PetscDualSpaceDestroy()
224: @*/
225: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
226: {
227: PetscBool iascii;
233: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
234: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
235: if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
236: return(0);
237: }
239: /*@
240: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
242: Collective on sp
244: Input Parameter:
245: . sp - the PetscDualSpace object to set options for
247: Options Database:
248: . -petscspace_degree the approximation order of the space
250: Level: intermediate
252: .seealso PetscDualSpaceView()
253: @*/
254: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
255: {
256: PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
257: PetscInt refDim = 0;
258: PetscBool flg;
259: const char *defaultType;
260: char name[256];
261: PetscErrorCode ierr;
265: if (!((PetscObject) sp)->type_name) {
266: defaultType = PETSCDUALSPACELAGRANGE;
267: } else {
268: defaultType = ((PetscObject) sp)->type_name;
269: }
270: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
272: PetscObjectOptionsBegin((PetscObject) sp);
273: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
274: if (flg) {
275: PetscDualSpaceSetType(sp, name);
276: } else if (!((PetscObject) sp)->type_name) {
277: PetscDualSpaceSetType(sp, defaultType);
278: }
279: PetscOptionsBoundedInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
280: PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
281: if (sp->ops->setfromoptions) {
282: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
283: }
284: PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
285: PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
286: if (flg) {
287: DM K;
289: if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
290: PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
291: PetscDualSpaceSetDM(sp, K);
292: DMDestroy(&K);
293: }
295: /* process any options handlers added with PetscObjectAddOptionsHandler() */
296: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
297: PetscOptionsEnd();
298: sp->setfromoptionscalled = PETSC_TRUE;
299: return(0);
300: }
302: /*@
303: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
305: Collective on sp
307: Input Parameter:
308: . sp - the PetscDualSpace object to setup
310: Level: intermediate
312: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
313: @*/
314: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
315: {
320: if (sp->setupcalled) return(0);
321: sp->setupcalled = PETSC_TRUE;
322: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
323: if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
324: return(0);
325: }
327: /*@
328: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
330: Collective on sp
332: Input Parameter:
333: . sp - the PetscDualSpace object to destroy
335: Level: beginner
337: .seealso PetscDualSpaceView()
338: @*/
339: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
340: {
341: PetscInt dim, f;
345: if (!*sp) return(0);
348: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
349: ((PetscObject) (*sp))->refct = 0;
351: PetscDualSpaceGetDimension(*sp, &dim);
352: for (f = 0; f < dim; ++f) {
353: PetscQuadratureDestroy(&(*sp)->functional[f]);
354: }
355: PetscFree((*sp)->functional);
356: PetscQuadratureDestroy(&(*sp)->allPoints);
357: DMDestroy(&(*sp)->dm);
359: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
360: PetscHeaderDestroy(sp);
361: return(0);
362: }
364: /*@
365: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
367: Collective
369: Input Parameter:
370: . comm - The communicator for the PetscDualSpace object
372: Output Parameter:
373: . sp - The PetscDualSpace object
375: Level: beginner
377: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
378: @*/
379: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
380: {
381: PetscDualSpace s;
386: PetscCitationsRegister(FECitation,&FEcite);
387: *sp = NULL;
388: PetscFEInitializePackage();
390: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
392: s->order = 0;
393: s->Nc = 1;
394: s->k = 0;
395: s->setupcalled = PETSC_FALSE;
397: *sp = s;
398: return(0);
399: }
401: /*@
402: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
404: Collective on sp
406: Input Parameter:
407: . sp - The original PetscDualSpace
409: Output Parameter:
410: . spNew - The duplicate PetscDualSpace
412: Level: beginner
414: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
415: @*/
416: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
417: {
423: (*sp->ops->duplicate)(sp, spNew);
424: return(0);
425: }
427: /*@
428: PetscDualSpaceGetDM - Get the DM representing the reference cell
430: Not collective
432: Input Parameter:
433: . sp - The PetscDualSpace
435: Output Parameter:
436: . dm - The reference cell
438: Level: intermediate
440: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
441: @*/
442: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
443: {
447: *dm = sp->dm;
448: return(0);
449: }
451: /*@
452: PetscDualSpaceSetDM - Get the DM representing the reference cell
454: Not collective
456: Input Parameters:
457: + sp - The PetscDualSpace
458: - dm - The reference cell
460: Level: intermediate
462: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
463: @*/
464: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
465: {
471: DMDestroy(&sp->dm);
472: PetscObjectReference((PetscObject) dm);
473: sp->dm = dm;
474: return(0);
475: }
477: /*@
478: PetscDualSpaceGetOrder - Get the order of the dual space
480: Not collective
482: Input Parameter:
483: . sp - The PetscDualSpace
485: Output Parameter:
486: . order - The order
488: Level: intermediate
490: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
491: @*/
492: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
493: {
497: *order = sp->order;
498: return(0);
499: }
501: /*@
502: PetscDualSpaceSetOrder - Set the order of the dual space
504: Not collective
506: Input Parameters:
507: + sp - The PetscDualSpace
508: - order - The order
510: Level: intermediate
512: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
513: @*/
514: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
515: {
518: sp->order = order;
519: return(0);
520: }
522: /*@
523: PetscDualSpaceGetNumComponents - Return the number of components for this space
525: Input Parameter:
526: . sp - The PetscDualSpace
528: Output Parameter:
529: . Nc - The number of components
531: Note: A vector space, for example, will have d components, where d is the spatial dimension
533: Level: intermediate
535: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
536: @*/
537: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
538: {
542: *Nc = sp->Nc;
543: return(0);
544: }
546: /*@
547: PetscDualSpaceSetNumComponents - Set the number of components for this space
549: Input Parameters:
550: + sp - The PetscDualSpace
551: - order - The number of components
553: Level: intermediate
555: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
556: @*/
557: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
558: {
561: sp->Nc = Nc;
562: return(0);
563: }
565: /*@
566: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
568: Not collective
570: Input Parameters:
571: + sp - The PetscDualSpace
572: - i - The basis number
574: Output Parameter:
575: . functional - The basis functional
577: Level: intermediate
579: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
580: @*/
581: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
582: {
583: PetscInt dim;
589: PetscDualSpaceGetDimension(sp, &dim);
590: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
591: *functional = sp->functional[i];
592: return(0);
593: }
595: /*@
596: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
598: Not collective
600: Input Parameter:
601: . sp - The PetscDualSpace
603: Output Parameter:
604: . dim - The dimension
606: Level: intermediate
608: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
609: @*/
610: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
611: {
617: *dim = 0;
618: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
619: return(0);
620: }
622: /*@C
623: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
625: Not collective
627: Input Parameter:
628: . sp - The PetscDualSpace
630: Output Parameter:
631: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
633: Level: intermediate
635: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
636: @*/
637: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
638: {
644: (*sp->ops->getnumdof)(sp, numDof);
645: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
646: return(0);
647: }
649: /*@
650: PetscDualSpaceCreateSection - Create a PetscSection over the reference cell with the layout from this space
652: Collective on sp
654: Input Parameters:
655: + sp - The PetscDualSpace
657: Output Parameter:
658: . section - The section
660: Level: advanced
662: .seealso: PetscDualSpaceCreate(), DMPLEX
663: @*/
664: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
665: {
666: DM dm;
667: PetscInt pStart, pEnd, depth, h, offset;
668: const PetscInt *numDof;
672: PetscDualSpaceGetDM(sp,&dm);
673: DMPlexGetChart(dm,&pStart,&pEnd);
674: PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
675: PetscSectionSetChart(*section,pStart,pEnd);
676: DMPlexGetDepth(dm,&depth);
677: PetscDualSpaceGetNumDof(sp,&numDof);
678: for (h = 0; h <= depth; h++) {
679: PetscInt hStart, hEnd, p, dof;
681: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
682: dof = numDof[depth - h];
683: for (p = hStart; p < hEnd; p++) {
684: PetscSectionSetDof(*section,p,dof);
685: }
686: }
687: PetscSectionSetUp(*section);
688: for (h = 0, offset = 0; h <= depth; h++) {
689: PetscInt hStart, hEnd, p, dof;
691: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
692: dof = numDof[depth - h];
693: for (p = hStart; p < hEnd; p++) {
694: PetscSectionGetDof(*section,p,&dof);
695: PetscSectionSetOffset(*section,p,offset);
696: offset += dof;
697: }
698: }
699: return(0);
700: }
702: /*@
703: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
705: Collective on sp
707: Input Parameters:
708: + sp - The PetscDualSpace
709: . dim - The spatial dimension
710: - simplex - Flag for simplex, otherwise use a tensor-product cell
712: Output Parameter:
713: . refdm - The reference cell
715: Level: intermediate
717: .seealso: PetscDualSpaceCreate(), DMPLEX
718: @*/
719: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
720: {
724: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
725: return(0);
726: }
728: /*@C
729: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
731: Input Parameters:
732: + sp - The PetscDualSpace object
733: . f - The basis functional index
734: . time - The time
735: . 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)
736: . numComp - The number of components for the function
737: . func - The input function
738: - ctx - A context for the function
740: Output Parameter:
741: . value - numComp output values
743: Note: The calling sequence for the callback func is given by:
745: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
746: $ PetscInt numComponents, PetscScalar values[], void *ctx)
748: Level: beginner
750: .seealso: PetscDualSpaceCreate()
751: @*/
752: 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)
753: {
760: (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
761: return(0);
762: }
764: /*@C
765: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
767: Input Parameters:
768: + sp - The PetscDualSpace object
769: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
771: Output Parameter:
772: . spValue - The values of all dual space functionals
774: Level: beginner
776: .seealso: PetscDualSpaceCreate()
777: @*/
778: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
779: {
784: (*sp->ops->applyall)(sp, pointEval, spValue);
785: return(0);
786: }
788: /*@C
789: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
791: Input Parameters:
792: + sp - The PetscDualSpace object
793: . f - The basis functional index
794: . time - The time
795: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
796: . Nc - The number of components for the function
797: . func - The input function
798: - ctx - A context for the function
800: Output Parameter:
801: . value - The output value
803: Note: The calling sequence for the callback func is given by:
805: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
806: $ PetscInt numComponents, PetscScalar values[], void *ctx)
808: and the idea is to evaluate the functional as an integral
810: $ n(f) = int dx n(x) . f(x)
812: where both n and f have Nc components.
814: Level: beginner
816: .seealso: PetscDualSpaceCreate()
817: @*/
818: 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)
819: {
820: DM dm;
821: PetscQuadrature n;
822: const PetscReal *points, *weights;
823: PetscReal x[3];
824: PetscScalar *val;
825: PetscInt dim, dE, qNc, c, Nq, q;
826: PetscBool isAffine;
827: PetscErrorCode ierr;
832: PetscDualSpaceGetDM(sp, &dm);
833: PetscDualSpaceGetFunctional(sp, f, &n);
834: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
835: if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
836: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
837: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
838: *value = 0.0;
839: isAffine = cgeom->isAffine;
840: dE = cgeom->dimEmbed;
841: for (q = 0; q < Nq; ++q) {
842: if (isAffine) {
843: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
844: (*func)(dE, time, x, Nc, val, ctx);
845: } else {
846: (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
847: }
848: for (c = 0; c < Nc; ++c) {
849: *value += val[c]*weights[q*Nc+c];
850: }
851: }
852: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
853: return(0);
854: }
856: /*@C
857: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
859: Input Parameters:
860: + sp - The PetscDualSpace object
861: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
863: Output Parameter:
864: . spValue - The values of all dual space functionals
866: Level: beginner
868: .seealso: PetscDualSpaceCreate()
869: @*/
870: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
871: {
872: PetscQuadrature n;
873: const PetscReal *points, *weights;
874: PetscInt qNc, c, Nq, q, f, spdim, Nc;
875: PetscInt offset;
876: PetscErrorCode ierr;
882: PetscDualSpaceGetDimension(sp, &spdim);
883: PetscDualSpaceGetNumComponents(sp, &Nc);
884: for (f = 0, offset = 0; f < spdim; f++) {
885: PetscDualSpaceGetFunctional(sp, f, &n);
886: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
887: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
888: spValue[f] = 0.0;
889: for (q = 0; q < Nq; ++q) {
890: for (c = 0; c < Nc; ++c) {
891: spValue[f] += pointEval[offset++]*weights[q*Nc+c];
892: }
893: }
894: }
895: return(0);
896: }
898: /*@
899: PetscDualSpaceGetAllPoints - Get all quadrature points from this space
901: Input Parameter:
902: . sp - The dualspace
904: Output Parameter:
905: . allPoints - A PetscQuadrature object containing all evaluation points
907: Level: advanced
909: .seealso: PetscDualSpaceCreate()
910: @*/
911: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
912: {
918: if (!sp->allPoints && sp->ops->createallpoints) {
919: (*sp->ops->createallpoints)(sp,&sp->allPoints);
920: }
921: *allPoints = sp->allPoints;
922: return(0);
923: }
925: /*@
926: PetscDualSpaceCreateAllPointsDefault - Create all evaluation points by examining functionals
928: Input Parameter:
929: . sp - The dualspace
931: Output Parameter:
932: . allPoints - A PetscQuadrature object containing all evaluation points
934: Level: advanced
936: .seealso: PetscDualSpaceCreate()
937: @*/
938: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
939: {
940: PetscInt spdim;
941: PetscInt numPoints, offset;
942: PetscReal *points;
943: PetscInt f, dim;
944: PetscQuadrature q;
945: PetscErrorCode ierr;
948: PetscDualSpaceGetDimension(sp,&spdim);
949: if (!spdim) {
950: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
951: PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
952: }
953: PetscDualSpaceGetFunctional(sp,0,&q);
954: PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
955: for (f = 1; f < spdim; f++) {
956: PetscInt Np;
958: PetscDualSpaceGetFunctional(sp,f,&q);
959: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
960: numPoints += Np;
961: }
962: PetscMalloc1(dim*numPoints,&points);
963: for (f = 0, offset = 0; f < spdim; f++) {
964: const PetscReal *p;
965: PetscInt Np, i;
967: PetscDualSpaceGetFunctional(sp,f,&q);
968: PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
969: for (i = 0; i < Np * dim; i++) {
970: points[offset + i] = p[i];
971: }
972: offset += Np * dim;
973: }
974: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
975: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
976: return(0);
977: }
979: /*@C
980: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
982: Input Parameters:
983: + sp - The PetscDualSpace object
984: . f - The basis functional index
985: . time - The time
986: . cgeom - A context with geometric information for this cell, we currently just use the centroid
987: . Nc - The number of components for the function
988: . func - The input function
989: - ctx - A context for the function
991: Output Parameter:
992: . value - The output value (scalar)
994: Note: The calling sequence for the callback func is given by:
996: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
997: $ PetscInt numComponents, PetscScalar values[], void *ctx)
999: and the idea is to evaluate the functional as an integral
1001: $ n(f) = int dx n(x) . f(x)
1003: where both n and f have Nc components.
1005: Level: beginner
1007: .seealso: PetscDualSpaceCreate()
1008: @*/
1009: 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)
1010: {
1011: DM dm;
1012: PetscQuadrature n;
1013: const PetscReal *points, *weights;
1014: PetscScalar *val;
1015: PetscInt dimEmbed, qNc, c, Nq, q;
1016: PetscErrorCode ierr;
1021: PetscDualSpaceGetDM(sp, &dm);
1022: DMGetCoordinateDim(dm, &dimEmbed);
1023: PetscDualSpaceGetFunctional(sp, f, &n);
1024: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1025: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1026: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1027: *value = 0.;
1028: for (q = 0; q < Nq; ++q) {
1029: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1030: for (c = 0; c < Nc; ++c) {
1031: *value += val[c]*weights[q*Nc+c];
1032: }
1033: }
1034: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1035: return(0);
1036: }
1038: /*@
1039: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
1040: given height. This assumes that the reference cell is symmetric over points of this height.
1042: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1043: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1044: support extracting subspaces, then NULL is returned.
1046: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1048: Not collective
1050: Input Parameters:
1051: + sp - the PetscDualSpace object
1052: - height - the height of the mesh point for which the subspace is desired
1054: Output Parameter:
1055: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1056: point, which will be of lesser dimension if height > 0.
1058: Level: advanced
1060: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1061: @*/
1062: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1063: {
1069: *subsp = NULL;
1070: if (sp->ops->getheightsubspace) {(*sp->ops->getheightsubspace)(sp, height, subsp);}
1071: return(0);
1072: }
1074: /*@
1075: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
1077: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
1078: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
1079: subspaces, then NULL is returned.
1081: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1083: Not collective
1085: Input Parameters:
1086: + sp - the PetscDualSpace object
1087: - point - the point (in the dual space's DM) for which the subspace is desired
1089: Output Parameters:
1090: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1091: point, which will be of lesser dimension if height > 0.
1093: Level: advanced
1095: .seealso: PetscDualSpace
1096: @*/
1097: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1098: {
1104: *bdsp = NULL;
1105: if (sp->ops->getpointsubspace) {
1106: (*sp->ops->getpointsubspace)(sp,point,bdsp);
1107: } else if (sp->ops->getheightsubspace) {
1108: DM dm;
1109: DMLabel label;
1110: PetscInt dim, depth, height;
1112: PetscDualSpaceGetDM(sp,&dm);
1113: DMPlexGetDepth(dm,&dim);
1114: DMPlexGetDepthLabel(dm,&label);
1115: DMLabelGetValue(label,point,&depth);
1116: height = dim - depth;
1117: (*sp->ops->getheightsubspace)(sp,height,bdsp);
1118: }
1119: return(0);
1120: }
1122: /*@C
1123: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
1125: Not collective
1127: Input Parameter:
1128: . sp - the PetscDualSpace object
1130: Output Parameters:
1131: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
1132: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation
1134: Note: The permutation and flip arrays are organized in the following way
1135: $ perms[p][ornt][dof # on point] = new local dof #
1136: $ flips[p][ornt][dof # on point] = reversal or not
1138: Level: developer
1140: .seealso: PetscDualSpaceSetSymmetries()
1141: @*/
1142: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1143: {
1150: if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1151: return(0);
1152: }
1154: /*@
1155: PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space
1157: Input Parameter:
1158: . dsp - The PetscDualSpace
1160: Output Parameter:
1161: . k - The simplex dimension
1163: Level: developer
1165: Note: Currently supported values are
1166: $ 0: These are H_1 methods that only transform coordinates
1167: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1168: $ 2: These are the same as 1
1169: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1171: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1172: @*/
1173: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1174: {
1178: *k = dsp->k;
1179: return(0);
1180: }
1182: /*@C
1183: PetscDualSpaceTransform - Transform the function values
1185: Input Parameters:
1186: + dsp - The PetscDualSpace
1187: . trans - The type of transform
1188: . isInverse - Flag to invert the transform
1189: . fegeom - The cell geometry
1190: . Nv - The number of function samples
1191: . Nc - The number of function components
1192: - vals - The function values
1194: Output Parameter:
1195: . vals - The transformed function values
1197: Level: intermediate
1199: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1200: @*/
1201: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1202: {
1203: PetscInt dim, v, c;
1209: dim = dsp->dm->dim;
1210: /* Assume its a vector, otherwise assume its a bunch of scalars */
1211: if (Nc == 1 || Nc != dim) return(0);
1212: switch (trans) {
1213: case IDENTITY_TRANSFORM: break;
1214: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1215: if (isInverse) {
1216: for (v = 0; v < Nv; ++v) {
1217: switch (dim)
1218: {
1219: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1220: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1221: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1222: }
1223: }
1224: } else {
1225: for (v = 0; v < Nv; ++v) {
1226: switch (dim)
1227: {
1228: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1229: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1230: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1231: }
1232: }
1233: }
1234: break;
1235: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1236: if (isInverse) {
1237: for (v = 0; v < Nv; ++v) {
1238: switch (dim)
1239: {
1240: case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1241: case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1242: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1243: }
1244: for (c = 0; c < Nc; ++c) vals[v*Nc+c] *= fegeom->detJ[0];
1245: }
1246: } else {
1247: for (v = 0; v < Nv; ++v) {
1248: switch (dim)
1249: {
1250: case 2: DMPlex_Mult2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1251: case 3: DMPlex_Mult3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1252: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1253: }
1254: for (c = 0; c < Nc; ++c) vals[v*Nc+c] /= fegeom->detJ[0];
1255: }
1256: }
1257: break;
1258: }
1259: return(0);
1260: }
1261: /*@C
1262: PetscDualSpaceTransformGradient - Transform the function gradient values
1264: Input Parameters:
1265: + dsp - The PetscDualSpace
1266: . trans - The type of transform
1267: . isInverse - Flag to invert the transform
1268: . fegeom - The cell geometry
1269: . Nv - The number of function gradient samples
1270: . Nc - The number of function components
1271: - vals - The function gradient values
1273: Output Parameter:
1274: . vals - The transformed function values
1276: Level: intermediate
1278: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1279: @*/
1280: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1281: {
1282: PetscInt dim, v, c, d;
1288: dim = dsp->dm->dim;
1289: /* Transform gradient */
1290: for (v = 0; v < Nv; ++v) {
1291: for (c = 0; c < Nc; ++c) {
1292: switch (dim)
1293: {
1294: case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];
1295: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1296: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1297: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1298: }
1299: }
1300: }
1301: /* Assume its a vector, otherwise assume its a bunch of scalars */
1302: if (Nc == 1 || Nc != dim) return(0);
1303: switch (trans) {
1304: case IDENTITY_TRANSFORM: break;
1305: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1306: if (isInverse) {
1307: for (v = 0; v < Nv; ++v) {
1308: for (d = 0; d < dim; ++d) {
1309: switch (dim)
1310: {
1311: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1312: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1313: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1314: }
1315: }
1316: }
1317: } else {
1318: for (v = 0; v < Nv; ++v) {
1319: for (d = 0; d < dim; ++d) {
1320: switch (dim)
1321: {
1322: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1323: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1324: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1325: }
1326: }
1327: }
1328: }
1329: break;
1330: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1331: if (isInverse) {
1332: for (v = 0; v < Nv; ++v) {
1333: for (d = 0; d < dim; ++d) {
1334: switch (dim)
1335: {
1336: case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1337: case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1338: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1339: }
1340: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1341: }
1342: }
1343: } else {
1344: for (v = 0; v < Nv; ++v) {
1345: for (d = 0; d < dim; ++d) {
1346: switch (dim)
1347: {
1348: case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1349: case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1350: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1351: }
1352: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1353: }
1354: }
1355: }
1356: break;
1357: }
1358: return(0);
1359: }
1361: /*@C
1362: 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.
1364: Input Parameters:
1365: + dsp - The PetscDualSpace
1366: . fegeom - The geometry for this cell
1367: . Nq - The number of function samples
1368: . Nc - The number of function components
1369: - pointEval - The function values
1371: Output Parameter:
1372: . pointEval - The transformed function values
1374: Level: advanced
1376: 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.
1378: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1379: @*/
1380: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1381: {
1382: PetscDualSpaceTransformType trans;
1383: PetscErrorCode ierr;
1389: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1390: This determines their transformation properties. */
1391: switch (dsp->k)
1392: {
1393: case 0: /* H^1 point evaluations */
1394: trans = IDENTITY_TRANSFORM;break;
1395: case 1: /* Hcurl preserves tangential edge traces */
1396: case 2:
1397: trans = COVARIANT_PIOLA_TRANSFORM;break;
1398: case 3: /* Hdiv preserve normal traces */
1399: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1400: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1401: }
1402: PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
1403: return(0);
1404: }
1406: /*@C
1407: 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.
1409: Input Parameters:
1410: + dsp - The PetscDualSpace
1411: . fegeom - The geometry for this cell
1412: . Nq - The number of function samples
1413: . Nc - The number of function components
1414: - pointEval - The function values
1416: Output Parameter:
1417: . pointEval - The transformed function values
1419: Level: advanced
1421: 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.
1423: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1424: @*/
1425: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1426: {
1427: PetscDualSpaceTransformType trans;
1428: PetscErrorCode ierr;
1434: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1435: This determines their transformation properties. */
1436: switch (dsp->k)
1437: {
1438: case 0: /* H^1 point evaluations */
1439: trans = IDENTITY_TRANSFORM;break;
1440: case 1: /* Hcurl preserves tangential edge traces */
1441: case 2:
1442: trans = COVARIANT_PIOLA_TRANSFORM;break;
1443: case 3: /* Hdiv preserve normal traces */
1444: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1445: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1446: }
1447: PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1448: return(0);
1449: }
1451: /*@C
1452: 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.
1454: Input Parameters:
1455: + dsp - The PetscDualSpace
1456: . fegeom - The geometry for this cell
1457: . Nq - The number of function gradient samples
1458: . Nc - The number of function components
1459: - pointEval - The function gradient values
1461: Output Parameter:
1462: . pointEval - The transformed function gradient values
1464: Level: advanced
1466: 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.
1468: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1469: @*/
1470: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1471: {
1472: PetscDualSpaceTransformType trans;
1473: PetscErrorCode ierr;
1479: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1480: This determines their transformation properties. */
1481: switch (dsp->k)
1482: {
1483: case 0: /* H^1 point evaluations */
1484: trans = IDENTITY_TRANSFORM;break;
1485: case 1: /* Hcurl preserves tangential edge traces */
1486: case 2:
1487: trans = COVARIANT_PIOLA_TRANSFORM;break;
1488: case 3: /* Hdiv preserve normal traces */
1489: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1490: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1491: }
1492: PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1493: return(0);
1494: }