1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
PetscSpaceSetFromOptions_Polynomial(PetscSpace sp,PetscOptionItems PetscOptionsObject)3 static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp, PetscOptionItems PetscOptionsObject)
4 {
5 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
6
7 PetscFunctionBegin;
8 PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace polynomial options");
9 PetscCall(PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL));
10 PetscOptionsHeadEnd();
11 PetscFunctionReturn(PETSC_SUCCESS);
12 }
13
PetscSpacePolynomialView_Ascii(PetscSpace sp,PetscViewer v)14 static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
15 {
16 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
17
18 PetscFunctionBegin;
19 PetscCall(PetscViewerASCIIPrintf(v, "%s space of degree %" PetscInt_FMT "\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
PetscSpaceView_Polynomial(PetscSpace sp,PetscViewer viewer)23 static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
24 {
25 PetscBool isascii;
26
27 PetscFunctionBegin;
28 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
29 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
30 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
31 if (isascii) PetscCall(PetscSpacePolynomialView_Ascii(sp, viewer));
32 PetscFunctionReturn(PETSC_SUCCESS);
33 }
34
PetscSpaceDestroy_Polynomial(PetscSpace sp)35 static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
36 {
37 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
38
39 PetscFunctionBegin;
40 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL));
41 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", NULL));
42 if (poly->subspaces) {
43 PetscInt d;
44
45 for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&poly->subspaces[d]));
46 }
47 PetscCall(PetscFree(poly->subspaces));
48 PetscCall(PetscFree(poly));
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
PetscSpaceSetUp_Polynomial(PetscSpace sp)52 static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
53 {
54 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
55
56 PetscFunctionBegin;
57 if (poly->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
58 if (sp->Nv <= 1) poly->tensor = PETSC_FALSE;
59 if (sp->Nc != 1) {
60 PetscInt Nc = sp->Nc;
61 PetscBool tensor = poly->tensor;
62 PetscInt Nv = sp->Nv;
63 PetscInt degree = sp->degree;
64 const char *prefix;
65 const char *name;
66 char subname[PETSC_MAX_PATH_LEN];
67 PetscSpace subsp;
68
69 PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
70 PetscCall(PetscSpaceSumSetNumSubspaces(sp, Nc));
71 PetscCall(PetscSpaceSumSetInterleave(sp, PETSC_TRUE, PETSC_FALSE));
72 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
73 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
74 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
75 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
76 if (((PetscObject)sp)->name) {
77 PetscCall(PetscObjectGetName((PetscObject)sp, &name));
78 PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
79 PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
80 } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
81 PetscCall(PetscSpaceSetType(subsp, PETSCSPACEPOLYNOMIAL));
82 PetscCall(PetscSpaceSetDegree(subsp, degree, PETSC_DETERMINE));
83 PetscCall(PetscSpaceSetNumComponents(subsp, 1));
84 PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
85 PetscCall(PetscSpacePolynomialSetTensor(subsp, tensor));
86 PetscCall(PetscSpaceSetUp(subsp));
87 for (PetscInt i = 0; i < Nc; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
88 PetscCall(PetscSpaceDestroy(&subsp));
89 PetscCall(PetscSpaceSetUp(sp));
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92 if (poly->tensor) {
93 sp->maxDegree = PETSC_DETERMINE;
94 PetscCall(PetscSpaceSetType(sp, PETSCSPACETENSOR));
95 PetscCall(PetscSpaceSetUp(sp));
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98 PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT " invalid", sp->degree);
99 sp->maxDegree = sp->degree;
100 poly->setupcalled = PETSC_TRUE;
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
PetscSpaceGetDimension_Polynomial(PetscSpace sp,PetscInt * dim)104 static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
105 {
106 PetscInt deg = sp->degree;
107 PetscInt n = sp->Nv;
108
109 PetscFunctionBegin;
110 PetscCall(PetscDTBinomialInt(n + deg, n, dim));
111 *dim *= sp->Nc;
112 PetscFunctionReturn(PETSC_SUCCESS);
113 }
114
CoordinateBasis(PetscInt dim,PetscInt npoints,const PetscReal points[],PetscInt jet,PetscInt Njet,PetscReal pScalar[])115 static PetscErrorCode CoordinateBasis(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt jet, PetscInt Njet, PetscReal pScalar[])
116 {
117 PetscFunctionBegin;
118 PetscCall(PetscArrayzero(pScalar, (1 + dim) * Njet * npoints));
119 for (PetscInt b = 0; b < 1 + dim; b++) {
120 for (PetscInt j = 0; j < PetscMin(1 + dim, Njet); j++) {
121 if (j == 0) {
122 if (b == 0) {
123 for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
124 } else {
125 for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = points[pt * dim + (b - 1)];
126 }
127 } else if (j == b) {
128 for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
129 }
130 }
131 }
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
PetscSpaceEvaluate_Polynomial(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])135 static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
136 {
137 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
138 DM dm = sp->dm;
139 PetscInt dim = sp->Nv;
140 PetscInt Nb, jet, Njet;
141 PetscReal *pScalar;
142
143 PetscFunctionBegin;
144 if (!poly->setupcalled) {
145 PetscCall(PetscSpaceSetUp(sp));
146 PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149 PetscCheck(!poly->tensor && sp->Nc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "tensor and multicomponent spaces should have been converted");
150 PetscCall(PetscDTBinomialInt(dim + sp->degree, dim, &Nb));
151 if (H) {
152 jet = 2;
153 } else if (D) {
154 jet = 1;
155 } else {
156 jet = 0;
157 }
158 PetscCall(PetscDTBinomialInt(dim + jet, dim, &Njet));
159 PetscCall(DMGetWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar));
160 // Why are we handling the case degree == 1 specially? Because we don't want numerical noise when we evaluate hat
161 // functions at the vertices of a simplex, which happens when we invert the Vandermonde matrix of the PKD basis.
162 // We don't make any promise about which basis is used.
163 if (sp->degree == 1) {
164 PetscCall(CoordinateBasis(dim, npoints, points, jet, Njet, pScalar));
165 } else {
166 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, sp->degree, jet, pScalar));
167 }
168 if (B) {
169 PetscInt p_strl = Nb;
170 PetscInt b_strl = 1;
171
172 PetscInt b_strr = Njet * npoints;
173 PetscInt p_strr = 1;
174
175 PetscCall(PetscArrayzero(B, npoints * Nb));
176 for (PetscInt b = 0; b < Nb; b++) {
177 for (PetscInt p = 0; p < npoints; p++) B[p * p_strl + b * b_strl] = pScalar[b * b_strr + p * p_strr];
178 }
179 }
180 if (D) {
181 PetscInt p_strl = dim * Nb;
182 PetscInt b_strl = dim;
183 PetscInt d_strl = 1;
184
185 PetscInt b_strr = Njet * npoints;
186 PetscInt d_strr = npoints;
187 PetscInt p_strr = 1;
188
189 PetscCall(PetscArrayzero(D, npoints * Nb * dim));
190 for (PetscInt d = 0; d < dim; d++) {
191 for (PetscInt b = 0; b < Nb; b++) {
192 for (PetscInt p = 0; p < npoints; p++) D[p * p_strl + b * b_strl + d * d_strl] = pScalar[b * b_strr + (1 + d) * d_strr + p * p_strr];
193 }
194 }
195 }
196 if (H) {
197 PetscInt p_strl = dim * dim * Nb;
198 PetscInt b_strl = dim * dim;
199 PetscInt d1_strl = dim;
200 PetscInt d2_strl = 1;
201
202 PetscInt b_strr = Njet * npoints;
203 PetscInt j_strr = npoints;
204 PetscInt p_strr = 1;
205
206 PetscInt *derivs;
207 PetscCall(PetscCalloc1(dim, &derivs));
208 PetscCall(PetscArrayzero(H, npoints * Nb * dim * dim));
209 for (PetscInt d1 = 0; d1 < dim; d1++) {
210 for (PetscInt d2 = 0; d2 < dim; d2++) {
211 PetscInt j;
212 derivs[d1]++;
213 derivs[d2]++;
214 PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
215 derivs[d1]--;
216 derivs[d2]--;
217 for (PetscInt b = 0; b < Nb; b++) {
218 for (PetscInt p = 0; p < npoints; p++) H[p * p_strl + b * b_strl + d1 * d1_strl + d2 * d2_strl] = pScalar[b * b_strr + j * j_strr + p * p_strr];
219 }
220 }
221 }
222 PetscCall(PetscFree(derivs));
223 }
224 PetscCall(DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar));
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
228 /*@
229 PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials.
230
231 Input Parameters:
232 + sp - the function space object
233 - tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space
234
235 Options Database Key:
236 . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension
237
238 Level: intermediate
239
240 Notes:
241 It is a tensor space if it is spanned by polynomials whose degree in each variable is
242 bounded by the given order, as opposed to the space spanned by polynomials
243 whose total degree---summing over all variables---is bounded by the given order.
244
245 .seealso: `PetscSpace`, `PetscSpacePolynomialGetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
246 @*/
PetscSpacePolynomialSetTensor(PetscSpace sp,PetscBool tensor)247 PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
248 {
249 PetscFunctionBegin;
250 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
251 PetscTryMethod(sp, "PetscSpacePolynomialSetTensor_C", (PetscSpace, PetscBool), (sp, tensor));
252 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254
255 /*@
256 PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor
257 polynomials.
258
259 Input Parameter:
260 . sp - the function space object
261
262 Output Parameter:
263 . tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space
264
265 Level: intermediate
266
267 Notes:
268 The space is a tensor space if it is spanned by polynomials whose degree in each variable is
269 bounded by the given order, as opposed to the space spanned by polynomials
270 whose total degree---summing over all variables---is bounded by the given order.
271
272 .seealso: `PetscSpace`, `PetscSpacePolynomialSetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
273 @*/
PetscSpacePolynomialGetTensor(PetscSpace sp,PetscBool * tensor)274 PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
275 {
276 PetscFunctionBegin;
277 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
278 PetscAssertPointer(tensor, 2);
279 PetscTryMethod(sp, "PetscSpacePolynomialGetTensor_C", (PetscSpace, PetscBool *), (sp, tensor));
280 PetscFunctionReturn(PETSC_SUCCESS);
281 }
282
PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp,PetscBool tensor)283 static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
284 {
285 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
286
287 PetscFunctionBegin;
288 poly->tensor = tensor;
289 PetscFunctionReturn(PETSC_SUCCESS);
290 }
291
PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp,PetscBool * tensor)292 static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
293 {
294 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
295
296 PetscFunctionBegin;
297 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
298 PetscAssertPointer(tensor, 2);
299 *tensor = poly->tensor;
300 PetscFunctionReturn(PETSC_SUCCESS);
301 }
302
PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp,PetscInt height,PetscSpace * subsp)303 static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
304 {
305 PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
306 PetscInt Nc, dim, order;
307 PetscBool tensor;
308
309 PetscFunctionBegin;
310 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
311 PetscCall(PetscSpaceGetNumVariables(sp, &dim));
312 PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
313 PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
314 PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
315 if (!poly->subspaces) PetscCall(PetscCalloc1(dim, &poly->subspaces));
316 if (height <= dim) {
317 if (!poly->subspaces[height - 1]) {
318 PetscSpace sub;
319 const char *name;
320
321 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
322 PetscCall(PetscObjectGetName((PetscObject)sp, &name));
323 PetscCall(PetscObjectSetName((PetscObject)sub, name));
324 PetscCall(PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL));
325 PetscCall(PetscSpaceSetNumComponents(sub, Nc));
326 PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
327 PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
328 PetscCall(PetscSpacePolynomialSetTensor(sub, tensor));
329 PetscCall(PetscSpaceSetUp(sub));
330 poly->subspaces[height - 1] = sub;
331 }
332 *subsp = poly->subspaces[height - 1];
333 } else {
334 *subsp = NULL;
335 }
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338
PetscSpaceInitialize_Polynomial(PetscSpace sp)339 static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
340 {
341 PetscFunctionBegin;
342 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
343 sp->ops->setup = PetscSpaceSetUp_Polynomial;
344 sp->ops->view = PetscSpaceView_Polynomial;
345 sp->ops->destroy = PetscSpaceDestroy_Polynomial;
346 sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
347 sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
348 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
349 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial));
350 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial));
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353
354 /*MC
355 PETSCSPACEPOLYNOMIAL = "poly" - A `PetscSpace` object that encapsulates a polynomial space, e.g. P1 is the space of
356 linear polynomials. The space is replicated for each component.
357
358 Level: intermediate
359
360 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
361 M*/
362
PetscSpaceCreate_Polynomial(PetscSpace sp)363 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
364 {
365 PetscSpace_Poly *poly;
366
367 PetscFunctionBegin;
368 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
369 PetscCall(PetscNew(&poly));
370 sp->data = poly;
371
372 poly->tensor = PETSC_FALSE;
373 poly->subspaces = NULL;
374
375 PetscCall(PetscSpaceInitialize_Polynomial(sp));
376 PetscFunctionReturn(PETSC_SUCCESS);
377 }
378