Actual source code: spacepoly.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/petscfeimpl.h>
3: const char *const PetscSpacePolynomialTypes[] = {"P", "PMINUS_HDIV", "PMINUS_HCURL", "PetscSpacePolynomialType", "PETSCSPACE_POLYNOMIALTYPE_", NULL};
5: static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
6: {
7: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
8: PetscErrorCode ierr;
11: PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
12: PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
13: PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
14: PetscOptionsEnum("-petscspace_poly_type", "Type of polynomial space", "PetscSpacePolynomialSetType", PetscSpacePolynomialTypes, (PetscEnum)poly->ptype, (PetscEnum*)&poly->ptype, NULL);
15: PetscOptionsTail();
16: return(0);
17: }
19: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
20: {
21: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
22: PetscErrorCode ierr;
25: PetscViewerASCIIPrintf(v, "%s%s%s space of degree %D\n", poly->ptype ? PetscSpacePolynomialTypes[poly->ptype] : "", poly->ptype ? " " : "", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);
26: return(0);
27: }
29: static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
30: {
31: PetscBool iascii;
37: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
38: if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
39: return(0);
40: }
42: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
43: {
44: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
45: PetscInt ndegree = sp->degree+1;
46: PetscInt deg;
47: PetscErrorCode ierr;
50: if (poly->setupCalled) return(0);
51: PetscMalloc1(ndegree, &poly->degrees);
52: for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
53: if (poly->tensor) {
54: sp->maxDegree = sp->degree + PetscMax(sp->Nv - 1,0);
55: } else {
56: sp->maxDegree = sp->degree;
57: }
58: poly->setupCalled = PETSC_TRUE;
59: return(0);
60: }
62: static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
63: {
64: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
65: PetscErrorCode ierr;
68: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
69: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);
70: PetscFree(poly->degrees);
71: if (poly->subspaces) {
72: PetscInt d;
74: for (d = 0; d < sp->Nv; ++d) {
75: PetscSpaceDestroy(&poly->subspaces[d]);
76: }
77: }
78: PetscFree(poly->subspaces);
79: PetscFree(poly);
80: return(0);
81: }
83: /* We treat the space as a tensor product of scalar polynomial spaces, so the dimension is multiplied by Nc */
84: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
85: {
86: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
87: PetscInt deg = sp->degree;
88: PetscInt n = sp->Nv, N, i;
89: PetscReal D = 1.0;
92: if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) --deg;
93: if (poly->tensor) {
94: N = 1;
95: for (i = 0; i < n; ++i) N *= (deg+1);
96: } else {
97: for (i = 1; i <= n; ++i) {
98: D *= ((PetscReal) (deg+i))/i;
99: }
100: N = (PetscInt) (D + 0.5);
101: }
102: if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) {
103: N *= sp->Nc + 1;
104: } else {
105: N *= sp->Nc;
106: }
107: *dim = N;
108: return(0);
109: }
111: /*
112: LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.
114: Input Parameters:
115: + len - The length of the tuple
116: . sum - The sum of all entries in the tuple
117: - ind - The current multi-index of the tuple, initialized to the 0 tuple
119: Output Parameter:
120: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
121: . tup - A tuple of len integers addig to sum
123: Level: developer
125: .seealso:
126: */
127: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
128: {
129: PetscInt i;
133: if (len == 1) {
134: ind[0] = -1;
135: tup[0] = sum;
136: } else if (sum == 0) {
137: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
138: } else {
139: tup[0] = sum - ind[0];
140: LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
141: if (ind[1] < 0) {
142: if (ind[0] == sum) {ind[0] = -1;}
143: else {ind[1] = 0; ++ind[0];}
144: }
145: }
146: return(0);
147: }
149: /*
150: TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.
152: Input Parameters:
153: + len - The length of the tuple
154: . max - The max for all entries in the tuple
155: - ind - The current multi-index of the tuple, initialized to the 0 tuple
157: Output Parameter:
158: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
159: . tup - A tuple of len integers less than max
161: Level: developer
163: .seealso:
164: */
165: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
166: {
167: PetscInt i;
171: if (len == 1) {
172: tup[0] = ind[0]++;
173: ind[0] = ind[0] >= max ? -1 : ind[0];
174: } else if (max == 0) {
175: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
176: } else {
177: tup[0] = ind[0];
178: TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
179: if (ind[1] < 0) {
180: ind[1] = 0;
181: if (ind[0] == max-1) {ind[0] = -1;}
182: else {++ind[0];}
183: }
184: }
185: return(0);
186: }
188: /*
189: p in [0, npoints), i in [0, pdim), c in [0, Nc)
191: B[p][i][c] = B[p][i_scalar][c][c]
192: */
193: static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
194: {
195: const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
196: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
197: DM dm = sp->dm;
198: PetscInt Nc = sp->Nc;
199: PetscInt ndegree = sp->degree+1;
200: PetscInt *degrees = poly->degrees;
201: PetscInt dim = sp->Nv;
202: PetscReal *lpoints, *tmp, *LB, *LD, *LH;
203: PetscInt *ind, *tup;
204: PetscInt c, pdim, pdimRed, d, e, der, der2, i, p, deg, o;
205: PetscErrorCode ierr;
208: PetscSpaceGetDimension(sp, &pdim);
209: pdim /= Nc;
210: DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);
211: DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
212: if (B || D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
213: if (D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
214: if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
215: for (d = 0; d < dim; ++d) {
216: for (p = 0; p < npoints; ++p) {
217: lpoints[p] = points[p*dim+d];
218: }
219: PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
220: /* LB, LD, LH (ndegree * dim x npoints) */
221: for (deg = 0; deg < ndegree; ++deg) {
222: for (p = 0; p < npoints; ++p) {
223: if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
224: if (D || H) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
225: if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
226: }
227: }
228: }
229: /* Multiply by A (pdim x ndegree * dim) */
230: PetscMalloc2(dim,&ind,dim,&tup);
231: if (B) {
232: PetscInt topDegree = sp->degree;
234: /* B (npoints x pdim x Nc) */
235: PetscArrayzero(B, npoints*pdim*Nc*Nc);
236: if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) topDegree--;
237: /* Make complete space portion */
238: if (poly->tensor) {
239: if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Tensor spaces not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
240: i = 0;
241: PetscArrayzero(ind, dim);
242: while (ind[0] >= 0) {
243: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
244: for (p = 0; p < npoints; ++p) {
245: B[(p*pdim + i)*Nc*Nc] = 1.0;
246: for (d = 0; d < dim; ++d) {
247: B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
248: }
249: }
250: ++i;
251: }
252: } else {
253: i = 0;
254: for (o = 0; o <= topDegree; ++o) {
255: PetscArrayzero(ind, dim);
256: while (ind[0] >= 0) {
257: LatticePoint_Internal(dim, o, ind, tup);
258: for (p = 0; p < npoints; ++p) {
259: B[(p*pdim + i)*Nc*Nc] = 1.0;
260: for (d = 0; d < dim; ++d) {
261: B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
262: }
263: }
264: ++i;
265: }
266: }
267: }
268: pdimRed = i;
269: /* Make direct sum basis for multicomponent space */
270: for (p = 0; p < npoints; ++p) {
271: for (i = 0; i < pdimRed; ++i) {
272: for (c = 1; c < Nc; ++c) {
273: B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
274: }
275: }
276: }
277: /* Make homogeneous part */
278: if (topDegree < sp->degree) {
279: if (poly->tensor) {
280: } else {
281: i = pdimRed;
282: PetscArrayzero(ind, dim);
283: while (ind[0] >= 0) {
284: LatticePoint_Internal(dim, topDegree, ind, tup);
285: for (p = 0; p < npoints; ++p) {
286: for (c = 0; c < Nc; ++c) {
287: B[(p*pdim*Nc + i*Nc + c)*Nc + c] = 1.0;
288: for (d = 0; d < dim; ++d) {
289: B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= LB[(tup[d]*dim + d)*npoints + p];
290: }
291: switch (poly->ptype) {
292: case PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV:
293: B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= LB[(c*dim + d)*npoints + p];break;
294: case PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL:
295: {
296: PetscReal sum = 0.0;
297: for (d = 0; d < dim; ++d) for (e = 0; e < dim; ++e) sum += eps[c][d][e]*LB[(d*dim + d)*npoints + p];
298: B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= sum;
299: break;
300: }
301: default: SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Invalid polynomial type %s", PetscSpacePolynomialTypes[poly->ptype]);
302: }
303: }
304: }
305: ++i;
306: }
307: }
308: }
309: }
310: if (D) {
311: if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Derivatives not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
312: /* D (npoints x pdim x Nc x dim) */
313: PetscArrayzero(D, npoints*pdim*Nc*Nc*dim);
314: if (poly->tensor) {
315: i = 0;
316: PetscArrayzero(ind, dim);
317: while (ind[0] >= 0) {
318: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
319: for (p = 0; p < npoints; ++p) {
320: for (der = 0; der < dim; ++der) {
321: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
322: for (d = 0; d < dim; ++d) {
323: if (d == der) {
324: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
325: } else {
326: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
327: }
328: }
329: }
330: }
331: ++i;
332: }
333: } else {
334: i = 0;
335: for (o = 0; o <= sp->degree; ++o) {
336: PetscArrayzero(ind, dim);
337: while (ind[0] >= 0) {
338: LatticePoint_Internal(dim, o, ind, tup);
339: for (p = 0; p < npoints; ++p) {
340: for (der = 0; der < dim; ++der) {
341: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
342: for (d = 0; d < dim; ++d) {
343: if (d == der) {
344: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
345: } else {
346: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
347: }
348: }
349: }
350: }
351: ++i;
352: }
353: }
354: }
355: /* Make direct sum basis for multicomponent space */
356: for (p = 0; p < npoints; ++p) {
357: for (i = 0; i < pdim; ++i) {
358: for (c = 1; c < Nc; ++c) {
359: for (d = 0; d < dim; ++d) {
360: D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d];
361: }
362: }
363: }
364: }
365: }
366: if (H) {
367: if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Hessians not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
368: /* H (npoints x pdim x Nc x Nc x dim x dim) */
369: PetscArrayzero(H, npoints*pdim*Nc*Nc*dim*dim);
370: if (poly->tensor) {
371: i = 0;
372: PetscArrayzero(ind, dim);
373: while (ind[0] >= 0) {
374: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
375: for (p = 0; p < npoints; ++p) {
376: for (der = 0; der < dim; ++der) {
377: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0;
378: for (d = 0; d < dim; ++d) {
379: if (d == der) {
380: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
381: } else {
382: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
383: }
384: }
385: for (der2 = der + 1; der2 < dim; ++der2) {
386: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
387: for (d = 0; d < dim; ++d) {
388: if (d == der || d == der2) {
389: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
390: } else {
391: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
392: }
393: }
394: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
395: }
396: }
397: }
398: ++i;
399: }
400: } else {
401: i = 0;
402: for (o = 0; o <= sp->degree; ++o) {
403: PetscArrayzero(ind, dim);
404: while (ind[0] >= 0) {
405: LatticePoint_Internal(dim, o, ind, tup);
406: for (p = 0; p < npoints; ++p) {
407: for (der = 0; der < dim; ++der) {
408: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0;
409: for (d = 0; d < dim; ++d) {
410: if (d == der) {
411: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
412: } else {
413: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
414: }
415: }
416: for (der2 = der + 1; der2 < dim; ++der2) {
417: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
418: for (d = 0; d < dim; ++d) {
419: if (d == der || d == der2) {
420: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
421: } else {
422: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
423: }
424: }
425: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
426: }
427: }
428: }
429: ++i;
430: }
431: }
432: }
433: /* Make direct sum basis for multicomponent space */
434: for (p = 0; p < npoints; ++p) {
435: for (i = 0; i < pdim; ++i) {
436: for (c = 1; c < Nc; ++c) {
437: for (d = 0; d < dim; ++d) {
438: for (e = 0; e < dim; ++e) {
439: H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e];
440: }
441: }
442: }
443: }
444: }
445: }
446: PetscFree2(ind,tup);
447: if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
448: if (D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
449: if (B || D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
450: DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
451: DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);
452: return(0);
453: }
455: /*@
456: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
457: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
458: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
460: Input Parameters:
461: + sp - the function space object
462: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
464: Options Database:
465: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension
467: Level: intermediate
469: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
470: @*/
471: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
472: {
477: PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
478: return(0);
479: }
481: /*@
482: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
483: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
484: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
486: Input Parameters:
487: . sp - the function space object
489: Output Parameters:
490: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
492: Level: intermediate
494: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
495: @*/
496: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
497: {
503: PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
504: return(0);
505: }
507: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
508: {
509: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
512: poly->tensor = tensor;
513: return(0);
514: }
516: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
517: {
518: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
523: *tensor = poly->tensor;
524: return(0);
525: }
527: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
528: {
529: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
530: PetscInt Nc, dim, order;
531: PetscBool tensor;
532: PetscErrorCode ierr;
535: PetscSpaceGetNumComponents(sp, &Nc);
536: PetscSpaceGetNumVariables(sp, &dim);
537: PetscSpaceGetDegree(sp, &order, NULL);
538: PetscSpacePolynomialGetTensor(sp, &tensor);
539: if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
540: if (!poly->subspaces) {PetscCalloc1(dim, &poly->subspaces);}
541: if (height <= dim) {
542: if (!poly->subspaces[height-1]) {
543: PetscSpace sub;
544: const char *name;
546: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
547: PetscObjectGetName((PetscObject) sp, &name);
548: PetscObjectSetName((PetscObject) sub, name);
549: PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
550: PetscSpaceSetNumComponents(sub, Nc);
551: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
552: PetscSpaceSetNumVariables(sub, dim-height);
553: PetscSpacePolynomialSetTensor(sub, tensor);
554: PetscSpaceSetUp(sub);
555: poly->subspaces[height-1] = sub;
556: }
557: *subsp = poly->subspaces[height-1];
558: } else {
559: *subsp = NULL;
560: }
561: return(0);
562: }
564: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
565: {
569: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
570: sp->ops->setup = PetscSpaceSetUp_Polynomial;
571: sp->ops->view = PetscSpaceView_Polynomial;
572: sp->ops->destroy = PetscSpaceDestroy_Polynomial;
573: sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
574: sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
575: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
576: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
577: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
578: return(0);
579: }
581: /*MC
582: PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of
583: linear polynomials. The space is replicated for each component.
585: Level: intermediate
587: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
588: M*/
590: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
591: {
592: PetscSpace_Poly *poly;
593: PetscErrorCode ierr;
597: PetscNewLog(sp,&poly);
598: sp->data = poly;
600: poly->symmetric = PETSC_FALSE;
601: poly->tensor = PETSC_FALSE;
602: poly->degrees = NULL;
603: poly->subspaces = NULL;
605: PetscSpaceInitialize_Polynomial(sp);
606: return(0);
607: }
609: /*@
610: PetscSpacePolynomialSetSymmetric - Set whether a function space is a space of symmetric polynomials
612: Input Parameters:
613: + sp - the function space object
614: - sym - flag for symmetric polynomials
616: Options Database:
617: . -petscspace_poly_sym <bool> - Whether to use symmetric polynomials
619: Level: intermediate
621: .seealso: PetscSpacePolynomialGetSymmetric(), PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
622: @*/
623: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
624: {
625: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
629: poly->symmetric = sym;
630: return(0);
631: }
633: /*@
634: PetscSpacePolynomialGetSymmetric - Get whether a function space is a space of symmetric polynomials
636: Input Parameter:
637: . sp - the function space object
639: Output Parameter:
640: . sym - flag for symmetric polynomials
642: Level: intermediate
644: .seealso: PetscSpacePolynomialSetSymmetric(), PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
645: @*/
646: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
647: {
648: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
653: *sym = poly->symmetric;
654: return(0);
655: }