Actual source code: spacepoly.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/petscfeimpl.h>
3: PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
4: {
5: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
6: PetscErrorCode ierr;
9: PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
10: PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
11: PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
12: PetscOptionsTail();
13: return(0);
14: }
16: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
17: {
18: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
19: PetscErrorCode ierr;
22: PetscViewerASCIIPrintf(v, "%s space of degree %D\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);
23: return(0);
24: }
26: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
27: {
28: PetscBool iascii;
34: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
35: if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
36: return(0);
37: }
39: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
40: {
41: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
42: PetscInt ndegree = sp->degree+1;
43: PetscInt deg;
44: PetscErrorCode ierr;
47: if (poly->setupCalled) return(0);
48: PetscMalloc1(ndegree, &poly->degrees);
49: for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
50: if (poly->tensor) {
51: sp->maxDegree = sp->degree + PetscMax(sp->Nv - 1,0);
52: } else {
53: sp->maxDegree = sp->degree;
54: }
55: poly->setupCalled = PETSC_TRUE;
56: return(0);
57: }
59: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
60: {
61: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
62: PetscErrorCode ierr;
65: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
66: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);
67: PetscFree(poly->degrees);
68: if (poly->subspaces) {
69: PetscInt d;
71: for (d = 0; d < sp->Nv; ++d) {
72: PetscSpaceDestroy(&poly->subspaces[d]);
73: }
74: }
75: PetscFree(poly->subspaces);
76: PetscFree(poly);
77: return(0);
78: }
80: /* We treat the space as a tensor product of scalar polynomial spaces, so the dimension is multiplied by Nc */
81: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
82: {
83: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
84: PetscInt deg = sp->degree;
85: PetscInt n = sp->Nv, i;
86: PetscReal D = 1.0;
89: if (poly->tensor) {
90: *dim = 1;
91: for (i = 0; i < n; ++i) *dim *= (deg+1);
92: } else {
93: for (i = 1; i <= n; ++i) {
94: D *= ((PetscReal) (deg+i))/i;
95: }
96: *dim = (PetscInt) (D + 0.5);
97: }
98: *dim *= sp->Nc;
99: return(0);
100: }
102: /*
103: LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.
105: Input Parameters:
106: + len - The length of the tuple
107: . sum - The sum of all entries in the tuple
108: - ind - The current multi-index of the tuple, initialized to the 0 tuple
110: Output Parameter:
111: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
112: . tup - A tuple of len integers addig to sum
114: Level: developer
116: .seealso:
117: */
118: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
119: {
120: PetscInt i;
124: if (len == 1) {
125: ind[0] = -1;
126: tup[0] = sum;
127: } else if (sum == 0) {
128: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
129: } else {
130: tup[0] = sum - ind[0];
131: LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
132: if (ind[1] < 0) {
133: if (ind[0] == sum) {ind[0] = -1;}
134: else {ind[1] = 0; ++ind[0];}
135: }
136: }
137: return(0);
138: }
140: /*
141: TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.
143: Input Parameters:
144: + len - The length of the tuple
145: . max - The max for all entries in the tuple
146: - ind - The current multi-index of the tuple, initialized to the 0 tuple
148: Output Parameter:
149: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
150: . tup - A tuple of len integers less than max
152: Level: developer
154: .seealso:
155: */
156: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
157: {
158: PetscInt i;
162: if (len == 1) {
163: tup[0] = ind[0]++;
164: ind[0] = ind[0] >= max ? -1 : ind[0];
165: } else if (max == 0) {
166: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
167: } else {
168: tup[0] = ind[0];
169: TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
170: if (ind[1] < 0) {
171: ind[1] = 0;
172: if (ind[0] == max-1) {ind[0] = -1;}
173: else {++ind[0];}
174: }
175: }
176: return(0);
177: }
179: /*
180: p in [0, npoints), i in [0, pdim), c in [0, Nc)
182: B[p][i][c] = B[p][i_scalar][c][c]
183: */
184: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
185: {
186: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
187: DM dm = sp->dm;
188: PetscInt Nc = sp->Nc;
189: PetscInt ndegree = sp->degree+1;
190: PetscInt *degrees = poly->degrees;
191: PetscInt dim = sp->Nv;
192: PetscReal *lpoints, *tmp, *LB, *LD, *LH;
193: PetscInt *ind, *tup;
194: PetscInt c, pdim, d, e, der, der2, i, p, deg, o;
195: PetscErrorCode ierr;
198: PetscSpaceGetDimension(sp, &pdim);
199: pdim /= Nc;
200: DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);
201: DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
202: if (B || D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
203: if (D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
204: if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
205: for (d = 0; d < dim; ++d) {
206: for (p = 0; p < npoints; ++p) {
207: lpoints[p] = points[p*dim+d];
208: }
209: PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
210: /* LB, LD, LH (ndegree * dim x npoints) */
211: for (deg = 0; deg < ndegree; ++deg) {
212: for (p = 0; p < npoints; ++p) {
213: if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
214: if (D || H) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
215: if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
216: }
217: }
218: }
219: /* Multiply by A (pdim x ndegree * dim) */
220: PetscMalloc2(dim,&ind,dim,&tup);
221: if (B) {
222: /* B (npoints x pdim x Nc) */
223: PetscMemzero(B, npoints*pdim*Nc*Nc * sizeof(PetscReal));
224: if (poly->tensor) {
225: i = 0;
226: PetscMemzero(ind, dim * sizeof(PetscInt));
227: while (ind[0] >= 0) {
228: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
229: for (p = 0; p < npoints; ++p) {
230: B[(p*pdim + i)*Nc*Nc] = 1.0;
231: for (d = 0; d < dim; ++d) {
232: B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
233: }
234: }
235: ++i;
236: }
237: } else {
238: i = 0;
239: for (o = 0; o <= sp->degree; ++o) {
240: PetscMemzero(ind, dim * sizeof(PetscInt));
241: while (ind[0] >= 0) {
242: LatticePoint_Internal(dim, o, ind, tup);
243: for (p = 0; p < npoints; ++p) {
244: B[(p*pdim + i)*Nc*Nc] = 1.0;
245: for (d = 0; d < dim; ++d) {
246: B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
247: }
248: }
249: ++i;
250: }
251: }
252: }
253: /* Make direct sum basis for multicomponent space */
254: for (p = 0; p < npoints; ++p) {
255: for (i = 0; i < pdim; ++i) {
256: for (c = 1; c < Nc; ++c) {
257: B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
258: }
259: }
260: }
261: }
262: if (D) {
263: /* D (npoints x pdim x Nc x dim) */
264: PetscMemzero(D, npoints*pdim*Nc*Nc*dim * sizeof(PetscReal));
265: if (poly->tensor) {
266: i = 0;
267: PetscMemzero(ind, dim * sizeof(PetscInt));
268: while (ind[0] >= 0) {
269: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
270: for (p = 0; p < npoints; ++p) {
271: for (der = 0; der < dim; ++der) {
272: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
273: for (d = 0; d < dim; ++d) {
274: if (d == der) {
275: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
276: } else {
277: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
278: }
279: }
280: }
281: }
282: ++i;
283: }
284: } else {
285: i = 0;
286: for (o = 0; o <= sp->degree; ++o) {
287: PetscMemzero(ind, dim * sizeof(PetscInt));
288: while (ind[0] >= 0) {
289: LatticePoint_Internal(dim, o, ind, tup);
290: for (p = 0; p < npoints; ++p) {
291: for (der = 0; der < dim; ++der) {
292: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
293: for (d = 0; d < dim; ++d) {
294: if (d == der) {
295: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
296: } else {
297: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
298: }
299: }
300: }
301: }
302: ++i;
303: }
304: }
305: }
306: /* Make direct sum basis for multicomponent space */
307: for (p = 0; p < npoints; ++p) {
308: for (i = 0; i < pdim; ++i) {
309: for (c = 1; c < Nc; ++c) {
310: for (d = 0; d < dim; ++d) {
311: D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d];
312: }
313: }
314: }
315: }
316: }
317: if (H) {
318: /* H (npoints x pdim x Nc x Nc x dim x dim) */
319: PetscMemzero(H, npoints*pdim*Nc*Nc*dim*dim * sizeof(PetscReal));
320: if (poly->tensor) {
321: i = 0;
322: PetscMemzero(ind, dim * sizeof(PetscInt));
323: while (ind[0] >= 0) {
324: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
325: for (p = 0; p < npoints; ++p) {
326: for (der = 0; der < dim; ++der) {
327: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0;
328: for (d = 0; d < dim; ++d) {
329: if (d == der) {
330: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
331: } else {
332: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
333: }
334: }
335: for (der2 = der + 1; der2 < dim; ++der2) {
336: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
337: for (d = 0; d < dim; ++d) {
338: if (d == der || d == der2) {
339: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
340: } else {
341: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
342: }
343: }
344: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
345: }
346: }
347: }
348: ++i;
349: }
350: } else {
351: i = 0;
352: for (o = 0; o <= sp->degree; ++o) {
353: PetscMemzero(ind, dim * sizeof(PetscInt));
354: while (ind[0] >= 0) {
355: LatticePoint_Internal(dim, o, ind, tup);
356: for (p = 0; p < npoints; ++p) {
357: for (der = 0; der < dim; ++der) {
358: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0;
359: for (d = 0; d < dim; ++d) {
360: if (d == der) {
361: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
362: } else {
363: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
364: }
365: }
366: for (der2 = der + 1; der2 < dim; ++der2) {
367: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
368: for (d = 0; d < dim; ++d) {
369: if (d == der || d == der2) {
370: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
371: } else {
372: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
373: }
374: }
375: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
376: }
377: }
378: }
379: ++i;
380: }
381: }
382: }
383: /* Make direct sum basis for multicomponent space */
384: for (p = 0; p < npoints; ++p) {
385: for (i = 0; i < pdim; ++i) {
386: for (c = 1; c < Nc; ++c) {
387: for (d = 0; d < dim; ++d) {
388: for (e = 0; e < dim; ++e) {
389: H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e];
390: }
391: }
392: }
393: }
394: }
395: }
396: PetscFree2(ind,tup);
397: if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
398: if (D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
399: if (B || D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
400: DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
401: DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);
402: return(0);
403: }
405: /*@
406: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
407: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
408: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
410: Input Parameters:
411: + sp - the function space object
412: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
414: Options Database:
415: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension
417: Level: beginner
419: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
420: @*/
421: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
422: {
427: PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
428: return(0);
429: }
431: /*@
432: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
433: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
434: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
436: Input Parameters:
437: . sp - the function space object
439: Output Parameters:
440: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
442: Level: beginner
444: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
445: @*/
446: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
447: {
453: PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
454: return(0);
455: }
457: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
458: {
459: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
462: poly->tensor = tensor;
463: return(0);
464: }
466: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
467: {
468: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
473: *tensor = poly->tensor;
474: return(0);
475: }
477: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
478: {
479: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
480: PetscInt Nc, dim, order;
481: PetscBool tensor;
482: PetscErrorCode ierr;
485: PetscSpaceGetNumComponents(sp, &Nc);
486: PetscSpaceGetNumVariables(sp, &dim);
487: PetscSpaceGetDegree(sp, &order, NULL);
488: PetscSpacePolynomialGetTensor(sp, &tensor);
489: 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);}
490: if (!poly->subspaces) {PetscCalloc1(dim, &poly->subspaces);}
491: if (height <= dim) {
492: if (!poly->subspaces[height-1]) {
493: PetscSpace sub;
494: const char *name;
496: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
497: PetscObjectGetName((PetscObject) sp, &name);
498: PetscObjectSetName((PetscObject) sub, name);
499: PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
500: PetscSpaceSetNumComponents(sub, Nc);
501: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
502: PetscSpaceSetNumVariables(sub, dim-height);
503: PetscSpacePolynomialSetTensor(sub, tensor);
504: PetscSpaceSetUp(sub);
505: poly->subspaces[height-1] = sub;
506: }
507: *subsp = poly->subspaces[height-1];
508: } else {
509: *subsp = NULL;
510: }
511: return(0);
512: }
514: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
515: {
519: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
520: sp->ops->setup = PetscSpaceSetUp_Polynomial;
521: sp->ops->view = PetscSpaceView_Polynomial;
522: sp->ops->destroy = PetscSpaceDestroy_Polynomial;
523: sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
524: sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
525: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
526: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
527: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
528: return(0);
529: }
531: /*MC
532: PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of
533: linear polynomials. The space is replicated for each component.
535: Level: intermediate
537: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
538: M*/
540: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
541: {
542: PetscSpace_Poly *poly;
543: PetscErrorCode ierr;
547: PetscNewLog(sp,&poly);
548: sp->data = poly;
550: poly->symmetric = PETSC_FALSE;
551: poly->tensor = PETSC_FALSE;
552: poly->degrees = NULL;
553: poly->subspaces = NULL;
555: PetscSpaceInitialize_Polynomial(sp);
556: return(0);
557: }
559: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
560: {
561: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
565: poly->symmetric = sym;
566: return(0);
567: }
569: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
570: {
571: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
576: *sym = poly->symmetric;
577: return(0);
578: }