Actual source code: dspacesimple.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: static PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
5: {
6: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
7: DM dm = sp->dm;
8: PetscInt dim, pStart, pEnd;
9: PetscSection section;
11: DMGetDimension(dm, &dim);
12: DMPlexGetChart(dm, &pStart, &pEnd);
13: PetscSectionCreate(PETSC_COMM_SELF, §ion);
14: PetscSectionSetChart(section, pStart, pEnd);
15: PetscSectionSetDof(section, pStart, s->dim);
16: PetscSectionSetUp(section);
17: sp->pointSection = section;
18: return 0;
19: }
21: static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
22: {
23: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
25: PetscFree(s->numDof);
26: PetscFree(s);
27: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
28: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
29: return 0;
30: }
32: static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
33: {
34: PetscInt dim, d;
36: PetscDualSpaceGetDimension(sp, &dim);
37: PetscDualSpaceSimpleSetDimension(spNew, dim);
38: for (d = 0; d < dim; ++d) {
39: PetscQuadrature q;
41: PetscDualSpaceGetFunctional(sp, d, &q);
42: PetscDualSpaceSimpleSetFunctional(spNew, d, q);
43: }
44: return 0;
45: }
47: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
48: {
49: return 0;
50: }
52: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
53: {
54: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
55: DM dm;
56: PetscInt spatialDim, f;
58: for (f = 0; f < s->dim; ++f) PetscQuadratureDestroy(&sp->functional[f]);
59: PetscFree(sp->functional);
60: s->dim = dim;
61: PetscCalloc1(s->dim, &sp->functional);
62: PetscFree(s->numDof);
63: PetscDualSpaceGetDM(sp, &dm);
64: DMGetCoordinateDim(dm, &spatialDim);
65: PetscCalloc1(spatialDim + 1, &s->numDof);
66: s->numDof[spatialDim] = dim;
67: return 0;
68: }
70: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
71: {
72: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
73: PetscReal *weights;
74: PetscInt Nc, c, Nq, p;
77: PetscQuadratureDuplicate(q, &sp->functional[f]);
78: /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
79: PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **)&weights);
80: for (c = 0; c < Nc; ++c) {
81: PetscReal vol = 0.0;
83: for (p = 0; p < Nq; ++p) vol += weights[p * Nc + c];
84: for (p = 0; p < Nq; ++p) weights[p * Nc + c] /= (vol == 0.0 ? 1.0 : vol);
85: }
86: return 0;
87: }
89: /*@
90: PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis
92: Logically Collective on sp
94: Input Parameters:
95: + sp - the `PetscDualSpace`
96: - dim - the basis dimension
98: Level: intermediate
100: .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`
101: @*/
102: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
103: {
107: PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace, PetscInt), (sp, dim));
108: return 0;
109: }
111: /*@
112: PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space
114: Not Collective
116: Input Parameters:
117: + sp - the `PetscDualSpace`
118: . f - the basis index
119: - q - the basis functional
121: Level: intermediate
123: Note:
124: The quadrature will be reweighted so that it has unit volume.
126: .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()`
127: @*/
128: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
129: {
131: PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q));
132: return 0;
133: }
135: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
136: {
137: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
138: sp->ops->setup = PetscDualSpaceSetUp_Simple;
139: sp->ops->view = NULL;
140: sp->ops->destroy = PetscDualSpaceDestroy_Simple;
141: sp->ops->duplicate = PetscDualSpaceDuplicate_Simple;
142: sp->ops->createheightsubspace = NULL;
143: sp->ops->createpointsubspace = NULL;
144: sp->ops->getsymmetries = NULL;
145: sp->ops->apply = PetscDualSpaceApplyDefault;
146: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
147: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
148: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
149: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
150: return 0;
151: }
153: /*MC
154: PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of arbitrary functionals
156: Level: intermediate
158: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
159: M*/
161: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
162: {
163: PetscDualSpace_Simple *s;
166: PetscNew(&s);
167: sp->data = s;
169: s->dim = 0;
170: s->numDof = NULL;
172: PetscDualSpaceInitialize_Simple(sp);
173: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
174: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
175: return 0;
176: }