Actual source code: dt.c
1: /* Discretization tools */
3: #include <petscdt.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/dtimpl.h>
7: #include <petscviewer.h>
8: #include <petscdmplex.h>
9: #include <petscdmshell.h>
11: #if defined(PETSC_HAVE_MPFR)
12: #include <mpfr.h>
13: #endif
15: const char *const PetscDTNodeTypes[] = {"gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
17: static PetscBool GolubWelschCite = PETSC_FALSE;
18: const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
19: " author = {Golub and Welsch},\n"
20: " title = {Calculation of Quadrature Rules},\n"
21: " journal = {Math. Comp.},\n"
22: " volume = {23},\n"
23: " number = {106},\n"
24: " pages = {221--230},\n"
25: " year = {1969}\n}\n";
27: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
28: quadrature rules:
30: - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
31: - in single precision, Newton's method starts producing incorrect roots around n = 15, but
32: the weights from Golub & Welsch become a problem before then: they produces errors
33: in computing the Jacobi-polynomial Gram matrix around n = 6.
35: So we default to Newton's method (required fewer dependencies) */
36: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
38: PetscClassId PETSCQUADRATURE_CLASSID = 0;
40: /*@
41: PetscQuadratureCreate - Create a PetscQuadrature object
43: Collective
45: Input Parameter:
46: . comm - The communicator for the PetscQuadrature object
48: Output Parameter:
49: . q - The PetscQuadrature object
51: Level: beginner
53: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
54: @*/
55: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
56: {
61: DMInitializePackage();
62: PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
63: (*q)->dim = -1;
64: (*q)->Nc = 1;
65: (*q)->order = -1;
66: (*q)->numPoints = 0;
67: (*q)->points = NULL;
68: (*q)->weights = NULL;
69: return(0);
70: }
72: /*@
73: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
75: Collective on q
77: Input Parameter:
78: . q - The PetscQuadrature object
80: Output Parameter:
81: . r - The new PetscQuadrature object
83: Level: beginner
85: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
86: @*/
87: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
88: {
89: PetscInt order, dim, Nc, Nq;
90: const PetscReal *points, *weights;
91: PetscReal *p, *w;
92: PetscErrorCode ierr;
96: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
97: PetscQuadratureGetOrder(q, &order);
98: PetscQuadratureSetOrder(*r, order);
99: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
100: PetscMalloc1(Nq*dim, &p);
101: PetscMalloc1(Nq*Nc, &w);
102: PetscArraycpy(p, points, Nq*dim);
103: PetscArraycpy(w, weights, Nc * Nq);
104: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
105: return(0);
106: }
108: /*@
109: PetscQuadratureDestroy - Destroys a PetscQuadrature object
111: Collective on q
113: Input Parameter:
114: . q - The PetscQuadrature object
116: Level: beginner
118: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
119: @*/
120: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
121: {
125: if (!*q) return(0);
127: if (--((PetscObject)(*q))->refct > 0) {
128: *q = NULL;
129: return(0);
130: }
131: PetscFree((*q)->points);
132: PetscFree((*q)->weights);
133: PetscHeaderDestroy(q);
134: return(0);
135: }
137: /*@
138: PetscQuadratureGetOrder - Return the order of the method
140: Not collective
142: Input Parameter:
143: . q - The PetscQuadrature object
145: Output Parameter:
146: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
148: Level: intermediate
150: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
151: @*/
152: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
153: {
157: *order = q->order;
158: return(0);
159: }
161: /*@
162: PetscQuadratureSetOrder - Return the order of the method
164: Not collective
166: Input Parameters:
167: + q - The PetscQuadrature object
168: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
170: Level: intermediate
172: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
173: @*/
174: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
175: {
178: q->order = order;
179: return(0);
180: }
182: /*@
183: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
185: Not collective
187: Input Parameter:
188: . q - The PetscQuadrature object
190: Output Parameter:
191: . Nc - The number of components
193: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
195: Level: intermediate
197: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
198: @*/
199: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
200: {
204: *Nc = q->Nc;
205: return(0);
206: }
208: /*@
209: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
211: Not collective
213: Input Parameters:
214: + q - The PetscQuadrature object
215: - Nc - The number of components
217: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
219: Level: intermediate
221: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
222: @*/
223: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
224: {
227: q->Nc = Nc;
228: return(0);
229: }
231: /*@C
232: PetscQuadratureGetData - Returns the data defining the quadrature
234: Not collective
236: Input Parameter:
237: . q - The PetscQuadrature object
239: Output Parameters:
240: + dim - The spatial dimension
241: . Nc - The number of components
242: . npoints - The number of quadrature points
243: . points - The coordinates of each quadrature point
244: - weights - The weight of each quadrature point
246: Level: intermediate
248: Fortran Notes:
249: From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
251: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
252: @*/
253: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
254: {
257: if (dim) {
259: *dim = q->dim;
260: }
261: if (Nc) {
263: *Nc = q->Nc;
264: }
265: if (npoints) {
267: *npoints = q->numPoints;
268: }
269: if (points) {
271: *points = q->points;
272: }
273: if (weights) {
275: *weights = q->weights;
276: }
277: return(0);
278: }
280: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
281: {
282: PetscScalar *Js, *Jinvs;
283: PetscInt i, j, k;
284: PetscBLASInt bm, bn, info;
288: if (!m || !n) return(0);
289: PetscBLASIntCast(m, &bm);
290: PetscBLASIntCast(n, &bn);
291: #if defined(PETSC_USE_COMPLEX)
292: PetscMalloc2(m*n, &Js, m*n, &Jinvs);
293: for (i = 0; i < m*n; i++) Js[i] = J[i];
294: #else
295: Js = (PetscReal *) J;
296: Jinvs = Jinv;
297: #endif
298: if (m == n) {
299: PetscBLASInt *pivots;
300: PetscScalar *W;
302: PetscMalloc2(m, &pivots, m, &W);
304: PetscArraycpy(Jinvs, Js, m * m);
305: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
306: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
307: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
308: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
309: PetscFree2(pivots, W);
310: } else if (m < n) {
311: PetscScalar *JJT;
312: PetscBLASInt *pivots;
313: PetscScalar *W;
315: PetscMalloc1(m*m, &JJT);
316: PetscMalloc2(m, &pivots, m, &W);
317: for (i = 0; i < m; i++) {
318: for (j = 0; j < m; j++) {
319: PetscScalar val = 0.;
321: for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
322: JJT[i * m + j] = val;
323: }
324: }
326: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
327: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
328: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
329: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
330: for (i = 0; i < n; i++) {
331: for (j = 0; j < m; j++) {
332: PetscScalar val = 0.;
334: for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
335: Jinvs[i * m + j] = val;
336: }
337: }
338: PetscFree2(pivots, W);
339: PetscFree(JJT);
340: } else {
341: PetscScalar *JTJ;
342: PetscBLASInt *pivots;
343: PetscScalar *W;
345: PetscMalloc1(n*n, &JTJ);
346: PetscMalloc2(n, &pivots, n, &W);
347: for (i = 0; i < n; i++) {
348: for (j = 0; j < n; j++) {
349: PetscScalar val = 0.;
351: for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
352: JTJ[i * n + j] = val;
353: }
354: }
356: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
357: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
358: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
359: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
360: for (i = 0; i < n; i++) {
361: for (j = 0; j < m; j++) {
362: PetscScalar val = 0.;
364: for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
365: Jinvs[i * m + j] = val;
366: }
367: }
368: PetscFree2(pivots, W);
369: PetscFree(JTJ);
370: }
371: #if defined(PETSC_USE_COMPLEX)
372: for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
373: PetscFree2(Js, Jinvs);
374: #endif
375: return(0);
376: }
378: /*@
379: PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
381: Collecive on PetscQuadrature
383: Input Parameters:
384: + q - the quadrature functional
385: . imageDim - the dimension of the image of the transformation
386: . origin - a point in the original space
387: . originImage - the image of the origin under the transformation
388: . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
389: - 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]
391: Output Parameters:
392: . 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.
394: 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.
396: Level: intermediate
398: .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
399: @*/
400: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
401: {
402: PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
403: const PetscReal *points;
404: const PetscReal *weights;
405: PetscReal *imagePoints, *imageWeights;
406: PetscReal *Jinv;
407: PetscReal *Jinvstar;
408: PetscErrorCode ierr;
412: if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim);
413: PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);
414: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);
415: if (Nc % formSize) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of formSize %D\n", Nc, formSize);
416: Ncopies = Nc / formSize;
417: PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);
418: imageNc = Ncopies * imageFormSize;
419: PetscMalloc1(Npoints * imageDim, &imagePoints);
420: PetscMalloc1(Npoints * imageNc, &imageWeights);
421: PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);
422: PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);
423: PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);
424: for (pt = 0; pt < Npoints; pt++) {
425: const PetscReal *point = &points[pt * dim];
426: PetscReal *imagePoint = &imagePoints[pt * imageDim];
428: for (i = 0; i < imageDim; i++) {
429: PetscReal val = originImage[i];
431: for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
432: imagePoint[i] = val;
433: }
434: for (c = 0; c < Ncopies; c++) {
435: const PetscReal *form = &weights[pt * Nc + c * formSize];
436: PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
438: for (i = 0; i < imageFormSize; i++) {
439: PetscReal val = 0.;
441: for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
442: imageForm[i] = val;
443: }
444: }
445: }
446: PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);
447: PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);
448: PetscFree2(Jinv, Jinvstar);
449: return(0);
450: }
452: /*@C
453: PetscQuadratureSetData - Sets the data defining the quadrature
455: Not collective
457: Input Parameters:
458: + q - The PetscQuadrature object
459: . dim - The spatial dimension
460: . Nc - The number of components
461: . npoints - The number of quadrature points
462: . points - The coordinates of each quadrature point
463: - weights - The weight of each quadrature point
465: Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
467: Level: intermediate
469: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
470: @*/
471: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
472: {
475: if (dim >= 0) q->dim = dim;
476: if (Nc >= 0) q->Nc = Nc;
477: if (npoints >= 0) q->numPoints = npoints;
478: if (points) {
480: q->points = points;
481: }
482: if (weights) {
484: q->weights = weights;
485: }
486: return(0);
487: }
489: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
490: {
491: PetscInt q, d, c;
492: PetscViewerFormat format;
493: PetscErrorCode ierr;
496: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);}
497: else {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);}
498: PetscViewerGetFormat(v, &format);
499: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
500: for (q = 0; q < quad->numPoints; ++q) {
501: PetscViewerASCIIPrintf(v, "p%D (", q);
502: PetscViewerASCIIUseTabs(v, PETSC_FALSE);
503: for (d = 0; d < quad->dim; ++d) {
504: if (d) {PetscViewerASCIIPrintf(v, ", ");}
505: PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
506: }
507: PetscViewerASCIIPrintf(v, ") ");
508: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "w%D (", q);}
509: for (c = 0; c < quad->Nc; ++c) {
510: if (c) {PetscViewerASCIIPrintf(v, ", ");}
511: PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
512: }
513: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
514: PetscViewerASCIIPrintf(v, "\n");
515: PetscViewerASCIIUseTabs(v, PETSC_TRUE);
516: }
517: return(0);
518: }
520: /*@C
521: PetscQuadratureView - Views a PetscQuadrature object
523: Collective on quad
525: Input Parameters:
526: + quad - The PetscQuadrature object
527: - viewer - The PetscViewer object
529: Level: beginner
531: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
532: @*/
533: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
534: {
535: PetscBool iascii;
541: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
542: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
543: PetscViewerASCIIPushTab(viewer);
544: if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
545: PetscViewerASCIIPopTab(viewer);
546: return(0);
547: }
549: /*@C
550: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
552: Not collective
554: Input Parameters:
555: + q - The original PetscQuadrature
556: . numSubelements - The number of subelements the original element is divided into
557: . v0 - An array of the initial points for each subelement
558: - jac - An array of the Jacobian mappings from the reference to each subelement
560: Output Parameters:
561: . dim - The dimension
563: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
565: Not available from Fortran
567: Level: intermediate
569: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
570: @*/
571: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
572: {
573: const PetscReal *points, *weights;
574: PetscReal *pointsRef, *weightsRef;
575: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
576: PetscErrorCode ierr;
583: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
584: PetscQuadratureGetOrder(q, &order);
585: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
586: npointsRef = npoints*numSubelements;
587: PetscMalloc1(npointsRef*dim,&pointsRef);
588: PetscMalloc1(npointsRef*Nc, &weightsRef);
589: for (c = 0; c < numSubelements; ++c) {
590: for (p = 0; p < npoints; ++p) {
591: for (d = 0; d < dim; ++d) {
592: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
593: for (e = 0; e < dim; ++e) {
594: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
595: }
596: }
597: /* Could also use detJ here */
598: for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
599: }
600: }
601: PetscQuadratureSetOrder(*qref, order);
602: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
603: return(0);
604: }
606: /* Compute the coefficients for the Jacobi polynomial recurrence,
607: *
608: * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
609: */
610: #define PetscDTJacobiRecurrence_Internal(n,a,b,cnm1,cnm1x,cnm2) \
611: do { \
612: PetscReal _a = (a); \
613: PetscReal _b = (b); \
614: PetscReal _n = (n); \
615: if (n == 1) { \
616: (cnm1) = (_a-_b) * 0.5; \
617: (cnm1x) = (_a+_b+2.)*0.5; \
618: (cnm2) = 0.; \
619: } else { \
620: PetscReal _2n = _n+_n; \
621: PetscReal _d = (_2n*(_n+_a+_b)*(_2n+_a+_b-2)); \
622: PetscReal _n1 = (_2n+_a+_b-1.)*(_a*_a-_b*_b); \
623: PetscReal _n1x = (_2n+_a+_b-1.)*(_2n+_a+_b)*(_2n+_a+_b-2); \
624: PetscReal _n2 = 2.*((_n+_a-1.)*(_n+_b-1.)*(_2n+_a+_b)); \
625: (cnm1) = _n1 / _d; \
626: (cnm1x) = _n1x / _d; \
627: (cnm2) = _n2 / _d; \
628: } \
629: } while (0)
631: /*@
632: PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
634: $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
636: Input Parameters:
637: - alpha - the left exponent > -1
638: . beta - the right exponent > -1
639: + n - the polynomial degree
641: Output Parameter:
642: . norm - the weighted L2 norm
644: Level: beginner
646: .seealso: PetscDTJacobiEval()
647: @*/
648: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
649: {
650: PetscReal twoab1;
651: PetscReal gr;
654: if (alpha <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid\n", (double) alpha);
655: if (beta <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid\n", (double) beta);
656: if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %D < 0 invalid\n", n);
657: twoab1 = PetscPowReal(2., alpha + beta + 1.);
658: #if defined(PETSC_HAVE_LGAMMA)
659: if (!n) {
660: gr = PetscExpReal(PetscLGamma(alpha+1.) + PetscLGamma(beta+1.) - PetscLGamma(alpha+beta+2.));
661: } else {
662: gr = PetscExpReal(PetscLGamma(n+alpha+1.) + PetscLGamma(n+beta+1.) - (PetscLGamma(n+1.) + PetscLGamma(n+alpha+beta+1.))) / (n+n+alpha+beta+1.);
663: }
664: #else
665: {
666: PetscInt alphai = (PetscInt) alpha;
667: PetscInt betai = (PetscInt) beta;
668: PetscInt i;
670: gr = n ? (1. / (n+n+alpha+beta+1.)) : 1.;
671: if ((PetscReal) alphai == alpha) {
672: if (!n) {
673: for (i = 0; i < alphai; i++) gr *= (i+1.) / (beta+i+1.);
674: gr /= (alpha+beta+1.);
675: } else {
676: for (i = 0; i < alphai; i++) gr *= (n+i+1.) / (n+beta+i+1.);
677: }
678: } else if ((PetscReal) betai == beta) {
679: if (!n) {
680: for (i = 0; i < betai; i++) gr *= (i+1.) / (alpha+i+2.);
681: gr /= (alpha+beta+1.);
682: } else {
683: for (i = 0; i < betai; i++) gr *= (n+i+1.) / (n+alpha+i+1.);
684: }
685: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
686: }
687: #endif
688: *norm = PetscSqrtReal(twoab1 * gr);
689: return(0);
690: }
692: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
693: {
694: PetscReal ak, bk;
695: PetscReal abk1;
696: PetscInt i,l,maxdegree;
699: maxdegree = degrees[ndegree-1] - k;
700: ak = a + k;
701: bk = b + k;
702: abk1 = a + b + k + 1.;
703: if (maxdegree < 0) {
704: for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.;
705: return(0);
706: }
707: for (i=0; i<npoints; i++) {
708: PetscReal pm1,pm2,x;
709: PetscReal cnm1, cnm1x, cnm2;
710: PetscInt j,m;
712: x = points[i];
713: pm2 = 1.;
714: PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2);
715: pm1 = (cnm1 + cnm1x*x);
716: l = 0;
717: while (l < ndegree && degrees[l] - k < 0) {
718: p[l++] = 0.;
719: }
720: while (l < ndegree && degrees[l] - k == 0) {
721: p[l] = pm2;
722: for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
723: l++;
724: }
725: while (l < ndegree && degrees[l] - k == 1) {
726: p[l] = pm1;
727: for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
728: l++;
729: }
730: for (j=2; j<=maxdegree; j++) {
731: PetscReal pp;
733: PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2);
734: pp = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2;
735: pm2 = pm1;
736: pm1 = pp;
737: while (l < ndegree && degrees[l] - k == j) {
738: p[l] = pp;
739: for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
740: l++;
741: }
742: }
743: p += ndegree;
744: }
745: return(0);
746: }
748: /*@
749: 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$.
751: Input Parameters:
752: + alpha - the left exponent of the weight
753: . beta - the right exponetn of the weight
754: . npoints - the number of points to evaluate the polynomials at
755: . points - [npoints] array of point coordinates
756: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
757: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
759: Output Argments:
760: - p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
761: (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
762: (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
763: varying) dimension is the index of the evaluation point.
765: Level: advanced
767: .seealso: PetscDTJacobiEval(), PetscDTPKDEvalJet()
768: @*/
769: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
770: {
771: PetscInt i, j, l;
772: PetscInt *degrees;
773: PetscReal *psingle;
774: PetscErrorCode ierr;
777: if (degree == 0) {
778: PetscInt zero = 0;
780: for (i = 0; i <= k; i++) {
781: PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i*npoints]);
782: }
783: return(0);
784: }
785: PetscMalloc1(degree + 1, °rees);
786: PetscMalloc1((degree + 1) * npoints, &psingle);
787: for (i = 0; i <= degree; i++) degrees[i] = i;
788: for (i = 0; i <= k; i++) {
789: PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle);
790: for (j = 0; j <= degree; j++) {
791: for (l = 0; l < npoints; l++) {
792: p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
793: }
794: }
795: }
796: PetscFree(psingle);
797: PetscFree(degrees);
798: return(0);
799: }
801: /*@
802: PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$
803: at points
805: Not Collective
807: Input Parameters:
808: + npoints - number of spatial points to evaluate at
809: . alpha - the left exponent > -1
810: . beta - the right exponent > -1
811: . points - array of locations to evaluate at
812: . ndegree - number of basis degrees to evaluate
813: - degrees - sorted array of degrees to evaluate
815: Output Parameters:
816: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
817: . D - row-oriented derivative evaluation matrix (or NULL)
818: - D2 - row-oriented second derivative evaluation matrix (or NULL)
820: Level: intermediate
822: .seealso: PetscDTGaussQuadrature()
823: @*/
824: PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
825: {
829: if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
830: if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
831: if (!npoints || !ndegree) return(0);
832: if (B) {PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);}
833: if (D) {PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);}
834: if (D2) {PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);}
835: return(0);
836: }
838: /*@
839: PetscDTLegendreEval - evaluate Legendre polynomials at points
841: Not Collective
843: Input Parameters:
844: + npoints - number of spatial points to evaluate at
845: . points - array of locations to evaluate at
846: . ndegree - number of basis degrees to evaluate
847: - degrees - sorted array of degrees to evaluate
849: Output Parameters:
850: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
851: . D - row-oriented derivative evaluation matrix (or NULL)
852: - D2 - row-oriented second derivative evaluation matrix (or NULL)
854: Level: intermediate
856: .seealso: PetscDTGaussQuadrature()
857: @*/
858: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
859: {
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;
891: if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
892: if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
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;
939: if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
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 Prioriol-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 Argments:
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^3,y^1,z^2, which as 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;
1009: PetscErrorCode ierr;
1012: PetscCitationsRegister(PKDCitation, &PKDCite);
1013: PetscDTBinomialInt(dim + k, k, &Nk);
1014: PetscDTBinomialInt(degree + dim, degree, &Ndeg);
1015: PetscMalloc2(dim, °tup, dim, &ktup);
1016: PetscMalloc1(Ndeg, &scales);
1017: initscale = 1.;
1018: if (dim > 1) {
1019: PetscDTBinomial(dim,2,&scaleexp);
1020: initscale = PetscPowReal(2.,scaleexp*0.5);
1021: }
1022: for (degidx = 0; degidx < Ndeg; degidx++) {
1023: PetscInt e, i;
1024: PetscInt m1idx = -1, m2idx = -1;
1025: PetscInt n;
1026: PetscInt degsum;
1027: PetscReal alpha;
1028: PetscReal cnm1, cnm1x, cnm2;
1029: PetscReal norm;
1031: PetscDTIndexToGradedOrder(dim, degidx, degtup);
1032: for (d = dim - 1; d >= 0; d--) if (degtup[d]) break;
1033: if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1034: scales[degidx] = initscale;
1035: for (e = 0; e < dim; e++) {
1036: PetscDTJacobiNorm(e,0.,0,&norm);
1037: scales[degidx] /= norm;
1038: }
1039: for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1040: for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1041: continue;
1042: }
1043: n = degtup[d];
1044: degtup[d]--;
1045: PetscDTGradedOrderToIndex(dim, degtup, &m1idx);
1046: if (degtup[d] > 0) {
1047: degtup[d]--;
1048: PetscDTGradedOrderToIndex(dim, degtup, &m2idx);
1049: degtup[d]++;
1050: }
1051: degtup[d]++;
1052: for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1053: alpha = 2 * degsum + d;
1054: PetscDTJacobiRecurrence_Internal(n,alpha,0.,cnm1,cnm1x,cnm2);
1056: scales[degidx] = initscale;
1057: for (e = 0, degsum = 0; e < dim; e++) {
1058: PetscInt f;
1059: PetscReal ealpha;
1060: PetscReal enorm;
1062: ealpha = 2 * degsum + e;
1063: for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1064: PetscDTJacobiNorm(ealpha,0.,degtup[e],&enorm);
1065: scales[degidx] /= enorm;
1066: degsum += degtup[e];
1067: }
1069: for (pt = 0; pt < npoints; pt++) {
1070: /* compute the multipliers */
1071: PetscReal thetanm1, thetanm1x, thetanm2;
1073: thetanm1x = dim - (d+1) + 2.*points[pt * dim + d];
1074: for (e = d+1; e < dim; e++) thetanm1x += points[pt * dim + e];
1075: thetanm1x *= 0.5;
1076: thetanm1 = (2. - (dim-(d+1)));
1077: for (e = d+1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1078: thetanm1 *= 0.5;
1079: thetanm2 = thetanm1 * thetanm1;
1081: for (kidx = 0; kidx < Nk; kidx++) {
1082: PetscInt f;
1084: PetscDTIndexToGradedOrder(dim, kidx, ktup);
1085: /* first sum in the same derivative terms */
1086: p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1087: if (m2idx >= 0) {
1088: p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];
1089: }
1091: for (f = d; f < dim; f++) {
1092: PetscInt km1idx, mplty = ktup[f];
1094: if (!mplty) continue;
1095: ktup[f]--;
1096: PetscDTGradedOrderToIndex(dim, ktup, &km1idx);
1098: /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1099: /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1100: /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1101: if (f > d) {
1102: PetscInt f2;
1104: p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1105: if (m2idx >= 0) {
1106: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1107: /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1108: for (f2 = f; f2 < dim; f2++) {
1109: PetscInt km2idx, mplty2 = ktup[f2];
1110: PetscInt factor;
1112: if (!mplty2) continue;
1113: ktup[f2]--;
1114: PetscDTGradedOrderToIndex(dim, ktup, &km2idx);
1116: factor = mplty * mplty2;
1117: if (f == f2) factor /= 2;
1118: p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1119: ktup[f2]++;
1120: }
1121: }
1122: } else {
1123: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1124: }
1125: ktup[f]++;
1126: }
1127: }
1128: }
1129: }
1130: for (degidx = 0; degidx < Ndeg; degidx++) {
1131: PetscReal scale = scales[degidx];
1132: PetscInt i;
1134: for (i = 0; i < Nk * npoints; i++) p[degidx*Nk*npoints + i] *= scale;
1135: }
1136: PetscFree(scales);
1137: PetscFree2(degtup, ktup);
1138: return(0);
1139: }
1141: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1142: * with lds n; diag and subdiag are overwritten */
1143: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[],
1144: PetscReal eigs[], PetscScalar V[])
1145: {
1146: char jobz = 'V'; /* eigenvalues and eigenvectors */
1147: char range = 'A'; /* all eigenvalues will be found */
1148: PetscReal VL = 0.; /* ignored because range is 'A' */
1149: PetscReal VU = 0.; /* ignored because range is 'A' */
1150: PetscBLASInt IL = 0; /* ignored because range is 'A' */
1151: PetscBLASInt IU = 0; /* ignored because range is 'A' */
1152: PetscReal abstol = 0.; /* unused */
1153: PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1154: PetscBLASInt *isuppz;
1155: PetscBLASInt lwork, liwork;
1156: PetscReal workquery;
1157: PetscBLASInt iworkquery;
1158: PetscBLASInt *iwork;
1159: PetscBLASInt info;
1160: PetscReal *work = NULL;
1164: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1165: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1166: #endif
1167: PetscBLASIntCast(n, &bn);
1168: PetscBLASIntCast(n, &ldz);
1169: #if !defined(PETSC_MISSING_LAPACK_STEGR)
1170: PetscMalloc1(2 * n, &isuppz);
1171: lwork = -1;
1172: liwork = -1;
1173: PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info));
1174: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
1175: lwork = (PetscBLASInt) workquery;
1176: liwork = (PetscBLASInt) iworkquery;
1177: PetscMalloc2(lwork, &work, liwork, &iwork);
1178: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1179: PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info));
1180: PetscFPTrapPop();
1181: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
1182: PetscFree2(work, iwork);
1183: PetscFree(isuppz);
1184: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1185: jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1186: tridiagonal matrix. Z is initialized to the identity
1187: matrix. */
1188: PetscMalloc1(PetscMax(1,2*n-2),&work);
1189: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info));
1190: PetscFPTrapPop();
1191: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
1192: PetscFree(work);
1193: PetscArraycpy(eigs,diag,n);
1194: #endif
1195: return(0);
1196: }
1198: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1199: * quadrature rules on the interval [-1, 1] */
1200: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1201: {
1202: PetscReal twoab1;
1203: PetscInt m = n - 2;
1204: PetscReal a = alpha + 1.;
1205: PetscReal b = beta + 1.;
1206: PetscReal gra, grb;
1209: twoab1 = PetscPowReal(2., a + b - 1.);
1210: #if defined(PETSC_HAVE_LGAMMA)
1211: grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) -
1212: (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.)));
1213: gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) -
1214: (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.)));
1215: #else
1216: {
1217: PetscInt alphai = (PetscInt) alpha;
1218: PetscInt betai = (PetscInt) beta;
1221: if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) {
1222: PetscReal binom1, binom2;
1224: PetscDTBinomial(m+b, b, &binom1);
1225: PetscDTBinomial(m+a+b, b, &binom2);
1226: grb = 1./ (binom1 * binom2);
1227: PetscDTBinomial(m+a, a, &binom1);
1228: PetscDTBinomial(m+a+b, a, &binom2);
1229: gra = 1./ (binom1 * binom2);
1230: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1231: }
1232: #endif
1233: *leftw = twoab1 * grb / b;
1234: *rightw = twoab1 * gra / a;
1235: return(0);
1236: }
1238: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1239: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1240: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1241: {
1242: PetscReal pn1, pn2;
1243: PetscReal cnm1, cnm1x, cnm2;
1244: PetscInt k;
1247: if (!n) {*P = 1.0; return(0);}
1248: PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2);
1249: pn2 = 1.;
1250: pn1 = cnm1 + cnm1x*x;
1251: if (n == 1) {*P = pn1; return(0);}
1252: *P = 0.0;
1253: for (k = 2; k < n+1; ++k) {
1254: PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2);
1256: *P = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2;
1257: pn2 = pn1;
1258: pn1 = *P;
1259: }
1260: return(0);
1261: }
1263: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1264: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1265: {
1266: PetscReal nP;
1267: PetscInt i;
1271: *P = 0.0;
1272: if (k > n) return(0);
1273: PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);
1274: for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1275: *P = nP;
1276: return(0);
1277: }
1279: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1280: {
1281: PetscInt maxIter = 100;
1282: PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1283: PetscReal a1, a6, gf;
1284: PetscInt k;
1289: a1 = PetscPowReal(2.0, a+b+1);
1290: #if defined(PETSC_HAVE_LGAMMA)
1291: {
1292: PetscReal a2, a3, a4, a5;
1293: a2 = PetscLGamma(a + npoints + 1);
1294: a3 = PetscLGamma(b + npoints + 1);
1295: a4 = PetscLGamma(a + b + npoints + 1);
1296: a5 = PetscLGamma(npoints + 1);
1297: gf = PetscExpReal(a2 + a3 - (a4 + a5));
1298: }
1299: #else
1300: {
1301: PetscInt ia, ib;
1303: ia = (PetscInt) a;
1304: ib = (PetscInt) b;
1305: gf = 1.;
1306: if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1307: for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1308: } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1309: for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1310: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1311: }
1312: #endif
1314: a6 = a1 * gf;
1315: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1316: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1317: for (k = 0; k < npoints; ++k) {
1318: PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP;
1319: PetscInt j;
1321: if (k > 0) r = 0.5 * (r + x[k-1]);
1322: for (j = 0; j < maxIter; ++j) {
1323: PetscReal s = 0.0, delta, f, fp;
1324: PetscInt i;
1326: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1327: PetscDTComputeJacobi(a, b, npoints, r, &f);
1328: PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);
1329: delta = f / (fp - f * s);
1330: r = r - delta;
1331: if (PetscAbsReal(delta) < eps) break;
1332: }
1333: x[k] = r;
1334: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);
1335: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1336: }
1337: return(0);
1338: }
1340: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1341: * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1342: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1343: {
1344: PetscInt i;
1347: for (i = 0; i < nPoints; i++) {
1348: PetscReal A, B, C;
1350: PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C);
1351: d[i] = -A / B;
1352: if (i) s[i-1] *= C / B;
1353: if (i < nPoints - 1) s[i] = 1. / B;
1354: }
1355: return(0);
1356: }
1358: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1359: {
1360: PetscReal mu0;
1361: PetscReal ga, gb, gab;
1362: PetscInt i;
1366: PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);
1368: #if defined(PETSC_HAVE_TGAMMA)
1369: ga = PetscTGamma(a + 1);
1370: gb = PetscTGamma(b + 1);
1371: gab = PetscTGamma(a + b + 2);
1372: #else
1373: {
1374: PetscInt ia, ib;
1376: ia = (PetscInt) a;
1377: ib = (PetscInt) b;
1378: if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1379: PetscDTFactorial(ia, &ga);
1380: PetscDTFactorial(ib, &gb);
1381: PetscDTFactorial(ia + ib + 1, &gb);
1382: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
1383: }
1384: #endif
1385: mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab;
1387: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1388: {
1389: PetscReal *diag, *subdiag;
1390: PetscScalar *V;
1392: PetscMalloc2(npoints, &diag, npoints, &subdiag);
1393: PetscMalloc1(npoints*npoints, &V);
1394: PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);
1395: for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1396: PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);
1397: for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1398: PetscFree(V);
1399: PetscFree2(diag, subdiag);
1400: }
1401: #else
1402: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1403: #endif
1404: { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1405: eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1406: the eigenvalues are sorted */
1407: PetscBool sorted;
1409: PetscSortedReal(npoints, x, &sorted);
1410: if (!sorted) {
1411: PetscInt *order, i;
1412: PetscReal *tmp;
1414: PetscMalloc2(npoints, &order, npoints, &tmp);
1415: for (i = 0; i < npoints; i++) order[i] = i;
1416: PetscSortRealWithPermutation(npoints, x, order);
1417: PetscArraycpy(tmp, x, npoints);
1418: for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1419: PetscArraycpy(tmp, w, npoints);
1420: for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1421: PetscFree2(order, tmp);
1422: }
1423: }
1424: return(0);
1425: }
1427: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1428: {
1432: if (npoints < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1433: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1434: if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1435: if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
1437: if (newton) {
1438: PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);
1439: } else {
1440: PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);
1441: }
1442: if (alpha == beta) { /* symmetrize */
1443: PetscInt i;
1444: for (i = 0; i < (npoints + 1) / 2; i++) {
1445: PetscInt j = npoints - 1 - i;
1446: PetscReal xi = x[i];
1447: PetscReal xj = x[j];
1448: PetscReal wi = w[i];
1449: PetscReal wj = w[j];
1451: x[i] = (xi - xj) / 2.;
1452: x[j] = (xj - xi) / 2.;
1453: w[i] = w[j] = (wi + wj) / 2.;
1454: }
1455: }
1456: return(0);
1457: }
1459: /*@
1460: PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1461: $(x-a)^\alpha (x-b)^\beta$.
1463: Not collective
1465: Input Parameters:
1466: + npoints - the number of points in the quadrature rule
1467: . a - the left endpoint of the interval
1468: . b - the right endpoint of the interval
1469: . alpha - the left exponent
1470: - beta - the right exponent
1472: Output Parameters:
1473: + x - array of length npoints, the locations of the quadrature points
1474: - w - array of length npoints, the weights of the quadrature points
1476: Level: intermediate
1478: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1479: @*/
1480: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1481: {
1482: PetscInt i;
1486: PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1487: if (a != -1. || b != 1.) { /* shift */
1488: for (i = 0; i < npoints; i++) {
1489: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1490: w[i] *= (b - a) / 2.;
1491: }
1492: }
1493: return(0);
1494: }
1496: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1497: {
1498: PetscInt i;
1502: if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1503: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1504: if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1505: if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
1507: x[0] = -1.;
1508: x[npoints-1] = 1.;
1509: if (npoints > 2) {
1510: PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);
1511: }
1512: for (i = 1; i < npoints - 1; i++) {
1513: w[i] /= (1. - x[i]*x[i]);
1514: }
1515: PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);
1516: return(0);
1517: }
1519: /*@
1520: PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1521: $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.
1523: Not collective
1525: Input Parameters:
1526: + npoints - the number of points in the quadrature rule
1527: . a - the left endpoint of the interval
1528: . b - the right endpoint of the interval
1529: . alpha - the left exponent
1530: - beta - the right exponent
1532: Output Parameters:
1533: + x - array of length npoints, the locations of the quadrature points
1534: - w - array of length npoints, the weights of the quadrature points
1536: Level: intermediate
1538: Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1539: @*/
1540: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1541: {
1542: PetscInt i;
1546: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1547: if (a != -1. || b != 1.) { /* shift */
1548: for (i = 0; i < npoints; i++) {
1549: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1550: w[i] *= (b - a) / 2.;
1551: }
1552: }
1553: return(0);
1554: }
1556: /*@
1557: PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1559: Not Collective
1561: Input Parameters:
1562: + npoints - number of points
1563: . a - left end of interval (often-1)
1564: - b - right end of interval (often +1)
1566: Output Parameters:
1567: + x - quadrature points
1568: - w - quadrature weights
1570: Level: intermediate
1572: References:
1573: . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
1575: .seealso: PetscDTLegendreEval()
1576: @*/
1577: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
1578: {
1579: PetscInt i;
1583: PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);
1584: if (a != -1. || b != 1.) { /* shift */
1585: for (i = 0; i < npoints; i++) {
1586: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1587: w[i] *= (b - a) / 2.;
1588: }
1589: }
1590: return(0);
1591: }
1593: /*@C
1594: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1595: nodes of a given size on the domain [-1,1]
1597: Not Collective
1599: Input Parameters:
1600: + n - number of grid nodes
1601: - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON
1603: Output Parameters:
1604: + x - quadrature points
1605: - w - quadrature weights
1607: Notes:
1608: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1609: close enough to the desired solution
1611: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1613: 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
1615: Level: intermediate
1617: .seealso: PetscDTGaussQuadrature()
1619: @*/
1620: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
1621: {
1622: PetscBool newton;
1626: if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
1627: newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1628: PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);
1629: return(0);
1630: }
1632: /*@
1633: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1635: Not Collective
1637: Input Parameters:
1638: + dim - The spatial dimension
1639: . Nc - The number of components
1640: . npoints - number of points in one dimension
1641: . a - left end of interval (often-1)
1642: - b - right end of interval (often +1)
1644: Output Parameter:
1645: . q - A PetscQuadrature object
1647: Level: intermediate
1649: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
1650: @*/
1651: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1652: {
1653: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
1654: PetscReal *x, *w, *xw, *ww;
1658: PetscMalloc1(totpoints*dim,&x);
1659: PetscMalloc1(totpoints*Nc,&w);
1660: /* Set up the Golub-Welsch system */
1661: switch (dim) {
1662: case 0:
1663: PetscFree(x);
1664: PetscFree(w);
1665: PetscMalloc1(1, &x);
1666: PetscMalloc1(Nc, &w);
1667: x[0] = 0.0;
1668: for (c = 0; c < Nc; ++c) w[c] = 1.0;
1669: break;
1670: case 1:
1671: PetscMalloc1(npoints,&ww);
1672: PetscDTGaussQuadrature(npoints, a, b, x, ww);
1673: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
1674: PetscFree(ww);
1675: break;
1676: case 2:
1677: PetscMalloc2(npoints,&xw,npoints,&ww);
1678: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1679: for (i = 0; i < npoints; ++i) {
1680: for (j = 0; j < npoints; ++j) {
1681: x[(i*npoints+j)*dim+0] = xw[i];
1682: x[(i*npoints+j)*dim+1] = xw[j];
1683: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
1684: }
1685: }
1686: PetscFree2(xw,ww);
1687: break;
1688: case 3:
1689: PetscMalloc2(npoints,&xw,npoints,&ww);
1690: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1691: for (i = 0; i < npoints; ++i) {
1692: for (j = 0; j < npoints; ++j) {
1693: for (k = 0; k < npoints; ++k) {
1694: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
1695: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
1696: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
1697: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
1698: }
1699: }
1700: }
1701: PetscFree2(xw,ww);
1702: break;
1703: default:
1704: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1705: }
1706: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1707: PetscQuadratureSetOrder(*q, 2*npoints-1);
1708: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1709: PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
1710: return(0);
1711: }
1713: /*@
1714: PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex
1716: Not Collective
1718: Input Parameters:
1719: + dim - The simplex dimension
1720: . Nc - The number of components
1721: . npoints - The number of points in one dimension
1722: . a - left end of interval (often-1)
1723: - b - right end of interval (often +1)
1725: Output Parameter:
1726: . q - A PetscQuadrature object
1728: Level: intermediate
1730: References:
1731: . 1. - Karniadakis and Sherwin. FIAT
1733: Note: For dim == 1, this is Gauss-Legendre quadrature
1735: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
1736: @*/
1737: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1738: {
1739: PetscInt totprev, totrem;
1740: PetscInt totpoints;
1741: PetscReal *p1, *w1;
1742: PetscReal *x, *w;
1743: PetscInt i, j, k, l, m, pt, c;
1747: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
1748: totpoints = 1;
1749: for (i = 0, totpoints = 1; i < dim; i++) totpoints *= npoints;
1750: PetscMalloc1(totpoints*dim, &x);
1751: PetscMalloc1(totpoints*Nc, &w);
1752: PetscMalloc2(npoints, &p1, npoints, &w1);
1753: for (i = 0; i < totpoints*Nc; i++) w[i] = 1.;
1754: for (i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; i++) {
1755: PetscReal mul;
1757: mul = PetscPowReal(2.,-i);
1758: PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1);
1759: for (pt = 0, l = 0; l < totprev; l++) {
1760: for (j = 0; j < npoints; j++) {
1761: for (m = 0; m < totrem; m++, pt++) {
1762: for (k = 0; k < i; k++) x[pt*dim+k] = (x[pt*dim+k]+1.)*(1.-p1[j])*0.5 - 1.;
1763: x[pt * dim + i] = p1[j];
1764: for (c = 0; c < Nc; c++) w[pt*Nc + c] *= mul * w1[j];
1765: }
1766: }
1767: }
1768: totprev *= npoints;
1769: totrem /= npoints;
1770: }
1771: PetscFree2(p1, w1);
1772: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1773: PetscQuadratureSetOrder(*q, 2*npoints-1);
1774: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1775: PetscObjectChangeTypeName((PetscObject)*q,"StroudConical");
1776: return(0);
1777: }
1779: /*@
1780: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
1782: Not Collective
1784: Input Parameters:
1785: + dim - The cell dimension
1786: . level - The number of points in one dimension, 2^l
1787: . a - left end of interval (often-1)
1788: - b - right end of interval (often +1)
1790: Output Parameter:
1791: . q - A PetscQuadrature object
1793: Level: intermediate
1795: .seealso: PetscDTGaussTensorQuadrature()
1796: @*/
1797: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1798: {
1799: const PetscInt p = 16; /* Digits of precision in the evaluation */
1800: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1801: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1802: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1803: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
1804: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
1805: PetscReal *x, *w;
1806: PetscInt K, k, npoints;
1807: PetscErrorCode ierr;
1810: if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1811: if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1812: /* Find K such that the weights are < 32 digits of precision */
1813: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1814: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1815: }
1816: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1817: PetscQuadratureSetOrder(*q, 2*K+1);
1818: npoints = 2*K-1;
1819: PetscMalloc1(npoints*dim, &x);
1820: PetscMalloc1(npoints, &w);
1821: /* Center term */
1822: x[0] = beta;
1823: w[0] = 0.5*alpha*PETSC_PI;
1824: for (k = 1; k < K; ++k) {
1825: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1826: xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1827: x[2*k-1] = -alpha*xk+beta;
1828: w[2*k-1] = wk;
1829: x[2*k+0] = alpha*xk+beta;
1830: w[2*k+0] = wk;
1831: }
1832: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1833: return(0);
1834: }
1836: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1837: {
1838: const PetscInt p = 16; /* Digits of precision in the evaluation */
1839: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1840: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1841: PetscReal h = 1.0; /* Step size, length between x_k */
1842: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1843: PetscReal osum = 0.0; /* Integral on last level */
1844: PetscReal psum = 0.0; /* Integral on the level before the last level */
1845: PetscReal sum; /* Integral on current level */
1846: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1847: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1848: PetscReal wk; /* Quadrature weight at x_k */
1849: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1850: PetscInt d; /* Digits of precision in the integral */
1853: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1854: /* Center term */
1855: func(beta, &lval);
1856: sum = 0.5*alpha*PETSC_PI*lval;
1857: /* */
1858: do {
1859: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1860: PetscInt k = 1;
1862: ++l;
1863: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1864: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1865: psum = osum;
1866: osum = sum;
1867: h *= 0.5;
1868: sum *= 0.5;
1869: do {
1870: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1871: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1872: lx = -alpha*(1.0 - yk)+beta;
1873: rx = alpha*(1.0 - yk)+beta;
1874: func(lx, &lval);
1875: func(rx, &rval);
1876: lterm = alpha*wk*lval;
1877: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1878: sum += lterm;
1879: rterm = alpha*wk*rval;
1880: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1881: sum += rterm;
1882: ++k;
1883: /* Only need to evaluate every other point on refined levels */
1884: if (l != 1) ++k;
1885: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
1887: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1888: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1889: d3 = PetscLog10Real(maxTerm) - p;
1890: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1891: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1892: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1893: } while (d < digits && l < 12);
1894: *sol = sum;
1896: return(0);
1897: }
1899: #if defined(PETSC_HAVE_MPFR)
1900: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1901: {
1902: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
1903: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1904: mpfr_t alpha; /* Half-width of the integration interval */
1905: mpfr_t beta; /* Center of the integration interval */
1906: mpfr_t h; /* Step size, length between x_k */
1907: mpfr_t osum; /* Integral on last level */
1908: mpfr_t psum; /* Integral on the level before the last level */
1909: mpfr_t sum; /* Integral on current level */
1910: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1911: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1912: mpfr_t wk; /* Quadrature weight at x_k */
1913: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1914: PetscInt d; /* Digits of precision in the integral */
1915: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
1918: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1919: /* Create high precision storage */
1920: 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);
1921: /* Initialization */
1922: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1923: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
1924: mpfr_set_d(osum, 0.0, MPFR_RNDN);
1925: mpfr_set_d(psum, 0.0, MPFR_RNDN);
1926: mpfr_set_d(h, 1.0, MPFR_RNDN);
1927: mpfr_const_pi(pi2, MPFR_RNDN);
1928: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1929: /* Center term */
1930: func(0.5*(b+a), &lval);
1931: mpfr_set(sum, pi2, MPFR_RNDN);
1932: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1933: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1934: /* */
1935: do {
1936: PetscReal d1, d2, d3, d4;
1937: PetscInt k = 1;
1939: ++l;
1940: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1941: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1942: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1943: mpfr_set(psum, osum, MPFR_RNDN);
1944: mpfr_set(osum, sum, MPFR_RNDN);
1945: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
1946: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1947: do {
1948: mpfr_set_si(kh, k, MPFR_RNDN);
1949: mpfr_mul(kh, kh, h, MPFR_RNDN);
1950: /* Weight */
1951: mpfr_set(wk, h, MPFR_RNDN);
1952: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1953: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1954: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1955: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1956: mpfr_sqr(tmp, tmp, MPFR_RNDN);
1957: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1958: mpfr_div(wk, wk, tmp, MPFR_RNDN);
1959: /* Abscissa */
1960: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1961: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1962: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1963: mpfr_exp(tmp, msinh, MPFR_RNDN);
1964: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1965: /* Quadrature points */
1966: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1967: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1968: mpfr_add(lx, lx, beta, MPFR_RNDU);
1969: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1970: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1971: mpfr_add(rx, rx, beta, MPFR_RNDD);
1972: /* Evaluation */
1973: func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1974: func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1975: /* Update */
1976: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1977: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1978: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1979: mpfr_abs(tmp, tmp, MPFR_RNDN);
1980: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1981: mpfr_set(curTerm, tmp, MPFR_RNDN);
1982: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1983: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1984: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1985: mpfr_abs(tmp, tmp, MPFR_RNDN);
1986: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1987: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1988: ++k;
1989: /* Only need to evaluate every other point on refined levels */
1990: if (l != 1) ++k;
1991: mpfr_log10(tmp, wk, MPFR_RNDN);
1992: mpfr_abs(tmp, tmp, MPFR_RNDN);
1993: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1994: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1995: mpfr_abs(tmp, tmp, MPFR_RNDN);
1996: mpfr_log10(tmp, tmp, MPFR_RNDN);
1997: d1 = mpfr_get_d(tmp, MPFR_RNDN);
1998: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1999: mpfr_abs(tmp, tmp, MPFR_RNDN);
2000: mpfr_log10(tmp, tmp, MPFR_RNDN);
2001: d2 = mpfr_get_d(tmp, MPFR_RNDN);
2002: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2003: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2004: mpfr_log10(tmp, curTerm, MPFR_RNDN);
2005: d4 = mpfr_get_d(tmp, MPFR_RNDN);
2006: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
2007: } while (d < digits && l < 8);
2008: *sol = mpfr_get_d(sum, MPFR_RNDN);
2009: /* Cleanup */
2010: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2011: return(0);
2012: }
2013: #else
2015: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
2016: {
2017: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2018: }
2019: #endif
2021: /* Overwrites A. Can only handle full-rank problems with m>=n
2022: * A in column-major format
2023: * Ainv in row-major format
2024: * tau has length m
2025: * worksize must be >= max(1,n)
2026: */
2027: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2028: {
2030: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2031: PetscScalar *A,*Ainv,*R,*Q,Alpha;
2034: #if defined(PETSC_USE_COMPLEX)
2035: {
2036: PetscInt i,j;
2037: PetscMalloc2(m*n,&A,m*n,&Ainv);
2038: for (j=0; j<n; j++) {
2039: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
2040: }
2041: mstride = m;
2042: }
2043: #else
2044: A = A_in;
2045: Ainv = Ainv_out;
2046: #endif
2048: PetscBLASIntCast(m,&M);
2049: PetscBLASIntCast(n,&N);
2050: PetscBLASIntCast(mstride,&lda);
2051: PetscBLASIntCast(worksize,&ldwork);
2052: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2053: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2054: PetscFPTrapPop();
2055: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2056: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2058: /* Extract an explicit representation of Q */
2059: Q = Ainv;
2060: PetscArraycpy(Q,A,mstride*n);
2061: K = N; /* full rank */
2062: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2063: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2065: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2066: Alpha = 1.0;
2067: ldb = lda;
2068: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
2069: /* Ainv is Q, overwritten with inverse */
2071: #if defined(PETSC_USE_COMPLEX)
2072: {
2073: PetscInt i;
2074: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2075: PetscFree2(A,Ainv);
2076: }
2077: #endif
2078: return(0);
2079: }
2081: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2082: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
2083: {
2085: PetscReal *Bv;
2086: PetscInt i,j;
2089: PetscMalloc1((ninterval+1)*ndegree,&Bv);
2090: /* Point evaluation of L_p on all the source vertices */
2091: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
2092: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2093: for (i=0; i<ninterval; i++) {
2094: for (j=0; j<ndegree; j++) {
2095: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2096: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2097: }
2098: }
2099: PetscFree(Bv);
2100: return(0);
2101: }
2103: /*@
2104: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2106: Not Collective
2108: Input Parameters:
2109: + degree - degree of reconstruction polynomial
2110: . nsource - number of source intervals
2111: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2112: . ntarget - number of target intervals
2113: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
2115: Output Parameter:
2116: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2118: Level: advanced
2120: .seealso: PetscDTLegendreEval()
2121: @*/
2122: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
2123: {
2125: PetscInt i,j,k,*bdegrees,worksize;
2126: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
2127: PetscScalar *tau,*work;
2133: if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource);
2134: if (PetscDefined(USE_DEBUG)) {
2135: for (i=0; i<nsource; i++) {
2136: if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%g,%g)",i,(double)sourcex[i],(double)sourcex[i+1]);
2137: }
2138: for (i=0; i<ntarget; i++) {
2139: if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%g,%g)",i,(double)targetx[i],(double)targetx[i+1]);
2140: }
2141: }
2142: xmin = PetscMin(sourcex[0],targetx[0]);
2143: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
2144: center = (xmin + xmax)/2;
2145: hscale = (xmax - xmin)/2;
2146: worksize = nsource;
2147: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
2148: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
2149: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
2150: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
2151: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
2152: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
2153: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
2154: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
2155: for (i=0; i<ntarget; i++) {
2156: PetscReal rowsum = 0;
2157: for (j=0; j<nsource; j++) {
2158: PetscReal sum = 0;
2159: for (k=0; k<degree+1; k++) {
2160: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
2161: }
2162: R[i*nsource+j] = sum;
2163: rowsum += sum;
2164: }
2165: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
2166: }
2167: PetscFree4(bdegrees,sourcey,Bsource,work);
2168: PetscFree4(tau,Bsinv,targety,Btarget);
2169: return(0);
2170: }
2172: /*@C
2173: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2175: Not Collective
2177: Input Parameters:
2178: + n - the number of GLL nodes
2179: . nodes - the GLL nodes
2180: . weights - the GLL weights
2181: - f - the function values at the nodes
2183: Output Parameter:
2184: . in - the value of the integral
2186: Level: beginner
2188: .seealso: PetscDTGaussLobattoLegendreQuadrature()
2190: @*/
2191: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
2192: {
2193: PetscInt i;
2196: *in = 0.;
2197: for (i=0; i<n; i++) {
2198: *in += f[i]*f[i]*weights[i];
2199: }
2200: return(0);
2201: }
2203: /*@C
2204: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2206: Not Collective
2208: Input Parameters:
2209: + n - the number of GLL nodes
2210: . nodes - the GLL nodes
2211: - weights - the GLL weights
2213: Output Parameter:
2214: . A - the stiffness element
2216: Level: beginner
2218: Notes:
2219: Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()
2221: 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)
2223: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2225: @*/
2226: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2227: {
2228: PetscReal **A;
2229: PetscErrorCode ierr;
2230: const PetscReal *gllnodes = nodes;
2231: const PetscInt p = n-1;
2232: PetscReal z0,z1,z2 = -1,x,Lpj,Lpr;
2233: PetscInt i,j,nn,r;
2236: PetscMalloc1(n,&A);
2237: PetscMalloc1(n*n,&A[0]);
2238: for (i=1; i<n; i++) A[i] = A[i-1]+n;
2240: for (j=1; j<p; j++) {
2241: x = gllnodes[j];
2242: z0 = 1.;
2243: z1 = x;
2244: for (nn=1; nn<p; nn++) {
2245: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2246: z0 = z1;
2247: z1 = z2;
2248: }
2249: Lpj=z2;
2250: for (r=1; r<p; r++) {
2251: if (r == j) {
2252: A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
2253: } else {
2254: x = gllnodes[r];
2255: z0 = 1.;
2256: z1 = x;
2257: for (nn=1; nn<p; nn++) {
2258: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2259: z0 = z1;
2260: z1 = z2;
2261: }
2262: Lpr = z2;
2263: A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
2264: }
2265: }
2266: }
2267: for (j=1; j<p+1; j++) {
2268: x = gllnodes[j];
2269: z0 = 1.;
2270: z1 = x;
2271: for (nn=1; nn<p; nn++) {
2272: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2273: z0 = z1;
2274: z1 = z2;
2275: }
2276: Lpj = z2;
2277: A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
2278: A[0][j] = A[j][0];
2279: }
2280: for (j=0; j<p; j++) {
2281: x = gllnodes[j];
2282: z0 = 1.;
2283: z1 = x;
2284: for (nn=1; nn<p; nn++) {
2285: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2286: z0 = z1;
2287: z1 = z2;
2288: }
2289: Lpj=z2;
2291: A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
2292: A[j][p] = A[p][j];
2293: }
2294: A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
2295: A[p][p]=A[0][0];
2296: *AA = A;
2297: return(0);
2298: }
2300: /*@C
2301: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
2303: Not Collective
2305: Input Parameters:
2306: + n - the number of GLL nodes
2307: . nodes - the GLL nodes
2308: . weights - the GLL weightss
2309: - A - the stiffness element
2311: Level: beginner
2313: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()
2315: @*/
2316: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2317: {
2321: PetscFree((*AA)[0]);
2322: PetscFree(*AA);
2323: *AA = NULL;
2324: return(0);
2325: }
2327: /*@C
2328: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
2330: Not Collective
2332: Input Parameter:
2333: + n - the number of GLL nodes
2334: . nodes - the GLL nodes
2335: . weights - the GLL weights
2337: Output Parameters:
2338: . AA - the stiffness element
2339: - AAT - the transpose of AA (pass in NULL if you do not need this array)
2341: Level: beginner
2343: Notes:
2344: Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()
2346: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2348: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2350: @*/
2351: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2352: {
2353: PetscReal **A, **AT = NULL;
2354: PetscErrorCode ierr;
2355: const PetscReal *gllnodes = nodes;
2356: const PetscInt p = n-1;
2357: PetscReal Li, Lj,d0;
2358: PetscInt i,j;
2361: PetscMalloc1(n,&A);
2362: PetscMalloc1(n*n,&A[0]);
2363: for (i=1; i<n; i++) A[i] = A[i-1]+n;
2365: if (AAT) {
2366: PetscMalloc1(n,&AT);
2367: PetscMalloc1(n*n,&AT[0]);
2368: for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
2369: }
2371: if (n==1) {A[0][0] = 0.;}
2372: d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
2373: for (i=0; i<n; i++) {
2374: for (j=0; j<n; j++) {
2375: A[i][j] = 0.;
2376: PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);
2377: PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);
2378: if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
2379: if ((j==i) && (i==0)) A[i][j] = -d0;
2380: if (j==i && i==p) A[i][j] = d0;
2381: if (AT) AT[j][i] = A[i][j];
2382: }
2383: }
2384: if (AAT) *AAT = AT;
2385: *AA = A;
2386: return(0);
2387: }
2389: /*@C
2390: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()
2392: Not Collective
2394: Input Parameters:
2395: + n - the number of GLL nodes
2396: . nodes - the GLL nodes
2397: . weights - the GLL weights
2398: . AA - the stiffness element
2399: - AAT - the transpose of the element
2401: Level: beginner
2403: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()
2405: @*/
2406: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2407: {
2411: PetscFree((*AA)[0]);
2412: PetscFree(*AA);
2413: *AA = NULL;
2414: if (*AAT) {
2415: PetscFree((*AAT)[0]);
2416: PetscFree(*AAT);
2417: *AAT = NULL;
2418: }
2419: return(0);
2420: }
2422: /*@C
2423: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
2425: Not Collective
2427: Input Parameters:
2428: + n - the number of GLL nodes
2429: . nodes - the GLL nodes
2430: - weights - the GLL weightss
2432: Output Parameter:
2433: . AA - the stiffness element
2435: Level: beginner
2437: Notes:
2438: Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()
2440: This is the same as the Gradient operator multiplied by the diagonal mass matrix
2442: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2444: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()
2446: @*/
2447: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2448: {
2449: PetscReal **D;
2450: PetscErrorCode ierr;
2451: const PetscReal *gllweights = weights;
2452: const PetscInt glln = n;
2453: PetscInt i,j;
2456: PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
2457: for (i=0; i<glln; i++) {
2458: for (j=0; j<glln; j++) {
2459: D[i][j] = gllweights[i]*D[i][j];
2460: }
2461: }
2462: *AA = D;
2463: return(0);
2464: }
2466: /*@C
2467: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
2469: Not Collective
2471: Input Parameters:
2472: + n - the number of GLL nodes
2473: . nodes - the GLL nodes
2474: . weights - the GLL weights
2475: - A - advection
2477: Level: beginner
2479: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()
2481: @*/
2482: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2483: {
2487: PetscFree((*AA)[0]);
2488: PetscFree(*AA);
2489: *AA = NULL;
2490: return(0);
2491: }
2493: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2494: {
2495: PetscReal **A;
2496: PetscErrorCode ierr;
2497: const PetscReal *gllweights = weights;
2498: const PetscInt glln = n;
2499: PetscInt i,j;
2502: PetscMalloc1(glln,&A);
2503: PetscMalloc1(glln*glln,&A[0]);
2504: for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
2505: if (glln==1) {A[0][0] = 0.;}
2506: for (i=0; i<glln; i++) {
2507: for (j=0; j<glln; j++) {
2508: A[i][j] = 0.;
2509: if (j==i) A[i][j] = gllweights[i];
2510: }
2511: }
2512: *AA = A;
2513: return(0);
2514: }
2516: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2517: {
2521: PetscFree((*AA)[0]);
2522: PetscFree(*AA);
2523: *AA = NULL;
2524: return(0);
2525: }
2527: /*@
2528: PetscDTIndexToBary - convert an index into a barycentric coordinate.
2530: Input Parameters:
2531: + 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)
2532: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2533: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
2535: Output Parameter:
2536: . coord - will be filled with the barycentric coordinate
2538: Level: beginner
2540: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2541: least significant and the last index is the most significant.
2543: .seealso: PetscDTBaryToIndex()
2544: @*/
2545: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
2546: {
2547: PetscInt c, d, s, total, subtotal, nexttotal;
2550: if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2551: if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
2552: if (!len) {
2553: if (!sum && !index) return(0);
2554: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2555: }
2556: for (c = 1, total = 1; c <= len; c++) {
2557: /* total is the number of ways to have a tuple of length c with sum */
2558: if (index < total) break;
2559: total = (total * (sum + c)) / c;
2560: }
2561: if (c > len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
2562: for (d = c; d < len; d++) coord[d] = 0;
2563: for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
2564: /* subtotal is the number of ways to have a tuple of length c with sum s */
2565: /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
2566: if ((index + subtotal) >= total) {
2567: coord[--c] = sum - s;
2568: index -= (total - subtotal);
2569: sum = s;
2570: total = nexttotal;
2571: subtotal = 1;
2572: nexttotal = 1;
2573: s = 0;
2574: } else {
2575: subtotal = (subtotal * (c + s)) / (s + 1);
2576: nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
2577: s++;
2578: }
2579: }
2580: return(0);
2581: }
2583: /*@
2584: PetscDTBaryToIndex - convert a barycentric coordinate to an index
2586: Input Parameters:
2587: + 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)
2588: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2589: - coord - a barycentric coordinate with the given length and sum
2591: Output Parameter:
2592: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
2594: Level: beginner
2596: Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2597: least significant and the last index is the most significant.
2599: .seealso: PetscDTIndexToBary
2600: @*/
2601: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
2602: {
2603: PetscInt c;
2604: PetscInt i;
2605: PetscInt total;
2608: if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2609: if (!len) {
2610: if (!sum) {
2611: *index = 0;
2612: return(0);
2613: }
2614: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2615: }
2616: for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
2617: i = total - 1;
2618: c = len - 1;
2619: sum -= coord[c];
2620: while (sum > 0) {
2621: PetscInt subtotal;
2622: PetscInt s;
2624: for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
2625: i -= subtotal;
2626: sum -= coord[--c];
2627: }
2628: *index = i;
2629: return(0);
2630: }