Actual source code: dt.c
1: /* Discretization tools */
3: #include <petscdt.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/dtimpl.h>
7: #include <petscviewer.h>
8: #include <petscdmplex.h>
9: #include <petscdmshell.h>
11: #if defined(PETSC_HAVE_MPFR)
12: #include <mpfr.h>
13: #endif
15: const char *const PetscDTNodeTypes[] = {"gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
17: static PetscBool GolubWelschCite = PETSC_FALSE;
18: const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
19: " author = {Golub and Welsch},\n"
20: " title = {Calculation of Quadrature Rules},\n"
21: " journal = {Math. Comp.},\n"
22: " volume = {23},\n"
23: " number = {106},\n"
24: " pages = {221--230},\n"
25: " year = {1969}\n}\n";
27: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
28: quadrature rules:
30: - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
31: - in single precision, Newton's method starts producing incorrect roots around n = 15, but
32: the weights from Golub & Welsch become a problem before then: they produces errors
33: in computing the Jacobi-polynomial Gram matrix around n = 6.
35: So we default to Newton's method (required fewer dependencies) */
36: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
38: PetscClassId PETSCQUADRATURE_CLASSID = 0;
40: /*@
41: PetscQuadratureCreate - Create a PetscQuadrature object
43: Collective
45: Input Parameter:
46: . comm - The communicator for the PetscQuadrature object
48: Output Parameter:
49: . q - The PetscQuadrature object
51: Level: beginner
53: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
54: @*/
55: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
56: {
58: DMInitializePackage();
59: PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
60: (*q)->dim = -1;
61: (*q)->Nc = 1;
62: (*q)->order = -1;
63: (*q)->numPoints = 0;
64: (*q)->points = NULL;
65: (*q)->weights = NULL;
66: return 0;
67: }
69: /*@
70: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
72: Collective on q
74: Input Parameter:
75: . q - The PetscQuadrature object
77: Output Parameter:
78: . r - The new PetscQuadrature object
80: Level: beginner
82: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
83: @*/
84: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
85: {
86: PetscInt order, dim, Nc, Nq;
87: const PetscReal *points, *weights;
88: PetscReal *p, *w;
91: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
92: PetscQuadratureGetOrder(q, &order);
93: PetscQuadratureSetOrder(*r, order);
94: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
95: PetscMalloc1(Nq*dim, &p);
96: PetscMalloc1(Nq*Nc, &w);
97: PetscArraycpy(p, points, Nq*dim);
98: PetscArraycpy(w, weights, Nc * Nq);
99: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
100: return 0;
101: }
103: /*@
104: PetscQuadratureDestroy - Destroys a PetscQuadrature object
106: Collective on q
108: Input Parameter:
109: . q - The PetscQuadrature object
111: Level: beginner
113: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
114: @*/
115: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
116: {
117: if (!*q) return 0;
119: if (--((PetscObject)(*q))->refct > 0) {
120: *q = NULL;
121: return 0;
122: }
123: PetscFree((*q)->points);
124: PetscFree((*q)->weights);
125: PetscHeaderDestroy(q);
126: return 0;
127: }
129: /*@
130: PetscQuadratureGetOrder - Return the order of the method
132: Not collective
134: Input Parameter:
135: . q - The PetscQuadrature object
137: Output Parameter:
138: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
140: Level: intermediate
142: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
143: @*/
144: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
145: {
148: *order = q->order;
149: return 0;
150: }
152: /*@
153: PetscQuadratureSetOrder - Return the order of the method
155: Not collective
157: Input Parameters:
158: + q - The PetscQuadrature object
159: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
161: Level: intermediate
163: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
164: @*/
165: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
166: {
168: q->order = order;
169: return 0;
170: }
172: /*@
173: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
175: Not collective
177: Input Parameter:
178: . q - The PetscQuadrature object
180: Output Parameter:
181: . Nc - The number of components
183: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
185: Level: intermediate
187: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
188: @*/
189: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
190: {
193: *Nc = q->Nc;
194: return 0;
195: }
197: /*@
198: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
200: Not collective
202: Input Parameters:
203: + q - The PetscQuadrature object
204: - Nc - The number of components
206: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
208: Level: intermediate
210: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
211: @*/
212: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
213: {
215: q->Nc = Nc;
216: return 0;
217: }
219: /*@C
220: PetscQuadratureGetData - Returns the data defining the quadrature
222: Not collective
224: Input Parameter:
225: . q - The PetscQuadrature object
227: Output Parameters:
228: + dim - The spatial dimension
229: . Nc - The number of components
230: . npoints - The number of quadrature points
231: . points - The coordinates of each quadrature point
232: - weights - The weight of each quadrature point
234: Level: intermediate
236: Fortran Notes:
237: From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
239: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
240: @*/
241: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
242: {
244: if (dim) {
246: *dim = q->dim;
247: }
248: if (Nc) {
250: *Nc = q->Nc;
251: }
252: if (npoints) {
254: *npoints = q->numPoints;
255: }
256: if (points) {
258: *points = q->points;
259: }
260: if (weights) {
262: *weights = q->weights;
263: }
264: return 0;
265: }
267: /*@
268: PetscQuadratureEqual - determine whether two quadratures are equivalent
270: Input Parameters:
271: + A - A PetscQuadrature object
272: - B - Another PetscQuadrature object
274: Output Parameters:
275: . equal - PETSC_TRUE if the quadratures are the same
277: Level: intermediate
279: .seealso: PetscQuadratureCreate()
280: @*/
281: PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal)
282: {
286: *equal = PETSC_FALSE;
287: if (A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) {
288: return 0;
289: }
290: for (PetscInt i=0; i<A->numPoints*A->dim; i++) {
291: if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) {
292: return 0;
293: }
294: }
295: if (!A->weights && !B->weights) {
296: *equal = PETSC_TRUE;
297: return 0;
298: }
299: if (A->weights && B->weights) {
300: for (PetscInt i=0; i<A->numPoints; i++) {
301: if (PetscAbsReal(A->weights[i] - B->weights[i]) > PETSC_SMALL) {
302: return 0;
303: }
304: }
305: *equal = PETSC_TRUE;
306: }
307: return 0;
308: }
310: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
311: {
312: PetscScalar *Js, *Jinvs;
313: PetscInt i, j, k;
314: PetscBLASInt bm, bn, info;
316: if (!m || !n) return 0;
317: PetscBLASIntCast(m, &bm);
318: PetscBLASIntCast(n, &bn);
319: #if defined(PETSC_USE_COMPLEX)
320: PetscMalloc2(m*n, &Js, m*n, &Jinvs);
321: for (i = 0; i < m*n; i++) Js[i] = J[i];
322: #else
323: Js = (PetscReal *) J;
324: Jinvs = Jinv;
325: #endif
326: if (m == n) {
327: PetscBLASInt *pivots;
328: PetscScalar *W;
330: PetscMalloc2(m, &pivots, m, &W);
332: PetscArraycpy(Jinvs, Js, m * m);
333: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
335: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
337: PetscFree2(pivots, W);
338: } else if (m < n) {
339: PetscScalar *JJT;
340: PetscBLASInt *pivots;
341: PetscScalar *W;
343: PetscMalloc1(m*m, &JJT);
344: PetscMalloc2(m, &pivots, m, &W);
345: for (i = 0; i < m; i++) {
346: for (j = 0; j < m; j++) {
347: PetscScalar val = 0.;
349: for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
350: JJT[i * m + j] = val;
351: }
352: }
354: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
356: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
358: for (i = 0; i < n; i++) {
359: for (j = 0; j < m; j++) {
360: PetscScalar val = 0.;
362: for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
363: Jinvs[i * m + j] = val;
364: }
365: }
366: PetscFree2(pivots, W);
367: PetscFree(JJT);
368: } else {
369: PetscScalar *JTJ;
370: PetscBLASInt *pivots;
371: PetscScalar *W;
373: PetscMalloc1(n*n, &JTJ);
374: PetscMalloc2(n, &pivots, n, &W);
375: for (i = 0; i < n; i++) {
376: for (j = 0; j < n; j++) {
377: PetscScalar val = 0.;
379: for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
380: JTJ[i * n + j] = val;
381: }
382: }
384: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
386: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
388: for (i = 0; i < n; i++) {
389: for (j = 0; j < m; j++) {
390: PetscScalar val = 0.;
392: for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
393: Jinvs[i * m + j] = val;
394: }
395: }
396: PetscFree2(pivots, W);
397: PetscFree(JTJ);
398: }
399: #if defined(PETSC_USE_COMPLEX)
400: for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
401: PetscFree2(Js, Jinvs);
402: #endif
403: return 0;
404: }
406: /*@
407: PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
409: Collecive on PetscQuadrature
411: Input Parameters:
412: + q - the quadrature functional
413: . imageDim - the dimension of the image of the transformation
414: . origin - a point in the original space
415: . originImage - the image of the origin under the transformation
416: . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
417: - formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose formDegree), it is assumed that they represent multiple k-forms) [see PetscDTAltVPullback() for interpretation of formDegree]
419: Output Parameters:
420: . Jinvstarq - a quadrature rule where each point is the image of a point in the original quadrature rule, and where the k-form weights have been pulled-back by the pseudoinverse of J to the k-form weights in the image space.
422: Note: the new quadrature rule will have a different number of components if spaces have different dimensions. For example, pushing a 2-form forward from a two dimensional space to a three dimensional space changes the number of components from 1 to 3.
424: Level: intermediate
426: .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
427: @*/
428: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
429: {
430: PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
431: const PetscReal *points;
432: const PetscReal *weights;
433: PetscReal *imagePoints, *imageWeights;
434: PetscReal *Jinv;
435: PetscReal *Jinvstar;
439: PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);
440: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);
442: Ncopies = Nc / formSize;
443: PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);
444: imageNc = Ncopies * imageFormSize;
445: PetscMalloc1(Npoints * imageDim, &imagePoints);
446: PetscMalloc1(Npoints * imageNc, &imageWeights);
447: PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);
448: PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);
449: PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);
450: for (pt = 0; pt < Npoints; pt++) {
451: const PetscReal *point = &points[pt * dim];
452: PetscReal *imagePoint = &imagePoints[pt * imageDim];
454: for (i = 0; i < imageDim; i++) {
455: PetscReal val = originImage[i];
457: for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
458: imagePoint[i] = val;
459: }
460: for (c = 0; c < Ncopies; c++) {
461: const PetscReal *form = &weights[pt * Nc + c * formSize];
462: PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
464: for (i = 0; i < imageFormSize; i++) {
465: PetscReal val = 0.;
467: for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
468: imageForm[i] = val;
469: }
470: }
471: }
472: PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);
473: PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);
474: PetscFree2(Jinv, Jinvstar);
475: return 0;
476: }
478: /*@C
479: PetscQuadratureSetData - Sets the data defining the quadrature
481: Not collective
483: Input Parameters:
484: + q - The PetscQuadrature object
485: . dim - The spatial dimension
486: . Nc - The number of components
487: . npoints - The number of quadrature points
488: . points - The coordinates of each quadrature point
489: - weights - The weight of each quadrature point
491: Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
493: Level: intermediate
495: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
496: @*/
497: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
498: {
500: if (dim >= 0) q->dim = dim;
501: if (Nc >= 0) q->Nc = Nc;
502: if (npoints >= 0) q->numPoints = npoints;
503: if (points) {
505: q->points = points;
506: }
507: if (weights) {
509: q->weights = weights;
510: }
511: return 0;
512: }
514: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
515: {
516: PetscInt q, d, c;
517: PetscViewerFormat format;
519: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);
520: else PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);
521: PetscViewerGetFormat(v, &format);
522: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
523: for (q = 0; q < quad->numPoints; ++q) {
524: PetscViewerASCIIPrintf(v, "p%D (", q);
525: PetscViewerASCIIUseTabs(v, PETSC_FALSE);
526: for (d = 0; d < quad->dim; ++d) {
527: if (d) PetscViewerASCIIPrintf(v, ", ");
528: PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
529: }
530: PetscViewerASCIIPrintf(v, ") ");
531: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, "w%D (", q);
532: for (c = 0; c < quad->Nc; ++c) {
533: if (c) PetscViewerASCIIPrintf(v, ", ");
534: PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
535: }
536: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, ")");
537: PetscViewerASCIIPrintf(v, "\n");
538: PetscViewerASCIIUseTabs(v, PETSC_TRUE);
539: }
540: return 0;
541: }
543: /*@C
544: PetscQuadratureView - Views a PetscQuadrature object
546: Collective on quad
548: Input Parameters:
549: + quad - The PetscQuadrature object
550: - viewer - The PetscViewer object
552: Level: beginner
554: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
555: @*/
556: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
557: {
558: PetscBool iascii;
562: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);
563: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
564: PetscViewerASCIIPushTab(viewer);
565: if (iascii) PetscQuadratureView_Ascii(quad, viewer);
566: PetscViewerASCIIPopTab(viewer);
567: return 0;
568: }
570: /*@C
571: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
573: Not collective
575: Input Parameters:
576: + q - The original PetscQuadrature
577: . numSubelements - The number of subelements the original element is divided into
578: . v0 - An array of the initial points for each subelement
579: - jac - An array of the Jacobian mappings from the reference to each subelement
581: Output Parameters:
582: . dim - The dimension
584: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
586: Not available from Fortran
588: Level: intermediate
590: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
591: @*/
592: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
593: {
594: const PetscReal *points, *weights;
595: PetscReal *pointsRef, *weightsRef;
596: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
602: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
603: PetscQuadratureGetOrder(q, &order);
604: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
605: npointsRef = npoints*numSubelements;
606: PetscMalloc1(npointsRef*dim,&pointsRef);
607: PetscMalloc1(npointsRef*Nc, &weightsRef);
608: for (c = 0; c < numSubelements; ++c) {
609: for (p = 0; p < npoints; ++p) {
610: for (d = 0; d < dim; ++d) {
611: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
612: for (e = 0; e < dim; ++e) {
613: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
614: }
615: }
616: /* Could also use detJ here */
617: for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
618: }
619: }
620: PetscQuadratureSetOrder(*qref, order);
621: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
622: return 0;
623: }
625: /* Compute the coefficients for the Jacobi polynomial recurrence,
626: *
627: * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
628: */
629: #define PetscDTJacobiRecurrence_Internal(n,a,b,cnm1,cnm1x,cnm2) \
630: do { \
631: PetscReal _a = (a); \
632: PetscReal _b = (b); \
633: PetscReal _n = (n); \
634: if (n == 1) { \
635: (cnm1) = (_a-_b) * 0.5; \
636: (cnm1x) = (_a+_b+2.)*0.5; \
637: (cnm2) = 0.; \
638: } else { \
639: PetscReal _2n = _n+_n; \
640: PetscReal _d = (_2n*(_n+_a+_b)*(_2n+_a+_b-2)); \
641: PetscReal _n1 = (_2n+_a+_b-1.)*(_a*_a-_b*_b); \
642: PetscReal _n1x = (_2n+_a+_b-1.)*(_2n+_a+_b)*(_2n+_a+_b-2); \
643: PetscReal _n2 = 2.*((_n+_a-1.)*(_n+_b-1.)*(_2n+_a+_b)); \
644: (cnm1) = _n1 / _d; \
645: (cnm1x) = _n1x / _d; \
646: (cnm2) = _n2 / _d; \
647: } \
648: } while (0)
650: /*@
651: PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
653: $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
655: Input Parameters:
656: - alpha - the left exponent > -1
657: . beta - the right exponent > -1
658: + n - the polynomial degree
660: Output Parameter:
661: . norm - the weighted L2 norm
663: Level: beginner
665: .seealso: PetscDTJacobiEval()
666: @*/
667: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
668: {
669: PetscReal twoab1;
670: PetscReal gr;
675: twoab1 = PetscPowReal(2., alpha + beta + 1.);
676: #if defined(PETSC_HAVE_LGAMMA)
677: if (!n) {
678: gr = PetscExpReal(PetscLGamma(alpha+1.) + PetscLGamma(beta+1.) - PetscLGamma(alpha+beta+2.));
679: } else {
680: gr = PetscExpReal(PetscLGamma(n+alpha+1.) + PetscLGamma(n+beta+1.) - (PetscLGamma(n+1.) + PetscLGamma(n+alpha+beta+1.))) / (n+n+alpha+beta+1.);
681: }
682: #else
683: {
684: PetscInt alphai = (PetscInt) alpha;
685: PetscInt betai = (PetscInt) beta;
686: PetscInt i;
688: gr = n ? (1. / (n+n+alpha+beta+1.)) : 1.;
689: if ((PetscReal) alphai == alpha) {
690: if (!n) {
691: for (i = 0; i < alphai; i++) gr *= (i+1.) / (beta+i+1.);
692: gr /= (alpha+beta+1.);
693: } else {
694: for (i = 0; i < alphai; i++) gr *= (n+i+1.) / (n+beta+i+1.);
695: }
696: } else if ((PetscReal) betai == beta) {
697: if (!n) {
698: for (i = 0; i < betai; i++) gr *= (i+1.) / (alpha+i+2.);
699: gr /= (alpha+beta+1.);
700: } else {
701: for (i = 0; i < betai; i++) gr *= (n+i+1.) / (n+alpha+i+1.);
702: }
703: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
704: }
705: #endif
706: *norm = PetscSqrtReal(twoab1 * gr);
707: return 0;
708: }
710: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
711: {
712: PetscReal ak, bk;
713: PetscReal abk1;
714: PetscInt i,l,maxdegree;
716: maxdegree = degrees[ndegree-1] - k;
717: ak = a + k;
718: bk = b + k;
719: abk1 = a + b + k + 1.;
720: if (maxdegree < 0) {
721: for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.;
722: return 0;
723: }
724: for (i=0; i<npoints; i++) {
725: PetscReal pm1,pm2,x;
726: PetscReal cnm1, cnm1x, cnm2;
727: PetscInt j,m;
729: x = points[i];
730: pm2 = 1.;
731: PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2);
732: pm1 = (cnm1 + cnm1x*x);
733: l = 0;
734: while (l < ndegree && degrees[l] - k < 0) {
735: p[l++] = 0.;
736: }
737: while (l < ndegree && degrees[l] - k == 0) {
738: p[l] = pm2;
739: for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
740: l++;
741: }
742: while (l < ndegree && degrees[l] - k == 1) {
743: p[l] = pm1;
744: for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
745: l++;
746: }
747: for (j=2; j<=maxdegree; j++) {
748: PetscReal pp;
750: PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2);
751: pp = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2;
752: pm2 = pm1;
753: pm1 = pp;
754: while (l < ndegree && degrees[l] - k == j) {
755: p[l] = pp;
756: for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
757: l++;
758: }
759: }
760: p += ndegree;
761: }
762: return 0;
763: }
765: /*@
766: PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree. The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta) f(x) g(x) dx$.
768: Input Parameters:
769: + alpha - the left exponent of the weight
770: . beta - the right exponetn of the weight
771: . npoints - the number of points to evaluate the polynomials at
772: . points - [npoints] array of point coordinates
773: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
774: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
776: Output Parameters:
777: - p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
778: (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
779: (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
780: varying) dimension is the index of the evaluation point.
782: Level: advanced
784: .seealso: PetscDTJacobiEval(), PetscDTPKDEvalJet()
785: @*/
786: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
787: {
788: PetscInt i, j, l;
789: PetscInt *degrees;
790: PetscReal *psingle;
792: if (degree == 0) {
793: PetscInt zero = 0;
795: for (i = 0; i <= k; i++) {
796: PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i*npoints]);
797: }
798: return 0;
799: }
800: PetscMalloc1(degree + 1, °rees);
801: PetscMalloc1((degree + 1) * npoints, &psingle);
802: for (i = 0; i <= degree; i++) degrees[i] = i;
803: for (i = 0; i <= k; i++) {
804: PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle);
805: for (j = 0; j <= degree; j++) {
806: for (l = 0; l < npoints; l++) {
807: p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
808: }
809: }
810: }
811: PetscFree(psingle);
812: PetscFree(degrees);
813: return 0;
814: }
816: /*@
817: PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$
818: at points
820: Not Collective
822: Input Parameters:
823: + npoints - number of spatial points to evaluate at
824: . alpha - the left exponent > -1
825: . beta - the right exponent > -1
826: . points - array of locations to evaluate at
827: . ndegree - number of basis degrees to evaluate
828: - degrees - sorted array of degrees to evaluate
830: Output Parameters:
831: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
832: . D - row-oriented derivative evaluation matrix (or NULL)
833: - D2 - row-oriented second derivative evaluation matrix (or NULL)
835: Level: intermediate
837: .seealso: PetscDTGaussQuadrature()
838: @*/
839: PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
840: {
843: if (!npoints || !ndegree) return 0;
844: if (B) PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);
845: if (D) PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);
846: if (D2) PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);
847: return 0;
848: }
850: /*@
851: PetscDTLegendreEval - evaluate Legendre polynomials at points
853: Not Collective
855: Input Parameters:
856: + npoints - number of spatial points to evaluate at
857: . points - array of locations to evaluate at
858: . ndegree - number of basis degrees to evaluate
859: - degrees - sorted array of degrees to evaluate
861: Output Parameters:
862: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
863: . D - row-oriented derivative evaluation matrix (or NULL)
864: - D2 - row-oriented second derivative evaluation matrix (or NULL)
866: Level: intermediate
868: .seealso: PetscDTGaussQuadrature()
869: @*/
870: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
871: {
872: PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);
873: return 0;
874: }
876: /*@
877: PetscDTIndexToGradedOrder - convert an index into a tuple of monomial degrees in a graded order (that is, if the degree sum of tuple x is less than the degree sum of tuple y, then the index of x is smaller than the index of y)
879: Input Parameters:
880: + len - the desired length of the degree tuple
881: - index - the index to convert: should be >= 0
883: Output Parameter:
884: . degtup - will be filled with a tuple of degrees
886: Level: beginner
888: Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
889: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
890: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
892: .seealso: PetscDTGradedOrderToIndex()
893: @*/
894: PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
895: {
896: PetscInt i, total;
897: PetscInt sum;
902: total = 1;
903: sum = 0;
904: while (index >= total) {
905: index -= total;
906: total = (total * (len + sum)) / (sum + 1);
907: sum++;
908: }
909: for (i = 0; i < len; i++) {
910: PetscInt c;
912: degtup[i] = sum;
913: for (c = 0, total = 1; c < sum; c++) {
914: /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
915: if (index < total) break;
916: index -= total;
917: total = (total * (len - 1 - i + c)) / (c + 1);
918: degtup[i]--;
919: }
920: sum -= degtup[i];
921: }
922: return 0;
923: }
925: /*@
926: PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of PetscDTIndexToGradedOrder().
928: Input Parameters:
929: + len - the length of the degree tuple
930: - degtup - tuple with this length
932: Output Parameter:
933: . index - index in graded order: >= 0
935: Level: Beginner
937: Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
938: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
939: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
941: .seealso: PetscDTIndexToGradedOrder()
942: @*/
943: PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
944: {
945: PetscInt i, idx, sum, total;
949: for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
950: idx = 0;
951: total = 1;
952: for (i = 0; i < sum; i++) {
953: idx += total;
954: total = (total * (len + i)) / (i + 1);
955: }
956: for (i = 0; i < len - 1; i++) {
957: PetscInt c;
959: total = 1;
960: sum -= degtup[i];
961: for (c = 0; c < sum; c++) {
962: idx += total;
963: total = (total * (len - 1 - i + c)) / (c + 1);
964: }
965: }
966: *index = idx;
967: return 0;
968: }
970: static PetscBool PKDCite = PETSC_FALSE;
971: const char PKDCitation[] = "@article{Kirby2010,\n"
972: " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
973: " author={Kirby, Robert C},\n"
974: " journal={ACM Transactions on Mathematical Software (TOMS)},\n"
975: " volume={37},\n"
976: " number={1},\n"
977: " pages={1--16},\n"
978: " year={2010},\n"
979: " publisher={ACM New York, NY, USA}\n}\n";
981: /*@
982: PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for
983: the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used
984: as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating
985: polynomials in that domain.
987: Input Parameters:
988: + dim - the number of variables in the multivariate polynomials
989: . npoints - the number of points to evaluate the polynomials at
990: . points - [npoints x dim] array of point coordinates
991: . degree - the degree (sum of degrees on the variables in a monomial) of the polynomial space to evaluate. There are ((dim + degree) choose dim) polynomials in this space.
992: - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives
993: in the jet. Choosing k = 0 means to evaluate just the function and no derivatives
995: Output Parameters:
996: - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree)
997: choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
998: three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
999: index; the third (fastest varying) dimension is the index of the evaluation point.
1001: Level: advanced
1003: Note: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
1004: ordering of PetscDTIndexToGradedOrder() and PetscDTGradedOrderToIndex(). For example, in 3D, the polynomial with
1005: leading monomial x^2,y^0,z^1, which has degree tuple (2,0,1), which by PetscDTGradedOrderToIndex() has index 12 (it is the 13th basis function in the space);
1006: the partial derivative $\partial_x \partial_z$ has order tuple (1,0,1), appears at index 6 in the jet (it is the 7th partial derivative in the jet).
1008: The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006.
1010: .seealso: PetscDTGradedOrderToIndex(), PetscDTIndexToGradedOrder(), PetscDTJacobiEvalJet()
1011: @*/
1012: PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1013: {
1014: PetscInt degidx, kidx, d, pt;
1015: PetscInt Nk, Ndeg;
1016: PetscInt *ktup, *degtup;
1017: PetscReal *scales, initscale, scaleexp;
1019: PetscCitationsRegister(PKDCitation, &PKDCite);
1020: PetscDTBinomialInt(dim + k, k, &Nk);
1021: PetscDTBinomialInt(degree + dim, degree, &Ndeg);
1022: PetscMalloc2(dim, °tup, dim, &ktup);
1023: PetscMalloc1(Ndeg, &scales);
1024: initscale = 1.;
1025: if (dim > 1) {
1026: PetscDTBinomial(dim,2,&scaleexp);
1027: initscale = PetscPowReal(2.,scaleexp*0.5);
1028: }
1029: for (degidx = 0; degidx < Ndeg; degidx++) {
1030: PetscInt e, i;
1031: PetscInt m1idx = -1, m2idx = -1;
1032: PetscInt n;
1033: PetscInt degsum;
1034: PetscReal alpha;
1035: PetscReal cnm1, cnm1x, cnm2;
1036: PetscReal norm;
1038: PetscDTIndexToGradedOrder(dim, degidx, degtup);
1039: for (d = dim - 1; d >= 0; d--) if (degtup[d]) break;
1040: if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1041: scales[degidx] = initscale;
1042: for (e = 0; e < dim; e++) {
1043: PetscDTJacobiNorm(e,0.,0,&norm);
1044: scales[degidx] /= norm;
1045: }
1046: for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1047: for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1048: continue;
1049: }
1050: n = degtup[d];
1051: degtup[d]--;
1052: PetscDTGradedOrderToIndex(dim, degtup, &m1idx);
1053: if (degtup[d] > 0) {
1054: degtup[d]--;
1055: PetscDTGradedOrderToIndex(dim, degtup, &m2idx);
1056: degtup[d]++;
1057: }
1058: degtup[d]++;
1059: for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1060: alpha = 2 * degsum + d;
1061: PetscDTJacobiRecurrence_Internal(n,alpha,0.,cnm1,cnm1x,cnm2);
1063: scales[degidx] = initscale;
1064: for (e = 0, degsum = 0; e < dim; e++) {
1065: PetscInt f;
1066: PetscReal ealpha;
1067: PetscReal enorm;
1069: ealpha = 2 * degsum + e;
1070: for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1071: PetscDTJacobiNorm(ealpha,0.,degtup[e],&enorm);
1072: scales[degidx] /= enorm;
1073: degsum += degtup[e];
1074: }
1076: for (pt = 0; pt < npoints; pt++) {
1077: /* compute the multipliers */
1078: PetscReal thetanm1, thetanm1x, thetanm2;
1080: thetanm1x = dim - (d+1) + 2.*points[pt * dim + d];
1081: for (e = d+1; e < dim; e++) thetanm1x += points[pt * dim + e];
1082: thetanm1x *= 0.5;
1083: thetanm1 = (2. - (dim-(d+1)));
1084: for (e = d+1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1085: thetanm1 *= 0.5;
1086: thetanm2 = thetanm1 * thetanm1;
1088: for (kidx = 0; kidx < Nk; kidx++) {
1089: PetscInt f;
1091: PetscDTIndexToGradedOrder(dim, kidx, ktup);
1092: /* first sum in the same derivative terms */
1093: p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1094: if (m2idx >= 0) {
1095: p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];
1096: }
1098: for (f = d; f < dim; f++) {
1099: PetscInt km1idx, mplty = ktup[f];
1101: if (!mplty) continue;
1102: ktup[f]--;
1103: PetscDTGradedOrderToIndex(dim, ktup, &km1idx);
1105: /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1106: /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1107: /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1108: if (f > d) {
1109: PetscInt f2;
1111: p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1112: if (m2idx >= 0) {
1113: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1114: /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1115: for (f2 = f; f2 < dim; f2++) {
1116: PetscInt km2idx, mplty2 = ktup[f2];
1117: PetscInt factor;
1119: if (!mplty2) continue;
1120: ktup[f2]--;
1121: PetscDTGradedOrderToIndex(dim, ktup, &km2idx);
1123: factor = mplty * mplty2;
1124: if (f == f2) factor /= 2;
1125: p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1126: ktup[f2]++;
1127: }
1128: }
1129: } else {
1130: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1131: }
1132: ktup[f]++;
1133: }
1134: }
1135: }
1136: }
1137: for (degidx = 0; degidx < Ndeg; degidx++) {
1138: PetscReal scale = scales[degidx];
1139: PetscInt i;
1141: for (i = 0; i < Nk * npoints; i++) p[degidx*Nk*npoints + i] *= scale;
1142: }
1143: PetscFree(scales);
1144: PetscFree2(degtup, ktup);
1145: return 0;
1146: }
1148: /*@
1149: PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree,
1150: which can be evaluated in PetscDTPTrimmedEvalJet().
1152: Input Parameters:
1153: + dim - the number of variables in the multivariate polynomials
1154: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space.
1155: - formDegree - the degree of the form
1157: Output Parameters:
1158: - size - The number ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree))
1160: Level: advanced
1162: .seealso: PetscDTPTrimmedEvalJet()
1163: @*/
1164: PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size)
1165: {
1166: PetscInt Nrk, Nbpt; // number of trimmed polynomials
1168: formDegree = PetscAbsInt(formDegree);
1169: PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt);
1170: PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk);
1171: Nbpt *= Nrk;
1172: *size = Nbpt;
1173: return 0;
1174: }
1176: /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it
1177: * was inferior to this implementation */
1178: static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1179: {
1180: PetscInt formDegreeOrig = formDegree;
1181: PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE;
1183: formDegree = PetscAbsInt(formDegreeOrig);
1184: if (formDegree == 0) {
1185: PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p);
1186: return 0;
1187: }
1188: if (formDegree == dim) {
1189: PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p);
1190: return 0;
1191: }
1192: PetscInt Nbpt;
1193: PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt);
1194: PetscInt Nf;
1195: PetscDTBinomialInt(dim, formDegree, &Nf);
1196: PetscInt Nk;
1197: PetscDTBinomialInt(dim + jetDegree, dim, &Nk);
1198: PetscArrayzero(p, Nbpt * Nf * Nk * npoints);
1200: PetscInt Nbpm1; // number of scalar polynomials up to degree - 1;
1201: PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1);
1202: PetscReal *p_scalar;
1203: PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar);
1204: PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar);
1205: PetscInt total = 0;
1206: // First add the full polynomials up to degree - 1 into the basis: take the scalar
1207: // and copy one for each form component
1208: for (PetscInt i = 0; i < Nbpm1; i++) {
1209: const PetscReal *src = &p_scalar[i * Nk * npoints];
1210: for (PetscInt f = 0; f < Nf; f++) {
1211: PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints];
1212: PetscArraycpy(dest, src, Nk * npoints);
1213: }
1214: }
1215: PetscInt *form_atoms;
1216: PetscMalloc1(formDegree + 1, &form_atoms);
1217: // construct the interior product pattern
1218: PetscInt (*pattern)[3];
1219: PetscInt Nf1; // number of formDegree + 1 forms
1220: PetscDTBinomialInt(dim, formDegree + 1, &Nf1);
1221: PetscInt nnz = Nf1 * (formDegree+1);
1222: PetscMalloc1(Nf1 * (formDegree+1), &pattern);
1223: PetscDTAltVInteriorPattern(dim, formDegree+1, pattern);
1224: PetscReal centroid = (1. - dim) / (dim + 1.);
1225: PetscInt *deriv;
1226: PetscMalloc1(dim, &deriv);
1227: for (PetscInt d = dim; d >= formDegree + 1; d--) {
1228: PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0
1229: // (equal to the number of formDegree forms in dimension d-1)
1230: PetscDTBinomialInt(d - 1, formDegree, &Nfd1);
1231: // The number of homogeneous (degree-1) scalar polynomials in d variables
1232: PetscInt Nh;
1233: PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh);
1234: const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints];
1235: for (PetscInt b = 0; b < Nh; b++) {
1236: const PetscReal *h_s = &h_scalar[b * Nk * npoints];
1237: for (PetscInt f = 0; f < Nfd1; f++) {
1238: // construct all formDegree+1 forms that start with dx_(dim - d) /\ ...
1239: form_atoms[0] = dim - d;
1240: PetscDTEnumSubset(d-1, formDegree, f, &form_atoms[1]);
1241: for (PetscInt i = 0; i < formDegree; i++) {
1242: form_atoms[1+i] += form_atoms[0] + 1;
1243: }
1244: PetscInt f_ind; // index of the resulting form
1245: PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind);
1246: PetscReal *p_f = &p[total++ * Nf * Nk * npoints];
1247: for (PetscInt nz = 0; nz < nnz; nz++) {
1248: PetscInt i = pattern[nz][0]; // formDegree component
1249: PetscInt j = pattern[nz][1]; // (formDegree + 1) component
1250: PetscInt v = pattern[nz][2]; // coordinate component
1251: PetscReal scale = v < 0 ? -1. : 1.;
1253: i = formNegative ? (Nf - 1 - i) : i;
1254: scale = (formNegative && (i & 1)) ? -scale : scale;
1255: v = v < 0 ? -(v + 1) : v;
1256: if (j != f_ind) {
1257: continue;
1258: }
1259: PetscReal *p_i = &p_f[i * Nk * npoints];
1260: for (PetscInt jet = 0; jet < Nk; jet++) {
1261: const PetscReal *h_jet = &h_s[jet * npoints];
1262: PetscReal *p_jet = &p_i[jet * npoints];
1264: for (PetscInt pt = 0; pt < npoints; pt++) {
1265: p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid);
1266: }
1267: PetscDTIndexToGradedOrder(dim, jet, deriv);
1268: deriv[v]++;
1269: PetscReal mult = deriv[v];
1270: PetscInt l;
1271: PetscDTGradedOrderToIndex(dim, deriv, &l);
1272: if (l >= Nk) {
1273: continue;
1274: }
1275: p_jet = &p_i[l * npoints];
1276: for (PetscInt pt = 0; pt < npoints; pt++) {
1277: p_jet[pt] += scale * mult * h_jet[pt];
1278: }
1279: deriv[v]--;
1280: }
1281: }
1282: }
1283: }
1284: }
1286: PetscFree(deriv);
1287: PetscFree(pattern);
1288: PetscFree(form_atoms);
1289: PetscFree(p_scalar);
1290: return 0;
1291: }
1293: /*@
1294: PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to
1295: a given degree.
1297: Input Parameters:
1298: + dim - the number of variables in the multivariate polynomials
1299: . npoints - the number of points to evaluate the polynomials at
1300: . points - [npoints x dim] array of point coordinates
1301: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate.
1302: There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space.
1303: (You can use PetscDTPTrimmedSize() to compute this size.)
1304: . formDegree - the degree of the form
1305: - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives
1306: in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives
1308: Output Parameters:
1309: - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is
1310: PetscDTPTrimmedSize() x ((dim + formDegree) choose dim) x ((dim + k) choose dim) x npoints,
1311: which also describes the order of the dimensions of this
1312: four-dimensional array:
1313: the first (slowest varying) dimension is basis function index;
1314: the second dimension is component of the form;
1315: the third dimension is jet index;
1316: the fourth (fastest varying) dimension is the index of the evaluation point.
1318: Level: advanced
1320: Note: The ordering of the basis functions is not graded, so the basis functions are not nested by degree like PetscDTPKDEvalJet().
1321: The basis functions are not an L2-orthonormal basis on any particular domain.
1323: The implementation is based on the description of the trimmed polynomials up to degree r as
1324: the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to
1325: homogeneous polynomials of degree (r-1).
1327: .seealso: PetscDTPKDEvalJet(), PetscDTPTrimmedSize()
1328: @*/
1329: PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1330: {
1331: PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p);
1332: return 0;
1333: }
1335: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1336: * with lds n; diag and subdiag are overwritten */
1337: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[],
1338: PetscReal eigs[], PetscScalar V[])
1339: {
1340: char jobz = 'V'; /* eigenvalues and eigenvectors */
1341: char range = 'A'; /* all eigenvalues will be found */
1342: PetscReal VL = 0.; /* ignored because range is 'A' */
1343: PetscReal VU = 0.; /* ignored because range is 'A' */
1344: PetscBLASInt IL = 0; /* ignored because range is 'A' */
1345: PetscBLASInt IU = 0; /* ignored because range is 'A' */
1346: PetscReal abstol = 0.; /* unused */
1347: PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1348: PetscBLASInt *isuppz;
1349: PetscBLASInt lwork, liwork;
1350: PetscReal workquery;
1351: PetscBLASInt iworkquery;
1352: PetscBLASInt *iwork;
1353: PetscBLASInt info;
1354: PetscReal *work = NULL;
1356: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1357: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1358: #endif
1359: PetscBLASIntCast(n, &bn);
1360: PetscBLASIntCast(n, &ldz);
1361: #if !defined(PETSC_MISSING_LAPACK_STEGR)
1362: PetscMalloc1(2 * n, &isuppz);
1363: lwork = -1;
1364: liwork = -1;
1365: PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info));
1367: lwork = (PetscBLASInt) workquery;
1368: liwork = (PetscBLASInt) iworkquery;
1369: PetscMalloc2(lwork, &work, liwork, &iwork);
1370: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1371: PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info));
1372: PetscFPTrapPop();
1374: PetscFree2(work, iwork);
1375: PetscFree(isuppz);
1376: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1377: jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1378: tridiagonal matrix. Z is initialized to the identity
1379: matrix. */
1380: PetscMalloc1(PetscMax(1,2*n-2),&work);
1381: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info));
1382: PetscFPTrapPop();
1384: PetscFree(work);
1385: PetscArraycpy(eigs,diag,n);
1386: #endif
1387: return 0;
1388: }
1390: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1391: * quadrature rules on the interval [-1, 1] */
1392: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1393: {
1394: PetscReal twoab1;
1395: PetscInt m = n - 2;
1396: PetscReal a = alpha + 1.;
1397: PetscReal b = beta + 1.;
1398: PetscReal gra, grb;
1400: twoab1 = PetscPowReal(2., a + b - 1.);
1401: #if defined(PETSC_HAVE_LGAMMA)
1402: grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) -
1403: (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.)));
1404: gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) -
1405: (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.)));
1406: #else
1407: {
1408: PetscInt alphai = (PetscInt) alpha;
1409: PetscInt betai = (PetscInt) beta;
1411: if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) {
1412: PetscReal binom1, binom2;
1414: PetscDTBinomial(m+b, b, &binom1);
1415: PetscDTBinomial(m+a+b, b, &binom2);
1416: grb = 1./ (binom1 * binom2);
1417: PetscDTBinomial(m+a, a, &binom1);
1418: PetscDTBinomial(m+a+b, a, &binom2);
1419: gra = 1./ (binom1 * binom2);
1420: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1421: }
1422: #endif
1423: *leftw = twoab1 * grb / b;
1424: *rightw = twoab1 * gra / a;
1425: return 0;
1426: }
1428: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1429: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1430: static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1431: {
1432: PetscReal pn1, pn2;
1433: PetscReal cnm1, cnm1x, cnm2;
1434: PetscInt k;
1436: if (!n) {*P = 1.0; return 0;}
1437: PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2);
1438: pn2 = 1.;
1439: pn1 = cnm1 + cnm1x*x;
1440: if (n == 1) {*P = pn1; return 0;}
1441: *P = 0.0;
1442: for (k = 2; k < n+1; ++k) {
1443: PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2);
1445: *P = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2;
1446: pn2 = pn1;
1447: pn1 = *P;
1448: }
1449: return 0;
1450: }
1452: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1453: static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1454: {
1455: PetscReal nP;
1456: PetscInt i;
1458: *P = 0.0;
1459: if (k > n) return 0;
1460: PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);
1461: for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1462: *P = nP;
1463: return 0;
1464: }
1466: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1467: {
1468: PetscInt maxIter = 100;
1469: PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1470: PetscReal a1, a6, gf;
1471: PetscInt k;
1474: a1 = PetscPowReal(2.0, a+b+1);
1475: #if defined(PETSC_HAVE_LGAMMA)
1476: {
1477: PetscReal a2, a3, a4, a5;
1478: a2 = PetscLGamma(a + npoints + 1);
1479: a3 = PetscLGamma(b + npoints + 1);
1480: a4 = PetscLGamma(a + b + npoints + 1);
1481: a5 = PetscLGamma(npoints + 1);
1482: gf = PetscExpReal(a2 + a3 - (a4 + a5));
1483: }
1484: #else
1485: {
1486: PetscInt ia, ib;
1488: ia = (PetscInt) a;
1489: ib = (PetscInt) b;
1490: gf = 1.;
1491: if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1492: for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1493: } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1494: for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1495: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1496: }
1497: #endif
1499: a6 = a1 * gf;
1500: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1501: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1502: for (k = 0; k < npoints; ++k) {
1503: PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP;
1504: PetscInt j;
1506: if (k > 0) r = 0.5 * (r + x[k-1]);
1507: for (j = 0; j < maxIter; ++j) {
1508: PetscReal s = 0.0, delta, f, fp;
1509: PetscInt i;
1511: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1512: PetscDTComputeJacobi(a, b, npoints, r, &f);
1513: PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);
1514: delta = f / (fp - f * s);
1515: r = r - delta;
1516: if (PetscAbsReal(delta) < eps) break;
1517: }
1518: x[k] = r;
1519: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);
1520: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1521: }
1522: return 0;
1523: }
1525: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1526: * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1527: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1528: {
1529: PetscInt i;
1531: for (i = 0; i < nPoints; i++) {
1532: PetscReal A, B, C;
1534: PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C);
1535: d[i] = -A / B;
1536: if (i) s[i-1] *= C / B;
1537: if (i < nPoints - 1) s[i] = 1. / B;
1538: }
1539: return 0;
1540: }
1542: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1543: {
1544: PetscReal mu0;
1545: PetscReal ga, gb, gab;
1546: PetscInt i;
1548: PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);
1550: #if defined(PETSC_HAVE_TGAMMA)
1551: ga = PetscTGamma(a + 1);
1552: gb = PetscTGamma(b + 1);
1553: gab = PetscTGamma(a + b + 2);
1554: #else
1555: {
1556: PetscInt ia, ib;
1558: ia = (PetscInt) a;
1559: ib = (PetscInt) b;
1560: if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1561: PetscDTFactorial(ia, &ga);
1562: PetscDTFactorial(ib, &gb);
1563: PetscDTFactorial(ia + ib + 1, &gb);
1564: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
1565: }
1566: #endif
1567: mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab;
1569: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1570: {
1571: PetscReal *diag, *subdiag;
1572: PetscScalar *V;
1574: PetscMalloc2(npoints, &diag, npoints, &subdiag);
1575: PetscMalloc1(npoints*npoints, &V);
1576: PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);
1577: for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1578: PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);
1579: for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1580: PetscFree(V);
1581: PetscFree2(diag, subdiag);
1582: }
1583: #else
1584: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1585: #endif
1586: { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1587: eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1588: the eigenvalues are sorted */
1589: PetscBool sorted;
1591: PetscSortedReal(npoints, x, &sorted);
1592: if (!sorted) {
1593: PetscInt *order, i;
1594: PetscReal *tmp;
1596: PetscMalloc2(npoints, &order, npoints, &tmp);
1597: for (i = 0; i < npoints; i++) order[i] = i;
1598: PetscSortRealWithPermutation(npoints, x, order);
1599: PetscArraycpy(tmp, x, npoints);
1600: for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1601: PetscArraycpy(tmp, w, npoints);
1602: for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1603: PetscFree2(order, tmp);
1604: }
1605: }
1606: return 0;
1607: }
1609: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1610: {
1612: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1616: if (newton) {
1617: PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);
1618: } else {
1619: PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);
1620: }
1621: if (alpha == beta) { /* symmetrize */
1622: PetscInt i;
1623: for (i = 0; i < (npoints + 1) / 2; i++) {
1624: PetscInt j = npoints - 1 - i;
1625: PetscReal xi = x[i];
1626: PetscReal xj = x[j];
1627: PetscReal wi = w[i];
1628: PetscReal wj = w[j];
1630: x[i] = (xi - xj) / 2.;
1631: x[j] = (xj - xi) / 2.;
1632: w[i] = w[j] = (wi + wj) / 2.;
1633: }
1634: }
1635: return 0;
1636: }
1638: /*@
1639: PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1640: $(x-a)^\alpha (x-b)^\beta$.
1642: Not collective
1644: Input Parameters:
1645: + npoints - the number of points in the quadrature rule
1646: . a - the left endpoint of the interval
1647: . b - the right endpoint of the interval
1648: . alpha - the left exponent
1649: - beta - the right exponent
1651: Output Parameters:
1652: + x - array of length npoints, the locations of the quadrature points
1653: - w - array of length npoints, the weights of the quadrature points
1655: Level: intermediate
1657: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1658: @*/
1659: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1660: {
1661: PetscInt i;
1663: PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1664: if (a != -1. || b != 1.) { /* shift */
1665: for (i = 0; i < npoints; i++) {
1666: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1667: w[i] *= (b - a) / 2.;
1668: }
1669: }
1670: return 0;
1671: }
1673: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1674: {
1675: PetscInt i;
1678: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1682: x[0] = -1.;
1683: x[npoints-1] = 1.;
1684: if (npoints > 2) {
1685: PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);
1686: }
1687: for (i = 1; i < npoints - 1; i++) {
1688: w[i] /= (1. - x[i]*x[i]);
1689: }
1690: PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);
1691: return 0;
1692: }
1694: /*@
1695: PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1696: $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.
1698: Not collective
1700: Input Parameters:
1701: + npoints - the number of points in the quadrature rule
1702: . a - the left endpoint of the interval
1703: . b - the right endpoint of the interval
1704: . alpha - the left exponent
1705: - beta - the right exponent
1707: Output Parameters:
1708: + x - array of length npoints, the locations of the quadrature points
1709: - w - array of length npoints, the weights of the quadrature points
1711: Level: intermediate
1713: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1714: @*/
1715: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1716: {
1717: PetscInt i;
1719: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1720: if (a != -1. || b != 1.) { /* shift */
1721: for (i = 0; i < npoints; i++) {
1722: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1723: w[i] *= (b - a) / 2.;
1724: }
1725: }
1726: return 0;
1727: }
1729: /*@
1730: PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1732: Not Collective
1734: Input Parameters:
1735: + npoints - number of points
1736: . a - left end of interval (often-1)
1737: - b - right end of interval (often +1)
1739: Output Parameters:
1740: + x - quadrature points
1741: - w - quadrature weights
1743: Level: intermediate
1745: References:
1746: . * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
1748: .seealso: PetscDTLegendreEval()
1749: @*/
1750: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
1751: {
1752: PetscInt i;
1754: PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);
1755: if (a != -1. || b != 1.) { /* shift */
1756: for (i = 0; i < npoints; i++) {
1757: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1758: w[i] *= (b - a) / 2.;
1759: }
1760: }
1761: return 0;
1762: }
1764: /*@C
1765: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1766: nodes of a given size on the domain [-1,1]
1768: Not Collective
1770: Input Parameters:
1771: + n - number of grid nodes
1772: - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON
1774: Output Parameters:
1775: + x - quadrature points
1776: - w - quadrature weights
1778: Notes:
1779: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1780: close enough to the desired solution
1782: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1784: See https://epubs.siam.org/doi/abs/10.1137/110855442 https://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes
1786: Level: intermediate
1788: .seealso: PetscDTGaussQuadrature()
1790: @*/
1791: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
1792: {
1793: PetscBool newton;
1796: newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1797: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);
1798: return 0;
1799: }
1801: /*@
1802: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1804: Not Collective
1806: Input Parameters:
1807: + dim - The spatial dimension
1808: . Nc - The number of components
1809: . npoints - number of points in one dimension
1810: . a - left end of interval (often-1)
1811: - b - right end of interval (often +1)
1813: Output Parameter:
1814: . q - A PetscQuadrature object
1816: Level: intermediate
1818: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
1819: @*/
1820: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1821: {
1822: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
1823: PetscReal *x, *w, *xw, *ww;
1825: PetscMalloc1(totpoints*dim,&x);
1826: PetscMalloc1(totpoints*Nc,&w);
1827: /* Set up the Golub-Welsch system */
1828: switch (dim) {
1829: case 0:
1830: PetscFree(x);
1831: PetscFree(w);
1832: PetscMalloc1(1, &x);
1833: PetscMalloc1(Nc, &w);
1834: x[0] = 0.0;
1835: for (c = 0; c < Nc; ++c) w[c] = 1.0;
1836: break;
1837: case 1:
1838: PetscMalloc1(npoints,&ww);
1839: PetscDTGaussQuadrature(npoints, a, b, x, ww);
1840: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
1841: PetscFree(ww);
1842: break;
1843: case 2:
1844: PetscMalloc2(npoints,&xw,npoints,&ww);
1845: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1846: for (i = 0; i < npoints; ++i) {
1847: for (j = 0; j < npoints; ++j) {
1848: x[(i*npoints+j)*dim+0] = xw[i];
1849: x[(i*npoints+j)*dim+1] = xw[j];
1850: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
1851: }
1852: }
1853: PetscFree2(xw,ww);
1854: break;
1855: case 3:
1856: PetscMalloc2(npoints,&xw,npoints,&ww);
1857: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1858: for (i = 0; i < npoints; ++i) {
1859: for (j = 0; j < npoints; ++j) {
1860: for (k = 0; k < npoints; ++k) {
1861: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
1862: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
1863: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
1864: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
1865: }
1866: }
1867: }
1868: PetscFree2(xw,ww);
1869: break;
1870: default:
1871: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1872: }
1873: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1874: PetscQuadratureSetOrder(*q, 2*npoints-1);
1875: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1876: PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
1877: return 0;
1878: }
1880: /*@
1881: PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex
1883: Not Collective
1885: Input Parameters:
1886: + dim - The simplex dimension
1887: . Nc - The number of components
1888: . npoints - The number of points in one dimension
1889: . a - left end of interval (often-1)
1890: - b - right end of interval (often +1)
1892: Output Parameter:
1893: . q - A PetscQuadrature object
1895: Level: intermediate
1897: References:
1898: . * - Karniadakis and Sherwin. FIAT
1900: Note: For dim == 1, this is Gauss-Legendre quadrature
1902: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
1903: @*/
1904: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1905: {
1906: PetscInt totprev, totrem;
1907: PetscInt totpoints;
1908: PetscReal *p1, *w1;
1909: PetscReal *x, *w;
1910: PetscInt i, j, k, l, m, pt, c;
1913: totpoints = 1;
1914: for (i = 0, totpoints = 1; i < dim; i++) totpoints *= npoints;
1915: PetscMalloc1(totpoints*dim, &x);
1916: PetscMalloc1(totpoints*Nc, &w);
1917: PetscMalloc2(npoints, &p1, npoints, &w1);
1918: for (i = 0; i < totpoints*Nc; i++) w[i] = 1.;
1919: for (i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; i++) {
1920: PetscReal mul;
1922: mul = PetscPowReal(2.,-i);
1923: PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1);
1924: for (pt = 0, l = 0; l < totprev; l++) {
1925: for (j = 0; j < npoints; j++) {
1926: for (m = 0; m < totrem; m++, pt++) {
1927: for (k = 0; k < i; k++) x[pt*dim+k] = (x[pt*dim+k]+1.)*(1.-p1[j])*0.5 - 1.;
1928: x[pt * dim + i] = p1[j];
1929: for (c = 0; c < Nc; c++) w[pt*Nc + c] *= mul * w1[j];
1930: }
1931: }
1932: }
1933: totprev *= npoints;
1934: totrem /= npoints;
1935: }
1936: PetscFree2(p1, w1);
1937: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1938: PetscQuadratureSetOrder(*q, 2*npoints-1);
1939: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1940: PetscObjectChangeTypeName((PetscObject)*q,"StroudConical");
1941: return 0;
1942: }
1944: /*@
1945: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
1947: Not Collective
1949: Input Parameters:
1950: + dim - The cell dimension
1951: . level - The number of points in one dimension, 2^l
1952: . a - left end of interval (often-1)
1953: - b - right end of interval (often +1)
1955: Output Parameter:
1956: . q - A PetscQuadrature object
1958: Level: intermediate
1960: .seealso: PetscDTGaussTensorQuadrature()
1961: @*/
1962: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1963: {
1964: const PetscInt p = 16; /* Digits of precision in the evaluation */
1965: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1966: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1967: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1968: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
1969: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
1970: PetscReal *x, *w;
1971: PetscInt K, k, npoints;
1975: /* Find K such that the weights are < 32 digits of precision */
1976: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1977: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1978: }
1979: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1980: PetscQuadratureSetOrder(*q, 2*K+1);
1981: npoints = 2*K-1;
1982: PetscMalloc1(npoints*dim, &x);
1983: PetscMalloc1(npoints, &w);
1984: /* Center term */
1985: x[0] = beta;
1986: w[0] = 0.5*alpha*PETSC_PI;
1987: for (k = 1; k < K; ++k) {
1988: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1989: xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1990: x[2*k-1] = -alpha*xk+beta;
1991: w[2*k-1] = wk;
1992: x[2*k+0] = alpha*xk+beta;
1993: w[2*k+0] = wk;
1994: }
1995: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1996: return 0;
1997: }
1999: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2000: {
2001: const PetscInt p = 16; /* Digits of precision in the evaluation */
2002: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
2003: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
2004: PetscReal h = 1.0; /* Step size, length between x_k */
2005: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2006: PetscReal osum = 0.0; /* Integral on last level */
2007: PetscReal psum = 0.0; /* Integral on the level before the last level */
2008: PetscReal sum; /* Integral on current level */
2009: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2010: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2011: PetscReal wk; /* Quadrature weight at x_k */
2012: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
2013: PetscInt d; /* Digits of precision in the integral */
2016: /* Center term */
2017: func(&beta, ctx, &lval);
2018: sum = 0.5*alpha*PETSC_PI*lval;
2019: /* */
2020: do {
2021: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
2022: PetscInt k = 1;
2024: ++l;
2025: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
2026: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2027: psum = osum;
2028: osum = sum;
2029: h *= 0.5;
2030: sum *= 0.5;
2031: do {
2032: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
2033: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
2034: lx = -alpha*(1.0 - yk)+beta;
2035: rx = alpha*(1.0 - yk)+beta;
2036: func(&lx, ctx, &lval);
2037: func(&rx, ctx, &rval);
2038: lterm = alpha*wk*lval;
2039: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
2040: sum += lterm;
2041: rterm = alpha*wk*rval;
2042: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
2043: sum += rterm;
2044: ++k;
2045: /* Only need to evaluate every other point on refined levels */
2046: if (l != 1) ++k;
2047: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
2049: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
2050: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
2051: d3 = PetscLog10Real(maxTerm) - p;
2052: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
2053: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
2054: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
2055: } while (d < digits && l < 12);
2056: *sol = sum;
2058: return 0;
2059: }
2061: #if defined(PETSC_HAVE_MPFR)
2062: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2063: {
2064: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
2065: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2066: mpfr_t alpha; /* Half-width of the integration interval */
2067: mpfr_t beta; /* Center of the integration interval */
2068: mpfr_t h; /* Step size, length between x_k */
2069: mpfr_t osum; /* Integral on last level */
2070: mpfr_t psum; /* Integral on the level before the last level */
2071: mpfr_t sum; /* Integral on current level */
2072: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2073: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2074: mpfr_t wk; /* Quadrature weight at x_k */
2075: PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */
2076: PetscInt d; /* Digits of precision in the integral */
2077: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
2080: /* Create high precision storage */
2081: mpfr_inits2(PetscCeilReal(safetyFactor*digits*PetscLogReal(10.)/PetscLogReal(2.)), alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2082: /* Initialization */
2083: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
2084: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
2085: mpfr_set_d(osum, 0.0, MPFR_RNDN);
2086: mpfr_set_d(psum, 0.0, MPFR_RNDN);
2087: mpfr_set_d(h, 1.0, MPFR_RNDN);
2088: mpfr_const_pi(pi2, MPFR_RNDN);
2089: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
2090: /* Center term */
2091: rtmp = 0.5*(b+a);
2092: func(&rtmp, ctx, &lval);
2093: mpfr_set(sum, pi2, MPFR_RNDN);
2094: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
2095: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
2096: /* */
2097: do {
2098: PetscReal d1, d2, d3, d4;
2099: PetscInt k = 1;
2101: ++l;
2102: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
2103: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
2104: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2105: mpfr_set(psum, osum, MPFR_RNDN);
2106: mpfr_set(osum, sum, MPFR_RNDN);
2107: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
2108: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
2109: do {
2110: mpfr_set_si(kh, k, MPFR_RNDN);
2111: mpfr_mul(kh, kh, h, MPFR_RNDN);
2112: /* Weight */
2113: mpfr_set(wk, h, MPFR_RNDN);
2114: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
2115: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
2116: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
2117: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2118: mpfr_sqr(tmp, tmp, MPFR_RNDN);
2119: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
2120: mpfr_div(wk, wk, tmp, MPFR_RNDN);
2121: /* Abscissa */
2122: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
2123: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2124: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2125: mpfr_exp(tmp, msinh, MPFR_RNDN);
2126: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2127: /* Quadrature points */
2128: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
2129: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
2130: mpfr_add(lx, lx, beta, MPFR_RNDU);
2131: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
2132: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
2133: mpfr_add(rx, rx, beta, MPFR_RNDD);
2134: /* Evaluation */
2135: rtmp = mpfr_get_d(lx, MPFR_RNDU);
2136: func(&rtmp, ctx, &lval);
2137: rtmp = mpfr_get_d(rx, MPFR_RNDD);
2138: func(&rtmp, ctx, &rval);
2139: /* Update */
2140: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2141: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
2142: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2143: mpfr_abs(tmp, tmp, MPFR_RNDN);
2144: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2145: mpfr_set(curTerm, tmp, MPFR_RNDN);
2146: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2147: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
2148: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2149: mpfr_abs(tmp, tmp, MPFR_RNDN);
2150: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2151: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
2152: ++k;
2153: /* Only need to evaluate every other point on refined levels */
2154: if (l != 1) ++k;
2155: mpfr_log10(tmp, wk, MPFR_RNDN);
2156: mpfr_abs(tmp, tmp, MPFR_RNDN);
2157: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
2158: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
2159: mpfr_abs(tmp, tmp, MPFR_RNDN);
2160: mpfr_log10(tmp, tmp, MPFR_RNDN);
2161: d1 = mpfr_get_d(tmp, MPFR_RNDN);
2162: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2163: mpfr_abs(tmp, tmp, MPFR_RNDN);
2164: mpfr_log10(tmp, tmp, MPFR_RNDN);
2165: d2 = mpfr_get_d(tmp, MPFR_RNDN);
2166: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2167: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2168: mpfr_log10(tmp, curTerm, MPFR_RNDN);
2169: d4 = mpfr_get_d(tmp, MPFR_RNDN);
2170: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
2171: } while (d < digits && l < 8);
2172: *sol = mpfr_get_d(sum, MPFR_RNDN);
2173: /* Cleanup */
2174: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2175: return 0;
2176: }
2177: #else
2179: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2180: {
2181: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2182: }
2183: #endif
2185: /*@
2186: PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures
2188: Not Collective
2190: Input Parameters:
2191: + q1 - The first quadrature
2192: - q2 - The second quadrature
2194: Output Parameter:
2195: . q - A PetscQuadrature object
2197: Level: intermediate
2199: .seealso: PetscDTGaussTensorQuadrature()
2200: @*/
2201: PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q)
2202: {
2203: const PetscReal *x1, *w1, *x2, *w2;
2204: PetscReal *x, *w;
2205: PetscInt dim1, Nc1, Np1, order1, qa, d1;
2206: PetscInt dim2, Nc2, Np2, order2, qb, d2;
2207: PetscInt dim, Nc, Np, order, qc, d;
2212: PetscQuadratureGetOrder(q1, &order1);
2213: PetscQuadratureGetOrder(q2, &order2);
2215: PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1);
2216: PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2);
2219: dim = dim1 + dim2;
2220: Nc = Nc1;
2221: Np = Np1 * Np2;
2222: order = order1;
2223: PetscQuadratureCreate(PETSC_COMM_SELF, q);
2224: PetscQuadratureSetOrder(*q, order);
2225: PetscMalloc1(Np*dim, &x);
2226: PetscMalloc1(Np, &w);
2227: for (qa = 0, qc = 0; qa < Np1; ++qa) {
2228: for (qb = 0; qb < Np2; ++qb, ++qc) {
2229: for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) {
2230: x[qc*dim+d] = x1[qa*dim1+d1];
2231: }
2232: for (d2 = 0; d2 < dim2; ++d2, ++d) {
2233: x[qc*dim+d] = x2[qb*dim2+d2];
2234: }
2235: w[qc] = w1[qa] * w2[qb];
2236: }
2237: }
2238: PetscQuadratureSetData(*q, dim, Nc, Np, x, w);
2239: return 0;
2240: }
2242: /* Overwrites A. Can only handle full-rank problems with m>=n
2243: * A in column-major format
2244: * Ainv in row-major format
2245: * tau has length m
2246: * worksize must be >= max(1,n)
2247: */
2248: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2249: {
2250: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2251: PetscScalar *A,*Ainv,*R,*Q,Alpha;
2253: #if defined(PETSC_USE_COMPLEX)
2254: {
2255: PetscInt i,j;
2256: PetscMalloc2(m*n,&A,m*n,&Ainv);
2257: for (j=0; j<n; j++) {
2258: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
2259: }
2260: mstride = m;
2261: }
2262: #else
2263: A = A_in;
2264: Ainv = Ainv_out;
2265: #endif
2267: PetscBLASIntCast(m,&M);
2268: PetscBLASIntCast(n,&N);
2269: PetscBLASIntCast(mstride,&lda);
2270: PetscBLASIntCast(worksize,&ldwork);
2271: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2272: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2273: PetscFPTrapPop();
2275: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2277: /* Extract an explicit representation of Q */
2278: Q = Ainv;
2279: PetscArraycpy(Q,A,mstride*n);
2280: K = N; /* full rank */
2281: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2284: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2285: Alpha = 1.0;
2286: ldb = lda;
2287: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
2288: /* Ainv is Q, overwritten with inverse */
2290: #if defined(PETSC_USE_COMPLEX)
2291: {
2292: PetscInt i;
2293: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2294: PetscFree2(A,Ainv);
2295: }
2296: #endif
2297: return 0;
2298: }
2300: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2301: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
2302: {
2303: PetscReal *Bv;
2304: PetscInt i,j;
2306: PetscMalloc1((ninterval+1)*ndegree,&Bv);
2307: /* Point evaluation of L_p on all the source vertices */
2308: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
2309: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2310: for (i=0; i<ninterval; i++) {
2311: for (j=0; j<ndegree; j++) {
2312: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2313: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2314: }
2315: }
2316: PetscFree(Bv);
2317: return 0;
2318: }
2320: /*@
2321: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2323: Not Collective
2325: Input Parameters:
2326: + degree - degree of reconstruction polynomial
2327: . nsource - number of source intervals
2328: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2329: . ntarget - number of target intervals
2330: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
2332: Output Parameter:
2333: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2335: Level: advanced
2337: .seealso: PetscDTLegendreEval()
2338: @*/
2339: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
2340: {
2341: PetscInt i,j,k,*bdegrees,worksize;
2342: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
2343: PetscScalar *tau,*work;
2349: if (PetscDefined(USE_DEBUG)) {
2350: for (i=0; i<nsource; i++) {
2352: }
2353: for (i=0; i<ntarget; i++) {
2355: }
2356: }
2357: xmin = PetscMin(sourcex[0],targetx[0]);
2358: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
2359: center = (xmin + xmax)/2;
2360: hscale = (xmax - xmin)/2;
2361: worksize = nsource;
2362: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
2363: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
2364: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
2365: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
2366: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
2367: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
2368: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
2369: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
2370: for (i=0; i<ntarget; i++) {
2371: PetscReal rowsum = 0;
2372: for (j=0; j<nsource; j++) {
2373: PetscReal sum = 0;
2374: for (k=0; k<degree+1; k++) {
2375: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
2376: }
2377: R[i*nsource+j] = sum;
2378: rowsum += sum;
2379: }
2380: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
2381: }
2382: PetscFree4(bdegrees,sourcey,Bsource,work);
2383: PetscFree4(tau,Bsinv,targety,Btarget);
2384: return 0;
2385: }
2387: /*@C
2388: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2390: Not Collective
2392: Input Parameters:
2393: + n - the number of GLL nodes
2394: . nodes - the GLL nodes
2395: . weights - the GLL weights
2396: - f - the function values at the nodes
2398: Output Parameter:
2399: . in - the value of the integral
2401: Level: beginner
2403: .seealso: PetscDTGaussLobattoLegendreQuadrature()
2405: @*/
2406: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
2407: {
2408: PetscInt i;
2410: *in = 0.;
2411: for (i=0; i<n; i++) {
2412: *in += f[i]*f[i]*weights[i];
2413: }
2414: return 0;
2415: }
2417: /*@C
2418: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2420: Not Collective
2422: Input Parameters:
2423: + n - the number of GLL nodes
2424: . nodes - the GLL nodes
2425: - weights - the GLL weights
2427: Output Parameter:
2428: . A - the stiffness element
2430: Level: beginner
2432: Notes:
2433: Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()
2435: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)
2437: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2439: @*/
2440: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2441: {
2442: PetscReal **A;
2443: const PetscReal *gllnodes = nodes;
2444: const PetscInt p = n-1;
2445: PetscReal z0,z1,z2 = -1,x,Lpj,Lpr;
2446: PetscInt i,j,nn,r;
2448: PetscMalloc1(n,&A);
2449: PetscMalloc1(n*n,&A[0]);
2450: for (i=1; i<n; i++) A[i] = A[i-1]+n;
2452: for (j=1; j<p; j++) {
2453: x = gllnodes[j];
2454: z0 = 1.;
2455: z1 = x;
2456: for (nn=1; nn<p; nn++) {
2457: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2458: z0 = z1;
2459: z1 = z2;
2460: }
2461: Lpj=z2;
2462: for (r=1; r<p; r++) {
2463: if (r == j) {
2464: A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
2465: } else {
2466: x = gllnodes[r];
2467: z0 = 1.;
2468: z1 = x;
2469: for (nn=1; nn<p; nn++) {
2470: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2471: z0 = z1;
2472: z1 = z2;
2473: }
2474: Lpr = z2;
2475: A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
2476: }
2477: }
2478: }
2479: for (j=1; j<p+1; j++) {
2480: x = gllnodes[j];
2481: z0 = 1.;
2482: z1 = x;
2483: for (nn=1; nn<p; nn++) {
2484: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2485: z0 = z1;
2486: z1 = z2;
2487: }
2488: Lpj = z2;
2489: A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
2490: A[0][j] = A[j][0];
2491: }
2492: for (j=0; j<p; j++) {
2493: x = gllnodes[j];
2494: z0 = 1.;
2495: z1 = x;
2496: for (nn=1; nn<p; nn++) {
2497: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2498: z0 = z1;
2499: z1 = z2;
2500: }
2501: Lpj=z2;
2503: A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
2504: A[j][p] = A[p][j];
2505: }
2506: A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
2507: A[p][p]=A[0][0];
2508: *AA = A;
2509: return 0;
2510: }
2512: /*@C
2513: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
2515: Not Collective
2517: Input Parameters:
2518: + n - the number of GLL nodes
2519: . nodes - the GLL nodes
2520: . weights - the GLL weightss
2521: - A - the stiffness element
2523: Level: beginner
2525: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()
2527: @*/
2528: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2529: {
2530: PetscFree((*AA)[0]);
2531: PetscFree(*AA);
2532: *AA = NULL;
2533: return 0;
2534: }
2536: /*@C
2537: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
2539: Not Collective
2541: Input Parameter:
2542: + n - the number of GLL nodes
2543: . nodes - the GLL nodes
2544: . weights - the GLL weights
2546: Output Parameters:
2547: . AA - the stiffness element
2548: - AAT - the transpose of AA (pass in NULL if you do not need this array)
2550: Level: beginner
2552: Notes:
2553: Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()
2555: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2557: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2559: @*/
2560: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2561: {
2562: PetscReal **A, **AT = NULL;
2563: const PetscReal *gllnodes = nodes;
2564: const PetscInt p = n-1;
2565: PetscReal Li, Lj,d0;
2566: PetscInt i,j;
2568: PetscMalloc1(n,&A);
2569: PetscMalloc1(n*n,&A[0]);
2570: for (i=1; i<n; i++) A[i] = A[i-1]+n;
2572: if (AAT) {
2573: PetscMalloc1(n,&AT);
2574: PetscMalloc1(n*n,&AT[0]);
2575: for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
2576: }
2578: if (n==1) {A[0][0] = 0.;}
2579: d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
2580: for (i=0; i<n; i++) {
2581: for (j=0; j<n; j++) {
2582: A[i][j] = 0.;
2583: PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);
2584: PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);
2585: if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
2586: if ((j==i) && (i==0)) A[i][j] = -d0;
2587: if (j==i && i==p) A[i][j] = d0;
2588: if (AT) AT[j][i] = A[i][j];
2589: }
2590: }
2591: if (AAT) *AAT = AT;
2592: *AA = A;
2593: return 0;
2594: }
2596: /*@C
2597: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()
2599: Not Collective
2601: Input Parameters:
2602: + n - the number of GLL nodes
2603: . nodes - the GLL nodes
2604: . weights - the GLL weights
2605: . AA - the stiffness element
2606: - AAT - the transpose of the element
2608: Level: beginner
2610: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()
2612: @*/
2613: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2614: {
2615: PetscFree((*AA)[0]);
2616: PetscFree(*AA);
2617: *AA = NULL;
2618: if (*AAT) {
2619: PetscFree((*AAT)[0]);
2620: PetscFree(*AAT);
2621: *AAT = NULL;
2622: }
2623: return 0;
2624: }
2626: /*@C
2627: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
2629: Not Collective
2631: Input Parameters:
2632: + n - the number of GLL nodes
2633: . nodes - the GLL nodes
2634: - weights - the GLL weightss
2636: Output Parameter:
2637: . AA - the stiffness element
2639: Level: beginner
2641: Notes:
2642: Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()
2644: This is the same as the Gradient operator multiplied by the diagonal mass matrix
2646: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2648: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()
2650: @*/
2651: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2652: {
2653: PetscReal **D;
2654: const PetscReal *gllweights = weights;
2655: const PetscInt glln = n;
2656: PetscInt i,j;
2658: PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
2659: for (i=0; i<glln; i++) {
2660: for (j=0; j<glln; j++) {
2661: D[i][j] = gllweights[i]*D[i][j];
2662: }
2663: }
2664: *AA = D;
2665: return 0;
2666: }
2668: /*@C
2669: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
2671: Not Collective
2673: Input Parameters:
2674: + n - the number of GLL nodes
2675: . nodes - the GLL nodes
2676: . weights - the GLL weights
2677: - A - advection
2679: Level: beginner
2681: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()
2683: @*/
2684: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2685: {
2686: PetscFree((*AA)[0]);
2687: PetscFree(*AA);
2688: *AA = NULL;
2689: return 0;
2690: }
2692: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2693: {
2694: PetscReal **A;
2695: const PetscReal *gllweights = weights;
2696: const PetscInt glln = n;
2697: PetscInt i,j;
2699: PetscMalloc1(glln,&A);
2700: PetscMalloc1(glln*glln,&A[0]);
2701: for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
2702: if (glln==1) {A[0][0] = 0.;}
2703: for (i=0; i<glln; i++) {
2704: for (j=0; j<glln; j++) {
2705: A[i][j] = 0.;
2706: if (j==i) A[i][j] = gllweights[i];
2707: }
2708: }
2709: *AA = A;
2710: return 0;
2711: }
2713: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2714: {
2715: PetscFree((*AA)[0]);
2716: PetscFree(*AA);
2717: *AA = NULL;
2718: return 0;
2719: }
2721: /*@
2722: PetscDTIndexToBary - convert an index into a barycentric coordinate.
2724: Input Parameters:
2725: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
2726: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2727: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
2729: Output Parameter:
2730: . coord - will be filled with the barycentric coordinate
2732: Level: beginner
2734: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2735: least significant and the last index is the most significant.
2737: .seealso: PetscDTBaryToIndex()
2738: @*/
2739: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
2740: {
2741: PetscInt c, d, s, total, subtotal, nexttotal;
2746: if (!len) {
2747: if (!sum && !index) return 0;
2748: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2749: }
2750: for (c = 1, total = 1; c <= len; c++) {
2751: /* total is the number of ways to have a tuple of length c with sum */
2752: if (index < total) break;
2753: total = (total * (sum + c)) / c;
2754: }
2756: for (d = c; d < len; d++) coord[d] = 0;
2757: for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
2758: /* subtotal is the number of ways to have a tuple of length c with sum s */
2759: /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
2760: if ((index + subtotal) >= total) {
2761: coord[--c] = sum - s;
2762: index -= (total - subtotal);
2763: sum = s;
2764: total = nexttotal;
2765: subtotal = 1;
2766: nexttotal = 1;
2767: s = 0;
2768: } else {
2769: subtotal = (subtotal * (c + s)) / (s + 1);
2770: nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
2771: s++;
2772: }
2773: }
2774: return 0;
2775: }
2777: /*@
2778: PetscDTBaryToIndex - convert a barycentric coordinate to an index
2780: Input Parameters:
2781: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
2782: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2783: - coord - a barycentric coordinate with the given length and sum
2785: Output Parameter:
2786: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
2788: Level: beginner
2790: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2791: least significant and the last index is the most significant.
2793: .seealso: PetscDTIndexToBary
2794: @*/
2795: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
2796: {
2797: PetscInt c;
2798: PetscInt i;
2799: PetscInt total;
2803: if (!len) {
2804: if (!sum) {
2805: *index = 0;
2806: return 0;
2807: }
2808: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2809: }
2810: for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
2811: i = total - 1;
2812: c = len - 1;
2813: sum -= coord[c];
2814: while (sum > 0) {
2815: PetscInt subtotal;
2816: PetscInt s;
2818: for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
2819: i -= subtotal;
2820: sum -= coord[--c];
2821: }
2822: *index = i;
2823: return 0;
2824: }