Actual source code: spacepoly.c
1: #include <petsc/private/petscfeimpl.h>
3: static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
4: {
5: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
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: }
14: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
15: {
16: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
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: }
23: static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
24: {
25: PetscBool iascii;
27: PetscFunctionBegin;
30: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
31: if (iascii) PetscCall(PetscSpacePolynomialView_Ascii(sp, viewer));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
36: {
37: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
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;
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: }
52: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
53: {
54: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
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;
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: }
104: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
105: {
106: PetscInt deg = sp->degree;
107: PetscInt n = sp->Nv;
109: PetscFunctionBegin;
110: PetscCall(PetscDTBinomialInt(n + deg, n, dim));
111: *dim *= sp->Nc;
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
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: }
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;
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;
172: PetscInt b_strr = Njet * npoints;
173: PetscInt p_strr = 1;
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;
185: PetscInt b_strr = Njet * npoints;
186: PetscInt d_strr = npoints;
187: PetscInt p_strr = 1;
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;
202: PetscInt b_strr = Njet * npoints;
203: PetscInt j_strr = npoints;
204: PetscInt p_strr = 1;
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: }
228: /*@
229: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials.
231: Input Parameters:
232: + sp - the function space object
233: - tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space
235: Options Database Key:
236: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension
238: Level: intermediate
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.
245: .seealso: `PetscSpace`, `PetscSpacePolynomialGetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
246: @*/
247: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
248: {
249: PetscFunctionBegin;
251: PetscTryMethod(sp, "PetscSpacePolynomialSetTensor_C", (PetscSpace, PetscBool), (sp, tensor));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor
257: polynomials.
259: Input Parameter:
260: . sp - the function space object
262: Output Parameter:
263: . tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space
265: Level: intermediate
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.
272: .seealso: `PetscSpace`, `PetscSpacePolynomialSetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
273: @*/
274: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
275: {
276: PetscFunctionBegin;
278: PetscAssertPointer(tensor, 2);
279: PetscTryMethod(sp, "PetscSpacePolynomialGetTensor_C", (PetscSpace, PetscBool *), (sp, tensor));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
284: {
285: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
287: PetscFunctionBegin;
288: poly->tensor = tensor;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
293: {
294: PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
296: PetscFunctionBegin;
298: PetscAssertPointer(tensor, 2);
299: *tensor = poly->tensor;
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
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;
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;
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: }
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: }
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.
358: Level: intermediate
360: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
361: M*/
363: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
364: {
365: PetscSpace_Poly *poly;
367: PetscFunctionBegin;
369: PetscCall(PetscNew(&poly));
370: sp->data = poly;
372: poly->tensor = PETSC_FALSE;
373: poly->subspaces = NULL;
375: PetscCall(PetscSpaceInitialize_Polynomial(sp));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }