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_shifted[] = {"default", "gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
16: const char *const *const PetscDTNodeTypes = PetscDTNodeTypes_shifted + 1;
18: const char *const PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "PETSCDTSIMPLEXQUAD_", NULL};
19: const char *const *const PetscDTSimplexQuadratureTypes = PetscDTSimplexQuadratureTypes_shifted + 1;
21: static PetscBool GolubWelschCite = PETSC_FALSE;
22: const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
23: " author = {Golub and Welsch},\n"
24: " title = {Calculation of Quadrature Rules},\n"
25: " journal = {Math. Comp.},\n"
26: " volume = {23},\n"
27: " number = {106},\n"
28: " pages = {221--230},\n"
29: " year = {1969}\n}\n";
31: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
32: quadrature rules:
34: - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
35: - in single precision, Newton's method starts producing incorrect roots around n = 15, but
36: the weights from Golub & Welsch become a problem before then: they produces errors
37: in computing the Jacobi-polynomial Gram matrix around n = 6.
39: So we default to Newton's method (required fewer dependencies) */
40: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
42: PetscClassId PETSCQUADRATURE_CLASSID = 0;
44: /*@
45: PetscQuadratureCreate - Create a PetscQuadrature object
47: Collective
49: Input Parameter:
50: . comm - The communicator for the PetscQuadrature object
52: Output Parameter:
53: . q - The PetscQuadrature object
55: Level: beginner
57: .seealso: `PetscQuadratureDestroy()`, `PetscQuadratureGetData()`
58: @*/
59: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
60: {
62: DMInitializePackage();
63: PetscHeaderCreate(*q, PETSCQUADRATURE_CLASSID, "PetscQuadrature", "Quadrature", "DT", comm, PetscQuadratureDestroy, PetscQuadratureView);
64: (*q)->dim = -1;
65: (*q)->Nc = 1;
66: (*q)->order = -1;
67: (*q)->numPoints = 0;
68: (*q)->points = NULL;
69: (*q)->weights = NULL;
70: return 0;
71: }
73: /*@
74: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
76: Collective on q
78: Input Parameter:
79: . q - The PetscQuadrature object
81: Output Parameter:
82: . r - The new PetscQuadrature object
84: Level: beginner
86: .seealso: `PetscQuadratureCreate()`, `PetscQuadratureDestroy()`, `PetscQuadratureGetData()`
87: @*/
88: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
89: {
90: PetscInt order, dim, Nc, Nq;
91: const PetscReal *points, *weights;
92: PetscReal *p, *w;
95: PetscQuadratureCreate(PetscObjectComm((PetscObject)q), r);
96: PetscQuadratureGetOrder(q, &order);
97: PetscQuadratureSetOrder(*r, order);
98: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
99: PetscMalloc1(Nq * dim, &p);
100: PetscMalloc1(Nq * Nc, &w);
101: PetscArraycpy(p, points, Nq * dim);
102: PetscArraycpy(w, weights, Nc * Nq);
103: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
104: return 0;
105: }
107: /*@
108: PetscQuadratureDestroy - Destroys a PetscQuadrature object
110: Collective on q
112: Input Parameter:
113: . q - The PetscQuadrature object
115: Level: beginner
117: .seealso: `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
118: @*/
119: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
120: {
121: if (!*q) return 0;
123: if (--((PetscObject)(*q))->refct > 0) {
124: *q = NULL;
125: return 0;
126: }
127: PetscFree((*q)->points);
128: PetscFree((*q)->weights);
129: PetscHeaderDestroy(q);
130: return 0;
131: }
133: /*@
134: PetscQuadratureGetOrder - Return the order of the method
136: Not collective
138: Input Parameter:
139: . q - The PetscQuadrature object
141: Output Parameter:
142: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
144: Level: intermediate
146: .seealso: `PetscQuadratureSetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
147: @*/
148: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
149: {
152: *order = q->order;
153: return 0;
154: }
156: /*@
157: PetscQuadratureSetOrder - Return the order of the method
159: Not collective
161: Input Parameters:
162: + q - The PetscQuadrature object
163: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
165: Level: intermediate
167: .seealso: `PetscQuadratureGetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
168: @*/
169: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
170: {
172: q->order = order;
173: return 0;
174: }
176: /*@
177: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
179: Not collective
181: Input Parameter:
182: . q - The PetscQuadrature object
184: Output Parameter:
185: . Nc - The number of components
187: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
189: Level: intermediate
191: .seealso: `PetscQuadratureSetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
192: @*/
193: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
194: {
197: *Nc = q->Nc;
198: return 0;
199: }
201: /*@
202: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
204: Not collective
206: Input Parameters:
207: + q - The PetscQuadrature object
208: - Nc - The number of components
210: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
212: Level: intermediate
214: .seealso: `PetscQuadratureGetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
215: @*/
216: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
217: {
219: q->Nc = Nc;
220: return 0;
221: }
223: /*@C
224: PetscQuadratureGetData - Returns the data defining the quadrature
226: Not collective
228: Input Parameter:
229: . q - The PetscQuadrature object
231: Output Parameters:
232: + dim - The spatial dimension
233: . Nc - The number of components
234: . npoints - The number of quadrature points
235: . points - The coordinates of each quadrature point
236: - weights - The weight of each quadrature point
238: Level: intermediate
240: Fortran Notes:
241: From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
243: .seealso: `PetscQuadratureCreate()`, `PetscQuadratureSetData()`
244: @*/
245: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
246: {
248: if (dim) {
250: *dim = q->dim;
251: }
252: if (Nc) {
254: *Nc = q->Nc;
255: }
256: if (npoints) {
258: *npoints = q->numPoints;
259: }
260: if (points) {
262: *points = q->points;
263: }
264: if (weights) {
266: *weights = q->weights;
267: }
268: return 0;
269: }
271: /*@
272: PetscQuadratureEqual - determine whether two quadratures are equivalent
274: Input Parameters:
275: + A - A PetscQuadrature object
276: - B - Another PetscQuadrature object
278: Output Parameters:
279: . equal - PETSC_TRUE if the quadratures are the same
281: Level: intermediate
283: .seealso: `PetscQuadratureCreate()`
284: @*/
285: PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal)
286: {
290: *equal = PETSC_FALSE;
291: if (A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) return 0;
292: for (PetscInt i = 0; i < A->numPoints * A->dim; i++) {
293: if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) return 0;
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) return 0;
302: }
303: *equal = PETSC_TRUE;
304: }
305: return 0;
306: }
308: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
309: {
310: PetscScalar *Js, *Jinvs;
311: PetscInt i, j, k;
312: PetscBLASInt bm, bn, info;
314: if (!m || !n) return 0;
315: PetscBLASIntCast(m, &bm);
316: PetscBLASIntCast(n, &bn);
317: #if defined(PETSC_USE_COMPLEX)
318: PetscMalloc2(m * n, &Js, m * n, &Jinvs);
319: for (i = 0; i < m * n; i++) Js[i] = J[i];
320: #else
321: Js = (PetscReal *)J;
322: Jinvs = Jinv;
323: #endif
324: if (m == n) {
325: PetscBLASInt *pivots;
326: PetscScalar *W;
328: PetscMalloc2(m, &pivots, m, &W);
330: PetscArraycpy(Jinvs, Js, m * m);
331: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
333: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
335: PetscFree2(pivots, W);
336: } else if (m < n) {
337: PetscScalar *JJT;
338: PetscBLASInt *pivots;
339: PetscScalar *W;
341: PetscMalloc1(m * m, &JJT);
342: PetscMalloc2(m, &pivots, m, &W);
343: for (i = 0; i < m; i++) {
344: for (j = 0; j < m; j++) {
345: PetscScalar val = 0.;
347: for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
348: JJT[i * m + j] = val;
349: }
350: }
352: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
354: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
356: for (i = 0; i < n; i++) {
357: for (j = 0; j < m; j++) {
358: PetscScalar val = 0.;
360: for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
361: Jinvs[i * m + j] = val;
362: }
363: }
364: PetscFree2(pivots, W);
365: PetscFree(JJT);
366: } else {
367: PetscScalar *JTJ;
368: PetscBLASInt *pivots;
369: PetscScalar *W;
371: PetscMalloc1(n * n, &JTJ);
372: PetscMalloc2(n, &pivots, n, &W);
373: for (i = 0; i < n; i++) {
374: for (j = 0; j < n; j++) {
375: PetscScalar val = 0.;
377: for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
378: JTJ[i * n + j] = val;
379: }
380: }
382: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
384: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
386: for (i = 0; i < n; i++) {
387: for (j = 0; j < m; j++) {
388: PetscScalar val = 0.;
390: for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
391: Jinvs[i * m + j] = val;
392: }
393: }
394: PetscFree2(pivots, W);
395: PetscFree(JTJ);
396: }
397: #if defined(PETSC_USE_COMPLEX)
398: for (i = 0; i < m * n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
399: PetscFree2(Js, Jinvs);
400: #endif
401: return 0;
402: }
404: /*@
405: PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
407: Collecive on PetscQuadrature
409: Input Parameters:
410: + q - the quadrature functional
411: . imageDim - the dimension of the image of the transformation
412: . origin - a point in the original space
413: . originImage - the image of the origin under the transformation
414: . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
415: - 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]
417: Output Parameters:
418: . 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.
420: 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.
422: Level: intermediate
424: .seealso: `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
425: @*/
426: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
427: {
428: PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
429: const PetscReal *points;
430: const PetscReal *weights;
431: PetscReal *imagePoints, *imageWeights;
432: PetscReal *Jinv;
433: PetscReal *Jinvstar;
437: PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);
438: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);
440: Ncopies = Nc / formSize;
441: PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);
442: imageNc = Ncopies * imageFormSize;
443: PetscMalloc1(Npoints * imageDim, &imagePoints);
444: PetscMalloc1(Npoints * imageNc, &imageWeights);
445: PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);
446: PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);
447: PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);
448: for (pt = 0; pt < Npoints; pt++) {
449: const PetscReal *point = &points[pt * dim];
450: PetscReal *imagePoint = &imagePoints[pt * imageDim];
452: for (i = 0; i < imageDim; i++) {
453: PetscReal val = originImage[i];
455: for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
456: imagePoint[i] = val;
457: }
458: for (c = 0; c < Ncopies; c++) {
459: const PetscReal *form = &weights[pt * Nc + c * formSize];
460: PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
462: for (i = 0; i < imageFormSize; i++) {
463: PetscReal val = 0.;
465: for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
466: imageForm[i] = val;
467: }
468: }
469: }
470: PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);
471: PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);
472: PetscFree2(Jinv, Jinvstar);
473: return 0;
474: }
476: /*@C
477: PetscQuadratureSetData - Sets the data defining the quadrature
479: Not collective
481: Input Parameters:
482: + q - The PetscQuadrature object
483: . dim - The spatial dimension
484: . Nc - The number of components
485: . npoints - The number of quadrature points
486: . points - The coordinates of each quadrature point
487: - weights - The weight of each quadrature point
489: Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
491: Level: intermediate
493: .seealso: `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
494: @*/
495: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
496: {
498: if (dim >= 0) q->dim = dim;
499: if (Nc >= 0) q->Nc = Nc;
500: if (npoints >= 0) q->numPoints = npoints;
501: if (points) {
503: q->points = points;
504: }
505: if (weights) {
507: q->weights = weights;
508: }
509: return 0;
510: }
512: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
513: {
514: PetscInt q, d, c;
515: PetscViewerFormat format;
517: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, "Quadrature of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ") with %" PetscInt_FMT " components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);
518: else PetscViewerASCIIPrintf(v, "Quadrature of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ")\n", quad->order, quad->numPoints, quad->dim);
519: PetscViewerGetFormat(v, &format);
520: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
521: for (q = 0; q < quad->numPoints; ++q) {
522: PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q);
523: PetscViewerASCIIUseTabs(v, PETSC_FALSE);
524: for (d = 0; d < quad->dim; ++d) {
525: if (d) PetscViewerASCIIPrintf(v, ", ");
526: PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d]);
527: }
528: PetscViewerASCIIPrintf(v, ") ");
529: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q);
530: for (c = 0; c < quad->Nc; ++c) {
531: if (c) PetscViewerASCIIPrintf(v, ", ");
532: PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c]);
533: }
534: if (quad->Nc > 1) PetscViewerASCIIPrintf(v, ")");
535: PetscViewerASCIIPrintf(v, "\n");
536: PetscViewerASCIIUseTabs(v, PETSC_TRUE);
537: }
538: return 0;
539: }
541: /*@C
542: PetscQuadratureView - Views a PetscQuadrature object
544: Collective on quad
546: Input Parameters:
547: + quad - The PetscQuadrature object
548: - viewer - The PetscViewer object
550: Level: beginner
552: .seealso: `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
553: @*/
554: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
555: {
556: PetscBool iascii;
560: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer);
561: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
562: PetscViewerASCIIPushTab(viewer);
563: if (iascii) PetscQuadratureView_Ascii(quad, viewer);
564: PetscViewerASCIIPopTab(viewer);
565: return 0;
566: }
568: /*@C
569: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
571: Not collective
573: Input Parameters:
574: + q - The original PetscQuadrature
575: . numSubelements - The number of subelements the original element is divided into
576: . v0 - An array of the initial points for each subelement
577: - jac - An array of the Jacobian mappings from the reference to each subelement
579: Output Parameters:
580: . dim - The dimension
582: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
584: Not available from Fortran
586: Level: intermediate
588: .seealso: `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
589: @*/
590: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
591: {
592: const PetscReal *points, *weights;
593: PetscReal *pointsRef, *weightsRef;
594: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
600: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
601: PetscQuadratureGetOrder(q, &order);
602: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
603: npointsRef = npoints * numSubelements;
604: PetscMalloc1(npointsRef * dim, &pointsRef);
605: PetscMalloc1(npointsRef * Nc, &weightsRef);
606: for (c = 0; c < numSubelements; ++c) {
607: for (p = 0; p < npoints; ++p) {
608: for (d = 0; d < dim; ++d) {
609: pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d];
610: for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0);
611: }
612: /* Could also use detJ here */
613: for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements;
614: }
615: }
616: PetscQuadratureSetOrder(*qref, order);
617: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
618: return 0;
619: }
621: /* Compute the coefficients for the Jacobi polynomial recurrence,
622: *
623: * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
624: */
625: #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \
626: do { \
627: PetscReal _a = (a); \
628: PetscReal _b = (b); \
629: PetscReal _n = (n); \
630: if (n == 1) { \
631: (cnm1) = (_a - _b) * 0.5; \
632: (cnm1x) = (_a + _b + 2.) * 0.5; \
633: (cnm2) = 0.; \
634: } else { \
635: PetscReal _2n = _n + _n; \
636: PetscReal _d = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \
637: PetscReal _n1 = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \
638: PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \
639: PetscReal _n2 = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \
640: (cnm1) = _n1 / _d; \
641: (cnm1x) = _n1x / _d; \
642: (cnm2) = _n2 / _d; \
643: } \
644: } while (0)
646: /*@
647: PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
649: $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
651: Input Parameters:
652: - alpha - the left exponent > -1
653: . beta - the right exponent > -1
654: + n - the polynomial degree
656: Output Parameter:
657: . norm - the weighted L2 norm
659: Level: beginner
661: .seealso: `PetscDTJacobiEval()`
662: @*/
663: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
664: {
665: PetscReal twoab1;
666: PetscReal gr;
671: twoab1 = PetscPowReal(2., alpha + beta + 1.);
672: #if defined(PETSC_HAVE_LGAMMA)
673: if (!n) {
674: gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.));
675: } else {
676: gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.);
677: }
678: #else
679: {
680: PetscInt alphai = (PetscInt)alpha;
681: PetscInt betai = (PetscInt)beta;
682: PetscInt i;
684: gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.;
685: if ((PetscReal)alphai == alpha) {
686: if (!n) {
687: for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.);
688: gr /= (alpha + beta + 1.);
689: } else {
690: for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.);
691: }
692: } else if ((PetscReal)betai == beta) {
693: if (!n) {
694: for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.);
695: gr /= (alpha + beta + 1.);
696: } else {
697: for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.);
698: }
699: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
700: }
701: #endif
702: *norm = PetscSqrtReal(twoab1 * gr);
703: return 0;
704: }
706: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
707: {
708: PetscReal ak, bk;
709: PetscReal abk1;
710: PetscInt i, l, maxdegree;
712: maxdegree = degrees[ndegree - 1] - k;
713: ak = a + k;
714: bk = b + k;
715: abk1 = a + b + k + 1.;
716: if (maxdegree < 0) {
717: for (i = 0; i < npoints; i++)
718: for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.;
719: return 0;
720: }
721: for (i = 0; i < npoints; i++) {
722: PetscReal pm1, pm2, x;
723: PetscReal cnm1, cnm1x, cnm2;
724: PetscInt j, m;
726: x = points[i];
727: pm2 = 1.;
728: PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2);
729: pm1 = (cnm1 + cnm1x * x);
730: l = 0;
731: while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.;
732: while (l < ndegree && degrees[l] - k == 0) {
733: p[l] = pm2;
734: for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
735: l++;
736: }
737: while (l < ndegree && degrees[l] - k == 1) {
738: p[l] = pm1;
739: for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
740: l++;
741: }
742: for (j = 2; j <= maxdegree; j++) {
743: PetscReal pp;
745: PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2);
746: pp = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2;
747: pm2 = pm1;
748: pm1 = pp;
749: while (l < ndegree && degrees[l] - k == j) {
750: p[l] = pp;
751: for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
752: l++;
753: }
754: }
755: p += ndegree;
756: }
757: return 0;
758: }
760: /*@
761: 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$.
763: Input Parameters:
764: + alpha - the left exponent of the weight
765: . beta - the right exponetn of the weight
766: . npoints - the number of points to evaluate the polynomials at
767: . points - [npoints] array of point coordinates
768: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
769: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
771: Output Parameters:
772: - p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
773: (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
774: (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
775: varying) dimension is the index of the evaluation point.
777: Level: advanced
779: .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()`
780: @*/
781: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
782: {
783: PetscInt i, j, l;
784: PetscInt *degrees;
785: PetscReal *psingle;
787: if (degree == 0) {
788: PetscInt zero = 0;
790: for (i = 0; i <= k; i++) PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints]);
791: return 0;
792: }
793: PetscMalloc1(degree + 1, °rees);
794: PetscMalloc1((degree + 1) * npoints, &psingle);
795: for (i = 0; i <= degree; i++) degrees[i] = i;
796: for (i = 0; i <= k; i++) {
797: PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle);
798: for (j = 0; j <= degree; j++) {
799: for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
800: }
801: }
802: PetscFree(psingle);
803: PetscFree(degrees);
804: return 0;
805: }
807: /*@
808: PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$
809: at points
811: Not Collective
813: Input Parameters:
814: + npoints - number of spatial points to evaluate at
815: . alpha - the left exponent > -1
816: . beta - the right exponent > -1
817: . points - array of locations to evaluate at
818: . ndegree - number of basis degrees to evaluate
819: - degrees - sorted array of degrees to evaluate
821: Output Parameters:
822: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
823: . D - row-oriented derivative evaluation matrix (or NULL)
824: - D2 - row-oriented second derivative evaluation matrix (or NULL)
826: Level: intermediate
828: .seealso: `PetscDTGaussQuadrature()`
829: @*/
830: PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
831: {
834: if (!npoints || !ndegree) return 0;
835: if (B) PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);
836: if (D) PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);
837: if (D2) PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);
838: return 0;
839: }
841: /*@
842: PetscDTLegendreEval - evaluate Legendre polynomials at points
844: Not Collective
846: Input Parameters:
847: + npoints - number of spatial points to evaluate at
848: . points - array of locations to evaluate at
849: . ndegree - number of basis degrees to evaluate
850: - degrees - sorted array of degrees to evaluate
852: Output Parameters:
853: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
854: . D - row-oriented derivative evaluation matrix (or NULL)
855: - D2 - row-oriented second derivative evaluation matrix (or NULL)
857: Level: intermediate
859: .seealso: `PetscDTGaussQuadrature()`
860: @*/
861: PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
862: {
863: PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);
864: return 0;
865: }
867: /*@
868: 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)
870: Input Parameters:
871: + len - the desired length of the degree tuple
872: - index - the index to convert: should be >= 0
874: Output Parameter:
875: . degtup - will be filled with a tuple of degrees
877: Level: beginner
879: Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
880: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
881: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
883: .seealso: `PetscDTGradedOrderToIndex()`
884: @*/
885: PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
886: {
887: PetscInt i, total;
888: PetscInt sum;
893: total = 1;
894: sum = 0;
895: while (index >= total) {
896: index -= total;
897: total = (total * (len + sum)) / (sum + 1);
898: sum++;
899: }
900: for (i = 0; i < len; i++) {
901: PetscInt c;
903: degtup[i] = sum;
904: for (c = 0, total = 1; c < sum; c++) {
905: /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
906: if (index < total) break;
907: index -= total;
908: total = (total * (len - 1 - i + c)) / (c + 1);
909: degtup[i]--;
910: }
911: sum -= degtup[i];
912: }
913: return 0;
914: }
916: /*@
917: PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of PetscDTIndexToGradedOrder().
919: Input Parameters:
920: + len - the length of the degree tuple
921: - degtup - tuple with this length
923: Output Parameter:
924: . index - index in graded order: >= 0
926: Level: Beginner
928: Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
929: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
930: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
932: .seealso: `PetscDTIndexToGradedOrder()`
933: @*/
934: PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
935: {
936: PetscInt i, idx, sum, total;
940: for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
941: idx = 0;
942: total = 1;
943: for (i = 0; i < sum; i++) {
944: idx += total;
945: total = (total * (len + i)) / (i + 1);
946: }
947: for (i = 0; i < len - 1; i++) {
948: PetscInt c;
950: total = 1;
951: sum -= degtup[i];
952: for (c = 0; c < sum; c++) {
953: idx += total;
954: total = (total * (len - 1 - i + c)) / (c + 1);
955: }
956: }
957: *index = idx;
958: return 0;
959: }
961: static PetscBool PKDCite = PETSC_FALSE;
962: const char PKDCitation[] = "@article{Kirby2010,\n"
963: " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
964: " author={Kirby, Robert C},\n"
965: " journal={ACM Transactions on Mathematical Software (TOMS)},\n"
966: " volume={37},\n"
967: " number={1},\n"
968: " pages={1--16},\n"
969: " year={2010},\n"
970: " publisher={ACM New York, NY, USA}\n}\n";
972: /*@
973: PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for
974: the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used
975: as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating
976: polynomials in that domain.
978: Input Parameters:
979: + dim - the number of variables in the multivariate polynomials
980: . npoints - the number of points to evaluate the polynomials at
981: . points - [npoints x dim] array of point coordinates
982: . 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.
983: - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives
984: in the jet. Choosing k = 0 means to evaluate just the function and no derivatives
986: Output Parameters:
987: - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree)
988: choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
989: three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
990: index; the third (fastest varying) dimension is the index of the evaluation point.
992: Level: advanced
994: Note: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
995: ordering of PetscDTIndexToGradedOrder() and PetscDTGradedOrderToIndex(). For example, in 3D, the polynomial with
996: 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);
997: 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).
999: The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006.
1001: .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()`
1002: @*/
1003: PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1004: {
1005: PetscInt degidx, kidx, d, pt;
1006: PetscInt Nk, Ndeg;
1007: PetscInt *ktup, *degtup;
1008: PetscReal *scales, initscale, scaleexp;
1010: PetscCitationsRegister(PKDCitation, &PKDCite);
1011: PetscDTBinomialInt(dim + k, k, &Nk);
1012: PetscDTBinomialInt(degree + dim, degree, &Ndeg);
1013: PetscMalloc2(dim, °tup, dim, &ktup);
1014: PetscMalloc1(Ndeg, &scales);
1015: initscale = 1.;
1016: if (dim > 1) {
1017: PetscDTBinomial(dim, 2, &scaleexp);
1018: initscale = PetscPowReal(2., scaleexp * 0.5);
1019: }
1020: for (degidx = 0; degidx < Ndeg; degidx++) {
1021: PetscInt e, i;
1022: PetscInt m1idx = -1, m2idx = -1;
1023: PetscInt n;
1024: PetscInt degsum;
1025: PetscReal alpha;
1026: PetscReal cnm1, cnm1x, cnm2;
1027: PetscReal norm;
1029: PetscDTIndexToGradedOrder(dim, degidx, degtup);
1030: for (d = dim - 1; d >= 0; d--)
1031: if (degtup[d]) break;
1032: if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1033: scales[degidx] = initscale;
1034: for (e = 0; e < dim; e++) {
1035: PetscDTJacobiNorm(e, 0., 0, &norm);
1036: scales[degidx] /= norm;
1037: }
1038: for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1039: for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1040: continue;
1041: }
1042: n = degtup[d];
1043: degtup[d]--;
1044: PetscDTGradedOrderToIndex(dim, degtup, &m1idx);
1045: if (degtup[d] > 0) {
1046: degtup[d]--;
1047: PetscDTGradedOrderToIndex(dim, degtup, &m2idx);
1048: degtup[d]++;
1049: }
1050: degtup[d]++;
1051: for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1052: alpha = 2 * degsum + d;
1053: PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2);
1055: scales[degidx] = initscale;
1056: for (e = 0, degsum = 0; e < dim; e++) {
1057: PetscInt f;
1058: PetscReal ealpha;
1059: PetscReal enorm;
1061: ealpha = 2 * degsum + e;
1062: for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1063: PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm);
1064: scales[degidx] /= enorm;
1065: degsum += degtup[e];
1066: }
1068: for (pt = 0; pt < npoints; pt++) {
1069: /* compute the multipliers */
1070: PetscReal thetanm1, thetanm1x, thetanm2;
1072: thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d];
1073: for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e];
1074: thetanm1x *= 0.5;
1075: thetanm1 = (2. - (dim - (d + 1)));
1076: for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1077: thetanm1 *= 0.5;
1078: thetanm2 = thetanm1 * thetanm1;
1080: for (kidx = 0; kidx < Nk; kidx++) {
1081: PetscInt f;
1083: PetscDTIndexToGradedOrder(dim, kidx, ktup);
1084: /* first sum in the same derivative terms */
1085: p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1086: if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];
1088: for (f = d; f < dim; f++) {
1089: PetscInt km1idx, mplty = ktup[f];
1091: if (!mplty) continue;
1092: ktup[f]--;
1093: PetscDTGradedOrderToIndex(dim, ktup, &km1idx);
1095: /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1096: /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1097: /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1098: if (f > d) {
1099: PetscInt f2;
1101: p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1102: if (m2idx >= 0) {
1103: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1104: /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1105: for (f2 = f; f2 < dim; f2++) {
1106: PetscInt km2idx, mplty2 = ktup[f2];
1107: PetscInt factor;
1109: if (!mplty2) continue;
1110: ktup[f2]--;
1111: PetscDTGradedOrderToIndex(dim, ktup, &km2idx);
1113: factor = mplty * mplty2;
1114: if (f == f2) factor /= 2;
1115: p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1116: ktup[f2]++;
1117: }
1118: }
1119: } else {
1120: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1121: }
1122: ktup[f]++;
1123: }
1124: }
1125: }
1126: }
1127: for (degidx = 0; degidx < Ndeg; degidx++) {
1128: PetscReal scale = scales[degidx];
1129: PetscInt i;
1131: for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale;
1132: }
1133: PetscFree(scales);
1134: PetscFree2(degtup, ktup);
1135: return 0;
1136: }
1138: /*@
1139: PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree,
1140: which can be evaluated in PetscDTPTrimmedEvalJet().
1142: Input Parameters:
1143: + dim - the number of variables in the multivariate polynomials
1144: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space.
1145: - formDegree - the degree of the form
1147: Output Parameters:
1148: - size - The number ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree))
1150: Level: advanced
1152: .seealso: `PetscDTPTrimmedEvalJet()`
1153: @*/
1154: PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size)
1155: {
1156: PetscInt Nrk, Nbpt; // number of trimmed polynomials
1158: formDegree = PetscAbsInt(formDegree);
1159: PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt);
1160: PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk);
1161: Nbpt *= Nrk;
1162: *size = Nbpt;
1163: return 0;
1164: }
1166: /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it
1167: * was inferior to this implementation */
1168: static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1169: {
1170: PetscInt formDegreeOrig = formDegree;
1171: PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE;
1173: formDegree = PetscAbsInt(formDegreeOrig);
1174: if (formDegree == 0) {
1175: PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p);
1176: return 0;
1177: }
1178: if (formDegree == dim) {
1179: PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p);
1180: return 0;
1181: }
1182: PetscInt Nbpt;
1183: PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt);
1184: PetscInt Nf;
1185: PetscDTBinomialInt(dim, formDegree, &Nf);
1186: PetscInt Nk;
1187: PetscDTBinomialInt(dim + jetDegree, dim, &Nk);
1188: PetscArrayzero(p, Nbpt * Nf * Nk * npoints);
1190: PetscInt Nbpm1; // number of scalar polynomials up to degree - 1;
1191: PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1);
1192: PetscReal *p_scalar;
1193: PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar);
1194: PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar);
1195: PetscInt total = 0;
1196: // First add the full polynomials up to degree - 1 into the basis: take the scalar
1197: // and copy one for each form component
1198: for (PetscInt i = 0; i < Nbpm1; i++) {
1199: const PetscReal *src = &p_scalar[i * Nk * npoints];
1200: for (PetscInt f = 0; f < Nf; f++) {
1201: PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints];
1202: PetscArraycpy(dest, src, Nk * npoints);
1203: }
1204: }
1205: PetscInt *form_atoms;
1206: PetscMalloc1(formDegree + 1, &form_atoms);
1207: // construct the interior product pattern
1208: PetscInt(*pattern)[3];
1209: PetscInt Nf1; // number of formDegree + 1 forms
1210: PetscDTBinomialInt(dim, formDegree + 1, &Nf1);
1211: PetscInt nnz = Nf1 * (formDegree + 1);
1212: PetscMalloc1(Nf1 * (formDegree + 1), &pattern);
1213: PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern);
1214: PetscReal centroid = (1. - dim) / (dim + 1.);
1215: PetscInt *deriv;
1216: PetscMalloc1(dim, &deriv);
1217: for (PetscInt d = dim; d >= formDegree + 1; d--) {
1218: PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0
1219: // (equal to the number of formDegree forms in dimension d-1)
1220: PetscDTBinomialInt(d - 1, formDegree, &Nfd1);
1221: // The number of homogeneous (degree-1) scalar polynomials in d variables
1222: PetscInt Nh;
1223: PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh);
1224: const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints];
1225: for (PetscInt b = 0; b < Nh; b++) {
1226: const PetscReal *h_s = &h_scalar[b * Nk * npoints];
1227: for (PetscInt f = 0; f < Nfd1; f++) {
1228: // construct all formDegree+1 forms that start with dx_(dim - d) /\ ...
1229: form_atoms[0] = dim - d;
1230: PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1]);
1231: for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1;
1232: PetscInt f_ind; // index of the resulting form
1233: PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind);
1234: PetscReal *p_f = &p[total++ * Nf * Nk * npoints];
1235: for (PetscInt nz = 0; nz < nnz; nz++) {
1236: PetscInt i = pattern[nz][0]; // formDegree component
1237: PetscInt j = pattern[nz][1]; // (formDegree + 1) component
1238: PetscInt v = pattern[nz][2]; // coordinate component
1239: PetscReal scale = v < 0 ? -1. : 1.;
1241: i = formNegative ? (Nf - 1 - i) : i;
1242: scale = (formNegative && (i & 1)) ? -scale : scale;
1243: v = v < 0 ? -(v + 1) : v;
1244: if (j != f_ind) continue;
1245: PetscReal *p_i = &p_f[i * Nk * npoints];
1246: for (PetscInt jet = 0; jet < Nk; jet++) {
1247: const PetscReal *h_jet = &h_s[jet * npoints];
1248: PetscReal *p_jet = &p_i[jet * npoints];
1250: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid);
1251: PetscDTIndexToGradedOrder(dim, jet, deriv);
1252: deriv[v]++;
1253: PetscReal mult = deriv[v];
1254: PetscInt l;
1255: PetscDTGradedOrderToIndex(dim, deriv, &l);
1256: if (l >= Nk) continue;
1257: p_jet = &p_i[l * npoints];
1258: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt];
1259: deriv[v]--;
1260: }
1261: }
1262: }
1263: }
1264: }
1266: PetscFree(deriv);
1267: PetscFree(pattern);
1268: PetscFree(form_atoms);
1269: PetscFree(p_scalar);
1270: return 0;
1271: }
1273: /*@
1274: PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to
1275: a given degree.
1277: Input Parameters:
1278: + dim - the number of variables in the multivariate polynomials
1279: . npoints - the number of points to evaluate the polynomials at
1280: . points - [npoints x dim] array of point coordinates
1281: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate.
1282: There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space.
1283: (You can use PetscDTPTrimmedSize() to compute this size.)
1284: . formDegree - the degree of the form
1285: - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives
1286: in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives
1288: Output Parameters:
1289: - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is
1290: PetscDTPTrimmedSize() x ((dim + formDegree) choose dim) x ((dim + k) choose dim) x npoints,
1291: which also describes the order of the dimensions of this
1292: four-dimensional array:
1293: the first (slowest varying) dimension is basis function index;
1294: the second dimension is component of the form;
1295: the third dimension is jet index;
1296: the fourth (fastest varying) dimension is the index of the evaluation point.
1298: Level: advanced
1300: Note: The ordering of the basis functions is not graded, so the basis functions are not nested by degree like PetscDTPKDEvalJet().
1301: The basis functions are not an L2-orthonormal basis on any particular domain.
1303: The implementation is based on the description of the trimmed polynomials up to degree r as
1304: the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to
1305: homogeneous polynomials of degree (r-1).
1307: .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()`
1308: @*/
1309: PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1310: {
1311: PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p);
1312: return 0;
1313: }
1315: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1316: * with lds n; diag and subdiag are overwritten */
1317: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[])
1318: {
1319: char jobz = 'V'; /* eigenvalues and eigenvectors */
1320: char range = 'A'; /* all eigenvalues will be found */
1321: PetscReal VL = 0.; /* ignored because range is 'A' */
1322: PetscReal VU = 0.; /* ignored because range is 'A' */
1323: PetscBLASInt IL = 0; /* ignored because range is 'A' */
1324: PetscBLASInt IU = 0; /* ignored because range is 'A' */
1325: PetscReal abstol = 0.; /* unused */
1326: PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1327: PetscBLASInt *isuppz;
1328: PetscBLASInt lwork, liwork;
1329: PetscReal workquery;
1330: PetscBLASInt iworkquery;
1331: PetscBLASInt *iwork;
1332: PetscBLASInt info;
1333: PetscReal *work = NULL;
1335: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1336: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1337: #endif
1338: PetscBLASIntCast(n, &bn);
1339: PetscBLASIntCast(n, &ldz);
1340: #if !defined(PETSC_MISSING_LAPACK_STEGR)
1341: PetscMalloc1(2 * n, &isuppz);
1342: lwork = -1;
1343: liwork = -1;
1344: PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info));
1346: lwork = (PetscBLASInt)workquery;
1347: liwork = (PetscBLASInt)iworkquery;
1348: PetscMalloc2(lwork, &work, liwork, &iwork);
1349: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1350: PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info));
1351: PetscFPTrapPop();
1353: PetscFree2(work, iwork);
1354: PetscFree(isuppz);
1355: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1356: jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1357: tridiagonal matrix. Z is initialized to the identity
1358: matrix. */
1359: PetscMalloc1(PetscMax(1, 2 * n - 2), &work);
1360: PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info));
1361: PetscFPTrapPop();
1363: PetscFree(work);
1364: PetscArraycpy(eigs, diag, n);
1365: #endif
1366: return 0;
1367: }
1369: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1370: * quadrature rules on the interval [-1, 1] */
1371: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1372: {
1373: PetscReal twoab1;
1374: PetscInt m = n - 2;
1375: PetscReal a = alpha + 1.;
1376: PetscReal b = beta + 1.;
1377: PetscReal gra, grb;
1379: twoab1 = PetscPowReal(2., a + b - 1.);
1380: #if defined(PETSC_HAVE_LGAMMA)
1381: grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.)));
1382: gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.)));
1383: #else
1384: {
1385: PetscInt alphai = (PetscInt)alpha;
1386: PetscInt betai = (PetscInt)beta;
1388: if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) {
1389: PetscReal binom1, binom2;
1391: PetscDTBinomial(m + b, b, &binom1);
1392: PetscDTBinomial(m + a + b, b, &binom2);
1393: grb = 1. / (binom1 * binom2);
1394: PetscDTBinomial(m + a, a, &binom1);
1395: PetscDTBinomial(m + a + b, a, &binom2);
1396: gra = 1. / (binom1 * binom2);
1397: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1398: }
1399: #endif
1400: *leftw = twoab1 * grb / b;
1401: *rightw = twoab1 * gra / a;
1402: return 0;
1403: }
1405: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1406: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1407: static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1408: {
1409: PetscReal pn1, pn2;
1410: PetscReal cnm1, cnm1x, cnm2;
1411: PetscInt k;
1413: if (!n) {
1414: *P = 1.0;
1415: return 0;
1416: }
1417: PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2);
1418: pn2 = 1.;
1419: pn1 = cnm1 + cnm1x * x;
1420: if (n == 1) {
1421: *P = pn1;
1422: return 0;
1423: }
1424: *P = 0.0;
1425: for (k = 2; k < n + 1; ++k) {
1426: PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2);
1428: *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2;
1429: pn2 = pn1;
1430: pn1 = *P;
1431: }
1432: return 0;
1433: }
1435: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1436: static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1437: {
1438: PetscReal nP;
1439: PetscInt i;
1441: *P = 0.0;
1442: if (k > n) return 0;
1443: PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP);
1444: for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1445: *P = nP;
1446: return 0;
1447: }
1449: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1450: {
1451: PetscInt maxIter = 100;
1452: PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1453: PetscReal a1, a6, gf;
1454: PetscInt k;
1457: a1 = PetscPowReal(2.0, a + b + 1);
1458: #if defined(PETSC_HAVE_LGAMMA)
1459: {
1460: PetscReal a2, a3, a4, a5;
1461: a2 = PetscLGamma(a + npoints + 1);
1462: a3 = PetscLGamma(b + npoints + 1);
1463: a4 = PetscLGamma(a + b + npoints + 1);
1464: a5 = PetscLGamma(npoints + 1);
1465: gf = PetscExpReal(a2 + a3 - (a4 + a5));
1466: }
1467: #else
1468: {
1469: PetscInt ia, ib;
1471: ia = (PetscInt)a;
1472: ib = (PetscInt)b;
1473: gf = 1.;
1474: if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1475: for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1476: } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1477: for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1478: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1479: }
1480: #endif
1482: a6 = a1 * gf;
1483: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1484: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1485: for (k = 0; k < npoints; ++k) {
1486: PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP;
1487: PetscInt j;
1489: if (k > 0) r = 0.5 * (r + x[k - 1]);
1490: for (j = 0; j < maxIter; ++j) {
1491: PetscReal s = 0.0, delta, f, fp;
1492: PetscInt i;
1494: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1495: PetscDTComputeJacobi(a, b, npoints, r, &f);
1496: PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);
1497: delta = f / (fp - f * s);
1498: r = r - delta;
1499: if (PetscAbsReal(delta) < eps) break;
1500: }
1501: x[k] = r;
1502: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);
1503: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1504: }
1505: return 0;
1506: }
1508: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1509: * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1510: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1511: {
1512: PetscInt i;
1514: for (i = 0; i < nPoints; i++) {
1515: PetscReal A, B, C;
1517: PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C);
1518: d[i] = -A / B;
1519: if (i) s[i - 1] *= C / B;
1520: if (i < nPoints - 1) s[i] = 1. / B;
1521: }
1522: return 0;
1523: }
1525: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1526: {
1527: PetscReal mu0;
1528: PetscReal ga, gb, gab;
1529: PetscInt i;
1531: PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);
1533: #if defined(PETSC_HAVE_TGAMMA)
1534: ga = PetscTGamma(a + 1);
1535: gb = PetscTGamma(b + 1);
1536: gab = PetscTGamma(a + b + 2);
1537: #else
1538: {
1539: PetscInt ia, ib;
1541: ia = (PetscInt)a;
1542: ib = (PetscInt)b;
1543: if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1544: PetscDTFactorial(ia, &ga);
1545: PetscDTFactorial(ib, &gb);
1546: PetscDTFactorial(ia + ib + 1, &gb);
1547: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable.");
1548: }
1549: #endif
1550: mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab;
1552: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1553: {
1554: PetscReal *diag, *subdiag;
1555: PetscScalar *V;
1557: PetscMalloc2(npoints, &diag, npoints, &subdiag);
1558: PetscMalloc1(npoints * npoints, &V);
1559: PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);
1560: for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1561: PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);
1562: for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1563: PetscFree(V);
1564: PetscFree2(diag, subdiag);
1565: }
1566: #else
1567: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1568: #endif
1569: { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1570: eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1571: the eigenvalues are sorted */
1572: PetscBool sorted;
1574: PetscSortedReal(npoints, x, &sorted);
1575: if (!sorted) {
1576: PetscInt *order, i;
1577: PetscReal *tmp;
1579: PetscMalloc2(npoints, &order, npoints, &tmp);
1580: for (i = 0; i < npoints; i++) order[i] = i;
1581: PetscSortRealWithPermutation(npoints, x, order);
1582: PetscArraycpy(tmp, x, npoints);
1583: for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1584: PetscArraycpy(tmp, w, npoints);
1585: for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1586: PetscFree2(order, tmp);
1587: }
1588: }
1589: return 0;
1590: }
1592: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1593: {
1595: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1599: if (newton) PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);
1600: else PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);
1601: if (alpha == beta) { /* symmetrize */
1602: PetscInt i;
1603: for (i = 0; i < (npoints + 1) / 2; i++) {
1604: PetscInt j = npoints - 1 - i;
1605: PetscReal xi = x[i];
1606: PetscReal xj = x[j];
1607: PetscReal wi = w[i];
1608: PetscReal wj = w[j];
1610: x[i] = (xi - xj) / 2.;
1611: x[j] = (xj - xi) / 2.;
1612: w[i] = w[j] = (wi + wj) / 2.;
1613: }
1614: }
1615: return 0;
1616: }
1618: /*@
1619: PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1620: $(x-a)^\alpha (x-b)^\beta$.
1622: Not collective
1624: Input Parameters:
1625: + npoints - the number of points in the quadrature rule
1626: . a - the left endpoint of the interval
1627: . b - the right endpoint of the interval
1628: . alpha - the left exponent
1629: - beta - the right exponent
1631: Output Parameters:
1632: + x - array of length npoints, the locations of the quadrature points
1633: - w - array of length npoints, the weights of the quadrature points
1635: Level: intermediate
1637: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1638: @*/
1639: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1640: {
1641: PetscInt i;
1643: PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1644: if (a != -1. || b != 1.) { /* shift */
1645: for (i = 0; i < npoints; i++) {
1646: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1647: w[i] *= (b - a) / 2.;
1648: }
1649: }
1650: return 0;
1651: }
1653: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1654: {
1655: PetscInt i;
1658: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1662: x[0] = -1.;
1663: x[npoints - 1] = 1.;
1664: if (npoints > 2) PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton);
1665: for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]);
1666: PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1]);
1667: return 0;
1668: }
1670: /*@
1671: PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1672: $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.
1674: Not collective
1676: Input Parameters:
1677: + npoints - the number of points in the quadrature rule
1678: . a - the left endpoint of the interval
1679: . b - the right endpoint of the interval
1680: . alpha - the left exponent
1681: - beta - the right exponent
1683: Output Parameters:
1684: + x - array of length npoints, the locations of the quadrature points
1685: - w - array of length npoints, the weights of the quadrature points
1687: Level: intermediate
1689: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1690: @*/
1691: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1692: {
1693: PetscInt i;
1695: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1696: if (a != -1. || b != 1.) { /* shift */
1697: for (i = 0; i < npoints; i++) {
1698: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1699: w[i] *= (b - a) / 2.;
1700: }
1701: }
1702: return 0;
1703: }
1705: /*@
1706: PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1708: Not Collective
1710: Input Parameters:
1711: + npoints - number of points
1712: . a - left end of interval (often-1)
1713: - b - right end of interval (often +1)
1715: Output Parameters:
1716: + x - quadrature points
1717: - w - quadrature weights
1719: Level: intermediate
1721: References:
1722: . * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
1724: .seealso: `PetscDTLegendreEval()`
1725: @*/
1726: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1727: {
1728: PetscInt i;
1730: PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);
1731: if (a != -1. || b != 1.) { /* shift */
1732: for (i = 0; i < npoints; i++) {
1733: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1734: w[i] *= (b - a) / 2.;
1735: }
1736: }
1737: return 0;
1738: }
1740: /*@C
1741: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1742: nodes of a given size on the domain [-1,1]
1744: Not Collective
1746: Input Parameters:
1747: + n - number of grid nodes
1748: - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON
1750: Output Parameters:
1751: + x - quadrature points
1752: - w - quadrature weights
1754: Notes:
1755: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1756: close enough to the desired solution
1758: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1760: 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
1762: Level: intermediate
1764: .seealso: `PetscDTGaussQuadrature()`
1766: @*/
1767: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal *x, PetscReal *w)
1768: {
1769: PetscBool newton;
1772: newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1773: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);
1774: return 0;
1775: }
1777: /*@
1778: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1780: Not Collective
1782: Input Parameters:
1783: + dim - The spatial dimension
1784: . Nc - The number of components
1785: . npoints - number of points in one dimension
1786: . a - left end of interval (often-1)
1787: - b - right end of interval (often +1)
1789: Output Parameter:
1790: . q - A PetscQuadrature object
1792: Level: intermediate
1794: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
1795: @*/
1796: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1797: {
1798: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
1799: PetscReal *x, *w, *xw, *ww;
1801: PetscMalloc1(totpoints * dim, &x);
1802: PetscMalloc1(totpoints * Nc, &w);
1803: /* Set up the Golub-Welsch system */
1804: switch (dim) {
1805: case 0:
1806: PetscFree(x);
1807: PetscFree(w);
1808: PetscMalloc1(1, &x);
1809: PetscMalloc1(Nc, &w);
1810: x[0] = 0.0;
1811: for (c = 0; c < Nc; ++c) w[c] = 1.0;
1812: break;
1813: case 1:
1814: PetscMalloc1(npoints, &ww);
1815: PetscDTGaussQuadrature(npoints, a, b, x, ww);
1816: for (i = 0; i < npoints; ++i)
1817: for (c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i];
1818: PetscFree(ww);
1819: break;
1820: case 2:
1821: PetscMalloc2(npoints, &xw, npoints, &ww);
1822: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1823: for (i = 0; i < npoints; ++i) {
1824: for (j = 0; j < npoints; ++j) {
1825: x[(i * npoints + j) * dim + 0] = xw[i];
1826: x[(i * npoints + j) * dim + 1] = xw[j];
1827: for (c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j];
1828: }
1829: }
1830: PetscFree2(xw, ww);
1831: break;
1832: case 3:
1833: PetscMalloc2(npoints, &xw, npoints, &ww);
1834: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1835: for (i = 0; i < npoints; ++i) {
1836: for (j = 0; j < npoints; ++j) {
1837: for (k = 0; k < npoints; ++k) {
1838: x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i];
1839: x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j];
1840: x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k];
1841: for (c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k];
1842: }
1843: }
1844: }
1845: PetscFree2(xw, ww);
1846: break;
1847: default:
1848: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim);
1849: }
1850: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1851: PetscQuadratureSetOrder(*q, 2 * npoints - 1);
1852: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1853: PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor");
1854: return 0;
1855: }
1857: /*@
1858: PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex
1860: Not Collective
1862: Input Parameters:
1863: + dim - The simplex dimension
1864: . Nc - The number of components
1865: . npoints - The number of points in one dimension
1866: . a - left end of interval (often-1)
1867: - b - right end of interval (often +1)
1869: Output Parameter:
1870: . q - A PetscQuadrature object
1872: Level: intermediate
1874: References:
1875: . * - Karniadakis and Sherwin. FIAT
1877: Note: For dim == 1, this is Gauss-Legendre quadrature
1879: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()`
1880: @*/
1881: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1882: {
1883: PetscInt totprev, totrem;
1884: PetscInt totpoints;
1885: PetscReal *p1, *w1;
1886: PetscReal *x, *w;
1887: PetscInt i, j, k, l, m, pt, c;
1890: totpoints = 1;
1891: for (i = 0, totpoints = 1; i < dim; i++) totpoints *= npoints;
1892: PetscMalloc1(totpoints * dim, &x);
1893: PetscMalloc1(totpoints * Nc, &w);
1894: PetscMalloc2(npoints, &p1, npoints, &w1);
1895: for (i = 0; i < totpoints * Nc; i++) w[i] = 1.;
1896: for (i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; i++) {
1897: PetscReal mul;
1899: mul = PetscPowReal(2., -i);
1900: PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1);
1901: for (pt = 0, l = 0; l < totprev; l++) {
1902: for (j = 0; j < npoints; j++) {
1903: for (m = 0; m < totrem; m++, pt++) {
1904: for (k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.;
1905: x[pt * dim + i] = p1[j];
1906: for (c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j];
1907: }
1908: }
1909: }
1910: totprev *= npoints;
1911: totrem /= npoints;
1912: }
1913: PetscFree2(p1, w1);
1914: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1915: PetscQuadratureSetOrder(*q, 2 * npoints - 1);
1916: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1917: PetscObjectChangeTypeName((PetscObject)*q, "StroudConical");
1918: return 0;
1919: }
1921: static PetscBool MinSymTriQuadCite = PETSC_FALSE;
1922: const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n"
1923: " title = {On the identification of symmetric quadrature rules for finite element methods},\n"
1924: " journal = {Computers & Mathematics with Applications},\n"
1925: " volume = {69},\n"
1926: " number = {10},\n"
1927: " pages = {1232-1241},\n"
1928: " year = {2015},\n"
1929: " issn = {0898-1221},\n"
1930: " doi = {10.1016/j.camwa.2015.03.017},\n"
1931: " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n"
1932: " author = {F.D. Witherden and P.E. Vincent},\n"
1933: "}\n";
1935: #include "petscdttriquadrules.h"
1937: static PetscBool MinSymTetQuadCite = PETSC_FALSE;
1938: const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n"
1939: " author = {Jaskowiec, Jan and Sukumar, N.},\n"
1940: " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n"
1941: " journal = {International Journal for Numerical Methods in Engineering},\n"
1942: " volume = {122},\n"
1943: " number = {1},\n"
1944: " pages = {148-171},\n"
1945: " doi = {10.1002/nme.6528},\n"
1946: " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n"
1947: " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n"
1948: " year = {2021}\n"
1949: "}\n";
1951: #include "petscdttetquadrules.h"
1953: // https://en.wikipedia.org/wiki/Partition_(number_theory)
1954: static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p)
1955: {
1956: // sequence A000041 in the OEIS
1957: const PetscInt partition[] = {1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627, 792, 1002, 1255, 1575, 1958, 2436, 3010, 3718, 4565, 5604};
1958: PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1;
1961: // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high
1963: *p = partition[n];
1964: return 0;
1965: }
1967: /*@
1968: PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree.
1970: Not Collective
1972: Input Parameters:
1973: + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron)
1974: . degree - The largest polynomial degree that is required to be integrated exactly
1975: - type - left end of interval (often-1)
1977: Output Parameter:
1978: . quad - A PetscQuadrature object for integration over the biunit simplex
1979: (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for
1980: polynomials up to the given degree
1982: Level: intermediate
1984: .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`
1985: @*/
1986: PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad)
1987: {
1988: PetscDTSimplexQuadratureType orig_type = type;
1992: if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM;
1993: if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) {
1994: PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2);
1995: PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad);
1996: } else {
1997: PetscInt n = dim + 1;
1998: PetscInt fact = 1;
1999: PetscInt *part, *perm;
2000: PetscInt p = 0;
2001: PetscInt max_degree;
2002: const PetscInt *nodes_per_type = NULL;
2003: const PetscInt *all_num_full_nodes = NULL;
2004: const PetscReal **weights_list = NULL;
2005: const PetscReal **compact_nodes_list = NULL;
2006: const char *citation = NULL;
2007: PetscBool *cited = NULL;
2009: switch (dim) {
2010: case 2:
2011: cited = &MinSymTriQuadCite;
2012: citation = MinSymTriQuadCitation;
2013: max_degree = PetscDTWVTriQuad_max_degree;
2014: nodes_per_type = PetscDTWVTriQuad_num_orbits;
2015: all_num_full_nodes = PetscDTWVTriQuad_num_nodes;
2016: weights_list = PetscDTWVTriQuad_weights;
2017: compact_nodes_list = PetscDTWVTriQuad_orbits;
2018: break;
2019: case 3:
2020: cited = &MinSymTetQuadCite;
2021: citation = MinSymTetQuadCitation;
2022: max_degree = PetscDTJSTetQuad_max_degree;
2023: nodes_per_type = PetscDTJSTetQuad_num_orbits;
2024: all_num_full_nodes = PetscDTJSTetQuad_num_nodes;
2025: weights_list = PetscDTJSTetQuad_weights;
2026: compact_nodes_list = PetscDTJSTetQuad_orbits;
2027: break;
2028: default:
2029: max_degree = -1;
2030: break;
2031: }
2033: if (degree > max_degree) {
2034: if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) {
2035: // fall back to conic
2036: PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad);
2037: return 0;
2038: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Minimal symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", dim, degree);
2039: }
2041: PetscCitationsRegister(citation, cited);
2043: PetscDTPartitionNumber(n, &p);
2044: for (PetscInt d = 2; d <= n; d++) fact *= d;
2046: PetscInt num_full_nodes = all_num_full_nodes[degree];
2047: const PetscReal *all_compact_nodes = compact_nodes_list[degree];
2048: const PetscReal *all_compact_weights = weights_list[degree];
2049: nodes_per_type = &nodes_per_type[p * degree];
2051: PetscReal *points;
2052: PetscReal *counts;
2053: PetscReal *weights;
2054: PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit
2055: PetscQuadrature q;
2057: // compute the transformation
2058: PetscMalloc1(n * dim, &bary_to_biunit);
2059: for (PetscInt d = 0; d < dim; d++) {
2060: for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0;
2061: }
2063: PetscMalloc3(n, &part, n, &perm, n, &counts);
2064: PetscCalloc1(num_full_nodes * dim, &points);
2065: PetscMalloc1(num_full_nodes, &weights);
2067: // (0, 0, ...) is the first partition lexicographically
2068: PetscArrayzero(part, n);
2069: PetscArrayzero(counts, n);
2070: counts[0] = n;
2072: // for each partition
2073: for (PetscInt s = 0, node_offset = 0; s < p; s++) {
2074: PetscInt num_compact_coords = part[n - 1] + 1;
2076: const PetscReal *compact_nodes = all_compact_nodes;
2077: const PetscReal *compact_weights = all_compact_weights;
2078: all_compact_nodes += num_compact_coords * nodes_per_type[s];
2079: all_compact_weights += nodes_per_type[s];
2081: // for every permutation of the vertices
2082: for (PetscInt f = 0; f < fact; f++) {
2083: PetscDTEnumPerm(n, f, perm, NULL);
2085: // check if it is a valid permutation
2086: PetscInt digit;
2087: for (digit = 1; digit < n; digit++) {
2088: // skip permutations that would duplicate a node because it has a smaller symmetry group
2089: if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break;
2090: }
2091: if (digit < n) continue;
2093: // create full nodes from this permutation of the compact nodes
2094: PetscReal *full_nodes = &points[node_offset * dim];
2095: PetscReal *full_weights = &weights[node_offset];
2097: PetscArraycpy(full_weights, compact_weights, nodes_per_type[s]);
2098: for (PetscInt b = 0; b < n; b++) {
2099: for (PetscInt d = 0; d < dim; d++) {
2100: for (PetscInt node = 0; node < nodes_per_type[s]; node++) full_nodes[node * dim + d] += bary_to_biunit[d * n + perm[b]] * compact_nodes[node * num_compact_coords + part[b]];
2101: }
2102: }
2103: node_offset += nodes_per_type[s];
2104: }
2106: if (s < p - 1) { // Generate the next partition
2107: /* A partition is described by the number of coordinates that are in
2108: * each set of duplicates (counts) and redundantly by mapping each
2109: * index to its set of duplicates (part)
2110: *
2111: * Counts should always be in nonincreasing order
2112: *
2113: * We want to generate the partitions lexically by part, which means
2114: * finding the last index where count > 1 and reducing by 1.
2115: *
2116: * For the new counts beyond that index, we eagerly assign the remaining
2117: * capacity of the partition to smaller indices (ensures lexical ordering),
2118: * while respecting the nonincreasing invariant of the counts
2119: */
2120: PetscInt last_digit = part[n - 1];
2121: PetscInt last_digit_with_extra = last_digit;
2122: while (counts[last_digit_with_extra] == 1) last_digit_with_extra--;
2123: PetscInt limit = --counts[last_digit_with_extra];
2124: PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1;
2125: for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) {
2126: counts[digit] = PetscMin(limit, total_to_distribute);
2127: total_to_distribute -= PetscMin(limit, total_to_distribute);
2128: }
2129: for (PetscInt digit = 0, offset = 0; digit < n; digit++) {
2130: PetscInt count = counts[digit];
2131: for (PetscInt c = 0; c < count; c++) part[offset++] = digit;
2132: }
2133: }
2134: }
2135: PetscFree3(part, perm, counts);
2136: PetscFree(bary_to_biunit);
2137: PetscQuadratureCreate(PETSC_COMM_SELF, &q);
2138: PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights);
2139: *quad = q;
2140: }
2141: return 0;
2142: }
2144: /*@
2145: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
2147: Not Collective
2149: Input Parameters:
2150: + dim - The cell dimension
2151: . level - The number of points in one dimension, 2^l
2152: . a - left end of interval (often-1)
2153: - b - right end of interval (often +1)
2155: Output Parameter:
2156: . q - A PetscQuadrature object
2158: Level: intermediate
2160: .seealso: `PetscDTGaussTensorQuadrature()`
2161: @*/
2162: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
2163: {
2164: const PetscInt p = 16; /* Digits of precision in the evaluation */
2165: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2166: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2167: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
2168: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
2169: PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */
2170: PetscReal *x, *w;
2171: PetscInt K, k, npoints;
2175: /* Find K such that the weights are < 32 digits of precision */
2176: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2 * p; ++K) wk = 0.5 * h * PETSC_PI * PetscCoshReal(K * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(K * h)));
2177: PetscQuadratureCreate(PETSC_COMM_SELF, q);
2178: PetscQuadratureSetOrder(*q, 2 * K + 1);
2179: npoints = 2 * K - 1;
2180: PetscMalloc1(npoints * dim, &x);
2181: PetscMalloc1(npoints, &w);
2182: /* Center term */
2183: x[0] = beta;
2184: w[0] = 0.5 * alpha * PETSC_PI;
2185: for (k = 1; k < K; ++k) {
2186: wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2187: xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h));
2188: x[2 * k - 1] = -alpha * xk + beta;
2189: w[2 * k - 1] = wk;
2190: x[2 * k + 0] = alpha * xk + beta;
2191: w[2 * k + 0] = wk;
2192: }
2193: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
2194: return 0;
2195: }
2197: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2198: {
2199: const PetscInt p = 16; /* Digits of precision in the evaluation */
2200: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2201: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2202: PetscReal h = 1.0; /* Step size, length between x_k */
2203: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2204: PetscReal osum = 0.0; /* Integral on last level */
2205: PetscReal psum = 0.0; /* Integral on the level before the last level */
2206: PetscReal sum; /* Integral on current level */
2207: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2208: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2209: PetscReal wk; /* Quadrature weight at x_k */
2210: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
2211: PetscInt d; /* Digits of precision in the integral */
2214: /* Center term */
2215: func(&beta, ctx, &lval);
2216: sum = 0.5 * alpha * PETSC_PI * lval;
2217: /* */
2218: do {
2219: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
2220: PetscInt k = 1;
2222: ++l;
2223: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2224: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2225: psum = osum;
2226: osum = sum;
2227: h *= 0.5;
2228: sum *= 0.5;
2229: do {
2230: wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2231: yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2232: lx = -alpha * (1.0 - yk) + beta;
2233: rx = alpha * (1.0 - yk) + beta;
2234: func(&lx, ctx, &lval);
2235: func(&rx, ctx, &rval);
2236: lterm = alpha * wk * lval;
2237: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
2238: sum += lterm;
2239: rterm = alpha * wk * rval;
2240: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
2241: sum += rterm;
2242: ++k;
2243: /* Only need to evaluate every other point on refined levels */
2244: if (l != 1) ++k;
2245: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
2247: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
2248: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
2249: d3 = PetscLog10Real(maxTerm) - p;
2250: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
2251: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
2252: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2253: } while (d < digits && l < 12);
2254: *sol = sum;
2256: return 0;
2257: }
2259: #if defined(PETSC_HAVE_MPFR)
2260: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2261: {
2262: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
2263: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2264: mpfr_t alpha; /* Half-width of the integration interval */
2265: mpfr_t beta; /* Center of the integration interval */
2266: mpfr_t h; /* Step size, length between x_k */
2267: mpfr_t osum; /* Integral on last level */
2268: mpfr_t psum; /* Integral on the level before the last level */
2269: mpfr_t sum; /* Integral on current level */
2270: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2271: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2272: mpfr_t wk; /* Quadrature weight at x_k */
2273: PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */
2274: PetscInt d; /* Digits of precision in the integral */
2275: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
2278: /* Create high precision storage */
2279: 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);
2280: /* Initialization */
2281: mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN);
2282: mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN);
2283: mpfr_set_d(osum, 0.0, MPFR_RNDN);
2284: mpfr_set_d(psum, 0.0, MPFR_RNDN);
2285: mpfr_set_d(h, 1.0, MPFR_RNDN);
2286: mpfr_const_pi(pi2, MPFR_RNDN);
2287: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
2288: /* Center term */
2289: rtmp = 0.5 * (b + a);
2290: func(&rtmp, ctx, &lval);
2291: mpfr_set(sum, pi2, MPFR_RNDN);
2292: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
2293: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
2294: /* */
2295: do {
2296: PetscReal d1, d2, d3, d4;
2297: PetscInt k = 1;
2299: ++l;
2300: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
2301: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2302: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2303: mpfr_set(psum, osum, MPFR_RNDN);
2304: mpfr_set(osum, sum, MPFR_RNDN);
2305: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
2306: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
2307: do {
2308: mpfr_set_si(kh, k, MPFR_RNDN);
2309: mpfr_mul(kh, kh, h, MPFR_RNDN);
2310: /* Weight */
2311: mpfr_set(wk, h, MPFR_RNDN);
2312: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
2313: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
2314: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
2315: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2316: mpfr_sqr(tmp, tmp, MPFR_RNDN);
2317: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
2318: mpfr_div(wk, wk, tmp, MPFR_RNDN);
2319: /* Abscissa */
2320: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
2321: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2322: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2323: mpfr_exp(tmp, msinh, MPFR_RNDN);
2324: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2325: /* Quadrature points */
2326: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
2327: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
2328: mpfr_add(lx, lx, beta, MPFR_RNDU);
2329: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
2330: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
2331: mpfr_add(rx, rx, beta, MPFR_RNDD);
2332: /* Evaluation */
2333: rtmp = mpfr_get_d(lx, MPFR_RNDU);
2334: func(&rtmp, ctx, &lval);
2335: rtmp = mpfr_get_d(rx, MPFR_RNDD);
2336: func(&rtmp, ctx, &rval);
2337: /* Update */
2338: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2339: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
2340: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2341: mpfr_abs(tmp, tmp, MPFR_RNDN);
2342: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2343: mpfr_set(curTerm, tmp, MPFR_RNDN);
2344: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2345: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
2346: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2347: mpfr_abs(tmp, tmp, MPFR_RNDN);
2348: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2349: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
2350: ++k;
2351: /* Only need to evaluate every other point on refined levels */
2352: if (l != 1) ++k;
2353: mpfr_log10(tmp, wk, MPFR_RNDN);
2354: mpfr_abs(tmp, tmp, MPFR_RNDN);
2355: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
2356: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
2357: mpfr_abs(tmp, tmp, MPFR_RNDN);
2358: mpfr_log10(tmp, tmp, MPFR_RNDN);
2359: d1 = mpfr_get_d(tmp, MPFR_RNDN);
2360: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2361: mpfr_abs(tmp, tmp, MPFR_RNDN);
2362: mpfr_log10(tmp, tmp, MPFR_RNDN);
2363: d2 = mpfr_get_d(tmp, MPFR_RNDN);
2364: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2365: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2366: mpfr_log10(tmp, curTerm, MPFR_RNDN);
2367: d4 = mpfr_get_d(tmp, MPFR_RNDN);
2368: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2369: } while (d < digits && l < 8);
2370: *sol = mpfr_get_d(sum, MPFR_RNDN);
2371: /* Cleanup */
2372: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2373: return 0;
2374: }
2375: #else
2377: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2378: {
2379: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2380: }
2381: #endif
2383: /*@
2384: PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures
2386: Not Collective
2388: Input Parameters:
2389: + q1 - The first quadrature
2390: - q2 - The second quadrature
2392: Output Parameter:
2393: . q - A PetscQuadrature object
2395: Level: intermediate
2397: .seealso: `PetscDTGaussTensorQuadrature()`
2398: @*/
2399: PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q)
2400: {
2401: const PetscReal *x1, *w1, *x2, *w2;
2402: PetscReal *x, *w;
2403: PetscInt dim1, Nc1, Np1, order1, qa, d1;
2404: PetscInt dim2, Nc2, Np2, order2, qb, d2;
2405: PetscInt dim, Nc, Np, order, qc, d;
2410: PetscQuadratureGetOrder(q1, &order1);
2411: PetscQuadratureGetOrder(q2, &order2);
2413: PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1);
2414: PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2);
2417: dim = dim1 + dim2;
2418: Nc = Nc1;
2419: Np = Np1 * Np2;
2420: order = order1;
2421: PetscQuadratureCreate(PETSC_COMM_SELF, q);
2422: PetscQuadratureSetOrder(*q, order);
2423: PetscMalloc1(Np * dim, &x);
2424: PetscMalloc1(Np, &w);
2425: for (qa = 0, qc = 0; qa < Np1; ++qa) {
2426: for (qb = 0; qb < Np2; ++qb, ++qc) {
2427: for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1];
2428: for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2];
2429: w[qc] = w1[qa] * w2[qb];
2430: }
2431: }
2432: PetscQuadratureSetData(*q, dim, Nc, Np, x, w);
2433: return 0;
2434: }
2436: /* Overwrites A. Can only handle full-rank problems with m>=n
2437: * A in column-major format
2438: * Ainv in row-major format
2439: * tau has length m
2440: * worksize must be >= max(1,n)
2441: */
2442: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2443: {
2444: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2445: PetscScalar *A, *Ainv, *R, *Q, Alpha;
2447: #if defined(PETSC_USE_COMPLEX)
2448: {
2449: PetscInt i, j;
2450: PetscMalloc2(m * n, &A, m * n, &Ainv);
2451: for (j = 0; j < n; j++) {
2452: for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j];
2453: }
2454: mstride = m;
2455: }
2456: #else
2457: A = A_in;
2458: Ainv = Ainv_out;
2459: #endif
2461: PetscBLASIntCast(m, &M);
2462: PetscBLASIntCast(n, &N);
2463: PetscBLASIntCast(mstride, &lda);
2464: PetscBLASIntCast(worksize, &ldwork);
2465: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2466: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2467: PetscFPTrapPop();
2469: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2471: /* Extract an explicit representation of Q */
2472: Q = Ainv;
2473: PetscArraycpy(Q, A, mstride * n);
2474: K = N; /* full rank */
2475: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2478: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2479: Alpha = 1.0;
2480: ldb = lda;
2481: PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb));
2482: /* Ainv is Q, overwritten with inverse */
2484: #if defined(PETSC_USE_COMPLEX)
2485: {
2486: PetscInt i;
2487: for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2488: PetscFree2(A, Ainv);
2489: }
2490: #endif
2491: return 0;
2492: }
2494: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2495: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B)
2496: {
2497: PetscReal *Bv;
2498: PetscInt i, j;
2500: PetscMalloc1((ninterval + 1) * ndegree, &Bv);
2501: /* Point evaluation of L_p on all the source vertices */
2502: PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL);
2503: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2504: for (i = 0; i < ninterval; i++) {
2505: for (j = 0; j < ndegree; j++) {
2506: if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2507: else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2508: }
2509: }
2510: PetscFree(Bv);
2511: return 0;
2512: }
2514: /*@
2515: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2517: Not Collective
2519: Input Parameters:
2520: + degree - degree of reconstruction polynomial
2521: . nsource - number of source intervals
2522: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2523: . ntarget - number of target intervals
2524: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
2526: Output Parameter:
2527: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2529: Level: advanced
2531: .seealso: `PetscDTLegendreEval()`
2532: @*/
2533: PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal *sourcex, PetscInt ntarget, const PetscReal *targetx, PetscReal *R)
2534: {
2535: PetscInt i, j, k, *bdegrees, worksize;
2536: PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget;
2537: PetscScalar *tau, *work;
2543: if (PetscDefined(USE_DEBUG)) {
2546: }
2547: xmin = PetscMin(sourcex[0], targetx[0]);
2548: xmax = PetscMax(sourcex[nsource], targetx[ntarget]);
2549: center = (xmin + xmax) / 2;
2550: hscale = (xmax - xmin) / 2;
2551: worksize = nsource;
2552: PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work);
2553: PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget);
2554: for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale;
2555: for (i = 0; i <= degree; i++) bdegrees[i] = i + 1;
2556: PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource);
2557: PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work);
2558: for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale;
2559: PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget);
2560: for (i = 0; i < ntarget; i++) {
2561: PetscReal rowsum = 0;
2562: for (j = 0; j < nsource; j++) {
2563: PetscReal sum = 0;
2564: for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j];
2565: R[i * nsource + j] = sum;
2566: rowsum += sum;
2567: }
2568: for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */
2569: }
2570: PetscFree4(bdegrees, sourcey, Bsource, work);
2571: PetscFree4(tau, Bsinv, targety, Btarget);
2572: return 0;
2573: }
2575: /*@C
2576: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2578: Not Collective
2580: Input Parameters:
2581: + n - the number of GLL nodes
2582: . nodes - the GLL nodes
2583: . weights - the GLL weights
2584: - f - the function values at the nodes
2586: Output Parameter:
2587: . in - the value of the integral
2589: Level: beginner
2591: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`
2593: @*/
2594: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal *nodes, PetscReal *weights, const PetscReal *f, PetscReal *in)
2595: {
2596: PetscInt i;
2598: *in = 0.;
2599: for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i];
2600: return 0;
2601: }
2603: /*@C
2604: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2606: Not Collective
2608: Input Parameters:
2609: + n - the number of GLL nodes
2610: . nodes - the GLL nodes
2611: - weights - the GLL weights
2613: Output Parameter:
2614: . A - the stiffness element
2616: Level: beginner
2618: Notes:
2619: Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()
2621: 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)
2623: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2625: @*/
2626: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2627: {
2628: PetscReal **A;
2629: const PetscReal *gllnodes = nodes;
2630: const PetscInt p = n - 1;
2631: PetscReal z0, z1, z2 = -1, x, Lpj, Lpr;
2632: PetscInt i, j, nn, r;
2634: PetscMalloc1(n, &A);
2635: PetscMalloc1(n * n, &A[0]);
2636: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
2638: for (j = 1; j < p; j++) {
2639: x = gllnodes[j];
2640: z0 = 1.;
2641: z1 = x;
2642: for (nn = 1; nn < p; nn++) {
2643: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2644: z0 = z1;
2645: z1 = z2;
2646: }
2647: Lpj = z2;
2648: for (r = 1; r < p; r++) {
2649: if (r == j) {
2650: A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj);
2651: } else {
2652: x = gllnodes[r];
2653: z0 = 1.;
2654: z1 = x;
2655: for (nn = 1; nn < p; nn++) {
2656: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2657: z0 = z1;
2658: z1 = z2;
2659: }
2660: Lpr = z2;
2661: A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r]));
2662: }
2663: }
2664: }
2665: for (j = 1; j < p + 1; j++) {
2666: x = gllnodes[j];
2667: z0 = 1.;
2668: z1 = x;
2669: for (nn = 1; nn < p; nn++) {
2670: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2671: z0 = z1;
2672: z1 = z2;
2673: }
2674: Lpj = z2;
2675: A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j]));
2676: A[0][j] = A[j][0];
2677: }
2678: for (j = 0; j < p; j++) {
2679: x = gllnodes[j];
2680: z0 = 1.;
2681: z1 = x;
2682: for (nn = 1; nn < p; nn++) {
2683: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2684: z0 = z1;
2685: z1 = z2;
2686: }
2687: Lpj = z2;
2689: A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j]));
2690: A[j][p] = A[p][j];
2691: }
2692: A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.;
2693: A[p][p] = A[0][0];
2694: *AA = A;
2695: return 0;
2696: }
2698: /*@C
2699: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
2701: Not Collective
2703: Input Parameters:
2704: + n - the number of GLL nodes
2705: . nodes - the GLL nodes
2706: . weights - the GLL weightss
2707: - A - the stiffness element
2709: Level: beginner
2711: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`
2713: @*/
2714: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2715: {
2716: PetscFree((*AA)[0]);
2717: PetscFree(*AA);
2718: *AA = NULL;
2719: return 0;
2720: }
2722: /*@C
2723: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
2725: Not Collective
2727: Input Parameter:
2728: + n - the number of GLL nodes
2729: . nodes - the GLL nodes
2730: . weights - the GLL weights
2732: Output Parameters:
2733: . AA - the stiffness element
2734: - AAT - the transpose of AA (pass in NULL if you do not need this array)
2736: Level: beginner
2738: Notes:
2739: Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()
2741: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2743: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2745: @*/
2746: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
2747: {
2748: PetscReal **A, **AT = NULL;
2749: const PetscReal *gllnodes = nodes;
2750: const PetscInt p = n - 1;
2751: PetscReal Li, Lj, d0;
2752: PetscInt i, j;
2754: PetscMalloc1(n, &A);
2755: PetscMalloc1(n * n, &A[0]);
2756: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
2758: if (AAT) {
2759: PetscMalloc1(n, &AT);
2760: PetscMalloc1(n * n, &AT[0]);
2761: for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n;
2762: }
2764: if (n == 1) A[0][0] = 0.;
2765: d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.;
2766: for (i = 0; i < n; i++) {
2767: for (j = 0; j < n; j++) {
2768: A[i][j] = 0.;
2769: PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);
2770: PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);
2771: if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j]));
2772: if ((j == i) && (i == 0)) A[i][j] = -d0;
2773: if (j == i && i == p) A[i][j] = d0;
2774: if (AT) AT[j][i] = A[i][j];
2775: }
2776: }
2777: if (AAT) *AAT = AT;
2778: *AA = A;
2779: return 0;
2780: }
2782: /*@C
2783: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()
2785: Not Collective
2787: Input Parameters:
2788: + n - the number of GLL nodes
2789: . nodes - the GLL nodes
2790: . weights - the GLL weights
2791: . AA - the stiffness element
2792: - AAT - the transpose of the element
2794: Level: beginner
2796: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
2798: @*/
2799: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
2800: {
2801: PetscFree((*AA)[0]);
2802: PetscFree(*AA);
2803: *AA = NULL;
2804: if (*AAT) {
2805: PetscFree((*AAT)[0]);
2806: PetscFree(*AAT);
2807: *AAT = NULL;
2808: }
2809: return 0;
2810: }
2812: /*@C
2813: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
2815: Not Collective
2817: Input Parameters:
2818: + n - the number of GLL nodes
2819: . nodes - the GLL nodes
2820: - weights - the GLL weightss
2822: Output Parameter:
2823: . AA - the stiffness element
2825: Level: beginner
2827: Notes:
2828: Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()
2830: This is the same as the Gradient operator multiplied by the diagonal mass matrix
2832: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2834: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()`
2836: @*/
2837: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2838: {
2839: PetscReal **D;
2840: const PetscReal *gllweights = weights;
2841: const PetscInt glln = n;
2842: PetscInt i, j;
2844: PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL);
2845: for (i = 0; i < glln; i++) {
2846: for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j];
2847: }
2848: *AA = D;
2849: return 0;
2850: }
2852: /*@C
2853: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
2855: Not Collective
2857: Input Parameters:
2858: + n - the number of GLL nodes
2859: . nodes - the GLL nodes
2860: . weights - the GLL weights
2861: - A - advection
2863: Level: beginner
2865: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
2867: @*/
2868: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2869: {
2870: PetscFree((*AA)[0]);
2871: PetscFree(*AA);
2872: *AA = NULL;
2873: return 0;
2874: }
2876: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2877: {
2878: PetscReal **A;
2879: const PetscReal *gllweights = weights;
2880: const PetscInt glln = n;
2881: PetscInt i, j;
2883: PetscMalloc1(glln, &A);
2884: PetscMalloc1(glln * glln, &A[0]);
2885: for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln;
2886: if (glln == 1) A[0][0] = 0.;
2887: for (i = 0; i < glln; i++) {
2888: for (j = 0; j < glln; j++) {
2889: A[i][j] = 0.;
2890: if (j == i) A[i][j] = gllweights[i];
2891: }
2892: }
2893: *AA = A;
2894: return 0;
2895: }
2897: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2898: {
2899: PetscFree((*AA)[0]);
2900: PetscFree(*AA);
2901: *AA = NULL;
2902: return 0;
2903: }
2905: /*@
2906: PetscDTIndexToBary - convert an index into a barycentric coordinate.
2908: Input Parameters:
2909: + 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)
2910: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2911: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
2913: Output Parameter:
2914: . coord - will be filled with the barycentric coordinate
2916: Level: beginner
2918: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2919: least significant and the last index is the most significant.
2921: .seealso: `PetscDTBaryToIndex()`
2922: @*/
2923: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
2924: {
2925: PetscInt c, d, s, total, subtotal, nexttotal;
2930: if (!len) {
2931: if (!sum && !index) return 0;
2932: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2933: }
2934: for (c = 1, total = 1; c <= len; c++) {
2935: /* total is the number of ways to have a tuple of length c with sum */
2936: if (index < total) break;
2937: total = (total * (sum + c)) / c;
2938: }
2940: for (d = c; d < len; d++) coord[d] = 0;
2941: for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
2942: /* subtotal is the number of ways to have a tuple of length c with sum s */
2943: /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
2944: if ((index + subtotal) >= total) {
2945: coord[--c] = sum - s;
2946: index -= (total - subtotal);
2947: sum = s;
2948: total = nexttotal;
2949: subtotal = 1;
2950: nexttotal = 1;
2951: s = 0;
2952: } else {
2953: subtotal = (subtotal * (c + s)) / (s + 1);
2954: nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
2955: s++;
2956: }
2957: }
2958: return 0;
2959: }
2961: /*@
2962: PetscDTBaryToIndex - convert a barycentric coordinate to an index
2964: Input Parameters:
2965: + 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)
2966: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2967: - coord - a barycentric coordinate with the given length and sum
2969: Output Parameter:
2970: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
2972: Level: beginner
2974: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2975: least significant and the last index is the most significant.
2977: .seealso: `PetscDTIndexToBary`
2978: @*/
2979: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
2980: {
2981: PetscInt c;
2982: PetscInt i;
2983: PetscInt total;
2987: if (!len) {
2988: if (!sum) {
2989: *index = 0;
2990: return 0;
2991: }
2992: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2993: }
2994: for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
2995: i = total - 1;
2996: c = len - 1;
2997: sum -= coord[c];
2998: while (sum > 0) {
2999: PetscInt subtotal;
3000: PetscInt s;
3002: for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
3003: i -= subtotal;
3004: sum -= coord[--c];
3005: }
3006: *index = i;
3007: return 0;
3008: }