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