1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petscdmplex.h>
3
PetscDualSpaceSetUp_Simple(PetscDualSpace sp)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;
10
11 PetscFunctionBegin;
12 PetscCall(DMGetDimension(dm, &dim));
13 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
14 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
15 PetscCall(PetscSectionSetChart(section, pStart, pEnd));
16 PetscCall(PetscSectionSetDof(section, pStart, s->dim));
17 PetscCall(PetscSectionSetUp(section));
18 sp->pointSection = section;
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
PetscDualSpaceDestroy_Simple(PetscDualSpace sp)22 static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
23 {
24 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
25
26 PetscFunctionBegin;
27 PetscCall(PetscFree(s->numDof));
28 PetscCall(PetscFree(s));
29 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", NULL));
30 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", NULL));
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
PetscDualSpaceDuplicate_Simple(PetscDualSpace sp,PetscDualSpace spNew)34 static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
35 {
36 PetscInt dim, d;
37
38 PetscFunctionBegin;
39 PetscCall(PetscDualSpaceGetDimension(sp, &dim));
40 PetscCall(PetscDualSpaceSimpleSetDimension(spNew, dim));
41 for (d = 0; d < dim; ++d) {
42 PetscQuadrature q;
43
44 PetscCall(PetscDualSpaceGetFunctional(sp, d, &q));
45 PetscCall(PetscDualSpaceSimpleSetFunctional(spNew, d, q));
46 }
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp,const PetscInt dim)50 static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
51 {
52 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
53 DM dm;
54 PetscInt spatialDim, f;
55
56 PetscFunctionBegin;
57 for (f = 0; f < s->dim; ++f) PetscCall(PetscQuadratureDestroy(&sp->functional[f]));
58 PetscCall(PetscFree(sp->functional));
59 s->dim = dim;
60 PetscCall(PetscCalloc1(s->dim, &sp->functional));
61 PetscCall(PetscFree(s->numDof));
62 PetscCall(PetscDualSpaceGetDM(sp, &dm));
63 PetscCall(DMGetCoordinateDim(dm, &spatialDim));
64 PetscCall(PetscCalloc1(spatialDim + 1, &s->numDof));
65 s->numDof[spatialDim] = dim;
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68
PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp,PetscInt f,PetscQuadrature q)69 static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
70 {
71 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
72 PetscReal *weights;
73 PetscInt Nc, c, Nq, p;
74
75 PetscFunctionBegin;
76 PetscCheck(!(f < 0) && !(f >= s->dim), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", f, s->dim);
77 PetscCall(PetscQuadratureDuplicate(q, &sp->functional[f]));
78 /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
79 PetscCall(PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **)&weights));
80 for (c = 0; c < Nc; ++c) {
81 PetscReal vol = 0.0;
82
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 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
89 /*@
90 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis
91
92 Logically Collective
93
94 Input Parameters:
95 + sp - the `PetscDualSpace`
96 - dim - the basis dimension
97
98 Level: intermediate
99
100 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`
101 @*/
PetscDualSpaceSimpleSetDimension(PetscDualSpace sp,PetscInt dim)102 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
103 {
104 PetscFunctionBegin;
105 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
106 PetscValidLogicalCollectiveInt(sp, dim, 2);
107 PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up");
108 PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace, PetscInt), (sp, dim));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
112 /*@
113 PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space
114
115 Not Collective
116
117 Input Parameters:
118 + sp - the `PetscDualSpace`
119 . func - the basis index
120 - q - the basis functional
121
122 Level: intermediate
123
124 Note:
125 The quadrature will be reweighted so that it has unit volume.
126
127 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()`
128 @*/
PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp,PetscInt func,PetscQuadrature q)129 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
130 {
131 PetscFunctionBegin;
132 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
133 PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q));
134 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136
PetscDualSpaceInitialize_Simple(PetscDualSpace sp)137 static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
138 {
139 PetscFunctionBegin;
140 sp->ops->setup = PetscDualSpaceSetUp_Simple;
141 sp->ops->view = NULL;
142 sp->ops->destroy = PetscDualSpaceDestroy_Simple;
143 sp->ops->duplicate = PetscDualSpaceDuplicate_Simple;
144 sp->ops->createheightsubspace = NULL;
145 sp->ops->createpointsubspace = NULL;
146 sp->ops->getsymmetries = NULL;
147 sp->ops->apply = PetscDualSpaceApplyDefault;
148 sp->ops->applyall = PetscDualSpaceApplyAllDefault;
149 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
150 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
151 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
155 /*MC
156 PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of functionals provided with `PetscDualSpaceSimpleSetFunctional()`
157
158 Level: intermediate
159
160 Developer Note:
161 It is not clear this has a good name
162
163 .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
164 M*/
165
PetscDualSpaceCreate_Simple(PetscDualSpace sp)166 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
167 {
168 PetscDualSpace_Simple *s;
169
170 PetscFunctionBegin;
171 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
172 PetscCall(PetscNew(&s));
173 sp->data = s;
174
175 s->dim = 0;
176 s->numDof = NULL;
177
178 PetscCall(PetscDualSpaceInitialize_Simple(sp));
179 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple));
180 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple));
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183