Actual source code: dualspace.c
petsc-3.11.4 2019-09-28
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: /*@C
10: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
12: Not Collective
14: Input Parameters:
15: + name - The name of a new user-defined creation routine
16: - create_func - The creation routine itself
18: Notes:
19: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
21: Sample usage:
22: .vb
23: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
24: .ve
26: Then, your PetscDualSpace type can be chosen with the procedural interface via
27: .vb
28: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
29: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
30: .ve
31: or at runtime via the option
32: .vb
33: -petscdualspace_type my_dual_space
34: .ve
36: Level: advanced
38: .keywords: PetscDualSpace, register
39: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()
41: @*/
42: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
43: {
47: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
48: return(0);
49: }
51: /*@C
52: PetscDualSpaceSetType - Builds a particular PetscDualSpace
54: Collective on PetscDualSpace
56: Input Parameters:
57: + sp - The PetscDualSpace object
58: - name - The kind of space
60: Options Database Key:
61: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
63: Level: intermediate
65: .keywords: PetscDualSpace, set, type
66: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
67: @*/
68: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
69: {
70: PetscErrorCode (*r)(PetscDualSpace);
71: PetscBool match;
76: PetscObjectTypeCompare((PetscObject) sp, name, &match);
77: if (match) return(0);
79: if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
80: PetscFunctionListFind(PetscDualSpaceList, name, &r);
81: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
83: if (sp->ops->destroy) {
84: (*sp->ops->destroy)(sp);
85: sp->ops->destroy = NULL;
86: }
87: (*r)(sp);
88: PetscObjectChangeTypeName((PetscObject) sp, name);
89: return(0);
90: }
92: /*@C
93: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
95: Not Collective
97: Input Parameter:
98: . sp - The PetscDualSpace
100: Output Parameter:
101: . name - The PetscDualSpace type name
103: Level: intermediate
105: .keywords: PetscDualSpace, get, type, name
106: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
107: @*/
108: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
109: {
115: if (!PetscDualSpaceRegisterAllCalled) {
116: PetscDualSpaceRegisterAll();
117: }
118: *name = ((PetscObject) sp)->type_name;
119: return(0);
120: }
122: /*@
123: PetscDualSpaceView - Views a PetscDualSpace
125: Collective on PetscDualSpace
127: Input Parameter:
128: + sp - the PetscDualSpace object to view
129: - v - the viewer
131: Level: developer
133: .seealso PetscDualSpaceDestroy()
134: @*/
135: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
136: {
137: PetscBool iascii;
143: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
144: PetscObjectPrintClassNamePrefixType((PetscObject)sp, v);
145: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
146: PetscViewerASCIIPushTab(v);
147: if (iascii) {PetscViewerASCIIPrintf(v, "Dual space of order %D with %D components\n", sp->order, sp->Nc);}
148: if (sp->ops->view) {(*sp->ops->view)(sp, v);}
149: PetscViewerASCIIPopTab(v);
150: return(0);
151: }
153: /*@
154: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
156: Collective on PetscDualSpace
158: Input Parameter:
159: . sp - the PetscDualSpace object to set options for
161: Options Database:
162: . -petscspace_degree the approximation order of the space
164: Level: developer
166: .seealso PetscDualSpaceView()
167: @*/
168: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
169: {
170: const char *defaultType;
171: char name[256];
172: PetscBool flg;
177: if (!((PetscObject) sp)->type_name) {
178: defaultType = PETSCDUALSPACELAGRANGE;
179: } else {
180: defaultType = ((PetscObject) sp)->type_name;
181: }
182: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
184: PetscObjectOptionsBegin((PetscObject) sp);
185: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
186: if (flg) {
187: PetscDualSpaceSetType(sp, name);
188: } else if (!((PetscObject) sp)->type_name) {
189: PetscDualSpaceSetType(sp, defaultType);
190: }
191: PetscOptionsInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
192: PetscOptionsInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL);
193: if (sp->ops->setfromoptions) {
194: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
195: }
196: /* process any options handlers added with PetscObjectAddOptionsHandler() */
197: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
198: PetscOptionsEnd();
199: PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
200: return(0);
201: }
203: /*@
204: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
206: Collective on PetscDualSpace
208: Input Parameter:
209: . sp - the PetscDualSpace object to setup
211: Level: developer
213: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
214: @*/
215: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
216: {
221: if (sp->setupcalled) return(0);
222: sp->setupcalled = PETSC_TRUE;
223: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
224: return(0);
225: }
227: /*@
228: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
230: Collective on PetscDualSpace
232: Input Parameter:
233: . sp - the PetscDualSpace object to destroy
235: Level: developer
237: .seealso PetscDualSpaceView()
238: @*/
239: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
240: {
241: PetscInt dim, f;
245: if (!*sp) return(0);
248: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
249: ((PetscObject) (*sp))->refct = 0;
251: PetscDualSpaceGetDimension(*sp, &dim);
252: for (f = 0; f < dim; ++f) {
253: PetscQuadratureDestroy(&(*sp)->functional[f]);
254: }
255: PetscFree((*sp)->functional);
256: PetscQuadratureDestroy(&(*sp)->allPoints);
257: DMDestroy(&(*sp)->dm);
259: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
260: PetscHeaderDestroy(sp);
261: return(0);
262: }
264: /*@
265: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
267: Collective on MPI_Comm
269: Input Parameter:
270: . comm - The communicator for the PetscDualSpace object
272: Output Parameter:
273: . sp - The PetscDualSpace object
275: Level: beginner
277: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
278: @*/
279: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
280: {
281: PetscDualSpace s;
286: PetscCitationsRegister(FECitation,&FEcite);
287: *sp = NULL;
288: PetscFEInitializePackage();
290: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
292: s->order = 0;
293: s->Nc = 1;
294: s->setupcalled = PETSC_FALSE;
296: *sp = s;
297: return(0);
298: }
300: /*@
301: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
303: Collective on PetscDualSpace
305: Input Parameter:
306: . sp - The original PetscDualSpace
308: Output Parameter:
309: . spNew - The duplicate PetscDualSpace
311: Level: beginner
313: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
314: @*/
315: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
316: {
322: (*sp->ops->duplicate)(sp, spNew);
323: return(0);
324: }
326: /*@
327: PetscDualSpaceGetDM - Get the DM representing the reference cell
329: Not collective
331: Input Parameter:
332: . sp - The PetscDualSpace
334: Output Parameter:
335: . dm - The reference cell
337: Level: intermediate
339: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
340: @*/
341: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
342: {
346: *dm = sp->dm;
347: return(0);
348: }
350: /*@
351: PetscDualSpaceSetDM - Get the DM representing the reference cell
353: Not collective
355: Input Parameters:
356: + sp - The PetscDualSpace
357: - dm - The reference cell
359: Level: intermediate
361: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
362: @*/
363: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
364: {
370: DMDestroy(&sp->dm);
371: PetscObjectReference((PetscObject) dm);
372: sp->dm = dm;
373: return(0);
374: }
376: /*@
377: PetscDualSpaceGetOrder - Get the order of the dual space
379: Not collective
381: Input Parameter:
382: . sp - The PetscDualSpace
384: Output Parameter:
385: . order - The order
387: Level: intermediate
389: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
390: @*/
391: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
392: {
396: *order = sp->order;
397: return(0);
398: }
400: /*@
401: PetscDualSpaceSetOrder - Set the order of the dual space
403: Not collective
405: Input Parameters:
406: + sp - The PetscDualSpace
407: - order - The order
409: Level: intermediate
411: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
412: @*/
413: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
414: {
417: sp->order = order;
418: return(0);
419: }
421: /*@
422: PetscDualSpaceGetNumComponents - Return the number of components for this space
424: Input Parameter:
425: . sp - The PetscDualSpace
427: Output Parameter:
428: . Nc - The number of components
430: Note: A vector space, for example, will have d components, where d is the spatial dimension
432: Level: intermediate
434: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
435: @*/
436: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
437: {
441: *Nc = sp->Nc;
442: return(0);
443: }
445: /*@
446: PetscDualSpaceSetNumComponents - Set the number of components for this space
448: Input Parameters:
449: + sp - The PetscDualSpace
450: - order - The number of components
452: Level: intermediate
454: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
455: @*/
456: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
457: {
460: sp->Nc = Nc;
461: return(0);
462: }
464: /*@
465: PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
467: Not collective
469: Input Parameter:
470: . sp - The PetscDualSpace
472: Output Parameter:
473: . tensor - Whether the dual space has tensor layout (vs. simplicial)
475: Level: intermediate
477: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
478: @*/
479: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
480: {
486: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
487: return(0);
488: }
490: /*@
491: PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
493: Not collective
495: Input Parameters:
496: + sp - The PetscDualSpace
497: - tensor - Whether the dual space has tensor layout (vs. simplicial)
499: Level: intermediate
501: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
502: @*/
503: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
504: {
509: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
510: return(0);
511: }
513: /*@
514: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
516: Not collective
518: Input Parameters:
519: + sp - The PetscDualSpace
520: - i - The basis number
522: Output Parameter:
523: . functional - The basis functional
525: Level: intermediate
527: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
528: @*/
529: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
530: {
531: PetscInt dim;
537: PetscDualSpaceGetDimension(sp, &dim);
538: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
539: *functional = sp->functional[i];
540: return(0);
541: }
543: /*@
544: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
546: Not collective
548: Input Parameter:
549: . sp - The PetscDualSpace
551: Output Parameter:
552: . dim - The dimension
554: Level: intermediate
556: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
557: @*/
558: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
559: {
565: *dim = 0;
566: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
567: return(0);
568: }
570: /*@C
571: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
573: Not collective
575: Input Parameter:
576: . sp - The PetscDualSpace
578: Output Parameter:
579: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
581: Level: intermediate
583: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
584: @*/
585: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
586: {
592: (*sp->ops->getnumdof)(sp, numDof);
593: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
594: return(0);
595: }
597: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
598: {
599: DM dm;
600: PetscInt pStart, pEnd, depth, h, offset;
601: const PetscInt *numDof;
605: PetscDualSpaceGetDM(sp,&dm);
606: DMPlexGetChart(dm,&pStart,&pEnd);
607: PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
608: PetscSectionSetChart(*section,pStart,pEnd);
609: DMPlexGetDepth(dm,&depth);
610: PetscDualSpaceGetNumDof(sp,&numDof);
611: for (h = 0; h <= depth; h++) {
612: PetscInt hStart, hEnd, p, dof;
614: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
615: dof = numDof[depth - h];
616: for (p = hStart; p < hEnd; p++) {
617: PetscSectionSetDof(*section,p,dof);
618: }
619: }
620: PetscSectionSetUp(*section);
621: for (h = 0, offset = 0; h <= depth; h++) {
622: PetscInt hStart, hEnd, p, dof;
624: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
625: dof = numDof[depth - h];
626: for (p = hStart; p < hEnd; p++) {
627: PetscSectionGetDof(*section,p,&dof);
628: PetscSectionSetOffset(*section,p,offset);
629: offset += dof;
630: }
631: }
632: return(0);
633: }
635: /*@
636: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
638: Collective on PetscDualSpace
640: Input Parameters:
641: + sp - The PetscDualSpace
642: . dim - The spatial dimension
643: - simplex - Flag for simplex, otherwise use a tensor-product cell
645: Output Parameter:
646: . refdm - The reference cell
648: Level: advanced
650: .keywords: PetscDualSpace, reference cell
651: .seealso: PetscDualSpaceCreate(), DMPLEX
652: @*/
653: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
654: {
658: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
659: return(0);
660: }
662: /*@C
663: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
665: Input Parameters:
666: + sp - The PetscDualSpace object
667: . f - The basis functional index
668: . time - The time
669: . 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)
670: . numComp - The number of components for the function
671: . func - The input function
672: - ctx - A context for the function
674: Output Parameter:
675: . value - numComp output values
677: Note: The calling sequence for the callback func is given by:
679: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
680: $ PetscInt numComponents, PetscScalar values[], void *ctx)
682: Level: developer
684: .seealso: PetscDualSpaceCreate()
685: @*/
686: 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)
687: {
694: (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
695: return(0);
696: }
698: /*@C
699: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
701: Input Parameters:
702: + sp - The PetscDualSpace object
703: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
705: Output Parameter:
706: . spValue - The values of all dual space functionals
708: Level: developer
710: .seealso: PetscDualSpaceCreate()
711: @*/
712: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
713: {
718: (*sp->ops->applyall)(sp, pointEval, spValue);
719: return(0);
720: }
722: /*@C
723: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
725: Input Parameters:
726: + sp - The PetscDualSpace object
727: . f - The basis functional index
728: . time - The time
729: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
730: . Nc - The number of components for the function
731: . func - The input function
732: - ctx - A context for the function
734: Output Parameter:
735: . value - The output value
737: Note: The calling sequence for the callback func is given by:
739: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
740: $ PetscInt numComponents, PetscScalar values[], void *ctx)
742: and the idea is to evaluate the functional as an integral
744: $ n(f) = int dx n(x) . f(x)
746: where both n and f have Nc components.
748: Level: developer
750: .seealso: PetscDualSpaceCreate()
751: @*/
752: 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)
753: {
754: DM dm;
755: PetscQuadrature n;
756: const PetscReal *points, *weights;
757: PetscReal x[3];
758: PetscScalar *val;
759: PetscInt dim, dE, qNc, c, Nq, q;
760: PetscBool isAffine;
761: PetscErrorCode ierr;
766: PetscDualSpaceGetDM(sp, &dm);
767: PetscDualSpaceGetFunctional(sp, f, &n);
768: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
769: if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
770: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
771: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
772: *value = 0.0;
773: isAffine = cgeom->isAffine;
774: dE = cgeom->dimEmbed;
775: for (q = 0; q < Nq; ++q) {
776: if (isAffine) {
777: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
778: (*func)(dE, time, x, Nc, val, ctx);
779: } else {
780: (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
781: }
782: for (c = 0; c < Nc; ++c) {
783: *value += val[c]*weights[q*Nc+c];
784: }
785: }
786: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
787: return(0);
788: }
790: /*@C
791: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
793: Input Parameters:
794: + sp - The PetscDualSpace object
795: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
797: Output Parameter:
798: . spValue - The values of all dual space functionals
800: Level: developer
802: .seealso: PetscDualSpaceCreate()
803: @*/
804: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
805: {
806: PetscQuadrature n;
807: const PetscReal *points, *weights;
808: PetscInt qNc, c, Nq, q, f, spdim, Nc;
809: PetscInt offset;
810: PetscErrorCode ierr;
816: PetscDualSpaceGetDimension(sp, &spdim);
817: PetscDualSpaceGetNumComponents(sp, &Nc);
818: for (f = 0, offset = 0; f < spdim; f++) {
819: PetscDualSpaceGetFunctional(sp, f, &n);
820: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
821: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
822: spValue[f] = 0.0;
823: for (q = 0; q < Nq; ++q) {
824: for (c = 0; c < Nc; ++c) {
825: spValue[f] += pointEval[offset++]*weights[q*Nc+c];
826: }
827: }
828: }
829: return(0);
830: }
832: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
833: {
839: if (!sp->allPoints && sp->ops->createallpoints) {
840: (*sp->ops->createallpoints)(sp,&sp->allPoints);
841: }
842: *allPoints = sp->allPoints;
843: return(0);
844: }
846: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
847: {
848: PetscInt spdim;
849: PetscInt numPoints, offset;
850: PetscReal *points;
851: PetscInt f, dim;
852: PetscQuadrature q;
853: PetscErrorCode ierr;
856: PetscDualSpaceGetDimension(sp,&spdim);
857: if (!spdim) {
858: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
859: PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
860: }
861: PetscDualSpaceGetFunctional(sp,0,&q);
862: PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
863: for (f = 1; f < spdim; f++) {
864: PetscInt Np;
866: PetscDualSpaceGetFunctional(sp,f,&q);
867: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
868: numPoints += Np;
869: }
870: PetscMalloc1(dim*numPoints,&points);
871: for (f = 0, offset = 0; f < spdim; f++) {
872: const PetscReal *p;
873: PetscInt Np, i;
875: PetscDualSpaceGetFunctional(sp,f,&q);
876: PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
877: for (i = 0; i < Np * dim; i++) {
878: points[offset + i] = p[i];
879: }
880: offset += Np * dim;
881: }
882: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
883: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
884: return(0);
885: }
887: /*@C
888: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
890: Input Parameters:
891: + sp - The PetscDualSpace object
892: . f - The basis functional index
893: . time - The time
894: . cgeom - A context with geometric information for this cell, we currently just use the centroid
895: . Nc - The number of components for the function
896: . func - The input function
897: - ctx - A context for the function
899: Output Parameter:
900: . value - The output value (scalar)
902: Note: The calling sequence for the callback func is given by:
904: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
905: $ PetscInt numComponents, PetscScalar values[], void *ctx)
907: and the idea is to evaluate the functional as an integral
909: $ n(f) = int dx n(x) . f(x)
911: where both n and f have Nc components.
913: Level: developer
915: .seealso: PetscDualSpaceCreate()
916: @*/
917: 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)
918: {
919: DM dm;
920: PetscQuadrature n;
921: const PetscReal *points, *weights;
922: PetscScalar *val;
923: PetscInt dimEmbed, qNc, c, Nq, q;
924: PetscErrorCode ierr;
929: PetscDualSpaceGetDM(sp, &dm);
930: DMGetCoordinateDim(dm, &dimEmbed);
931: PetscDualSpaceGetFunctional(sp, f, &n);
932: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
933: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
934: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
935: *value = 0.;
936: for (q = 0; q < Nq; ++q) {
937: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
938: for (c = 0; c < Nc; ++c) {
939: *value += val[c]*weights[q*Nc+c];
940: }
941: }
942: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
943: return(0);
944: }
946: /*@
947: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
948: given height. This assumes that the reference cell is symmetric over points of this height.
950: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
951: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
952: support extracting subspaces, then NULL is returned.
954: This does not increment the reference count on the returned dual space, and the user should not destroy it.
956: Not collective
958: Input Parameters:
959: + sp - the PetscDualSpace object
960: - height - the height of the mesh point for which the subspace is desired
962: Output Parameter:
963: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
964: point, which will be of lesser dimension if height > 0.
966: Level: advanced
968: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
969: @*/
970: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
971: {
977: *subsp = NULL;
978: if (sp->ops->getheightsubspace) {
979: (*sp->ops->getheightsubspace)(sp, height, subsp);
980: }
981: return(0);
982: }
984: /*@
985: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
987: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
988: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
989: subspaces, then NULL is returned.
991: This does not increment the reference count on the returned dual space, and the user should not destroy it.
993: Not collective
995: Input Parameters:
996: + sp - the PetscDualSpace object
997: - point - the point (in the dual space's DM) for which the subspace is desired
999: Output Parameters:
1000: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1001: point, which will be of lesser dimension if height > 0.
1003: Level: advanced
1005: .seealso: PetscDualSpace
1006: @*/
1007: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1008: {
1014: *bdsp = NULL;
1015: if (sp->ops->getpointsubspace) {
1016: (*sp->ops->getpointsubspace)(sp,point,bdsp);
1017: } else if (sp->ops->getheightsubspace) {
1018: DM dm;
1019: DMLabel label;
1020: PetscInt dim, depth, height;
1022: PetscDualSpaceGetDM(sp,&dm);
1023: DMPlexGetDepth(dm,&dim);
1024: DMPlexGetDepthLabel(dm,&label);
1025: DMLabelGetValue(label,point,&depth);
1026: height = dim - depth;
1027: (*sp->ops->getheightsubspace)(sp,height,bdsp);
1028: }
1029: return(0);
1030: }