Actual source code: dualspace.c
petsc-3.13.6 2020-09-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}, {0,2}.
16: Input Parameters:
17: + len - The length of the tuple
18: . max - The maximum sum
19: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
21: Output Parameter:
22: . tup - A tuple of len integers whos sum is at most 'max'
24: Level: developer
26: .seealso: PetscDualSpaceTensorPointLexicographic_Internal()
27: */
28: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
29: {
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: if (sp->k) {
196: PetscViewerASCIIPrintf(v, "Dual space for %D-forms %swith %D components, size %D\n", PetscAbsInt(sp->k), sp->k < 0 ? "(stored in dual form) ": "", sp->Nc, pdim);
197: } else {
198: PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
199: }
200: if (sp->ops->view) {(*sp->ops->view)(sp, v);}
201: PetscViewerGetFormat(v, &format);
202: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
203: PetscViewerASCIIPushTab(v);
204: for (f = 0; f < pdim; ++f) {
205: PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
206: PetscViewerASCIIPushTab(v);
207: PetscQuadratureView(sp->functional[f], v);
208: PetscViewerASCIIPopTab(v);
209: }
210: PetscViewerASCIIPopTab(v);
211: }
212: PetscViewerASCIIPopTab(v);
213: return(0);
214: }
216: /*@C
217: PetscDualSpaceViewFromOptions - View from Options
219: Collective on PetscDualSpace
221: Input Parameters:
222: + A - the PetscDualSpace object
223: . obj - Optional object, proivides prefix
224: - name - command line option
226: Level: intermediate
227: .seealso: PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate()
228: @*/
229: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[])
230: {
235: PetscObjectViewFromOptions((PetscObject)A,obj,name);
236: return(0);
237: }
239: /*@
240: PetscDualSpaceView - Views a PetscDualSpace
242: Collective on sp
244: Input Parameter:
245: + sp - the PetscDualSpace object to view
246: - v - the viewer
248: Level: beginner
250: .seealso PetscDualSpaceDestroy(), PetscDualSpace
251: @*/
252: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
253: {
254: PetscBool iascii;
260: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
261: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
262: if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
263: return(0);
264: }
266: /*@
267: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
269: Collective on sp
271: Input Parameter:
272: . sp - the PetscDualSpace object to set options for
274: Options Database:
275: . -petscspace_degree the approximation order of the space
277: Level: intermediate
279: .seealso PetscDualSpaceView(), PetscDualSpace, PetscObjectSetFromOptions()
280: @*/
281: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
282: {
283: PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
284: PetscInt refDim = 0;
285: PetscBool flg;
286: const char *defaultType;
287: char name[256];
288: PetscErrorCode ierr;
292: if (!((PetscObject) sp)->type_name) {
293: defaultType = PETSCDUALSPACELAGRANGE;
294: } else {
295: defaultType = ((PetscObject) sp)->type_name;
296: }
297: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
299: PetscObjectOptionsBegin((PetscObject) sp);
300: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
301: if (flg) {
302: PetscDualSpaceSetType(sp, name);
303: } else if (!((PetscObject) sp)->type_name) {
304: PetscDualSpaceSetType(sp, defaultType);
305: }
306: PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
307: PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);
308: PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
309: if (sp->ops->setfromoptions) {
310: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
311: }
312: PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
313: PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
314: if (flg) {
315: DM K;
317: if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
318: PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
319: PetscDualSpaceSetDM(sp, K);
320: DMDestroy(&K);
321: }
323: /* process any options handlers added with PetscObjectAddOptionsHandler() */
324: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
325: PetscOptionsEnd();
326: sp->setfromoptionscalled = PETSC_TRUE;
327: return(0);
328: }
330: /*@
331: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
333: Collective on sp
335: Input Parameter:
336: . sp - the PetscDualSpace object to setup
338: Level: intermediate
340: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace
341: @*/
342: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
343: {
348: if (sp->setupcalled) return(0);
349: sp->setupcalled = PETSC_TRUE;
350: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
351: if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
352: return(0);
353: }
355: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
356: {
357: PetscInt pStart = -1, pEnd = -1, depth = -1;
361: if (!dm) return(0);
362: DMPlexGetChart(dm, &pStart, &pEnd);
363: DMPlexGetDepth(dm, &depth);
365: if (sp->pointSpaces) {
366: PetscInt i;
368: for (i = 0; i < pEnd - pStart; i++) {
369: PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
370: }
371: }
372: PetscFree(sp->pointSpaces);
374: if (sp->heightSpaces) {
375: PetscInt i;
377: for (i = 0; i <= depth; i++) {
378: PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
379: }
380: }
381: PetscFree(sp->heightSpaces);
383: PetscSectionDestroy(&(sp->pointSection));
384: PetscQuadratureDestroy(&(sp->intNodes));
385: VecDestroy(&(sp->intDofValues));
386: VecDestroy(&(sp->intNodeValues));
387: MatDestroy(&(sp->intMat));
388: PetscQuadratureDestroy(&(sp->allNodes));
389: VecDestroy(&(sp->allDofValues));
390: VecDestroy(&(sp->allNodeValues));
391: MatDestroy(&(sp->allMat));
392: PetscFree(sp->numDof);
393: return(0);
394: }
397: /*@
398: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
400: Collective on sp
402: Input Parameter:
403: . sp - the PetscDualSpace object to destroy
405: Level: beginner
407: .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate()
408: @*/
409: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
410: {
411: PetscInt dim, f;
412: DM dm;
416: if (!*sp) return(0);
419: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
420: ((PetscObject) (*sp))->refct = 0;
422: PetscDualSpaceGetDimension(*sp, &dim);
423: dm = (*sp)->dm;
425: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
426: PetscDualSpaceClearDMData_Internal(*sp, dm);
428: for (f = 0; f < dim; ++f) {
429: PetscQuadratureDestroy(&(*sp)->functional[f]);
430: }
431: PetscFree((*sp)->functional);
432: DMDestroy(&(*sp)->dm);
433: PetscHeaderDestroy(sp);
434: return(0);
435: }
437: /*@
438: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
440: Collective
442: Input Parameter:
443: . comm - The communicator for the PetscDualSpace object
445: Output Parameter:
446: . sp - The PetscDualSpace object
448: Level: beginner
450: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
451: @*/
452: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
453: {
454: PetscDualSpace s;
459: PetscCitationsRegister(FECitation,&FEcite);
460: *sp = NULL;
461: PetscFEInitializePackage();
463: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
465: s->order = 0;
466: s->Nc = 1;
467: s->k = 0;
468: s->spdim = -1;
469: s->spintdim = -1;
470: s->uniform = PETSC_TRUE;
471: s->setupcalled = PETSC_FALSE;
473: *sp = s;
474: return(0);
475: }
477: /*@
478: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
480: Collective on sp
482: Input Parameter:
483: . sp - The original PetscDualSpace
485: Output Parameter:
486: . spNew - The duplicate PetscDualSpace
488: Level: beginner
490: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
491: @*/
492: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
493: {
494: DM dm;
495: PetscDualSpaceType type;
496: const char *name;
502: PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
503: PetscObjectGetName((PetscObject) sp, &name);
504: PetscObjectSetName((PetscObject) *spNew, name);
505: PetscDualSpaceGetType(sp, &type);
506: PetscDualSpaceSetType(*spNew, type);
507: PetscDualSpaceGetDM(sp, &dm);
508: PetscDualSpaceSetDM(*spNew, dm);
510: (*spNew)->order = sp->order;
511: (*spNew)->k = sp->k;
512: (*spNew)->Nc = sp->Nc;
513: (*spNew)->uniform = sp->uniform;
514: if (sp->ops->duplicate) {(*sp->ops->duplicate)(sp, *spNew);}
515: return(0);
516: }
518: /*@
519: PetscDualSpaceGetDM - Get the DM representing the reference cell
521: Not collective
523: Input Parameter:
524: . sp - The PetscDualSpace
526: Output Parameter:
527: . dm - The reference cell
529: Level: intermediate
531: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
532: @*/
533: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
534: {
538: *dm = sp->dm;
539: return(0);
540: }
542: /*@
543: PetscDualSpaceSetDM - Get the DM representing the reference cell
545: Not collective
547: Input Parameters:
548: + sp - The PetscDualSpace
549: - dm - The reference cell
551: Level: intermediate
553: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
554: @*/
555: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
556: {
562: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up");
563: PetscObjectReference((PetscObject) dm);
564: if (sp->dm && sp->dm != dm) {
565: PetscDualSpaceClearDMData_Internal(sp, sp->dm);
566: }
567: DMDestroy(&sp->dm);
568: sp->dm = dm;
569: return(0);
570: }
572: /*@
573: PetscDualSpaceGetOrder - Get the order of the dual space
575: Not collective
577: Input Parameter:
578: . sp - The PetscDualSpace
580: Output Parameter:
581: . order - The order
583: Level: intermediate
585: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
586: @*/
587: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
588: {
592: *order = sp->order;
593: return(0);
594: }
596: /*@
597: PetscDualSpaceSetOrder - Set the order of the dual space
599: Not collective
601: Input Parameters:
602: + sp - The PetscDualSpace
603: - order - The order
605: Level: intermediate
607: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
608: @*/
609: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
610: {
613: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up");
614: sp->order = order;
615: return(0);
616: }
618: /*@
619: PetscDualSpaceGetNumComponents - Return the number of components for this space
621: Input Parameter:
622: . sp - The PetscDualSpace
624: Output Parameter:
625: . Nc - The number of components
627: Note: A vector space, for example, will have d components, where d is the spatial dimension
629: Level: intermediate
631: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
632: @*/
633: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
634: {
638: *Nc = sp->Nc;
639: return(0);
640: }
642: /*@
643: PetscDualSpaceSetNumComponents - Set the number of components for this space
645: Input Parameters:
646: + sp - The PetscDualSpace
647: - order - The number of components
649: Level: intermediate
651: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
652: @*/
653: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
654: {
657: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
658: sp->Nc = Nc;
659: return(0);
660: }
662: /*@
663: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
665: Not collective
667: Input Parameters:
668: + sp - The PetscDualSpace
669: - i - The basis number
671: Output Parameter:
672: . functional - The basis functional
674: Level: intermediate
676: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
677: @*/
678: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
679: {
680: PetscInt dim;
686: PetscDualSpaceGetDimension(sp, &dim);
687: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
688: *functional = sp->functional[i];
689: return(0);
690: }
692: /*@
693: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
695: Not collective
697: Input Parameter:
698: . sp - The PetscDualSpace
700: Output Parameter:
701: . dim - The dimension
703: Level: intermediate
705: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
706: @*/
707: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
708: {
714: if (sp->spdim < 0) {
715: PetscSection section;
717: PetscDualSpaceGetSection(sp, §ion);
718: if (section) {
719: PetscSectionGetStorageSize(section, &(sp->spdim));
720: } else sp->spdim = 0;
721: }
722: *dim = sp->spdim;
723: return(0);
724: }
726: /*@
727: PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain
729: Not collective
731: Input Parameter:
732: . sp - The PetscDualSpace
734: Output Parameter:
735: . dim - The dimension
737: Level: intermediate
739: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
740: @*/
741: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
742: {
748: if (sp->spintdim < 0) {
749: PetscSection section;
751: PetscDualSpaceGetSection(sp, §ion);
752: if (section) {
753: PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
754: } else sp->spintdim = 0;
755: }
756: *intdim = sp->spintdim;
757: return(0);
758: }
760: /*@
761: PetscDualSpaceGetUniform - Whether this dual space is uniform
763: Not collective
765: Input Parameters:
766: . sp - A dual space
768: Output Parameters:
769: . uniform - PETSC_TRUE if (a) the dual space is the same for each point in a stratum of the reference DMPlex, and
770: (b) every symmetry of each point in the reference DMPlex is also a symmetry of the point's dual space.
773: Level: advanced
775: Note: all of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells
776: with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform.
778: .seealso: PetscDualSpaceGetPointSubspace(), PetscDualSpaceGetSymmetries()
779: @*/
780: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
781: {
785: *uniform = sp->uniform;
786: return(0);
787: }
790: /*@C
791: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
793: Not collective
795: Input Parameter:
796: . sp - The PetscDualSpace
798: Output Parameter:
799: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
801: Level: intermediate
803: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
804: @*/
805: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
806: {
812: if (!sp->uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform space does not have a fixed number of dofs for each height");
813: if (!sp->numDof) {
814: DM dm;
815: PetscInt depth, d;
816: PetscSection section;
818: PetscDualSpaceGetDM(sp, &dm);
819: DMPlexGetDepth(dm, &depth);
820: PetscCalloc1(depth+1,&(sp->numDof));
821: PetscDualSpaceGetSection(sp, §ion);
822: for (d = 0; d <= depth; d++) {
823: PetscInt dStart, dEnd;
825: DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
826: if (dEnd <= dStart) continue;
827: PetscSectionGetDof(section, dStart, &(sp->numDof[d]));
829: }
830: }
831: *numDof = sp->numDof;
832: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
833: return(0);
834: }
836: /* create the section of the right size and set a permutation for topological ordering */
837: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
838: {
839: DM dm;
840: PetscInt pStart, pEnd, cStart, cEnd, c, depth, count, i;
841: PetscInt *seen, *perm;
842: PetscSection section;
846: dm = sp->dm;
847: PetscSectionCreate(PETSC_COMM_SELF, §ion);
848: DMPlexGetChart(dm, &pStart, &pEnd);
849: PetscSectionSetChart(section, pStart, pEnd);
850: PetscCalloc1(pEnd - pStart, &seen);
851: PetscMalloc1(pEnd - pStart, &perm);
852: DMPlexGetDepth(dm, &depth);
853: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
854: for (c = cStart, count = 0; c < cEnd; c++) {
855: PetscInt closureSize = -1, e;
856: PetscInt *closure = NULL;
858: perm[count++] = c;
859: seen[c-pStart] = 1;
860: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
861: for (e = 0; e < closureSize; e++) {
862: PetscInt point = closure[2*e];
864: if (seen[point-pStart]) continue;
865: perm[count++] = point;
866: seen[point-pStart] = 1;
867: }
868: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
869: }
870: if (count != pEnd - pStart) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering");
871: for (i = 0; i < pEnd - pStart; i++) if (perm[i] != i) break;
872: if (i < pEnd - pStart) {
873: IS permIS;
875: ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
876: ISSetPermutation(permIS);
877: PetscSectionSetPermutation(section, permIS);
878: ISDestroy(&permIS);
879: } else {
880: PetscFree(perm);
881: }
882: PetscFree(seen);
883: *topSection = section;
884: return(0);
885: }
887: /* mark boundary points and set up */
888: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
889: {
890: DM dm;
891: DMLabel boundary;
892: PetscInt pStart, pEnd, p;
896: dm = sp->dm;
897: DMLabelCreate(PETSC_COMM_SELF,"boundary",&boundary);
898: PetscDualSpaceGetDM(sp,&dm);
899: DMPlexMarkBoundaryFaces(dm,1,boundary);
900: DMPlexLabelComplete(dm,boundary);
901: DMPlexGetChart(dm, &pStart, &pEnd);
902: for (p = pStart; p < pEnd; p++) {
903: PetscInt bval;
905: DMLabelGetValue(boundary, p, &bval);
906: if (bval == 1) {
907: PetscInt dof;
909: PetscSectionGetDof(section, p, &dof);
910: PetscSectionSetConstraintDof(section, p, dof);
911: }
912: }
913: DMLabelDestroy(&boundary);
914: PetscSectionSetUp(section);
915: return(0);
916: }
918: /*@
919: PetscDualSpaceGetSection - Create a PetscSection over the reference cell with the layout from this space
921: Collective on sp
923: Input Parameters:
924: . sp - The PetscDualSpace
926: Output Parameter:
927: . section - The section
929: Level: advanced
931: .seealso: PetscDualSpaceCreate(), DMPLEX
932: @*/
933: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
934: {
935: PetscInt pStart, pEnd, p;
939: if (!sp->pointSection) {
940: /* mark the boundary */
941: PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
942: DMPlexGetChart(sp->dm,&pStart,&pEnd);
943: for (p = pStart; p < pEnd; p++) {
944: PetscDualSpace psp;
946: PetscDualSpaceGetPointSubspace(sp, p, &psp);
947: if (psp) {
948: PetscInt dof;
950: PetscDualSpaceGetInteriorDimension(psp, &dof);
951: PetscSectionSetDof(sp->pointSection,p,dof);
952: }
953: }
954: PetscDualSpaceSectionSetUp_Internal(sp,sp->pointSection);
955: }
956: *section = sp->pointSection;
957: return(0);
958: }
960: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
961: * have one cell */
962: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
963: {
964: PetscReal *sv0, *v0, *J;
965: PetscSection section;
966: PetscInt dim, s, k;
967: DM dm;
971: PetscDualSpaceGetDM(sp, &dm);
972: DMGetDimension(dm, &dim);
973: PetscDualSpaceGetSection(sp, §ion);
974: PetscMalloc3(dim, &v0, dim, &sv0, dim*dim, &J);
975: PetscDualSpaceGetFormDegree(sp, &k);
976: for (s = sStart; s < sEnd; s++) {
977: PetscReal detJ, hdetJ;
978: PetscDualSpace ssp;
979: PetscInt dof, off, f, sdim;
980: PetscInt i, j;
981: DM sdm;
983: PetscDualSpaceGetPointSubspace(sp, s, &ssp);
984: if (!ssp) continue;
985: PetscSectionGetDof(section, s, &dof);
986: PetscSectionGetOffset(section, s, &off);
987: /* get the first vertex of the reference cell */
988: PetscDualSpaceGetDM(ssp, &sdm);
989: DMGetDimension(sdm, &sdim);
990: DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);
991: DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);
992: /* compactify Jacobian */
993: for (i = 0; i < dim; i++) for (j = 0; j < sdim; j++) J[i* sdim + j] = J[i * dim + j];
994: for (f = 0; f < dof; f++) {
995: PetscQuadrature fn;
997: PetscDualSpaceGetFunctional(ssp, f, &fn);
998: PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off+f]));
999: }
1000: }
1001: PetscFree3(v0, sv0, J);
1002: return(0);
1003: }
1005: /*@
1006: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1008: Collective on sp
1010: Input Parameters:
1011: + sp - The PetscDualSpace
1012: . dim - The spatial dimension
1013: - simplex - Flag for simplex, otherwise use a tensor-product cell
1015: Output Parameter:
1016: . refdm - The reference cell
1018: Level: intermediate
1020: .seealso: PetscDualSpaceCreate(), DMPLEX
1021: @*/
1022: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1023: {
1027: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1028: return(0);
1029: }
1031: /*@C
1032: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1034: Input Parameters:
1035: + sp - The PetscDualSpace object
1036: . f - The basis functional index
1037: . time - The time
1038: . 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)
1039: . numComp - The number of components for the function
1040: . func - The input function
1041: - ctx - A context for the function
1043: Output Parameter:
1044: . value - numComp output values
1046: Note: The calling sequence for the callback func is given by:
1048: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1049: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1051: Level: beginner
1053: .seealso: PetscDualSpaceCreate()
1054: @*/
1055: 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)
1056: {
1063: (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
1064: return(0);
1065: }
1067: /*@C
1068: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
1070: Input Parameters:
1071: + sp - The PetscDualSpace object
1072: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
1074: Output Parameter:
1075: . spValue - The values of all dual space functionals
1077: Level: beginner
1079: .seealso: PetscDualSpaceCreate()
1080: @*/
1081: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1082: {
1087: (*sp->ops->applyall)(sp, pointEval, spValue);
1088: return(0);
1089: }
1091: /*@C
1092: PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
1094: Input Parameters:
1095: + sp - The PetscDualSpace object
1096: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1098: Output Parameter:
1099: . spValue - The values of interior dual space functionals
1101: Level: beginner
1103: .seealso: PetscDualSpaceCreate()
1104: @*/
1105: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1106: {
1111: (*sp->ops->applyint)(sp, pointEval, spValue);
1112: return(0);
1113: }
1115: /*@C
1116: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
1118: Input Parameters:
1119: + sp - The PetscDualSpace object
1120: . f - The basis functional index
1121: . time - The time
1122: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1123: . Nc - The number of components for the function
1124: . func - The input function
1125: - ctx - A context for the function
1127: Output Parameter:
1128: . value - The output value
1130: Note: The calling sequence for the callback func is given by:
1132: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1133: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1135: and the idea is to evaluate the functional as an integral
1137: $ n(f) = int dx n(x) . f(x)
1139: where both n and f have Nc components.
1141: Level: beginner
1143: .seealso: PetscDualSpaceCreate()
1144: @*/
1145: 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)
1146: {
1147: DM dm;
1148: PetscQuadrature n;
1149: const PetscReal *points, *weights;
1150: PetscReal x[3];
1151: PetscScalar *val;
1152: PetscInt dim, dE, qNc, c, Nq, q;
1153: PetscBool isAffine;
1154: PetscErrorCode ierr;
1159: PetscDualSpaceGetDM(sp, &dm);
1160: PetscDualSpaceGetFunctional(sp, f, &n);
1161: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
1162: if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
1163: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1164: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1165: *value = 0.0;
1166: isAffine = cgeom->isAffine;
1167: dE = cgeom->dimEmbed;
1168: for (q = 0; q < Nq; ++q) {
1169: if (isAffine) {
1170: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
1171: (*func)(dE, time, x, Nc, val, ctx);
1172: } else {
1173: (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
1174: }
1175: for (c = 0; c < Nc; ++c) {
1176: *value += val[c]*weights[q*Nc+c];
1177: }
1178: }
1179: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1180: return(0);
1181: }
1183: /*@C
1184: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
1186: Input Parameters:
1187: + sp - The PetscDualSpace object
1188: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
1190: Output Parameter:
1191: . spValue - The values of all dual space functionals
1193: Level: beginner
1195: .seealso: PetscDualSpaceCreate()
1196: @*/
1197: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1198: {
1199: Vec pointValues, dofValues;
1200: Mat allMat;
1201: PetscErrorCode ierr;
1207: PetscDualSpaceGetAllData(sp, NULL, &allMat);
1208: if (!(sp->allNodeValues)) {
1209: MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1210: }
1211: pointValues = sp->allNodeValues;
1212: if (!(sp->allDofValues)) {
1213: MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1214: }
1215: dofValues = sp->allDofValues;
1216: VecPlaceArray(pointValues, pointEval);
1217: VecPlaceArray(dofValues, spValue);
1218: MatMult(allMat, pointValues, dofValues);
1219: VecResetArray(dofValues);
1220: VecResetArray(pointValues);
1221: return(0);
1222: }
1224: /*@C
1225: PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
1227: Input Parameters:
1228: + sp - The PetscDualSpace object
1229: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1231: Output Parameter:
1232: . spValue - The values of interior dual space functionals
1234: Level: beginner
1236: .seealso: PetscDualSpaceCreate()
1237: @*/
1238: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1239: {
1240: Vec pointValues, dofValues;
1241: Mat intMat;
1242: PetscErrorCode ierr;
1248: PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1249: if (!(sp->intNodeValues)) {
1250: MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1251: }
1252: pointValues = sp->intNodeValues;
1253: if (!(sp->intDofValues)) {
1254: MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1255: }
1256: dofValues = sp->intDofValues;
1257: VecPlaceArray(pointValues, pointEval);
1258: VecPlaceArray(dofValues, spValue);
1259: MatMult(intMat, pointValues, dofValues);
1260: VecResetArray(dofValues);
1261: VecResetArray(pointValues);
1262: return(0);
1263: }
1265: /*@
1266: PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values
1268: Input Parameter:
1269: . sp - The dualspace
1271: Output Parameter:
1272: + allNodes - A PetscQuadrature object containing all evaluation nodes
1273: - allMat - A Mat for the node-to-dof transformation
1275: Level: advanced
1277: .seealso: PetscDualSpaceCreate()
1278: @*/
1279: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1280: {
1287: if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1288: PetscQuadrature qpoints;
1289: Mat amat;
1291: (*sp->ops->createalldata)(sp,&qpoints,&amat);
1292: PetscQuadratureDestroy(&(sp->allNodes));
1293: MatDestroy(&(sp->allMat));
1294: sp->allNodes = qpoints;
1295: sp->allMat = amat;
1296: }
1297: if (allNodes) *allNodes = sp->allNodes;
1298: if (allMat) *allMat = sp->allMat;
1299: return(0);
1300: }
1302: /*@
1303: PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals
1305: Input Parameter:
1306: . sp - The dualspace
1308: Output Parameter:
1309: + allNodes - A PetscQuadrature object containing all evaluation nodes
1310: - allMat - A Mat for the node-to-dof transformation
1312: Level: advanced
1314: .seealso: PetscDualSpaceCreate()
1315: @*/
1316: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1317: {
1318: PetscInt spdim;
1319: PetscInt numPoints, offset;
1320: PetscReal *points;
1321: PetscInt f, dim;
1322: PetscInt Nc, nrows, ncols;
1323: PetscInt maxNumPoints;
1324: PetscQuadrature q;
1325: Mat A;
1326: PetscErrorCode ierr;
1329: PetscDualSpaceGetNumComponents(sp, &Nc);
1330: PetscDualSpaceGetDimension(sp,&spdim);
1331: if (!spdim) {
1332: PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1333: PetscQuadratureSetData(*allNodes,0,0,0,NULL,NULL);
1334: }
1335: nrows = spdim;
1336: PetscDualSpaceGetFunctional(sp,0,&q);
1337: PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
1338: maxNumPoints = numPoints;
1339: for (f = 1; f < spdim; f++) {
1340: PetscInt Np;
1342: PetscDualSpaceGetFunctional(sp,f,&q);
1343: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1344: numPoints += Np;
1345: maxNumPoints = PetscMax(maxNumPoints,Np);
1346: }
1347: ncols = numPoints * Nc;
1348: PetscMalloc1(dim*numPoints,&points);
1349: MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1350: for (f = 0, offset = 0; f < spdim; f++) {
1351: const PetscReal *p, *w;
1352: PetscInt Np, i;
1353: PetscInt fnc;
1355: PetscDualSpaceGetFunctional(sp,f,&q);
1356: PetscQuadratureGetData(q,NULL,&fnc,&Np,&p,&w);
1357: if (fnc != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch");
1358: for (i = 0; i < Np * dim; i++) {
1359: points[offset* dim + i] = p[i];
1360: }
1361: for (i = 0; i < Np * Nc; i++) {
1362: MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1363: }
1364: offset += Np;
1365: }
1366: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1367: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1368: PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1369: PetscQuadratureSetData(*allNodes,dim,0,numPoints,points,NULL);
1370: *allMat = A;
1371: return(0);
1372: }
1374: /*@
1375: PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1376: this space, as well as the matrix that computes the degrees of freedom from the quadrature values. Degrees of
1377: freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the
1378: reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by
1379: PetscDualSpaceGetSection()).
1381: Input Parameter:
1382: . sp - The dualspace
1384: Output Parameter:
1385: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1386: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1387: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1388: npoints is the number of points in intNodes and nc is PetscDualSpaceGetNumComponents().
1390: Level: advanced
1392: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetDimension(), PetscDualSpaceGetNumComponents(), PetscQuadratureGetData()
1393: @*/
1394: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1395: {
1402: if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1403: PetscQuadrature qpoints;
1404: Mat imat;
1406: (*sp->ops->createintdata)(sp,&qpoints,&imat);
1407: PetscQuadratureDestroy(&(sp->intNodes));
1408: MatDestroy(&(sp->intMat));
1409: sp->intNodes = qpoints;
1410: sp->intMat = imat;
1411: }
1412: if (intNodes) *intNodes = sp->intNodes;
1413: if (intMat) *intMat = sp->intMat;
1414: return(0);
1415: }
1417: /*@
1418: PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values
1420: Input Parameter:
1421: . sp - The dualspace
1423: Output Parameter:
1424: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1425: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1426: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1427: npoints is the number of points in allNodes and nc is PetscDualSpaceGetNumComponents().
1429: Level: advanced
1431: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetInteriorData()
1432: @*/
1433: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1434: {
1435: DM dm;
1436: PetscInt spdim0;
1437: PetscInt Nc;
1438: PetscInt pStart, pEnd, p, f;
1439: PetscSection section;
1440: PetscInt numPoints, offset, matoffset;
1441: PetscReal *points;
1442: PetscInt dim;
1443: PetscInt *nnz;
1444: PetscQuadrature q;
1445: Mat imat;
1446: PetscErrorCode ierr;
1450: PetscDualSpaceGetSection(sp, §ion);
1451: PetscSectionGetConstrainedStorageSize(section, &spdim0);
1452: if (!spdim0) {
1453: *intNodes = NULL;
1454: *intMat = NULL;
1455: return(0);
1456: }
1457: PetscDualSpaceGetNumComponents(sp, &Nc);
1458: PetscSectionGetChart(section, &pStart, &pEnd);
1459: PetscDualSpaceGetDM(sp, &dm);
1460: DMGetDimension(dm, &dim);
1461: PetscMalloc1(spdim0, &nnz);
1462: for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1463: PetscInt dof, cdof, off, d;
1465: PetscSectionGetDof(section, p, &dof);
1466: PetscSectionGetConstraintDof(section, p, &cdof);
1467: if (!(dof - cdof)) continue;
1468: PetscSectionGetOffset(section, p, &off);
1469: for (d = 0; d < dof; d++, off++, f++) {
1470: PetscInt Np;
1472: PetscDualSpaceGetFunctional(sp,off,&q);
1473: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1474: nnz[f] = Np * Nc;
1475: numPoints += Np;
1476: }
1477: }
1478: MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1479: PetscFree(nnz);
1480: PetscMalloc1(dim*numPoints,&points);
1481: for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1482: PetscInt dof, cdof, off, d;
1484: PetscSectionGetDof(section, p, &dof);
1485: PetscSectionGetConstraintDof(section, p, &cdof);
1486: if (!(dof - cdof)) continue;
1487: PetscSectionGetOffset(section, p, &off);
1488: for (d = 0; d < dof; d++, off++, f++) {
1489: const PetscReal *p;
1490: const PetscReal *w;
1491: PetscInt Np, i;
1493: PetscDualSpaceGetFunctional(sp,off,&q);
1494: PetscQuadratureGetData(q,NULL,NULL,&Np,&p,&w);
1495: for (i = 0; i < Np * dim; i++) {
1496: points[offset + i] = p[i];
1497: }
1498: for (i = 0; i < Np * Nc; i++) {
1499: MatSetValue(imat, f, matoffset + i, w[i],INSERT_VALUES);
1500: }
1501: offset += Np * dim;
1502: matoffset += Np * Nc;
1503: }
1504: }
1505: PetscQuadratureCreate(PETSC_COMM_SELF,intNodes);
1506: PetscQuadratureSetData(*intNodes,dim,0,numPoints,points,NULL);
1507: MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1508: MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1509: *intMat = imat;
1510: return(0);
1511: }
1513: /*@C
1514: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
1516: Input Parameters:
1517: + sp - The PetscDualSpace object
1518: . f - The basis functional index
1519: . time - The time
1520: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1521: . Nc - The number of components for the function
1522: . func - The input function
1523: - ctx - A context for the function
1525: Output Parameter:
1526: . value - The output value (scalar)
1528: Note: The calling sequence for the callback func is given by:
1530: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1531: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1533: and the idea is to evaluate the functional as an integral
1535: $ n(f) = int dx n(x) . f(x)
1537: where both n and f have Nc components.
1539: Level: beginner
1541: .seealso: PetscDualSpaceCreate()
1542: @*/
1543: 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)
1544: {
1545: DM dm;
1546: PetscQuadrature n;
1547: const PetscReal *points, *weights;
1548: PetscScalar *val;
1549: PetscInt dimEmbed, qNc, c, Nq, q;
1550: PetscErrorCode ierr;
1555: PetscDualSpaceGetDM(sp, &dm);
1556: DMGetCoordinateDim(dm, &dimEmbed);
1557: PetscDualSpaceGetFunctional(sp, f, &n);
1558: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1559: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1560: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1561: *value = 0.;
1562: for (q = 0; q < Nq; ++q) {
1563: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1564: for (c = 0; c < Nc; ++c) {
1565: *value += val[c]*weights[q*Nc+c];
1566: }
1567: }
1568: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1569: return(0);
1570: }
1572: /*@
1573: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
1574: given height. This assumes that the reference cell is symmetric over points of this height.
1576: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1577: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1578: support extracting subspaces, then NULL is returned.
1580: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1582: Not collective
1584: Input Parameters:
1585: + sp - the PetscDualSpace object
1586: - height - the height of the mesh point for which the subspace is desired
1588: Output Parameter:
1589: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1590: point, which will be of lesser dimension if height > 0.
1592: Level: advanced
1594: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1595: @*/
1596: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1597: {
1598: PetscInt depth = -1, cStart, cEnd;
1599: DM dm;
1605: if (!(sp->uniform)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform dual space does not have a single dual space at each height");
1606: *subsp = NULL;
1607: dm = sp->dm;
1608: DMPlexGetDepth(dm, &depth);
1609: if (height < 0 || height > depth) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height");
1610: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1611: if (height == 0 && cEnd == cStart + 1) {
1612: *subsp = sp;
1613: return(0);
1614: }
1615: if (!sp->heightSpaces) {
1616: PetscInt h;
1617: PetscCalloc1(depth+1, &(sp->heightSpaces));
1619: for (h = 0; h <= depth; h++) {
1620: if (h == 0 && cEnd == cStart + 1) continue;
1621: if (sp->ops->createheightsubspace) {(*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));}
1622: else if (sp->pointSpaces) {
1623: PetscInt hStart, hEnd;
1625: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
1626: if (hEnd > hStart) {
1627: PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1628: sp->heightSpaces[h] = sp->pointSpaces[hStart];
1629: }
1630: }
1631: }
1632: }
1633: *subsp = sp->heightSpaces[height];
1634: return(0);
1635: }
1637: /*@
1638: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
1640: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
1641: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
1642: subspaces, then NULL is returned.
1644: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1646: Not collective
1648: Input Parameters:
1649: + sp - the PetscDualSpace object
1650: - point - the point (in the dual space's DM) for which the subspace is desired
1652: Output Parameters:
1653: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1654: point, which will be of lesser dimension if height > 0.
1656: Level: advanced
1658: .seealso: PetscDualSpace
1659: @*/
1660: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1661: {
1662: PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1663: DM dm;
1669: *bdsp = NULL;
1670: dm = sp->dm;
1671: DMPlexGetChart(dm, &pStart, &pEnd);
1672: if (point < pStart || point > pEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point");
1673: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1674: 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 */
1675: *bdsp = sp;
1676: return(0);
1677: }
1678: if (!sp->pointSpaces) {
1679: PetscInt p;
1680: PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));
1682: for (p = 0; p < pEnd - pStart; p++) {
1683: if (p + pStart == cStart && cEnd == cStart + 1) continue;
1684: if (sp->ops->createpointsubspace) {(*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));}
1685: else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1686: PetscInt dim, depth, height;
1687: DMLabel label;
1689: DMPlexGetDepth(dm,&dim);
1690: DMPlexGetDepthLabel(dm,&label);
1691: DMLabelGetValue(label,p+pStart,&depth);
1692: height = dim - depth;
1693: PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1694: PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1695: }
1696: }
1697: }
1698: *bdsp = sp->pointSpaces[point - pStart];
1699: return(0);
1700: }
1702: /*@C
1703: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
1705: Not collective
1707: Input Parameter:
1708: . sp - the PetscDualSpace object
1710: Output Parameters:
1711: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1712: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation
1714: Note: The permutation and flip arrays are organized in the following way
1715: $ perms[p][ornt][dof # on point] = new local dof #
1716: $ flips[p][ornt][dof # on point] = reversal or not
1718: Level: developer
1720: @*/
1721: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1722: {
1729: if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1730: return(0);
1731: }
1733: /*@
1734: PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1735: dual space's functionals.
1737: Input Parameter:
1738: . dsp - The PetscDualSpace
1740: Output Parameter:
1741: . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1742: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1743: 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).
1744: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1745: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1746: but are stored as 1-forms.
1748: Level: developer
1750: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1751: @*/
1752: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1753: {
1757: *k = dsp->k;
1758: return(0);
1759: }
1761: /*@
1762: PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1763: dual space's functionals.
1765: Input Parameter:
1766: + dsp - The PetscDualSpace
1767: - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1768: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1769: 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).
1770: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1771: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1772: but are stored as 1-forms.
1774: Level: developer
1776: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1777: @*/
1778: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1779: {
1780: PetscInt dim;
1784: if (dsp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
1785: dim = dsp->dm->dim;
1786: if (k < -dim || k > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %D-form on %D-dimensional reference cell", PetscAbsInt(k), dim);
1787: dsp->k = k;
1788: return(0);
1789: }
1791: /*@
1792: PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space
1794: Input Parameter:
1795: . dsp - The PetscDualSpace
1797: Output Parameter:
1798: . k - The simplex dimension
1800: Level: developer
1802: Note: Currently supported values are
1803: $ 0: These are H_1 methods that only transform coordinates
1804: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1805: $ 2: These are the same as 1
1806: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1808: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1809: @*/
1810: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1811: {
1812: PetscInt dim;
1817: dim = dsp->dm->dim;
1818: if (!dsp->k) *k = IDENTITY_TRANSFORM;
1819: else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1820: else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1821: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1822: return(0);
1823: }
1825: /*@C
1826: PetscDualSpaceTransform - Transform the function values
1828: Input Parameters:
1829: + dsp - The PetscDualSpace
1830: . trans - The type of transform
1831: . isInverse - Flag to invert the transform
1832: . fegeom - The cell geometry
1833: . Nv - The number of function samples
1834: . Nc - The number of function components
1835: - vals - The function values
1837: Output Parameter:
1838: . vals - The transformed function values
1840: Level: intermediate
1842: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1844: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1845: @*/
1846: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1847: {
1848: PetscReal Jstar[9] = {0};
1849: PetscInt dim, v, c, Nk;
1856: /* TODO: not handling dimEmbed != dim right now */
1857: dim = dsp->dm->dim;
1858: /* No change needed for 0-forms */
1859: if (!dsp->k) return(0);
1860: PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1861: /* TODO: use fegeom->isAffine */
1862: PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1863: for (v = 0; v < Nv; ++v) {
1864: switch (Nk) {
1865: case 1:
1866: for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0];
1867: break;
1868: case 2:
1869: for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1870: break;
1871: case 3:
1872: for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1873: break;
1874: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk);
1875: }
1876: }
1877: return(0);
1878: }
1880: /*@C
1881: PetscDualSpaceTransformGradient - Transform the function gradient values
1883: Input Parameters:
1884: + dsp - The PetscDualSpace
1885: . trans - The type of transform
1886: . isInverse - Flag to invert the transform
1887: . fegeom - The cell geometry
1888: . Nv - The number of function gradient samples
1889: . Nc - The number of function components
1890: - vals - The function gradient values
1892: Output Parameter:
1893: . vals - The transformed function values
1895: Level: intermediate
1897: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1899: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1900: @*/
1901: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1902: {
1903: PetscInt dim, v, c, d;
1909: dim = dsp->dm->dim;
1910: /* Transform gradient */
1911: for (v = 0; v < Nv; ++v) {
1912: for (c = 0; c < Nc; ++c) {
1913: switch (dim)
1914: {
1915: case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];break;
1916: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1917: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1918: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1919: }
1920: }
1921: }
1922: /* Assume its a vector, otherwise assume its a bunch of scalars */
1923: if (Nc == 1 || Nc != dim) return(0);
1924: switch (trans) {
1925: case IDENTITY_TRANSFORM: break;
1926: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1927: if (isInverse) {
1928: for (v = 0; v < Nv; ++v) {
1929: for (d = 0; d < dim; ++d) {
1930: switch (dim)
1931: {
1932: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1933: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1934: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1935: }
1936: }
1937: }
1938: } else {
1939: for (v = 0; v < Nv; ++v) {
1940: for (d = 0; d < dim; ++d) {
1941: switch (dim)
1942: {
1943: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1944: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1945: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1946: }
1947: }
1948: }
1949: }
1950: break;
1951: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1952: if (isInverse) {
1953: for (v = 0; v < Nv; ++v) {
1954: for (d = 0; d < dim; ++d) {
1955: switch (dim)
1956: {
1957: case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1958: case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1959: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1960: }
1961: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1962: }
1963: }
1964: } else {
1965: for (v = 0; v < Nv; ++v) {
1966: for (d = 0; d < dim; ++d) {
1967: switch (dim)
1968: {
1969: case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1970: case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1971: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1972: }
1973: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1974: }
1975: }
1976: }
1977: break;
1978: }
1979: return(0);
1980: }
1982: /*@C
1983: 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.
1985: Input Parameters:
1986: + dsp - The PetscDualSpace
1987: . fegeom - The geometry for this cell
1988: . Nq - The number of function samples
1989: . Nc - The number of function components
1990: - pointEval - The function values
1992: Output Parameter:
1993: . pointEval - The transformed function values
1995: Level: advanced
1997: 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.
1999: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2001: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2002: @*/
2003: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2004: {
2005: PetscDualSpaceTransformType trans;
2006: PetscInt k;
2007: PetscErrorCode ierr;
2013: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2014: This determines their transformation properties. */
2015: PetscDualSpaceGetDeRahm(dsp, &k);
2016: switch (k)
2017: {
2018: case 0: /* H^1 point evaluations */
2019: trans = IDENTITY_TRANSFORM;break;
2020: case 1: /* Hcurl preserves tangential edge traces */
2021: trans = COVARIANT_PIOLA_TRANSFORM;break;
2022: case 2:
2023: case 3: /* Hdiv preserve normal traces */
2024: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2025: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2026: }
2027: PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2028: return(0);
2029: }
2031: /*@C
2032: 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.
2034: Input Parameters:
2035: + dsp - The PetscDualSpace
2036: . fegeom - The geometry for this cell
2037: . Nq - The number of function samples
2038: . Nc - The number of function components
2039: - pointEval - The function values
2041: Output Parameter:
2042: . pointEval - The transformed function values
2044: Level: advanced
2046: 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.
2048: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2050: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2051: @*/
2052: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2053: {
2054: PetscDualSpaceTransformType trans;
2055: PetscInt k;
2056: PetscErrorCode ierr;
2062: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2063: This determines their transformation properties. */
2064: PetscDualSpaceGetDeRahm(dsp, &k);
2065: switch (k)
2066: {
2067: case 0: /* H^1 point evaluations */
2068: trans = IDENTITY_TRANSFORM;break;
2069: case 1: /* Hcurl preserves tangential edge traces */
2070: trans = COVARIANT_PIOLA_TRANSFORM;break;
2071: case 2:
2072: case 3: /* Hdiv preserve normal traces */
2073: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2074: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2075: }
2076: PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2077: return(0);
2078: }
2080: /*@C
2081: 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.
2083: Input Parameters:
2084: + dsp - The PetscDualSpace
2085: . fegeom - The geometry for this cell
2086: . Nq - The number of function gradient samples
2087: . Nc - The number of function components
2088: - pointEval - The function gradient values
2090: Output Parameter:
2091: . pointEval - The transformed function gradient values
2093: Level: advanced
2095: 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.
2097: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2099: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2100: @*/
2101: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2102: {
2103: PetscDualSpaceTransformType trans;
2104: PetscInt k;
2105: PetscErrorCode ierr;
2111: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2112: This determines their transformation properties. */
2113: PetscDualSpaceGetDeRahm(dsp, &k);
2114: switch (k)
2115: {
2116: case 0: /* H^1 point evaluations */
2117: trans = IDENTITY_TRANSFORM;break;
2118: case 1: /* Hcurl preserves tangential edge traces */
2119: trans = COVARIANT_PIOLA_TRANSFORM;break;
2120: case 2:
2121: case 3: /* Hdiv preserve normal traces */
2122: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2123: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2124: }
2125: PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2126: return(0);
2127: }