Actual source code: dualspace.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: PetscClassId PETSCDUALSPACE_CLASSID = 0;
6: PetscLogEvent PETSCDUALSPACE_SetUp;
8: PetscFunctionList PetscDualSpaceList = NULL;
9: PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE;
11: const char *const PetscDualSpaceReferenceCells[] = {"SIMPLEX", "TENSOR", "PetscDualSpaceReferenceCell", "PETSCDUALSPACE_REFCELL_", NULL};
13: /*
14: PetscDualSpaceLatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
15: Ordering is lexicographic with lowest index as least significant in ordering.
16: 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}.
18: Input Parameters:
19: + len - The length of the tuple
20: . max - The maximum sum
21: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
23: Output Parameter:
24: . tup - A tuple of len integers whos sum is at most 'max'
26: Level: developer
28: .seealso: PetscDualSpaceTensorPointLexicographic_Internal()
29: */
30: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
31: {
33: while (len--) {
34: max -= tup[len];
35: if (!max) {
36: tup[len] = 0;
37: break;
38: }
39: }
40: tup[++len]++;
41: return(0);
42: }
44: /*
45: PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
46: Ordering is lexicographic with lowest index as least significant in ordering.
47: 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}.
49: Input Parameters:
50: + len - The length of the tuple
51: . max - The maximum value
52: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
54: Output Parameter:
55: . tup - A tuple of len integers whos sum is at most 'max'
57: Level: developer
59: .seealso: PetscDualSpaceLatticePointLexicographic_Internal()
60: */
61: PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
62: {
63: PetscInt i;
66: for (i = 0; i < len; i++) {
67: if (tup[i] < max) {
68: break;
69: } else {
70: tup[i] = 0;
71: }
72: }
73: tup[i]++;
74: return(0);
75: }
77: /*@C
78: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
80: Not Collective
82: Input Parameters:
83: + name - The name of a new user-defined creation routine
84: - create_func - The creation routine itself
86: Notes:
87: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
89: Sample usage:
90: .vb
91: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
92: .ve
94: Then, your PetscDualSpace type can be chosen with the procedural interface via
95: .vb
96: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
97: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
98: .ve
99: or at runtime via the option
100: .vb
101: -petscdualspace_type my_dual_space
102: .ve
104: Level: advanced
106: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()
108: @*/
109: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
110: {
114: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
115: return(0);
116: }
118: /*@C
119: PetscDualSpaceSetType - Builds a particular PetscDualSpace
121: Collective on sp
123: Input Parameters:
124: + sp - The PetscDualSpace object
125: - name - The kind of space
127: Options Database Key:
128: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
130: Level: intermediate
132: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
133: @*/
134: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
135: {
136: PetscErrorCode (*r)(PetscDualSpace);
137: PetscBool match;
142: PetscObjectTypeCompare((PetscObject) sp, name, &match);
143: if (match) return(0);
145: if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
146: PetscFunctionListFind(PetscDualSpaceList, name, &r);
147: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
149: if (sp->ops->destroy) {
150: (*sp->ops->destroy)(sp);
151: sp->ops->destroy = NULL;
152: }
153: (*r)(sp);
154: PetscObjectChangeTypeName((PetscObject) sp, name);
155: return(0);
156: }
158: /*@C
159: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
161: Not Collective
163: Input Parameter:
164: . sp - The PetscDualSpace
166: Output Parameter:
167: . name - The PetscDualSpace type name
169: Level: intermediate
171: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
172: @*/
173: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
174: {
180: if (!PetscDualSpaceRegisterAllCalled) {
181: PetscDualSpaceRegisterAll();
182: }
183: *name = ((PetscObject) sp)->type_name;
184: return(0);
185: }
187: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
188: {
189: PetscViewerFormat format;
190: PetscInt pdim, f;
191: PetscErrorCode ierr;
194: PetscDualSpaceGetDimension(sp, &pdim);
195: PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);
196: PetscViewerASCIIPushTab(v);
197: if (sp->k) {
198: 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);
199: } else {
200: PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
201: }
202: if (sp->ops->view) {(*sp->ops->view)(sp, v);}
203: PetscViewerGetFormat(v, &format);
204: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
205: PetscViewerASCIIPushTab(v);
206: for (f = 0; f < pdim; ++f) {
207: PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
208: PetscViewerASCIIPushTab(v);
209: PetscQuadratureView(sp->functional[f], v);
210: PetscViewerASCIIPopTab(v);
211: }
212: PetscViewerASCIIPopTab(v);
213: }
214: PetscViewerASCIIPopTab(v);
215: return(0);
216: }
218: /*@C
219: PetscDualSpaceViewFromOptions - View from Options
221: Collective on PetscDualSpace
223: Input Parameters:
224: + A - the PetscDualSpace object
225: . obj - Optional object, proivides prefix
226: - name - command line option
228: Level: intermediate
229: .seealso: PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate()
230: @*/
231: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[])
232: {
237: PetscObjectViewFromOptions((PetscObject)A,obj,name);
238: return(0);
239: }
241: /*@
242: PetscDualSpaceView - Views a PetscDualSpace
244: Collective on sp
246: Input Parameters:
247: + sp - the PetscDualSpace object to view
248: - v - the viewer
250: Level: beginner
252: .seealso PetscDualSpaceDestroy(), PetscDualSpace
253: @*/
254: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
255: {
256: PetscBool iascii;
262: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
263: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
264: if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
265: return(0);
266: }
268: /*@
269: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
271: Collective on sp
273: Input Parameter:
274: . sp - the PetscDualSpace object to set options for
276: Options Database:
277: + -petscdualspace_order <order> - the approximation order of the space
278: . -petscdualspace_form_degree <deg> - the form degree, say 0 for point evaluations, or 2 for area integrals
279: . -petscdualspace_components <c> - the number of components, say d for a vector field
280: . -petscdualspace_refdim <d> - The spatial dimension of the reference cell
281: - -petscdualspace_refcell <celltype> - Reference cell type name
283: Level: intermediate
285: .seealso PetscDualSpaceView(), PetscDualSpace, PetscObjectSetFromOptions()
286: @*/
287: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
288: {
289: PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
290: PetscInt refDim = 0;
291: PetscBool flg;
292: const char *defaultType;
293: char name[256];
294: PetscErrorCode ierr;
298: if (!((PetscObject) sp)->type_name) {
299: defaultType = PETSCDUALSPACELAGRANGE;
300: } else {
301: defaultType = ((PetscObject) sp)->type_name;
302: }
303: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
305: PetscObjectOptionsBegin((PetscObject) sp);
306: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
307: if (flg) {
308: PetscDualSpaceSetType(sp, name);
309: } else if (!((PetscObject) sp)->type_name) {
310: PetscDualSpaceSetType(sp, defaultType);
311: }
312: PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
313: PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);
314: PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
315: if (sp->ops->setfromoptions) {
316: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
317: }
318: PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
319: PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
320: if (flg) {
321: DM K;
323: if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
324: PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
325: PetscDualSpaceSetDM(sp, K);
326: DMDestroy(&K);
327: }
329: /* process any options handlers added with PetscObjectAddOptionsHandler() */
330: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
331: PetscOptionsEnd();
332: sp->setfromoptionscalled = PETSC_TRUE;
333: return(0);
334: }
336: /*@
337: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
339: Collective on sp
341: Input Parameter:
342: . sp - the PetscDualSpace object to setup
344: Level: intermediate
346: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace
347: @*/
348: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
349: {
354: if (sp->setupcalled) return(0);
355: PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
356: sp->setupcalled = PETSC_TRUE;
357: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
358: PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
359: if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
360: return(0);
361: }
363: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
364: {
365: PetscInt pStart = -1, pEnd = -1, depth = -1;
369: if (!dm) return(0);
370: DMPlexGetChart(dm, &pStart, &pEnd);
371: DMPlexGetDepth(dm, &depth);
373: if (sp->pointSpaces) {
374: PetscInt i;
376: for (i = 0; i < pEnd - pStart; i++) {
377: PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
378: }
379: }
380: PetscFree(sp->pointSpaces);
382: if (sp->heightSpaces) {
383: PetscInt i;
385: for (i = 0; i <= depth; i++) {
386: PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
387: }
388: }
389: PetscFree(sp->heightSpaces);
391: PetscSectionDestroy(&(sp->pointSection));
392: PetscQuadratureDestroy(&(sp->intNodes));
393: VecDestroy(&(sp->intDofValues));
394: VecDestroy(&(sp->intNodeValues));
395: MatDestroy(&(sp->intMat));
396: PetscQuadratureDestroy(&(sp->allNodes));
397: VecDestroy(&(sp->allDofValues));
398: VecDestroy(&(sp->allNodeValues));
399: MatDestroy(&(sp->allMat));
400: PetscFree(sp->numDof);
401: return(0);
402: }
404: /*@
405: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
407: Collective on sp
409: Input Parameter:
410: . sp - the PetscDualSpace object to destroy
412: Level: beginner
414: .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate()
415: @*/
416: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
417: {
418: PetscInt dim, f;
419: DM dm;
423: if (!*sp) return(0);
426: if (--((PetscObject)(*sp))->refct > 0) {*sp = NULL; return(0);}
427: ((PetscObject) (*sp))->refct = 0;
429: PetscDualSpaceGetDimension(*sp, &dim);
430: dm = (*sp)->dm;
432: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
433: PetscDualSpaceClearDMData_Internal(*sp, dm);
435: for (f = 0; f < dim; ++f) {
436: PetscQuadratureDestroy(&(*sp)->functional[f]);
437: }
438: PetscFree((*sp)->functional);
439: DMDestroy(&(*sp)->dm);
440: PetscHeaderDestroy(sp);
441: return(0);
442: }
444: /*@
445: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
447: Collective
449: Input Parameter:
450: . comm - The communicator for the PetscDualSpace object
452: Output Parameter:
453: . sp - The PetscDualSpace object
455: Level: beginner
457: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
458: @*/
459: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
460: {
461: PetscDualSpace s;
466: PetscCitationsRegister(FECitation,&FEcite);
467: *sp = NULL;
468: PetscFEInitializePackage();
470: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
472: s->order = 0;
473: s->Nc = 1;
474: s->k = 0;
475: s->spdim = -1;
476: s->spintdim = -1;
477: s->uniform = PETSC_TRUE;
478: s->setupcalled = PETSC_FALSE;
480: *sp = s;
481: return(0);
482: }
484: /*@
485: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
487: Collective on sp
489: Input Parameter:
490: . sp - The original PetscDualSpace
492: Output Parameter:
493: . spNew - The duplicate PetscDualSpace
495: Level: beginner
497: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
498: @*/
499: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
500: {
501: DM dm;
502: PetscDualSpaceType type;
503: const char *name;
509: PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
510: PetscObjectGetName((PetscObject) sp, &name);
511: PetscObjectSetName((PetscObject) *spNew, name);
512: PetscDualSpaceGetType(sp, &type);
513: PetscDualSpaceSetType(*spNew, type);
514: PetscDualSpaceGetDM(sp, &dm);
515: PetscDualSpaceSetDM(*spNew, dm);
517: (*spNew)->order = sp->order;
518: (*spNew)->k = sp->k;
519: (*spNew)->Nc = sp->Nc;
520: (*spNew)->uniform = sp->uniform;
521: if (sp->ops->duplicate) {(*sp->ops->duplicate)(sp, *spNew);}
522: return(0);
523: }
525: /*@
526: PetscDualSpaceGetDM - Get the DM representing the reference cell
528: Not collective
530: Input Parameter:
531: . sp - The PetscDualSpace
533: Output Parameter:
534: . dm - The reference cell
536: Level: intermediate
538: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
539: @*/
540: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
541: {
545: *dm = sp->dm;
546: return(0);
547: }
549: /*@
550: PetscDualSpaceSetDM - Get the DM representing the reference cell
552: Not collective
554: Input Parameters:
555: + sp - The PetscDualSpace
556: - dm - The reference cell
558: Level: intermediate
560: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
561: @*/
562: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
563: {
569: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up");
570: PetscObjectReference((PetscObject) dm);
571: if (sp->dm && sp->dm != dm) {
572: PetscDualSpaceClearDMData_Internal(sp, sp->dm);
573: }
574: DMDestroy(&sp->dm);
575: sp->dm = dm;
576: return(0);
577: }
579: /*@
580: PetscDualSpaceGetOrder - Get the order of the dual space
582: Not collective
584: Input Parameter:
585: . sp - The PetscDualSpace
587: Output Parameter:
588: . order - The order
590: Level: intermediate
592: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
593: @*/
594: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
595: {
599: *order = sp->order;
600: return(0);
601: }
603: /*@
604: PetscDualSpaceSetOrder - Set the order of the dual space
606: Not collective
608: Input Parameters:
609: + sp - The PetscDualSpace
610: - order - The order
612: Level: intermediate
614: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
615: @*/
616: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
617: {
620: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up");
621: sp->order = order;
622: return(0);
623: }
625: /*@
626: PetscDualSpaceGetNumComponents - Return the number of components for this space
628: Input Parameter:
629: . sp - The PetscDualSpace
631: Output Parameter:
632: . Nc - The number of components
634: Note: A vector space, for example, will have d components, where d is the spatial dimension
636: Level: intermediate
638: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
639: @*/
640: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
641: {
645: *Nc = sp->Nc;
646: return(0);
647: }
649: /*@
650: PetscDualSpaceSetNumComponents - Set the number of components for this space
652: Input Parameters:
653: + sp - The PetscDualSpace
654: - order - The number of components
656: Level: intermediate
658: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
659: @*/
660: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
661: {
664: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
665: sp->Nc = Nc;
666: return(0);
667: }
669: /*@
670: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
672: Not collective
674: Input Parameters:
675: + sp - The PetscDualSpace
676: - i - The basis number
678: Output Parameter:
679: . functional - The basis functional
681: Level: intermediate
683: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
684: @*/
685: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
686: {
687: PetscInt dim;
693: PetscDualSpaceGetDimension(sp, &dim);
694: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
695: *functional = sp->functional[i];
696: return(0);
697: }
699: /*@
700: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
702: Not collective
704: Input Parameter:
705: . sp - The PetscDualSpace
707: Output Parameter:
708: . dim - The dimension
710: Level: intermediate
712: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
713: @*/
714: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
715: {
721: if (sp->spdim < 0) {
722: PetscSection section;
724: PetscDualSpaceGetSection(sp, §ion);
725: if (section) {
726: PetscSectionGetStorageSize(section, &(sp->spdim));
727: } else sp->spdim = 0;
728: }
729: *dim = sp->spdim;
730: return(0);
731: }
733: /*@
734: PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain
736: Not collective
738: Input Parameter:
739: . sp - The PetscDualSpace
741: Output Parameter:
742: . dim - The dimension
744: Level: intermediate
746: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
747: @*/
748: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
749: {
755: if (sp->spintdim < 0) {
756: PetscSection section;
758: PetscDualSpaceGetSection(sp, §ion);
759: if (section) {
760: PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
761: } else sp->spintdim = 0;
762: }
763: *intdim = sp->spintdim;
764: return(0);
765: }
767: /*@
768: PetscDualSpaceGetUniform - Whether this dual space is uniform
770: Not collective
772: Input Parameters:
773: . sp - A dual space
775: Output Parameters:
776: . uniform - PETSC_TRUE if (a) the dual space is the same for each point in a stratum of the reference DMPlex, and
777: (b) every symmetry of each point in the reference DMPlex is also a symmetry of the point's dual space.
779: Level: advanced
781: Note: all of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells
782: with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform.
784: .seealso: PetscDualSpaceGetPointSubspace(), PetscDualSpaceGetSymmetries()
785: @*/
786: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
787: {
791: *uniform = sp->uniform;
792: return(0);
793: }
795: /*@C
796: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
798: Not collective
800: Input Parameter:
801: . sp - The PetscDualSpace
803: Output Parameter:
804: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
806: Level: intermediate
808: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
809: @*/
810: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
811: {
817: 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");
818: if (!sp->numDof) {
819: DM dm;
820: PetscInt depth, d;
821: PetscSection section;
823: PetscDualSpaceGetDM(sp, &dm);
824: DMPlexGetDepth(dm, &depth);
825: PetscCalloc1(depth+1,&(sp->numDof));
826: PetscDualSpaceGetSection(sp, §ion);
827: for (d = 0; d <= depth; d++) {
828: PetscInt dStart, dEnd;
830: DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
831: if (dEnd <= dStart) continue;
832: PetscSectionGetDof(section, dStart, &(sp->numDof[d]));
834: }
835: }
836: *numDof = sp->numDof;
837: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
838: return(0);
839: }
841: /* create the section of the right size and set a permutation for topological ordering */
842: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
843: {
844: DM dm;
845: PetscInt pStart, pEnd, cStart, cEnd, c, depth, count, i;
846: PetscInt *seen, *perm;
847: PetscSection section;
851: dm = sp->dm;
852: PetscSectionCreate(PETSC_COMM_SELF, §ion);
853: DMPlexGetChart(dm, &pStart, &pEnd);
854: PetscSectionSetChart(section, pStart, pEnd);
855: PetscCalloc1(pEnd - pStart, &seen);
856: PetscMalloc1(pEnd - pStart, &perm);
857: DMPlexGetDepth(dm, &depth);
858: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
859: for (c = cStart, count = 0; c < cEnd; c++) {
860: PetscInt closureSize = -1, e;
861: PetscInt *closure = NULL;
863: perm[count++] = c;
864: seen[c-pStart] = 1;
865: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
866: for (e = 0; e < closureSize; e++) {
867: PetscInt point = closure[2*e];
869: if (seen[point-pStart]) continue;
870: perm[count++] = point;
871: seen[point-pStart] = 1;
872: }
873: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
874: }
875: if (count != pEnd - pStart) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering");
876: for (i = 0; i < pEnd - pStart; i++) if (perm[i] != i) break;
877: if (i < pEnd - pStart) {
878: IS permIS;
880: ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
881: ISSetPermutation(permIS);
882: PetscSectionSetPermutation(section, permIS);
883: ISDestroy(&permIS);
884: } else {
885: PetscFree(perm);
886: }
887: PetscFree(seen);
888: *topSection = section;
889: return(0);
890: }
892: /* mark boundary points and set up */
893: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
894: {
895: DM dm;
896: DMLabel boundary;
897: PetscInt pStart, pEnd, p;
901: dm = sp->dm;
902: DMLabelCreate(PETSC_COMM_SELF,"boundary",&boundary);
903: PetscDualSpaceGetDM(sp,&dm);
904: DMPlexMarkBoundaryFaces(dm,1,boundary);
905: DMPlexLabelComplete(dm,boundary);
906: DMPlexGetChart(dm, &pStart, &pEnd);
907: for (p = pStart; p < pEnd; p++) {
908: PetscInt bval;
910: DMLabelGetValue(boundary, p, &bval);
911: if (bval == 1) {
912: PetscInt dof;
914: PetscSectionGetDof(section, p, &dof);
915: PetscSectionSetConstraintDof(section, p, dof);
916: }
917: }
918: DMLabelDestroy(&boundary);
919: PetscSectionSetUp(section);
920: return(0);
921: }
923: /*@
924: PetscDualSpaceGetSection - Create a PetscSection over the reference cell with the layout from this space
926: Collective on sp
928: Input Parameters:
929: . sp - The PetscDualSpace
931: Output Parameter:
932: . section - The section
934: Level: advanced
936: .seealso: PetscDualSpaceCreate(), DMPLEX
937: @*/
938: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
939: {
940: PetscInt pStart, pEnd, p;
944: if (!sp->pointSection) {
945: /* mark the boundary */
946: PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
947: DMPlexGetChart(sp->dm,&pStart,&pEnd);
948: for (p = pStart; p < pEnd; p++) {
949: PetscDualSpace psp;
951: PetscDualSpaceGetPointSubspace(sp, p, &psp);
952: if (psp) {
953: PetscInt dof;
955: PetscDualSpaceGetInteriorDimension(psp, &dof);
956: PetscSectionSetDof(sp->pointSection,p,dof);
957: }
958: }
959: PetscDualSpaceSectionSetUp_Internal(sp,sp->pointSection);
960: }
961: *section = sp->pointSection;
962: return(0);
963: }
965: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
966: * have one cell */
967: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
968: {
969: PetscReal *sv0, *v0, *J;
970: PetscSection section;
971: PetscInt dim, s, k;
972: DM dm;
976: PetscDualSpaceGetDM(sp, &dm);
977: DMGetDimension(dm, &dim);
978: PetscDualSpaceGetSection(sp, §ion);
979: PetscMalloc3(dim, &v0, dim, &sv0, dim*dim, &J);
980: PetscDualSpaceGetFormDegree(sp, &k);
981: for (s = sStart; s < sEnd; s++) {
982: PetscReal detJ, hdetJ;
983: PetscDualSpace ssp;
984: PetscInt dof, off, f, sdim;
985: PetscInt i, j;
986: DM sdm;
988: PetscDualSpaceGetPointSubspace(sp, s, &ssp);
989: if (!ssp) continue;
990: PetscSectionGetDof(section, s, &dof);
991: PetscSectionGetOffset(section, s, &off);
992: /* get the first vertex of the reference cell */
993: PetscDualSpaceGetDM(ssp, &sdm);
994: DMGetDimension(sdm, &sdim);
995: DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);
996: DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);
997: /* compactify Jacobian */
998: for (i = 0; i < dim; i++) for (j = 0; j < sdim; j++) J[i* sdim + j] = J[i * dim + j];
999: for (f = 0; f < dof; f++) {
1000: PetscQuadrature fn;
1002: PetscDualSpaceGetFunctional(ssp, f, &fn);
1003: PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off+f]));
1004: }
1005: }
1006: PetscFree3(v0, sv0, J);
1007: return(0);
1008: }
1010: /*@
1011: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1013: Collective on sp
1015: Input Parameters:
1016: + sp - The PetscDualSpace
1017: . dim - The spatial dimension
1018: - simplex - Flag for simplex, otherwise use a tensor-product cell
1020: Output Parameter:
1021: . refdm - The reference cell
1023: Note: This DM is on PETSC_COMM_SELF.
1025: Level: intermediate
1027: .seealso: PetscDualSpaceCreate(), DMPLEX
1028: @*/
1029: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1030: {
1034: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, simplex), refdm);
1035: return(0);
1036: }
1038: /*@C
1039: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1041: Input Parameters:
1042: + sp - The PetscDualSpace object
1043: . f - The basis functional index
1044: . time - The time
1045: . 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)
1046: . numComp - The number of components for the function
1047: . func - The input function
1048: - ctx - A context for the function
1050: Output Parameter:
1051: . value - numComp output values
1053: Note: The calling sequence for the callback func is given by:
1055: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1056: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1058: Level: beginner
1060: .seealso: PetscDualSpaceCreate()
1061: @*/
1062: 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)
1063: {
1070: (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
1071: return(0);
1072: }
1074: /*@C
1075: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
1077: Input Parameters:
1078: + sp - The PetscDualSpace object
1079: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
1081: Output Parameter:
1082: . spValue - The values of all dual space functionals
1084: Level: beginner
1086: .seealso: PetscDualSpaceCreate()
1087: @*/
1088: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1089: {
1094: (*sp->ops->applyall)(sp, pointEval, spValue);
1095: return(0);
1096: }
1098: /*@C
1099: PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
1101: Input Parameters:
1102: + sp - The PetscDualSpace object
1103: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1105: Output Parameter:
1106: . spValue - The values of interior dual space functionals
1108: Level: beginner
1110: .seealso: PetscDualSpaceCreate()
1111: @*/
1112: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1113: {
1118: (*sp->ops->applyint)(sp, pointEval, spValue);
1119: return(0);
1120: }
1122: /*@C
1123: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
1125: Input Parameters:
1126: + sp - The PetscDualSpace object
1127: . f - The basis functional index
1128: . time - The time
1129: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1130: . Nc - The number of components for the function
1131: . func - The input function
1132: - ctx - A context for the function
1134: Output Parameter:
1135: . value - The output value
1137: Note: The calling sequence for the callback func is given by:
1139: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1140: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1142: and the idea is to evaluate the functional as an integral
1144: $ n(f) = int dx n(x) . f(x)
1146: where both n and f have Nc components.
1148: Level: beginner
1150: .seealso: PetscDualSpaceCreate()
1151: @*/
1152: 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)
1153: {
1154: DM dm;
1155: PetscQuadrature n;
1156: const PetscReal *points, *weights;
1157: PetscReal x[3];
1158: PetscScalar *val;
1159: PetscInt dim, dE, qNc, c, Nq, q;
1160: PetscBool isAffine;
1161: PetscErrorCode ierr;
1166: PetscDualSpaceGetDM(sp, &dm);
1167: PetscDualSpaceGetFunctional(sp, f, &n);
1168: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
1169: if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
1170: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1171: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1172: *value = 0.0;
1173: isAffine = cgeom->isAffine;
1174: dE = cgeom->dimEmbed;
1175: for (q = 0; q < Nq; ++q) {
1176: if (isAffine) {
1177: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
1178: (*func)(dE, time, x, Nc, val, ctx);
1179: } else {
1180: (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
1181: }
1182: for (c = 0; c < Nc; ++c) {
1183: *value += val[c]*weights[q*Nc+c];
1184: }
1185: }
1186: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1187: return(0);
1188: }
1190: /*@C
1191: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData()
1193: Input Parameters:
1194: + sp - The PetscDualSpace object
1195: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()
1197: Output Parameter:
1198: . spValue - The values of all dual space functionals
1200: Level: beginner
1202: .seealso: PetscDualSpaceCreate()
1203: @*/
1204: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1205: {
1206: Vec pointValues, dofValues;
1207: Mat allMat;
1208: PetscErrorCode ierr;
1214: PetscDualSpaceGetAllData(sp, NULL, &allMat);
1215: if (!(sp->allNodeValues)) {
1216: MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1217: }
1218: pointValues = sp->allNodeValues;
1219: if (!(sp->allDofValues)) {
1220: MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1221: }
1222: dofValues = sp->allDofValues;
1223: VecPlaceArray(pointValues, pointEval);
1224: VecPlaceArray(dofValues, spValue);
1225: MatMult(allMat, pointValues, dofValues);
1226: VecResetArray(dofValues);
1227: VecResetArray(pointValues);
1228: return(0);
1229: }
1231: /*@C
1232: PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()
1234: Input Parameters:
1235: + sp - The PetscDualSpace object
1236: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()
1238: Output Parameter:
1239: . spValue - The values of interior dual space functionals
1241: Level: beginner
1243: .seealso: PetscDualSpaceCreate()
1244: @*/
1245: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1246: {
1247: Vec pointValues, dofValues;
1248: Mat intMat;
1249: PetscErrorCode ierr;
1255: PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1256: if (!(sp->intNodeValues)) {
1257: MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1258: }
1259: pointValues = sp->intNodeValues;
1260: if (!(sp->intDofValues)) {
1261: MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1262: }
1263: dofValues = sp->intDofValues;
1264: VecPlaceArray(pointValues, pointEval);
1265: VecPlaceArray(dofValues, spValue);
1266: MatMult(intMat, pointValues, dofValues);
1267: VecResetArray(dofValues);
1268: VecResetArray(pointValues);
1269: return(0);
1270: }
1272: /*@
1273: PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values
1275: Input Parameter:
1276: . sp - The dualspace
1278: Output Parameters:
1279: + allNodes - A PetscQuadrature object containing all evaluation nodes
1280: - allMat - A Mat for the node-to-dof transformation
1282: Level: advanced
1284: .seealso: PetscDualSpaceCreate()
1285: @*/
1286: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1287: {
1294: if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1295: PetscQuadrature qpoints;
1296: Mat amat;
1298: (*sp->ops->createalldata)(sp,&qpoints,&amat);
1299: PetscQuadratureDestroy(&(sp->allNodes));
1300: MatDestroy(&(sp->allMat));
1301: sp->allNodes = qpoints;
1302: sp->allMat = amat;
1303: }
1304: if (allNodes) *allNodes = sp->allNodes;
1305: if (allMat) *allMat = sp->allMat;
1306: return(0);
1307: }
1309: /*@
1310: PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals
1312: Input Parameter:
1313: . sp - The dualspace
1315: Output Parameters:
1316: + allNodes - A PetscQuadrature object containing all evaluation nodes
1317: - allMat - A Mat for the node-to-dof transformation
1319: Level: advanced
1321: .seealso: PetscDualSpaceCreate()
1322: @*/
1323: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1324: {
1325: PetscInt spdim;
1326: PetscInt numPoints, offset;
1327: PetscReal *points;
1328: PetscInt f, dim;
1329: PetscInt Nc, nrows, ncols;
1330: PetscInt maxNumPoints;
1331: PetscQuadrature q;
1332: Mat A;
1333: PetscErrorCode ierr;
1336: PetscDualSpaceGetNumComponents(sp, &Nc);
1337: PetscDualSpaceGetDimension(sp,&spdim);
1338: if (!spdim) {
1339: PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1340: PetscQuadratureSetData(*allNodes,0,0,0,NULL,NULL);
1341: }
1342: nrows = spdim;
1343: PetscDualSpaceGetFunctional(sp,0,&q);
1344: PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
1345: maxNumPoints = numPoints;
1346: for (f = 1; f < spdim; f++) {
1347: PetscInt Np;
1349: PetscDualSpaceGetFunctional(sp,f,&q);
1350: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1351: numPoints += Np;
1352: maxNumPoints = PetscMax(maxNumPoints,Np);
1353: }
1354: ncols = numPoints * Nc;
1355: PetscMalloc1(dim*numPoints,&points);
1356: MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1357: for (f = 0, offset = 0; f < spdim; f++) {
1358: const PetscReal *p, *w;
1359: PetscInt Np, i;
1360: PetscInt fnc;
1362: PetscDualSpaceGetFunctional(sp,f,&q);
1363: PetscQuadratureGetData(q,NULL,&fnc,&Np,&p,&w);
1364: if (fnc != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch");
1365: for (i = 0; i < Np * dim; i++) {
1366: points[offset* dim + i] = p[i];
1367: }
1368: for (i = 0; i < Np * Nc; i++) {
1369: MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1370: }
1371: offset += Np;
1372: }
1373: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1374: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1375: PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1376: PetscQuadratureSetData(*allNodes,dim,0,numPoints,points,NULL);
1377: *allMat = A;
1378: return(0);
1379: }
1381: /*@
1382: PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1383: this space, as well as the matrix that computes the degrees of freedom from the quadrature values. Degrees of
1384: freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the
1385: reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by
1386: PetscDualSpaceGetSection()).
1388: Input Parameter:
1389: . sp - The dualspace
1391: Output Parameters:
1392: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1393: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1394: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1395: npoints is the number of points in intNodes and nc is PetscDualSpaceGetNumComponents().
1397: Level: advanced
1399: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetDimension(), PetscDualSpaceGetNumComponents(), PetscQuadratureGetData()
1400: @*/
1401: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1402: {
1409: if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1410: PetscQuadrature qpoints;
1411: Mat imat;
1413: (*sp->ops->createintdata)(sp,&qpoints,&imat);
1414: PetscQuadratureDestroy(&(sp->intNodes));
1415: MatDestroy(&(sp->intMat));
1416: sp->intNodes = qpoints;
1417: sp->intMat = imat;
1418: }
1419: if (intNodes) *intNodes = sp->intNodes;
1420: if (intMat) *intMat = sp->intMat;
1421: return(0);
1422: }
1424: /*@
1425: PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values
1427: Input Parameter:
1428: . sp - The dualspace
1430: Output Parameters:
1431: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1432: - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1433: the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1434: npoints is the number of points in allNodes and nc is PetscDualSpaceGetNumComponents().
1436: Level: advanced
1438: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetInteriorData()
1439: @*/
1440: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1441: {
1442: DM dm;
1443: PetscInt spdim0;
1444: PetscInt Nc;
1445: PetscInt pStart, pEnd, p, f;
1446: PetscSection section;
1447: PetscInt numPoints, offset, matoffset;
1448: PetscReal *points;
1449: PetscInt dim;
1450: PetscInt *nnz;
1451: PetscQuadrature q;
1452: Mat imat;
1453: PetscErrorCode ierr;
1457: PetscDualSpaceGetSection(sp, §ion);
1458: PetscSectionGetConstrainedStorageSize(section, &spdim0);
1459: if (!spdim0) {
1460: *intNodes = NULL;
1461: *intMat = NULL;
1462: return(0);
1463: }
1464: PetscDualSpaceGetNumComponents(sp, &Nc);
1465: PetscSectionGetChart(section, &pStart, &pEnd);
1466: PetscDualSpaceGetDM(sp, &dm);
1467: DMGetDimension(dm, &dim);
1468: PetscMalloc1(spdim0, &nnz);
1469: for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1470: PetscInt dof, cdof, off, d;
1472: PetscSectionGetDof(section, p, &dof);
1473: PetscSectionGetConstraintDof(section, p, &cdof);
1474: if (!(dof - cdof)) continue;
1475: PetscSectionGetOffset(section, p, &off);
1476: for (d = 0; d < dof; d++, off++, f++) {
1477: PetscInt Np;
1479: PetscDualSpaceGetFunctional(sp,off,&q);
1480: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1481: nnz[f] = Np * Nc;
1482: numPoints += Np;
1483: }
1484: }
1485: MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1486: PetscFree(nnz);
1487: PetscMalloc1(dim*numPoints,&points);
1488: for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1489: PetscInt dof, cdof, off, d;
1491: PetscSectionGetDof(section, p, &dof);
1492: PetscSectionGetConstraintDof(section, p, &cdof);
1493: if (!(dof - cdof)) continue;
1494: PetscSectionGetOffset(section, p, &off);
1495: for (d = 0; d < dof; d++, off++, f++) {
1496: const PetscReal *p;
1497: const PetscReal *w;
1498: PetscInt Np, i;
1500: PetscDualSpaceGetFunctional(sp,off,&q);
1501: PetscQuadratureGetData(q,NULL,NULL,&Np,&p,&w);
1502: for (i = 0; i < Np * dim; i++) {
1503: points[offset + i] = p[i];
1504: }
1505: for (i = 0; i < Np * Nc; i++) {
1506: MatSetValue(imat, f, matoffset + i, w[i],INSERT_VALUES);
1507: }
1508: offset += Np * dim;
1509: matoffset += Np * Nc;
1510: }
1511: }
1512: PetscQuadratureCreate(PETSC_COMM_SELF,intNodes);
1513: PetscQuadratureSetData(*intNodes,dim,0,numPoints,points,NULL);
1514: MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1515: MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1516: *intMat = imat;
1517: return(0);
1518: }
1520: /*@C
1521: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
1523: Input Parameters:
1524: + sp - The PetscDualSpace object
1525: . f - The basis functional index
1526: . time - The time
1527: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1528: . Nc - The number of components for the function
1529: . func - The input function
1530: - ctx - A context for the function
1532: Output Parameter:
1533: . value - The output value (scalar)
1535: Note: The calling sequence for the callback func is given by:
1537: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1538: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1540: and the idea is to evaluate the functional as an integral
1542: $ n(f) = int dx n(x) . f(x)
1544: where both n and f have Nc components.
1546: Level: beginner
1548: .seealso: PetscDualSpaceCreate()
1549: @*/
1550: 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)
1551: {
1552: DM dm;
1553: PetscQuadrature n;
1554: const PetscReal *points, *weights;
1555: PetscScalar *val;
1556: PetscInt dimEmbed, qNc, c, Nq, q;
1557: PetscErrorCode ierr;
1562: PetscDualSpaceGetDM(sp, &dm);
1563: DMGetCoordinateDim(dm, &dimEmbed);
1564: PetscDualSpaceGetFunctional(sp, f, &n);
1565: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1566: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1567: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1568: *value = 0.;
1569: for (q = 0; q < Nq; ++q) {
1570: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1571: for (c = 0; c < Nc; ++c) {
1572: *value += val[c]*weights[q*Nc+c];
1573: }
1574: }
1575: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1576: return(0);
1577: }
1579: /*@
1580: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
1581: given height. This assumes that the reference cell is symmetric over points of this height.
1583: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1584: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1585: support extracting subspaces, then NULL is returned.
1587: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1589: Not collective
1591: Input Parameters:
1592: + sp - the PetscDualSpace object
1593: - height - the height of the mesh point for which the subspace is desired
1595: Output Parameter:
1596: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1597: point, which will be of lesser dimension if height > 0.
1599: Level: advanced
1601: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1602: @*/
1603: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1604: {
1605: PetscInt depth = -1, cStart, cEnd;
1606: DM dm;
1612: 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");
1613: *subsp = NULL;
1614: dm = sp->dm;
1615: DMPlexGetDepth(dm, &depth);
1616: if (height < 0 || height > depth) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height");
1617: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1618: if (height == 0 && cEnd == cStart + 1) {
1619: *subsp = sp;
1620: return(0);
1621: }
1622: if (!sp->heightSpaces) {
1623: PetscInt h;
1624: PetscCalloc1(depth+1, &(sp->heightSpaces));
1626: for (h = 0; h <= depth; h++) {
1627: if (h == 0 && cEnd == cStart + 1) continue;
1628: if (sp->ops->createheightsubspace) {(*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));}
1629: else if (sp->pointSpaces) {
1630: PetscInt hStart, hEnd;
1632: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
1633: if (hEnd > hStart) {
1634: const char *name;
1636: PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1637: if (sp->pointSpaces[hStart]) {
1638: PetscObjectGetName((PetscObject) sp, &name);
1639: PetscObjectSetName((PetscObject) sp->pointSpaces[hStart], name);
1640: }
1641: sp->heightSpaces[h] = sp->pointSpaces[hStart];
1642: }
1643: }
1644: }
1645: }
1646: *subsp = sp->heightSpaces[height];
1647: return(0);
1648: }
1650: /*@
1651: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
1653: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
1654: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
1655: subspaces, then NULL is returned.
1657: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1659: Not collective
1661: Input Parameters:
1662: + sp - the PetscDualSpace object
1663: - point - the point (in the dual space's DM) for which the subspace is desired
1665: Output Parameters:
1666: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1667: point, which will be of lesser dimension if height > 0.
1669: Level: advanced
1671: .seealso: PetscDualSpace
1672: @*/
1673: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1674: {
1675: PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1676: DM dm;
1682: *bdsp = NULL;
1683: dm = sp->dm;
1684: DMPlexGetChart(dm, &pStart, &pEnd);
1685: if (point < pStart || point > pEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point");
1686: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1687: 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 */
1688: *bdsp = sp;
1689: return(0);
1690: }
1691: if (!sp->pointSpaces) {
1692: PetscInt p;
1693: PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));
1695: for (p = 0; p < pEnd - pStart; p++) {
1696: if (p + pStart == cStart && cEnd == cStart + 1) continue;
1697: if (sp->ops->createpointsubspace) {(*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));}
1698: else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1699: PetscInt dim, depth, height;
1700: DMLabel label;
1702: DMPlexGetDepth(dm,&dim);
1703: DMPlexGetDepthLabel(dm,&label);
1704: DMLabelGetValue(label,p+pStart,&depth);
1705: height = dim - depth;
1706: PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1707: PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1708: }
1709: }
1710: }
1711: *bdsp = sp->pointSpaces[point - pStart];
1712: return(0);
1713: }
1715: /*@C
1716: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
1718: Not collective
1720: Input Parameter:
1721: . sp - the PetscDualSpace object
1723: Output Parameters:
1724: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1725: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation
1727: Note: The permutation and flip arrays are organized in the following way
1728: $ perms[p][ornt][dof # on point] = new local dof #
1729: $ flips[p][ornt][dof # on point] = reversal or not
1731: Level: developer
1733: @*/
1734: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1735: {
1742: if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1743: return(0);
1744: }
1746: /*@
1747: PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1748: dual space's functionals.
1750: Input Parameter:
1751: . dsp - The PetscDualSpace
1753: Output Parameter:
1754: . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1755: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1756: 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).
1757: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1758: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1759: but are stored as 1-forms.
1761: Level: developer
1763: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1764: @*/
1765: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1766: {
1770: *k = dsp->k;
1771: return(0);
1772: }
1774: /*@
1775: PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1776: dual space's functionals.
1778: Input Parameters:
1779: + dsp - The PetscDualSpace
1780: - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1781: in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example,
1782: 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).
1783: If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1784: Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1785: but are stored as 1-forms.
1787: Level: developer
1789: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1790: @*/
1791: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1792: {
1793: PetscInt dim;
1797: if (dsp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
1798: dim = dsp->dm->dim;
1799: if (k < -dim || k > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %D-form on %D-dimensional reference cell", PetscAbsInt(k), dim);
1800: dsp->k = k;
1801: return(0);
1802: }
1804: /*@
1805: PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space
1807: Input Parameter:
1808: . dsp - The PetscDualSpace
1810: Output Parameter:
1811: . k - The simplex dimension
1813: Level: developer
1815: Note: Currently supported values are
1816: $ 0: These are H_1 methods that only transform coordinates
1817: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1818: $ 2: These are the same as 1
1819: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1821: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1822: @*/
1823: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1824: {
1825: PetscInt dim;
1830: dim = dsp->dm->dim;
1831: if (!dsp->k) *k = IDENTITY_TRANSFORM;
1832: else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1833: else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1834: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1835: return(0);
1836: }
1838: /*@C
1839: PetscDualSpaceTransform - Transform the function values
1841: Input Parameters:
1842: + dsp - The PetscDualSpace
1843: . trans - The type of transform
1844: . isInverse - Flag to invert the transform
1845: . fegeom - The cell geometry
1846: . Nv - The number of function samples
1847: . Nc - The number of function components
1848: - vals - The function values
1850: Output Parameter:
1851: . vals - The transformed function values
1853: Level: intermediate
1855: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1857: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpaceTransformHessian(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1858: @*/
1859: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1860: {
1861: PetscReal Jstar[9] = {0};
1862: PetscInt dim, v, c, Nk;
1869: /* TODO: not handling dimEmbed != dim right now */
1870: dim = dsp->dm->dim;
1871: /* No change needed for 0-forms */
1872: if (!dsp->k) return(0);
1873: PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1874: /* TODO: use fegeom->isAffine */
1875: PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1876: for (v = 0; v < Nv; ++v) {
1877: switch (Nk) {
1878: case 1:
1879: for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0];
1880: break;
1881: case 2:
1882: for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1883: break;
1884: case 3:
1885: for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1886: break;
1887: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk);
1888: }
1889: }
1890: return(0);
1891: }
1893: /*@C
1894: PetscDualSpaceTransformGradient - Transform the function gradient values
1896: Input Parameters:
1897: + dsp - The PetscDualSpace
1898: . trans - The type of transform
1899: . isInverse - Flag to invert the transform
1900: . fegeom - The cell geometry
1901: . Nv - The number of function gradient samples
1902: . Nc - The number of function components
1903: - vals - The function gradient values
1905: Output Parameter:
1906: . vals - The transformed function gradient values
1908: Level: intermediate
1910: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1912: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1913: @*/
1914: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1915: {
1916: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1917: PetscInt v, c, d;
1923: #ifdef PETSC_USE_DEBUG
1924: if (dE <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE);
1925: #endif
1926: /* Transform gradient */
1927: if (dim == dE) {
1928: for (v = 0; v < Nv; ++v) {
1929: for (c = 0; c < Nc; ++c) {
1930: switch (dim)
1931: {
1932: case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];break;
1933: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1934: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1935: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1936: }
1937: }
1938: }
1939: } else {
1940: for (v = 0; v < Nv; ++v) {
1941: for (c = 0; c < Nc; ++c) {
1942: DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v*Nc+c)*dE], &vals[(v*Nc+c)*dE]);
1943: }
1944: }
1945: }
1946: /* Assume its a vector, otherwise assume its a bunch of scalars */
1947: if (Nc == 1 || Nc != dim) return(0);
1948: switch (trans) {
1949: case IDENTITY_TRANSFORM: break;
1950: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1951: if (isInverse) {
1952: for (v = 0; v < Nv; ++v) {
1953: for (d = 0; d < dim; ++d) {
1954: switch (dim)
1955: {
1956: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1957: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1958: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1959: }
1960: }
1961: }
1962: } else {
1963: for (v = 0; v < Nv; ++v) {
1964: for (d = 0; d < dim; ++d) {
1965: switch (dim)
1966: {
1967: case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1968: case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1969: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1970: }
1971: }
1972: }
1973: }
1974: break;
1975: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1976: if (isInverse) {
1977: for (v = 0; v < Nv; ++v) {
1978: for (d = 0; d < dim; ++d) {
1979: switch (dim)
1980: {
1981: case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1982: case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1983: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1984: }
1985: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1986: }
1987: }
1988: } else {
1989: for (v = 0; v < Nv; ++v) {
1990: for (d = 0; d < dim; ++d) {
1991: switch (dim)
1992: {
1993: case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1994: case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1995: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1996: }
1997: for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1998: }
1999: }
2000: }
2001: break;
2002: }
2003: return(0);
2004: }
2006: /*@C
2007: PetscDualSpaceTransformHessian - Transform the function Hessian values
2009: Input Parameters:
2010: + dsp - The PetscDualSpace
2011: . trans - The type of transform
2012: . isInverse - Flag to invert the transform
2013: . fegeom - The cell geometry
2014: . Nv - The number of function Hessian samples
2015: . Nc - The number of function components
2016: - vals - The function gradient values
2018: Output Parameter:
2019: . vals - The transformed function Hessian values
2021: Level: intermediate
2023: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2025: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
2026: @*/
2027: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
2028: {
2029: const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
2030: PetscInt v, c;
2036: #ifdef PETSC_USE_DEBUG
2037: if (dE <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE);
2038: #endif
2039: /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
2040: if (dim == dE) {
2041: for (v = 0; v < Nv; ++v) {
2042: for (c = 0; c < Nc; ++c) {
2043: switch (dim)
2044: {
2045: case 1: vals[(v*Nc+c)*dim*dim] *= PetscSqr(fegeom->invJ[0]);break;
2046: case 2: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break;
2047: case 3: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break;
2048: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
2049: }
2050: }
2051: }
2052: } else {
2053: for (v = 0; v < Nv; ++v) {
2054: for (c = 0; c < Nc; ++c) {
2055: DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v*Nc+c)*dE*dE], &vals[(v*Nc+c)*dE*dE]);
2056: }
2057: }
2058: }
2059: /* Assume its a vector, otherwise assume its a bunch of scalars */
2060: if (Nc == 1 || Nc != dim) return(0);
2061: switch (trans) {
2062: case IDENTITY_TRANSFORM: break;
2063: case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
2064: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2065: case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2066: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2067: }
2068: return(0);
2069: }
2071: /*@C
2072: 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.
2074: Input Parameters:
2075: + dsp - The PetscDualSpace
2076: . fegeom - The geometry for this cell
2077: . Nq - The number of function samples
2078: . Nc - The number of function components
2079: - pointEval - The function values
2081: Output Parameter:
2082: . pointEval - The transformed function values
2084: Level: advanced
2086: 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.
2088: Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2090: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2091: @*/
2092: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2093: {
2094: PetscDualSpaceTransformType trans;
2095: PetscInt k;
2096: PetscErrorCode ierr;
2102: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2103: This determines their transformation properties. */
2104: PetscDualSpaceGetDeRahm(dsp, &k);
2105: switch (k)
2106: {
2107: case 0: /* H^1 point evaluations */
2108: trans = IDENTITY_TRANSFORM;break;
2109: case 1: /* Hcurl preserves tangential edge traces */
2110: trans = COVARIANT_PIOLA_TRANSFORM;break;
2111: case 2:
2112: case 3: /* Hdiv preserve normal traces */
2113: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2114: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2115: }
2116: PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2117: return(0);
2118: }
2120: /*@C
2121: 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.
2123: Input Parameters:
2124: + dsp - The PetscDualSpace
2125: . fegeom - The geometry for this cell
2126: . Nq - The number of function samples
2127: . Nc - The number of function components
2128: - pointEval - The function values
2130: Output Parameter:
2131: . pointEval - The transformed function values
2133: Level: advanced
2135: 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.
2137: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2139: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2140: @*/
2141: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2142: {
2143: PetscDualSpaceTransformType trans;
2144: PetscInt k;
2145: PetscErrorCode ierr;
2151: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2152: This determines their transformation properties. */
2153: PetscDualSpaceGetDeRahm(dsp, &k);
2154: switch (k)
2155: {
2156: case 0: /* H^1 point evaluations */
2157: trans = IDENTITY_TRANSFORM;break;
2158: case 1: /* Hcurl preserves tangential edge traces */
2159: trans = COVARIANT_PIOLA_TRANSFORM;break;
2160: case 2:
2161: case 3: /* Hdiv preserve normal traces */
2162: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2163: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2164: }
2165: PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2166: return(0);
2167: }
2169: /*@C
2170: 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.
2172: Input Parameters:
2173: + dsp - The PetscDualSpace
2174: . fegeom - The geometry for this cell
2175: . Nq - The number of function gradient samples
2176: . Nc - The number of function components
2177: - pointEval - The function gradient values
2179: Output Parameter:
2180: . pointEval - The transformed function gradient values
2182: Level: advanced
2184: 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.
2186: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2188: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2189: @*/
2190: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2191: {
2192: PetscDualSpaceTransformType trans;
2193: PetscInt k;
2194: PetscErrorCode ierr;
2200: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2201: This determines their transformation properties. */
2202: PetscDualSpaceGetDeRahm(dsp, &k);
2203: switch (k)
2204: {
2205: case 0: /* H^1 point evaluations */
2206: trans = IDENTITY_TRANSFORM;break;
2207: case 1: /* Hcurl preserves tangential edge traces */
2208: trans = COVARIANT_PIOLA_TRANSFORM;break;
2209: case 2:
2210: case 3: /* Hdiv preserve normal traces */
2211: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2212: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2213: }
2214: PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2215: return(0);
2216: }
2218: /*@C
2219: PetscDualSpacePushforwardHessian - Transform the given function Hessian so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.
2221: Input Parameters:
2222: + dsp - The PetscDualSpace
2223: . fegeom - The geometry for this cell
2224: . Nq - The number of function Hessian samples
2225: . Nc - The number of function components
2226: - pointEval - The function gradient values
2228: Output Parameter:
2229: . pointEval - The transformed function Hessian values
2231: Level: advanced
2233: 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.
2235: Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
2237: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2238: @*/
2239: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2240: {
2241: PetscDualSpaceTransformType trans;
2242: PetscInt k;
2243: PetscErrorCode ierr;
2249: /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2250: This determines their transformation properties. */
2251: PetscDualSpaceGetDeRahm(dsp, &k);
2252: switch (k)
2253: {
2254: case 0: /* H^1 point evaluations */
2255: trans = IDENTITY_TRANSFORM;break;
2256: case 1: /* Hcurl preserves tangential edge traces */
2257: trans = COVARIANT_PIOLA_TRANSFORM;break;
2258: case 2:
2259: case 3: /* Hdiv preserve normal traces */
2260: trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2261: default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2262: }
2263: PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2264: return(0);
2265: }