Actual source code: spacepoly.c
1: #include <petsc/private/petscfeimpl.h>
3: static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
4: {
5: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
7: PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
8: PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
9: PetscOptionsTail();
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 %D\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) {
42: PetscSpaceDestroy(&poly->subspaces[d]);
43: }
44: }
45: PetscFree(poly->subspaces);
46: PetscFree(poly);
47: return 0;
48: }
50: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
51: {
52: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
54: if (poly->setupCalled) return 0;
55: if (sp->Nv <= 1) {
56: poly->tensor = PETSC_FALSE;
57: }
58: if (sp->Nc != 1) {
59: PetscInt Nc = sp->Nc;
60: PetscBool tensor = poly->tensor;
61: PetscInt Nv = sp->Nv;
62: PetscInt degree = sp->degree;
63: const char *prefix;
64: const char *name;
65: char subname[PETSC_MAX_PATH_LEN];
66: PetscSpace subsp;
68: PetscSpaceSetType(sp, PETSCSPACESUM);
69: PetscSpaceSumSetNumSubspaces(sp, Nc);
70: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
71: PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
72: PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
73: PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
74: if (((PetscObject)sp)->name) {
75: PetscObjectGetName((PetscObject)sp, &name);
76: PetscSNPrintf(subname, PETSC_MAX_PATH_LEN-1, "%s sum component", name);
77: PetscObjectSetName((PetscObject)subsp, subname);
78: } else {
79: PetscObjectSetName((PetscObject)subsp, "sum component");
80: }
81: PetscSpaceSetType(subsp, PETSCSPACEPOLYNOMIAL);
82: PetscSpaceSetDegree(subsp, degree, PETSC_DETERMINE);
83: PetscSpaceSetNumComponents(subsp, 1);
84: PetscSpaceSetNumVariables(subsp, Nv);
85: PetscSpacePolynomialSetTensor(subsp, tensor);
86: PetscSpaceSetUp(subsp);
87: for (PetscInt i = 0; i < Nc; i++) {
88: PetscSpaceSumSetSubspace(sp, i, subsp);
89: }
90: PetscSpaceDestroy(&subsp);
91: PetscSpaceSetUp(sp);
92: return 0;
93: }
94: if (poly->tensor) {
95: sp->maxDegree = PETSC_DETERMINE;
96: PetscSpaceSetType(sp, PETSCSPACETENSOR);
97: PetscSpaceSetUp(sp);
98: return 0;
99: }
101: sp->maxDegree = sp->degree;
102: poly->setupCalled = PETSC_TRUE;
103: return 0;
104: }
106: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
107: {
108: PetscInt deg = sp->degree;
109: PetscInt n = sp->Nv;
111: PetscDTBinomialInt(n + deg, n, dim);
112: *dim *= sp->Nc;
113: return 0;
114: }
116: static PetscErrorCode CoordinateBasis(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt jet, PetscInt Njet, PetscReal pScalar[])
117: {
118: 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++) {
124: pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
125: }
126: } else {
127: for (PetscInt pt = 0; pt < npoints; pt++) {
128: pScalar[b * Njet * npoints + j * npoints + pt] = points[pt * dim + (b-1)];
129: }
130: }
131: } else if (j == b) {
132: for (PetscInt pt = 0; pt < npoints; pt++) {
133: pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
134: }
135: }
136: }
137: }
138: return 0;
139: }
141: static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
142: {
143: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
144: DM dm = sp->dm;
145: PetscInt dim = sp->Nv;
146: PetscInt Nb, jet, Njet;
147: PetscReal *pScalar;
149: if (!poly->setupCalled) {
150: PetscSpaceSetUp(sp);
151: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
152: return 0;
153: }
155: PetscDTBinomialInt(dim + sp->degree, dim, &Nb);
156: if (H) {
157: jet = 2;
158: } else if (D) {
159: jet = 1;
160: } else {
161: jet = 0;
162: }
163: PetscDTBinomialInt(dim + jet, dim, &Njet);
164: DMGetWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);
165: // Why are we handling the case degree == 1 specially? Because we don't want numerical noise when we evaluate hat
166: // functions at the vertices of a simplex, which happens when we invert the Vandermonde matrix of the PKD basis.
167: // We don't make any promise about which basis is used.
168: if (sp->degree == 1) {
169: CoordinateBasis(dim, npoints, points, jet, Njet, pScalar);
170: } else {
171: PetscDTPKDEvalJet(dim, npoints, points, sp->degree, jet, pScalar);
172: }
173: if (B) {
174: PetscInt p_strl = Nb;
175: PetscInt b_strl = 1;
177: PetscInt b_strr = Njet*npoints;
178: PetscInt p_strr = 1;
180: PetscArrayzero(B, npoints * Nb);
181: for (PetscInt b = 0; b < Nb; b++) {
182: for (PetscInt p = 0; p < npoints; p++) {
183: B[p*p_strl + b*b_strl] = pScalar[b*b_strr + p*p_strr];
184: }
185: }
186: }
187: if (D) {
188: PetscInt p_strl = dim*Nb;
189: PetscInt b_strl = dim;
190: PetscInt d_strl = 1;
192: PetscInt b_strr = Njet*npoints;
193: PetscInt d_strr = npoints;
194: PetscInt p_strr = 1;
196: PetscArrayzero(D, npoints * Nb * dim);
197: for (PetscInt d = 0; d < dim; d++) {
198: for (PetscInt b = 0; b < Nb; b++) {
199: for (PetscInt p = 0; p < npoints; p++) {
200: D[p*p_strl + b*b_strl + d*d_strl] = pScalar[b*b_strr + (1+d)*d_strr + p*p_strr];
201: }
202: }
203: }
204: }
205: if (H) {
206: PetscInt p_strl = dim*dim*Nb;
207: PetscInt b_strl = dim*dim;
208: PetscInt d1_strl = dim;
209: PetscInt d2_strl = 1;
211: PetscInt b_strr = Njet*npoints;
212: PetscInt j_strr = npoints;
213: PetscInt p_strr = 1;
215: PetscInt *derivs;
216: PetscCalloc1(dim, &derivs);
217: PetscArrayzero(H, npoints * Nb * dim * dim);
218: for (PetscInt d1 = 0; d1 < dim; d1++) {
219: for (PetscInt d2 = 0; d2 < dim; d2++) {
220: PetscInt j;
221: derivs[d1]++;
222: derivs[d2]++;
223: PetscDTGradedOrderToIndex(dim, derivs, &j);
224: derivs[d1]--;
225: derivs[d2]--;
226: for (PetscInt b = 0; b < Nb; b++) {
227: for (PetscInt p = 0; p < npoints; p++) {
228: H[p*p_strl + b*b_strl + d1*d1_strl + d2*d2_strl] = pScalar[b*b_strr + j*j_strr + p*p_strr];
229: }
230: }
231: }
232: }
233: PetscFree(derivs);
234: }
235: DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);
236: return 0;
237: }
239: /*@
240: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
241: by polynomials whose degree in each variable is bounded by the given order), as opposed to polynomials (the space is
242: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
244: Input Parameters:
245: + sp - the function space object
246: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
248: Options Database:
249: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension
251: Level: intermediate
253: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
254: @*/
255: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
256: {
258: PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
259: return 0;
260: }
262: /*@
263: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
264: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
265: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
267: Input Parameters:
268: . sp - the function space object
270: Output Parameters:
271: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
273: Level: intermediate
275: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
276: @*/
277: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
278: {
281: PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
282: return 0;
283: }
285: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
286: {
287: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
289: poly->tensor = tensor;
290: return 0;
291: }
293: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
294: {
295: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
299: *tensor = poly->tensor;
300: return 0;
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: PetscSpaceGetNumComponents(sp, &Nc);
310: PetscSpaceGetNumVariables(sp, &dim);
311: PetscSpaceGetDegree(sp, &order, NULL);
312: PetscSpacePolynomialGetTensor(sp, &tensor);
314: if (!poly->subspaces) PetscCalloc1(dim, &poly->subspaces);
315: if (height <= dim) {
316: if (!poly->subspaces[height-1]) {
317: PetscSpace sub;
318: const char *name;
320: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
321: PetscObjectGetName((PetscObject) sp, &name);
322: PetscObjectSetName((PetscObject) sub, name);
323: PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
324: PetscSpaceSetNumComponents(sub, Nc);
325: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
326: PetscSpaceSetNumVariables(sub, dim-height);
327: PetscSpacePolynomialSetTensor(sub, tensor);
328: PetscSpaceSetUp(sub);
329: poly->subspaces[height-1] = sub;
330: }
331: *subsp = poly->subspaces[height-1];
332: } else {
333: *subsp = NULL;
334: }
335: return 0;
336: }
338: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
339: {
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: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
348: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
349: return 0;
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: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
359: M*/
361: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
362: {
363: PetscSpace_Poly *poly;
366: PetscNewLog(sp,&poly);
367: sp->data = poly;
369: poly->tensor = PETSC_FALSE;
370: poly->subspaces = NULL;
372: PetscSpaceInitialize_Polynomial(sp);
373: return 0;
374: }