Actual source code: dualspace.c
petsc-3.10.5 2019-03-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: {
141: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
142: if (sp->ops->view) {(*sp->ops->view)(sp, v);}
143: return(0);
144: }
146: /*@
147: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
149: Collective on PetscDualSpace
151: Input Parameter:
152: . sp - the PetscDualSpace object to set options for
154: Options Database:
155: . -petscspace_degree the approximation order of the space
157: Level: developer
159: .seealso PetscDualSpaceView()
160: @*/
161: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
162: {
163: const char *defaultType;
164: char name[256];
165: PetscBool flg;
170: if (!((PetscObject) sp)->type_name) {
171: defaultType = PETSCDUALSPACELAGRANGE;
172: } else {
173: defaultType = ((PetscObject) sp)->type_name;
174: }
175: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
177: PetscObjectOptionsBegin((PetscObject) sp);
178: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
179: if (flg) {
180: PetscDualSpaceSetType(sp, name);
181: } else if (!((PetscObject) sp)->type_name) {
182: PetscDualSpaceSetType(sp, defaultType);
183: }
184: PetscOptionsInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
185: PetscOptionsInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL);
186: if (sp->ops->setfromoptions) {
187: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
188: }
189: /* process any options handlers added with PetscObjectAddOptionsHandler() */
190: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
191: PetscOptionsEnd();
192: PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
193: return(0);
194: }
196: /*@
197: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
199: Collective on PetscDualSpace
201: Input Parameter:
202: . sp - the PetscDualSpace object to setup
204: Level: developer
206: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
207: @*/
208: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
209: {
214: if (sp->setupcalled) return(0);
215: sp->setupcalled = PETSC_TRUE;
216: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
217: return(0);
218: }
220: /*@
221: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
223: Collective on PetscDualSpace
225: Input Parameter:
226: . sp - the PetscDualSpace object to destroy
228: Level: developer
230: .seealso PetscDualSpaceView()
231: @*/
232: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
233: {
234: PetscInt dim, f;
238: if (!*sp) return(0);
241: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
242: ((PetscObject) (*sp))->refct = 0;
244: PetscDualSpaceGetDimension(*sp, &dim);
245: for (f = 0; f < dim; ++f) {
246: PetscQuadratureDestroy(&(*sp)->functional[f]);
247: }
248: PetscFree((*sp)->functional);
249: PetscQuadratureDestroy(&(*sp)->allPoints);
250: DMDestroy(&(*sp)->dm);
252: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
253: PetscHeaderDestroy(sp);
254: return(0);
255: }
257: /*@
258: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
260: Collective on MPI_Comm
262: Input Parameter:
263: . comm - The communicator for the PetscDualSpace object
265: Output Parameter:
266: . sp - The PetscDualSpace object
268: Level: beginner
270: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
271: @*/
272: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
273: {
274: PetscDualSpace s;
279: PetscCitationsRegister(FECitation,&FEcite);
280: *sp = NULL;
281: PetscFEInitializePackage();
283: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
285: s->order = 0;
286: s->Nc = 1;
287: s->setupcalled = PETSC_FALSE;
289: *sp = s;
290: return(0);
291: }
293: /*@
294: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
296: Collective on PetscDualSpace
298: Input Parameter:
299: . sp - The original PetscDualSpace
301: Output Parameter:
302: . spNew - The duplicate PetscDualSpace
304: Level: beginner
306: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
307: @*/
308: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
309: {
315: (*sp->ops->duplicate)(sp, spNew);
316: return(0);
317: }
319: /*@
320: PetscDualSpaceGetDM - Get the DM representing the reference cell
322: Not collective
324: Input Parameter:
325: . sp - The PetscDualSpace
327: Output Parameter:
328: . dm - The reference cell
330: Level: intermediate
332: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
333: @*/
334: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
335: {
339: *dm = sp->dm;
340: return(0);
341: }
343: /*@
344: PetscDualSpaceSetDM - Get the DM representing the reference cell
346: Not collective
348: Input Parameters:
349: + sp - The PetscDualSpace
350: - dm - The reference cell
352: Level: intermediate
354: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
355: @*/
356: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
357: {
363: DMDestroy(&sp->dm);
364: PetscObjectReference((PetscObject) dm);
365: sp->dm = dm;
366: return(0);
367: }
369: /*@
370: PetscDualSpaceGetOrder - Get the order of the dual space
372: Not collective
374: Input Parameter:
375: . sp - The PetscDualSpace
377: Output Parameter:
378: . order - The order
380: Level: intermediate
382: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
383: @*/
384: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
385: {
389: *order = sp->order;
390: return(0);
391: }
393: /*@
394: PetscDualSpaceSetOrder - Set the order of the dual space
396: Not collective
398: Input Parameters:
399: + sp - The PetscDualSpace
400: - order - The order
402: Level: intermediate
404: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
405: @*/
406: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
407: {
410: sp->order = order;
411: return(0);
412: }
414: /*@
415: PetscDualSpaceGetNumComponents - Return the number of components for this space
417: Input Parameter:
418: . sp - The PetscDualSpace
420: Output Parameter:
421: . Nc - The number of components
423: Note: A vector space, for example, will have d components, where d is the spatial dimension
425: Level: intermediate
427: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
428: @*/
429: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
430: {
434: *Nc = sp->Nc;
435: return(0);
436: }
438: /*@
439: PetscDualSpaceSetNumComponents - Set the number of components for this space
441: Input Parameters:
442: + sp - The PetscDualSpace
443: - order - The number of components
445: Level: intermediate
447: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
448: @*/
449: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
450: {
453: sp->Nc = Nc;
454: return(0);
455: }
457: /*@
458: PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
460: Not collective
462: Input Parameter:
463: . sp - The PetscDualSpace
465: Output Parameter:
466: . tensor - Whether the dual space has tensor layout (vs. simplicial)
468: Level: intermediate
470: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
471: @*/
472: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
473: {
479: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
480: return(0);
481: }
483: /*@
484: PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
486: Not collective
488: Input Parameters:
489: + sp - The PetscDualSpace
490: - tensor - Whether the dual space has tensor layout (vs. simplicial)
492: Level: intermediate
494: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
495: @*/
496: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
497: {
502: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
503: return(0);
504: }
506: /*@
507: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
509: Not collective
511: Input Parameters:
512: + sp - The PetscDualSpace
513: - i - The basis number
515: Output Parameter:
516: . functional - The basis functional
518: Level: intermediate
520: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
521: @*/
522: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
523: {
524: PetscInt dim;
530: PetscDualSpaceGetDimension(sp, &dim);
531: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
532: *functional = sp->functional[i];
533: return(0);
534: }
536: /*@
537: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
539: Not collective
541: Input Parameter:
542: . sp - The PetscDualSpace
544: Output Parameter:
545: . dim - The dimension
547: Level: intermediate
549: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
550: @*/
551: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
552: {
558: *dim = 0;
559: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
560: return(0);
561: }
563: /*@C
564: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
566: Not collective
568: Input Parameter:
569: . sp - The PetscDualSpace
571: Output Parameter:
572: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
574: Level: intermediate
576: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
577: @*/
578: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
579: {
585: (*sp->ops->getnumdof)(sp, numDof);
586: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
587: return(0);
588: }
590: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
591: {
592: DM dm;
593: PetscInt pStart, pEnd, depth, h, offset;
594: const PetscInt *numDof;
598: PetscDualSpaceGetDM(sp,&dm);
599: DMPlexGetChart(dm,&pStart,&pEnd);
600: PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
601: PetscSectionSetChart(*section,pStart,pEnd);
602: DMPlexGetDepth(dm,&depth);
603: PetscDualSpaceGetNumDof(sp,&numDof);
604: for (h = 0; h <= depth; h++) {
605: PetscInt hStart, hEnd, p, dof;
607: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
608: dof = numDof[depth - h];
609: for (p = hStart; p < hEnd; p++) {
610: PetscSectionSetDof(*section,p,dof);
611: }
612: }
613: PetscSectionSetUp(*section);
614: for (h = 0, offset = 0; h <= depth; h++) {
615: PetscInt hStart, hEnd, p, dof;
617: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
618: dof = numDof[depth - h];
619: for (p = hStart; p < hEnd; p++) {
620: PetscSectionGetDof(*section,p,&dof);
621: PetscSectionSetOffset(*section,p,offset);
622: offset += dof;
623: }
624: }
625: return(0);
626: }
628: /*@
629: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
631: Collective on PetscDualSpace
633: Input Parameters:
634: + sp - The PetscDualSpace
635: . dim - The spatial dimension
636: - simplex - Flag for simplex, otherwise use a tensor-product cell
638: Output Parameter:
639: . refdm - The reference cell
641: Level: advanced
643: .keywords: PetscDualSpace, reference cell
644: .seealso: PetscDualSpaceCreate(), DMPLEX
645: @*/
646: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
647: {
651: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
652: return(0);
653: }
655: /*@C
656: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
658: Input Parameters:
659: + sp - The PetscDualSpace object
660: . f - The basis functional index
661: . time - The time
662: . 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)
663: . numComp - The number of components for the function
664: . func - The input function
665: - ctx - A context for the function
667: Output Parameter:
668: . value - numComp output values
670: Note: The calling sequence for the callback func is given by:
672: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
673: $ PetscInt numComponents, PetscScalar values[], void *ctx)
675: Level: developer
677: .seealso: PetscDualSpaceCreate()
678: @*/
679: 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)
680: {
687: (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
688: return(0);
689: }
691: /*@C
692: PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
694: Input Parameters:
695: + sp - The PetscDualSpace object
696: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
698: Output Parameter:
699: . spValue - The values of all dual space functionals
701: Level: developer
703: .seealso: PetscDualSpaceCreate()
704: @*/
705: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
706: {
711: (*sp->ops->applyall)(sp, pointEval, spValue);
712: return(0);
713: }
715: /*@C
716: PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.
718: Input Parameters:
719: + sp - The PetscDualSpace object
720: . f - The basis functional index
721: . time - The time
722: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
723: . Nc - The number of components for the function
724: . func - The input function
725: - ctx - A context for the function
727: Output Parameter:
728: . value - The output value
730: Note: The calling sequence for the callback func is given by:
732: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
733: $ PetscInt numComponents, PetscScalar values[], void *ctx)
735: and the idea is to evaluate the functional as an integral
737: $ n(f) = int dx n(x) . f(x)
739: where both n and f have Nc components.
741: Level: developer
743: .seealso: PetscDualSpaceCreate()
744: @*/
745: 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)
746: {
747: DM dm;
748: PetscQuadrature n;
749: const PetscReal *points, *weights;
750: PetscReal x[3];
751: PetscScalar *val;
752: PetscInt dim, dE, qNc, c, Nq, q;
753: PetscBool isAffine;
754: PetscErrorCode ierr;
759: PetscDualSpaceGetDM(sp, &dm);
760: PetscDualSpaceGetFunctional(sp, f, &n);
761: PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
762: if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
763: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
764: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
765: *value = 0.0;
766: isAffine = cgeom->isAffine;
767: dE = cgeom->dimEmbed;
768: for (q = 0; q < Nq; ++q) {
769: if (isAffine) {
770: CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
771: (*func)(dE, time, x, Nc, val, ctx);
772: } else {
773: (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
774: }
775: for (c = 0; c < Nc; ++c) {
776: *value += val[c]*weights[q*Nc+c];
777: }
778: }
779: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
780: return(0);
781: }
783: /*@C
784: PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()
786: Input Parameters:
787: + sp - The PetscDualSpace object
788: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()
790: Output Parameter:
791: . spValue - The values of all dual space functionals
793: Level: developer
795: .seealso: PetscDualSpaceCreate()
796: @*/
797: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
798: {
799: PetscQuadrature n;
800: const PetscReal *points, *weights;
801: PetscInt qNc, c, Nq, q, f, spdim, Nc;
802: PetscInt offset;
803: PetscErrorCode ierr;
809: PetscDualSpaceGetDimension(sp, &spdim);
810: PetscDualSpaceGetNumComponents(sp, &Nc);
811: for (f = 0, offset = 0; f < spdim; f++) {
812: PetscDualSpaceGetFunctional(sp, f, &n);
813: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
814: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
815: spValue[f] = 0.0;
816: for (q = 0; q < Nq; ++q) {
817: for (c = 0; c < Nc; ++c) {
818: spValue[f] += pointEval[offset++]*weights[q*Nc+c];
819: }
820: }
821: }
822: return(0);
823: }
825: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
826: {
832: if (!sp->allPoints && sp->ops->createallpoints) {
833: (*sp->ops->createallpoints)(sp,&sp->allPoints);
834: }
835: *allPoints = sp->allPoints;
836: return(0);
837: }
839: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
840: {
841: PetscInt spdim;
842: PetscInt numPoints, offset;
843: PetscReal *points;
844: PetscInt f, dim;
845: PetscQuadrature q;
846: PetscErrorCode ierr;
849: PetscDualSpaceGetDimension(sp,&spdim);
850: if (!spdim) {
851: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
852: PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
853: }
854: PetscDualSpaceGetFunctional(sp,0,&q);
855: PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
856: for (f = 1; f < spdim; f++) {
857: PetscInt Np;
859: PetscDualSpaceGetFunctional(sp,f,&q);
860: PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
861: numPoints += Np;
862: }
863: PetscMalloc1(dim*numPoints,&points);
864: for (f = 0, offset = 0; f < spdim; f++) {
865: const PetscReal *p;
866: PetscInt Np, i;
868: PetscDualSpaceGetFunctional(sp,f,&q);
869: PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
870: for (i = 0; i < Np * dim; i++) {
871: points[offset + i] = p[i];
872: }
873: offset += Np * dim;
874: }
875: PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
876: PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
877: return(0);
878: }
880: /*@C
881: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.
883: Input Parameters:
884: + sp - The PetscDualSpace object
885: . f - The basis functional index
886: . time - The time
887: . cgeom - A context with geometric information for this cell, we currently just use the centroid
888: . Nc - The number of components for the function
889: . func - The input function
890: - ctx - A context for the function
892: Output Parameter:
893: . value - The output value (scalar)
895: Note: The calling sequence for the callback func is given by:
897: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
898: $ PetscInt numComponents, PetscScalar values[], void *ctx)
900: and the idea is to evaluate the functional as an integral
902: $ n(f) = int dx n(x) . f(x)
904: where both n and f have Nc components.
906: Level: developer
908: .seealso: PetscDualSpaceCreate()
909: @*/
910: 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)
911: {
912: DM dm;
913: PetscQuadrature n;
914: const PetscReal *points, *weights;
915: PetscScalar *val;
916: PetscInt dimEmbed, qNc, c, Nq, q;
917: PetscErrorCode ierr;
922: PetscDualSpaceGetDM(sp, &dm);
923: DMGetCoordinateDim(dm, &dimEmbed);
924: PetscDualSpaceGetFunctional(sp, f, &n);
925: PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
926: if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
927: DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
928: *value = 0.;
929: for (q = 0; q < Nq; ++q) {
930: (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
931: for (c = 0; c < Nc; ++c) {
932: *value += val[c]*weights[q*Nc+c];
933: }
934: }
935: DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
936: return(0);
937: }
939: /*@
940: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
941: given height. This assumes that the reference cell is symmetric over points of this height.
943: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
944: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
945: support extracting subspaces, then NULL is returned.
947: This does not increment the reference count on the returned dual space, and the user should not destroy it.
949: Not collective
951: Input Parameters:
952: + sp - the PetscDualSpace object
953: - height - the height of the mesh point for which the subspace is desired
955: Output Parameter:
956: . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
957: point, which will be of lesser dimension if height > 0.
959: Level: advanced
961: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
962: @*/
963: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
964: {
970: *subsp = NULL;
971: if (sp->ops->getheightsubspace) {
972: (*sp->ops->getheightsubspace)(sp, height, subsp);
973: }
974: return(0);
975: }
977: /*@
978: PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.
980: If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
981: defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
982: subspaces, then NULL is returned.
984: This does not increment the reference count on the returned dual space, and the user should not destroy it.
986: Not collective
988: Input Parameters:
989: + sp - the PetscDualSpace object
990: - point - the point (in the dual space's DM) for which the subspace is desired
992: Output Parameters:
993: bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the
994: point, which will be of lesser dimension if height > 0.
996: Level: advanced
998: .seealso: PetscDualSpace
999: @*/
1000: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1001: {
1007: *bdsp = NULL;
1008: if (sp->ops->getpointsubspace) {
1009: (*sp->ops->getpointsubspace)(sp,point,bdsp);
1010: } else if (sp->ops->getheightsubspace) {
1011: DM dm;
1012: DMLabel label;
1013: PetscInt dim, depth, height;
1015: PetscDualSpaceGetDM(sp,&dm);
1016: DMPlexGetDepth(dm,&dim);
1017: DMPlexGetDepthLabel(dm,&label);
1018: DMLabelGetValue(label,point,&depth);
1019: height = dim - depth;
1020: (*sp->ops->getheightsubspace)(sp,height,bdsp);
1021: }
1022: return(0);
1023: }