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