Actual source code: spacepoly.c
petsc-3.10.5 2019-03-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 viewer)
17: {
18: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
19: PetscErrorCode ierr;
22: if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space of degree %D\n", sp->degree);}
23: else {PetscViewerASCIIPrintf(viewer, "Polynomial space of degree %D\n", sp->degree);}
24: return(0);
25: }
27: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
28: {
29: PetscBool iascii;
35: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
36: if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
37: return(0);
38: }
40: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
41: {
42: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
43: PetscInt ndegree = sp->degree+1;
44: PetscInt deg;
45: PetscErrorCode ierr;
48: if (poly->setupCalled) return(0);
49: PetscMalloc1(ndegree, &poly->degrees);
50: for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
51: if (poly->tensor) {
52: sp->maxDegree = sp->degree + PetscMax(sp->Nv - 1,0);
53: } else {
54: sp->maxDegree = sp->degree;
55: }
56: poly->setupCalled = PETSC_TRUE;
57: return(0);
58: }
60: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
61: {
62: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
63: PetscErrorCode ierr;
66: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
67: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);
68: PetscFree(poly->degrees);
69: if (poly->subspaces) {
70: PetscInt d;
72: for (d = 0; d < sp->Nv; ++d) {
73: PetscSpaceDestroy(&poly->subspaces[d]);
74: }
75: }
76: PetscFree(poly->subspaces);
77: PetscFree(poly);
78: return(0);
79: }
81: /* We treat the space as a tensor product of scalar polynomial spaces, so the dimension is multiplied by Nc */
82: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
83: {
84: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
85: PetscInt deg = sp->degree;
86: PetscInt n = sp->Nv, i;
87: PetscReal D = 1.0;
90: if (poly->tensor) {
91: *dim = 1;
92: for (i = 0; i < n; ++i) *dim *= (deg+1);
93: } else {
94: for (i = 1; i <= n; ++i) {
95: D *= ((PetscReal) (deg+i))/i;
96: }
97: *dim = (PetscInt) (D + 0.5);
98: }
99: *dim *= sp->Nc;
100: return(0);
101: }
103: /*
104: LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.
106: Input Parameters:
107: + len - The length of the tuple
108: . sum - The sum of all entries in the tuple
109: - ind - The current multi-index of the tuple, initialized to the 0 tuple
111: Output Parameter:
112: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
113: . tup - A tuple of len integers addig to sum
115: Level: developer
117: .seealso:
118: */
119: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
120: {
121: PetscInt i;
125: if (len == 1) {
126: ind[0] = -1;
127: tup[0] = sum;
128: } else if (sum == 0) {
129: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
130: } else {
131: tup[0] = sum - ind[0];
132: LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
133: if (ind[1] < 0) {
134: if (ind[0] == sum) {ind[0] = -1;}
135: else {ind[1] = 0; ++ind[0];}
136: }
137: }
138: return(0);
139: }
141: /*
142: TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.
144: Input Parameters:
145: + len - The length of the tuple
146: . max - The max for all entries in the tuple
147: - ind - The current multi-index of the tuple, initialized to the 0 tuple
149: Output Parameter:
150: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
151: . tup - A tuple of len integers less than max
153: Level: developer
155: .seealso:
156: */
157: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
158: {
159: PetscInt i;
163: if (len == 1) {
164: tup[0] = ind[0]++;
165: ind[0] = ind[0] >= max ? -1 : ind[0];
166: } else if (max == 0) {
167: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
168: } else {
169: tup[0] = ind[0];
170: TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
171: if (ind[1] < 0) {
172: ind[1] = 0;
173: if (ind[0] == max-1) {ind[0] = -1;}
174: else {++ind[0];}
175: }
176: }
177: return(0);
178: }
180: /*
181: p in [0, npoints), i in [0, pdim), c in [0, Nc)
183: B[p][i][c] = B[p][i_scalar][c][c]
184: */
185: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
186: {
187: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
188: DM dm = sp->dm;
189: PetscInt Nc = sp->Nc;
190: PetscInt ndegree = sp->degree+1;
191: PetscInt *degrees = poly->degrees;
192: PetscInt dim = sp->Nv;
193: PetscReal *lpoints, *tmp, *LB, *LD, *LH;
194: PetscInt *ind, *tup;
195: PetscInt c, pdim, d, e, der, der2, i, p, deg, o;
196: PetscErrorCode ierr;
199: PetscSpaceGetDimension(sp, &pdim);
200: pdim /= Nc;
201: DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);
202: DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
203: if (B || D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
204: if (D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
205: if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
206: for (d = 0; d < dim; ++d) {
207: for (p = 0; p < npoints; ++p) {
208: lpoints[p] = points[p*dim+d];
209: }
210: PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
211: /* LB, LD, LH (ndegree * dim x npoints) */
212: for (deg = 0; deg < ndegree; ++deg) {
213: for (p = 0; p < npoints; ++p) {
214: if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
215: if (D || H) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
216: if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
217: }
218: }
219: }
220: /* Multiply by A (pdim x ndegree * dim) */
221: PetscMalloc2(dim,&ind,dim,&tup);
222: if (B) {
223: /* B (npoints x pdim x Nc) */
224: PetscMemzero(B, npoints*pdim*Nc*Nc * sizeof(PetscReal));
225: if (poly->tensor) {
226: i = 0;
227: PetscMemzero(ind, dim * sizeof(PetscInt));
228: while (ind[0] >= 0) {
229: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
230: for (p = 0; p < npoints; ++p) {
231: B[(p*pdim + i)*Nc*Nc] = 1.0;
232: for (d = 0; d < dim; ++d) {
233: B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
234: }
235: }
236: ++i;
237: }
238: } else {
239: i = 0;
240: for (o = 0; o <= sp->degree; ++o) {
241: PetscMemzero(ind, dim * sizeof(PetscInt));
242: while (ind[0] >= 0) {
243: LatticePoint_Internal(dim, o, 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: }
253: }
254: /* Make direct sum basis for multicomponent space */
255: for (p = 0; p < npoints; ++p) {
256: for (i = 0; i < pdim; ++i) {
257: for (c = 1; c < Nc; ++c) {
258: B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
259: }
260: }
261: }
262: }
263: if (D) {
264: /* D (npoints x pdim x Nc x dim) */
265: PetscMemzero(D, npoints*pdim*Nc*Nc*dim * sizeof(PetscReal));
266: if (poly->tensor) {
267: i = 0;
268: PetscMemzero(ind, dim * sizeof(PetscInt));
269: while (ind[0] >= 0) {
270: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
271: for (p = 0; p < npoints; ++p) {
272: for (der = 0; der < dim; ++der) {
273: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
274: for (d = 0; d < dim; ++d) {
275: if (d == der) {
276: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
277: } else {
278: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
279: }
280: }
281: }
282: }
283: ++i;
284: }
285: } else {
286: i = 0;
287: for (o = 0; o <= sp->degree; ++o) {
288: PetscMemzero(ind, dim * sizeof(PetscInt));
289: while (ind[0] >= 0) {
290: LatticePoint_Internal(dim, o, ind, tup);
291: for (p = 0; p < npoints; ++p) {
292: for (der = 0; der < dim; ++der) {
293: D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
294: for (d = 0; d < dim; ++d) {
295: if (d == der) {
296: D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
297: } else {
298: D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
299: }
300: }
301: }
302: }
303: ++i;
304: }
305: }
306: }
307: /* Make direct sum basis for multicomponent space */
308: for (p = 0; p < npoints; ++p) {
309: for (i = 0; i < pdim; ++i) {
310: for (c = 1; c < Nc; ++c) {
311: for (d = 0; d < dim; ++d) {
312: D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d];
313: }
314: }
315: }
316: }
317: }
318: if (H) {
319: /* H (npoints x pdim x Nc x Nc x dim x dim) */
320: PetscMemzero(H, npoints*pdim*Nc*Nc*dim*dim * sizeof(PetscReal));
321: if (poly->tensor) {
322: i = 0;
323: PetscMemzero(ind, dim * sizeof(PetscInt));
324: while (ind[0] >= 0) {
325: TensorPoint_Internal(dim, sp->degree+1, ind, tup);
326: for (p = 0; p < npoints; ++p) {
327: for (der = 0; der < dim; ++der) {
328: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0;
329: for (d = 0; d < dim; ++d) {
330: if (d == der) {
331: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
332: } else {
333: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
334: }
335: }
336: for (der2 = der + 1; der2 < dim; ++der2) {
337: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
338: for (d = 0; d < dim; ++d) {
339: if (d == der || d == der2) {
340: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
341: } else {
342: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
343: }
344: }
345: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
346: }
347: }
348: }
349: ++i;
350: }
351: } else {
352: i = 0;
353: for (o = 0; o <= sp->degree; ++o) {
354: PetscMemzero(ind, dim * sizeof(PetscInt));
355: while (ind[0] >= 0) {
356: LatticePoint_Internal(dim, o, ind, tup);
357: for (p = 0; p < npoints; ++p) {
358: for (der = 0; der < dim; ++der) {
359: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0;
360: for (d = 0; d < dim; ++d) {
361: if (d == der) {
362: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
363: } else {
364: H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
365: }
366: }
367: for (der2 = der + 1; der2 < dim; ++der2) {
368: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
369: for (d = 0; d < dim; ++d) {
370: if (d == der || d == der2) {
371: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
372: } else {
373: H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
374: }
375: }
376: H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
377: }
378: }
379: }
380: ++i;
381: }
382: }
383: }
384: /* Make direct sum basis for multicomponent space */
385: for (p = 0; p < npoints; ++p) {
386: for (i = 0; i < pdim; ++i) {
387: for (c = 1; c < Nc; ++c) {
388: for (d = 0; d < dim; ++d) {
389: for (e = 0; e < dim; ++e) {
390: H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e];
391: }
392: }
393: }
394: }
395: }
396: }
397: PetscFree2(ind,tup);
398: if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
399: if (D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
400: if (B || D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
401: DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
402: DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);
403: return(0);
404: }
406: /*@
407: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
408: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
409: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
411: Input Parameters:
412: + sp - the function space object
413: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
415: Level: beginner
417: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
418: @*/
419: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
420: {
425: PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
426: return(0);
427: }
429: /*@
430: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
431: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
432: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
434: Input Parameters:
435: . sp - the function space object
437: Output Parameters:
438: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
440: Level: beginner
442: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
443: @*/
444: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
445: {
451: PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
452: return(0);
453: }
455: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
456: {
457: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
460: poly->tensor = tensor;
461: return(0);
462: }
464: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
465: {
466: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
471: *tensor = poly->tensor;
472: return(0);
473: }
475: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
476: {
477: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
478: PetscInt Nc, dim, order;
479: PetscBool tensor;
480: PetscErrorCode ierr;
483: PetscSpaceGetNumComponents(sp, &Nc);
484: PetscSpaceGetNumVariables(sp, &dim);
485: PetscSpaceGetDegree(sp, &order, NULL);
486: PetscSpacePolynomialGetTensor(sp, &tensor);
487: 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);}
488: if (!poly->subspaces) {PetscCalloc1(dim, &poly->subspaces);}
489: if (height <= dim) {
490: if (!poly->subspaces[height-1]) {
491: PetscSpace sub;
493: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
494: PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
495: PetscSpaceSetNumComponents(sub, Nc);
496: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
497: PetscSpaceSetNumVariables(sub, dim-height);
498: PetscSpacePolynomialSetTensor(sub, tensor);
499: PetscSpaceSetUp(sub);
500: poly->subspaces[height-1] = sub;
501: }
502: *subsp = poly->subspaces[height-1];
503: } else {
504: *subsp = NULL;
505: }
506: return(0);
507: }
509: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
510: {
514: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
515: sp->ops->setup = PetscSpaceSetUp_Polynomial;
516: sp->ops->view = PetscSpaceView_Polynomial;
517: sp->ops->destroy = PetscSpaceDestroy_Polynomial;
518: sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
519: sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
520: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
521: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
522: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
523: return(0);
524: }
526: /*MC
527: PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of
528: linear polynomials. The space is replicated for each component.
530: Level: intermediate
532: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
533: M*/
535: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
536: {
537: PetscSpace_Poly *poly;
538: PetscErrorCode ierr;
542: PetscNewLog(sp,&poly);
543: sp->data = poly;
545: poly->symmetric = PETSC_FALSE;
546: poly->tensor = PETSC_FALSE;
547: poly->degrees = NULL;
548: poly->subspaces = NULL;
550: PetscSpaceInitialize_Polynomial(sp);
551: return(0);
552: }
554: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
555: {
556: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
560: poly->symmetric = sym;
561: return(0);
562: }
564: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
565: {
566: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
571: *sym = poly->symmetric;
572: return(0);
573: }