Actual source code: dt.c
petsc-3.7.5 2017-01-01
1: /* Discretization tools */
3: #include <petscconf.h>
4: #if defined(PETSC_HAVE_MATHIMF_H)
5: #include <mathimf.h> /* this needs to be included before math.h */
6: #endif
7: #ifdef PETSC_HAVE_MPFR
8: #include <mpfr.h>
9: #endif
11: #include <petscdt.h> /*I "petscdt.h" I*/
12: #include <petscblaslapack.h>
13: #include <petsc/private/petscimpl.h>
14: #include <petsc/private/dtimpl.h>
15: #include <petscviewer.h>
16: #include <petscdmplex.h>
17: #include <petscdmshell.h>
19: static PetscBool GaussCite = PETSC_FALSE;
20: const char GaussCitation[] = "@article{GolubWelsch1969,\n"
21: " author = {Golub and Welsch},\n"
22: " title = {Calculation of Quadrature Rules},\n"
23: " journal = {Math. Comp.},\n"
24: " volume = {23},\n"
25: " number = {106},\n"
26: " pages = {221--230},\n"
27: " year = {1969}\n}\n";
31: /*@
32: PetscQuadratureCreate - Create a PetscQuadrature object
34: Collective on MPI_Comm
36: Input Parameter:
37: . comm - The communicator for the PetscQuadrature object
39: Output Parameter:
40: . q - The PetscQuadrature object
42: Level: beginner
44: .keywords: PetscQuadrature, quadrature, create
45: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
46: @*/
47: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
48: {
53: DMInitializePackage();
54: PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
55: (*q)->dim = -1;
56: (*q)->order = -1;
57: (*q)->numPoints = 0;
58: (*q)->points = NULL;
59: (*q)->weights = NULL;
60: return(0);
61: }
65: /*@
66: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
68: Collective on PetscQuadrature
70: Input Parameter:
71: . q - The PetscQuadrature object
73: Output Parameter:
74: . r - The new PetscQuadrature object
76: Level: beginner
78: .keywords: PetscQuadrature, quadrature, clone
79: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
80: @*/
81: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
82: {
83: PetscInt order, dim, Nq;
84: const PetscReal *points, *weights;
85: PetscReal *p, *w;
86: PetscErrorCode ierr;
90: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
91: PetscQuadratureGetOrder(q, &order);
92: PetscQuadratureSetOrder(*r, order);
93: PetscQuadratureGetData(q, &dim, &Nq, &points, &weights);
94: PetscMalloc1(Nq*dim, &p);
95: PetscMalloc1(Nq, &w);
96: PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
97: PetscMemcpy(w, weights, Nq * sizeof(PetscReal));
98: PetscQuadratureSetData(*r, dim, Nq, p, w);
99: return(0);
100: }
104: /*@
105: PetscQuadratureDestroy - Destroys a PetscQuadrature object
107: Collective on PetscQuadrature
109: Input Parameter:
110: . q - The PetscQuadrature object
112: Level: beginner
114: .keywords: PetscQuadrature, quadrature, destroy
115: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
116: @*/
117: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
118: {
122: if (!*q) return(0);
124: if (--((PetscObject)(*q))->refct > 0) {
125: *q = NULL;
126: return(0);
127: }
128: PetscFree((*q)->points);
129: PetscFree((*q)->weights);
130: PetscHeaderDestroy(q);
131: return(0);
132: }
136: /*@
137: PetscQuadratureGetOrder - Return the quadrature information
139: Not collective
141: Input Parameter:
142: . q - The PetscQuadrature object
144: Output Parameter:
145: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
147: Output Parameter:
149: Level: intermediate
151: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
152: @*/
153: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
154: {
158: *order = q->order;
159: return(0);
160: }
164: /*@
165: PetscQuadratureSetOrder - Return the quadrature information
167: Not collective
169: Input Parameters:
170: + q - The PetscQuadrature object
171: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
173: Level: intermediate
175: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
176: @*/
177: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
178: {
181: q->order = order;
182: return(0);
183: }
187: /*@C
188: PetscQuadratureGetData - Returns the data defining the quadrature
190: Not collective
192: Input Parameter:
193: . q - The PetscQuadrature object
195: Output Parameters:
196: + dim - The spatial dimension
197: . npoints - The number of quadrature points
198: . points - The coordinates of each quadrature point
199: - weights - The weight of each quadrature point
201: Level: intermediate
203: .keywords: PetscQuadrature, quadrature
204: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
205: @*/
206: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
207: {
210: if (dim) {
212: *dim = q->dim;
213: }
214: if (npoints) {
216: *npoints = q->numPoints;
217: }
218: if (points) {
220: *points = q->points;
221: }
222: if (weights) {
224: *weights = q->weights;
225: }
226: return(0);
227: }
231: /*@C
232: PetscQuadratureSetData - Sets the data defining the quadrature
234: Not collective
236: Input Parameters:
237: + q - The PetscQuadrature object
238: . dim - The spatial dimension
239: . npoints - The number of quadrature points
240: . points - The coordinates of each quadrature point
241: - weights - The weight of each quadrature point
243: Level: intermediate
245: .keywords: PetscQuadrature, quadrature
246: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
247: @*/
248: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
249: {
252: if (dim >= 0) q->dim = dim;
253: if (npoints >= 0) q->numPoints = npoints;
254: if (points) {
256: q->points = points;
257: }
258: if (weights) {
260: q->weights = weights;
261: }
262: return(0);
263: }
267: /*@C
268: PetscQuadratureView - Views a PetscQuadrature object
270: Collective on PetscQuadrature
272: Input Parameters:
273: + q - The PetscQuadrature object
274: - viewer - The PetscViewer object
276: Level: beginner
278: .keywords: PetscQuadrature, quadrature, view
279: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
280: @*/
281: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
282: {
283: PetscInt q, d;
287: PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
288: PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n (", quad->numPoints);
289: for (q = 0; q < quad->numPoints; ++q) {
290: for (d = 0; d < quad->dim; ++d) {
291: if (d) PetscViewerASCIIPrintf(viewer, ", ");
292: PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
293: }
294: PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
295: }
296: return(0);
297: }
301: /*@C
302: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
304: Not collective
306: Input Parameter:
307: + q - The original PetscQuadrature
308: . numSubelements - The number of subelements the original element is divided into
309: . v0 - An array of the initial points for each subelement
310: - jac - An array of the Jacobian mappings from the reference to each subelement
312: Output Parameters:
313: . dim - The dimension
315: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
317: Level: intermediate
319: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
320: @*/
321: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
322: {
323: const PetscReal *points, *weights;
324: PetscReal *pointsRef, *weightsRef;
325: PetscInt dim, order, npoints, npointsRef, c, p, d, e;
326: PetscErrorCode ierr;
333: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
334: PetscQuadratureGetOrder(q, &order);
335: PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
336: npointsRef = npoints*numSubelements;
337: PetscMalloc1(npointsRef*dim,&pointsRef);
338: PetscMalloc1(npointsRef,&weightsRef);
339: for (c = 0; c < numSubelements; ++c) {
340: for (p = 0; p < npoints; ++p) {
341: for (d = 0; d < dim; ++d) {
342: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
343: for (e = 0; e < dim; ++e) {
344: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
345: }
346: }
347: /* Could also use detJ here */
348: weightsRef[c*npoints+p] = weights[p]/numSubelements;
349: }
350: }
351: PetscQuadratureSetOrder(*qref, order);
352: PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
353: return(0);
354: }
358: /*@
359: PetscDTLegendreEval - evaluate Legendre polynomial at points
361: Not Collective
363: Input Arguments:
364: + npoints - number of spatial points to evaluate at
365: . points - array of locations to evaluate at
366: . ndegree - number of basis degrees to evaluate
367: - degrees - sorted array of degrees to evaluate
369: Output Arguments:
370: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
371: . D - row-oriented derivative evaluation matrix (or NULL)
372: - D2 - row-oriented second derivative evaluation matrix (or NULL)
374: Level: intermediate
376: .seealso: PetscDTGaussQuadrature()
377: @*/
378: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
379: {
380: PetscInt i,maxdegree;
383: if (!npoints || !ndegree) return(0);
384: maxdegree = degrees[ndegree-1];
385: for (i=0; i<npoints; i++) {
386: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
387: PetscInt j,k;
388: x = points[i];
389: pm2 = 0;
390: pm1 = 1;
391: pd2 = 0;
392: pd1 = 0;
393: pdd2 = 0;
394: pdd1 = 0;
395: k = 0;
396: if (degrees[k] == 0) {
397: if (B) B[i*ndegree+k] = pm1;
398: if (D) D[i*ndegree+k] = pd1;
399: if (D2) D2[i*ndegree+k] = pdd1;
400: k++;
401: }
402: for (j=1; j<=maxdegree; j++,k++) {
403: PetscReal p,d,dd;
404: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
405: d = pd2 + (2*j-1)*pm1;
406: dd = pdd2 + (2*j-1)*pd1;
407: pm2 = pm1;
408: pm1 = p;
409: pd2 = pd1;
410: pd1 = d;
411: pdd2 = pdd1;
412: pdd1 = dd;
413: if (degrees[k] == j) {
414: if (B) B[i*ndegree+k] = p;
415: if (D) D[i*ndegree+k] = d;
416: if (D2) D2[i*ndegree+k] = dd;
417: }
418: }
419: }
420: return(0);
421: }
425: /*@
426: PetscDTGaussQuadrature - create Gauss quadrature
428: Not Collective
430: Input Arguments:
431: + npoints - number of points
432: . a - left end of interval (often-1)
433: - b - right end of interval (often +1)
435: Output Arguments:
436: + x - quadrature points
437: - w - quadrature weights
439: Level: intermediate
441: References:
442: . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
444: .seealso: PetscDTLegendreEval()
445: @*/
446: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
447: {
449: PetscInt i;
450: PetscReal *work;
451: PetscScalar *Z;
452: PetscBLASInt N,LDZ,info;
455: PetscCitationsRegister(GaussCitation, &GaussCite);
456: /* Set up the Golub-Welsch system */
457: for (i=0; i<npoints; i++) {
458: x[i] = 0; /* diagonal is 0 */
459: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
460: }
461: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
462: PetscBLASIntCast(npoints,&N);
463: LDZ = N;
464: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
465: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
466: PetscFPTrapPop();
467: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
469: for (i=0; i<(npoints+1)/2; i++) {
470: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
471: x[i] = (a+b)/2 - y*(b-a)/2;
472: if (x[i] == -0.0) x[i] = 0.0;
473: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
475: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
476: }
477: PetscFree2(Z,work);
478: return(0);
479: }
483: /*@
484: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
486: Not Collective
488: Input Arguments:
489: + dim - The spatial dimension
490: . npoints - number of points in one dimension
491: . a - left end of interval (often-1)
492: - b - right end of interval (often +1)
494: Output Argument:
495: . q - A PetscQuadrature object
497: Level: intermediate
499: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
500: @*/
501: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
502: {
503: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
504: PetscReal *x, *w, *xw, *ww;
508: PetscMalloc1(totpoints*dim,&x);
509: PetscMalloc1(totpoints,&w);
510: /* Set up the Golub-Welsch system */
511: switch (dim) {
512: case 0:
513: PetscFree(x);
514: PetscFree(w);
515: PetscMalloc1(1, &x);
516: PetscMalloc1(1, &w);
517: x[0] = 0.0;
518: w[0] = 1.0;
519: break;
520: case 1:
521: PetscDTGaussQuadrature(npoints, a, b, x, w);
522: break;
523: case 2:
524: PetscMalloc2(npoints,&xw,npoints,&ww);
525: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
526: for (i = 0; i < npoints; ++i) {
527: for (j = 0; j < npoints; ++j) {
528: x[(i*npoints+j)*dim+0] = xw[i];
529: x[(i*npoints+j)*dim+1] = xw[j];
530: w[i*npoints+j] = ww[i] * ww[j];
531: }
532: }
533: PetscFree2(xw,ww);
534: break;
535: case 3:
536: PetscMalloc2(npoints,&xw,npoints,&ww);
537: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
538: for (i = 0; i < npoints; ++i) {
539: for (j = 0; j < npoints; ++j) {
540: for (k = 0; k < npoints; ++k) {
541: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
542: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
543: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
544: w[(i*npoints+j)*npoints+k] = ww[i] * ww[j] * ww[k];
545: }
546: }
547: }
548: PetscFree2(xw,ww);
549: break;
550: default:
551: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
552: }
553: PetscQuadratureCreate(PETSC_COMM_SELF, q);
554: PetscQuadratureSetOrder(*q, npoints);
555: PetscQuadratureSetData(*q, dim, totpoints, x, w);
556: return(0);
557: }
561: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
562: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
563: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
564: {
565: PetscReal f = 1.0;
566: PetscInt i;
569: for (i = 1; i < n+1; ++i) f *= i;
570: *factorial = f;
571: return(0);
572: }
576: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
577: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
578: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
579: {
580: PetscReal apb, pn1, pn2;
581: PetscInt k;
584: if (!n) {*P = 1.0; return(0);}
585: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
586: apb = a + b;
587: pn2 = 1.0;
588: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
589: *P = 0.0;
590: for (k = 2; k < n+1; ++k) {
591: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
592: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
593: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
594: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
596: a2 = a2 / a1;
597: a3 = a3 / a1;
598: a4 = a4 / a1;
599: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
600: pn2 = pn1;
601: pn1 = *P;
602: }
603: return(0);
604: }
608: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
609: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
610: {
611: PetscReal nP;
615: if (!n) {*P = 0.0; return(0);}
616: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
617: *P = 0.5 * (a + b + n + 1) * nP;
618: return(0);
619: }
623: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
624: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
625: {
627: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
628: *eta = y;
629: return(0);
630: }
634: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
635: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
636: {
638: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
639: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
640: *zeta = z;
641: return(0);
642: }
646: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
647: {
648: PetscInt maxIter = 100;
649: PetscReal eps = 1.0e-8;
650: PetscReal a1, a2, a3, a4, a5, a6;
651: PetscInt k;
656: a1 = PetscPowReal(2.0, a+b+1);
657: #if defined(PETSC_HAVE_TGAMMA)
658: a2 = PetscTGamma(a + npoints + 1);
659: a3 = PetscTGamma(b + npoints + 1);
660: a4 = PetscTGamma(a + b + npoints + 1);
661: #else
662: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
663: #endif
665: PetscDTFactorial_Internal(npoints, &a5);
666: a6 = a1 * a2 * a3 / a4 / a5;
667: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
668: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
669: for (k = 0; k < npoints; ++k) {
670: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
671: PetscInt j;
673: if (k > 0) r = 0.5 * (r + x[k-1]);
674: for (j = 0; j < maxIter; ++j) {
675: PetscReal s = 0.0, delta, f, fp;
676: PetscInt i;
678: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
679: PetscDTComputeJacobi(a, b, npoints, r, &f);
680: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
681: delta = f / (fp - f * s);
682: r = r - delta;
683: if (PetscAbsReal(delta) < eps) break;
684: }
685: x[k] = r;
686: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
687: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
688: }
689: return(0);
690: }
694: /*@C
695: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
697: Not Collective
699: Input Arguments:
700: + dim - The simplex dimension
701: . order - The number of points in one dimension
702: . a - left end of interval (often-1)
703: - b - right end of interval (often +1)
705: Output Argument:
706: . q - A PetscQuadrature object
708: Level: intermediate
710: References:
711: . 1. - Karniadakis and Sherwin. FIAT
713: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
714: @*/
715: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
716: {
717: PetscInt npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
718: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
719: PetscInt i, j, k;
723: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
724: PetscMalloc1(npoints*dim, &x);
725: PetscMalloc1(npoints, &w);
726: switch (dim) {
727: case 0:
728: PetscFree(x);
729: PetscFree(w);
730: PetscMalloc1(1, &x);
731: PetscMalloc1(1, &w);
732: x[0] = 0.0;
733: w[0] = 1.0;
734: break;
735: case 1:
736: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
737: break;
738: case 2:
739: PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
740: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
741: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
742: for (i = 0; i < order; ++i) {
743: for (j = 0; j < order; ++j) {
744: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
745: w[i*order+j] = 0.5 * wx[i] * wy[j];
746: }
747: }
748: PetscFree4(px,wx,py,wy);
749: break;
750: case 3:
751: PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
752: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
753: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
754: PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
755: for (i = 0; i < order; ++i) {
756: for (j = 0; j < order; ++j) {
757: for (k = 0; k < order; ++k) {
758: PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*order+j)*order+k)*3+0], &x[((i*order+j)*order+k)*3+1], &x[((i*order+j)*order+k)*3+2]);
759: w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
760: }
761: }
762: }
763: PetscFree6(px,wx,py,wy,pz,wz);
764: break;
765: default:
766: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
767: }
768: PetscQuadratureCreate(PETSC_COMM_SELF, q);
769: PetscQuadratureSetOrder(*q, order);
770: PetscQuadratureSetData(*q, dim, npoints, x, w);
771: return(0);
772: }
776: /*@C
777: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
779: Not Collective
781: Input Arguments:
782: + dim - The cell dimension
783: . level - The number of points in one dimension, 2^l
784: . a - left end of interval (often-1)
785: - b - right end of interval (often +1)
787: Output Argument:
788: . q - A PetscQuadrature object
790: Level: intermediate
792: .seealso: PetscDTGaussTensorQuadrature()
793: @*/
794: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
795: {
796: const PetscInt p = 16; /* Digits of precision in the evaluation */
797: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
798: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
799: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
800: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
801: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
802: PetscReal *x, *w;
803: PetscInt K, k, npoints;
804: PetscErrorCode ierr;
807: if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
808: if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
809: /* Find K such that the weights are < 32 digits of precision */
810: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
811: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
812: }
813: PetscQuadratureCreate(PETSC_COMM_SELF, q);
814: PetscQuadratureSetOrder(*q, 2*K+1);
815: npoints = 2*K-1;
816: PetscMalloc1(npoints*dim, &x);
817: PetscMalloc1(npoints, &w);
818: /* Center term */
819: x[0] = beta;
820: w[0] = 0.5*alpha*PETSC_PI;
821: for (k = 1; k < K; ++k) {
822: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
823: xk = tanh(0.5*PETSC_PI*PetscSinhReal(k*h));
824: x[2*k-1] = -alpha*xk+beta;
825: w[2*k-1] = wk;
826: x[2*k+0] = alpha*xk+beta;
827: w[2*k+0] = wk;
828: }
829: PetscQuadratureSetData(*q, dim, npoints, x, w);
830: return(0);
831: }
835: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
836: {
837: const PetscInt p = 16; /* Digits of precision in the evaluation */
838: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
839: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
840: PetscReal h = 1.0; /* Step size, length between x_k */
841: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
842: PetscReal osum = 0.0; /* Integral on last level */
843: PetscReal psum = 0.0; /* Integral on the level before the last level */
844: PetscReal sum; /* Integral on current level */
845: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
846: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
847: PetscReal wk; /* Quadrature weight at x_k */
848: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
849: PetscInt d; /* Digits of precision in the integral */
852: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
853: /* Center term */
854: func(beta, &lval);
855: sum = 0.5*alpha*PETSC_PI*lval;
856: /* */
857: do {
858: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
859: PetscInt k = 1;
861: ++l;
862: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
863: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
864: psum = osum;
865: osum = sum;
866: h *= 0.5;
867: sum *= 0.5;
868: do {
869: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
870: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
871: lx = -alpha*(1.0 - yk)+beta;
872: rx = alpha*(1.0 - yk)+beta;
873: func(lx, &lval);
874: func(rx, &rval);
875: lterm = alpha*wk*lval;
876: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
877: sum += lterm;
878: rterm = alpha*wk*rval;
879: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
880: sum += rterm;
881: /* if (l == 1) printf("k is %d and sum is %15.15f and wk is %15.15f\n", k, sum, wk); */
882: ++k;
883: /* Only need to evaluate every other point on refined levels */
884: if (l != 1) ++k;
885: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
887: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
888: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
889: d3 = PetscLog10Real(maxTerm) - p;
890: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
891: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
892: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
893: } while (d < digits && l < 12);
894: *sol = sum;
896: return(0);
897: }
899: #ifdef PETSC_HAVE_MPFR
902: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
903: {
904: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
905: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
906: mpfr_t alpha; /* Half-width of the integration interval */
907: mpfr_t beta; /* Center of the integration interval */
908: mpfr_t h; /* Step size, length between x_k */
909: mpfr_t osum; /* Integral on last level */
910: mpfr_t psum; /* Integral on the level before the last level */
911: mpfr_t sum; /* Integral on current level */
912: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
913: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
914: mpfr_t wk; /* Quadrature weight at x_k */
915: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
916: PetscInt d; /* Digits of precision in the integral */
917: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
920: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
921: /* Create high precision storage */
922: 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);
923: /* Initialization */
924: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
925: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
926: mpfr_set_d(osum, 0.0, MPFR_RNDN);
927: mpfr_set_d(psum, 0.0, MPFR_RNDN);
928: mpfr_set_d(h, 1.0, MPFR_RNDN);
929: mpfr_const_pi(pi2, MPFR_RNDN);
930: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
931: /* Center term */
932: func(0.5*(b+a), &lval);
933: mpfr_set(sum, pi2, MPFR_RNDN);
934: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
935: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
936: /* */
937: do {
938: PetscReal d1, d2, d3, d4;
939: PetscInt k = 1;
941: ++l;
942: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
943: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
944: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
945: mpfr_set(psum, osum, MPFR_RNDN);
946: mpfr_set(osum, sum, MPFR_RNDN);
947: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
948: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
949: do {
950: mpfr_set_si(kh, k, MPFR_RNDN);
951: mpfr_mul(kh, kh, h, MPFR_RNDN);
952: /* Weight */
953: mpfr_set(wk, h, MPFR_RNDN);
954: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
955: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
956: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
957: mpfr_cosh(tmp, msinh, MPFR_RNDN);
958: mpfr_sqr(tmp, tmp, MPFR_RNDN);
959: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
960: mpfr_div(wk, wk, tmp, MPFR_RNDN);
961: /* Abscissa */
962: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
963: mpfr_cosh(tmp, msinh, MPFR_RNDN);
964: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
965: mpfr_exp(tmp, msinh, MPFR_RNDN);
966: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
967: /* Quadrature points */
968: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
969: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
970: mpfr_add(lx, lx, beta, MPFR_RNDU);
971: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
972: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
973: mpfr_add(rx, rx, beta, MPFR_RNDD);
974: /* Evaluation */
975: func(mpfr_get_d(lx, MPFR_RNDU), &lval);
976: func(mpfr_get_d(rx, MPFR_RNDD), &rval);
977: /* Update */
978: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
979: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
980: mpfr_add(sum, sum, tmp, MPFR_RNDN);
981: mpfr_abs(tmp, tmp, MPFR_RNDN);
982: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
983: mpfr_set(curTerm, tmp, MPFR_RNDN);
984: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
985: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
986: mpfr_add(sum, sum, tmp, MPFR_RNDN);
987: mpfr_abs(tmp, tmp, MPFR_RNDN);
988: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
989: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
990: /* if (l == 1) printf("k is %d and sum is %15.15f and wk is %15.15f\n", k, sum, wk); */
991: ++k;
992: /* Only need to evaluate every other point on refined levels */
993: if (l != 1) ++k;
994: mpfr_log10(tmp, wk, MPFR_RNDN);
995: mpfr_abs(tmp, tmp, MPFR_RNDN);
996: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
997: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
998: mpfr_abs(tmp, tmp, MPFR_RNDN);
999: mpfr_log10(tmp, tmp, MPFR_RNDN);
1000: d1 = mpfr_get_d(tmp, MPFR_RNDN);
1001: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1002: mpfr_abs(tmp, tmp, MPFR_RNDN);
1003: mpfr_log10(tmp, tmp, MPFR_RNDN);
1004: d2 = mpfr_get_d(tmp, MPFR_RNDN);
1005: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1006: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1007: mpfr_log10(tmp, curTerm, MPFR_RNDN);
1008: d4 = mpfr_get_d(tmp, MPFR_RNDN);
1009: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1010: } while (d < digits && l < 8);
1011: *sol = mpfr_get_d(sum, MPFR_RNDN);
1012: /* Cleanup */
1013: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1014: return(0);
1015: }
1016: #else
1019: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1020: {
1021: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1022: }
1023: #endif
1027: /* Overwrites A. Can only handle full-rank problems with m>=n
1028: * A in column-major format
1029: * Ainv in row-major format
1030: * tau has length m
1031: * worksize must be >= max(1,n)
1032: */
1033: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1034: {
1036: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1037: PetscScalar *A,*Ainv,*R,*Q,Alpha;
1040: #if defined(PETSC_USE_COMPLEX)
1041: {
1042: PetscInt i,j;
1043: PetscMalloc2(m*n,&A,m*n,&Ainv);
1044: for (j=0; j<n; j++) {
1045: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1046: }
1047: mstride = m;
1048: }
1049: #else
1050: A = A_in;
1051: Ainv = Ainv_out;
1052: #endif
1054: PetscBLASIntCast(m,&M);
1055: PetscBLASIntCast(n,&N);
1056: PetscBLASIntCast(mstride,&lda);
1057: PetscBLASIntCast(worksize,&ldwork);
1058: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1059: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1060: PetscFPTrapPop();
1061: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1062: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1064: /* Extract an explicit representation of Q */
1065: Q = Ainv;
1066: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1067: K = N; /* full rank */
1068: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1069: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1071: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1072: Alpha = 1.0;
1073: ldb = lda;
1074: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1075: /* Ainv is Q, overwritten with inverse */
1077: #if defined(PETSC_USE_COMPLEX)
1078: {
1079: PetscInt i;
1080: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1081: PetscFree2(A,Ainv);
1082: }
1083: #endif
1084: return(0);
1085: }
1089: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1090: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1091: {
1093: PetscReal *Bv;
1094: PetscInt i,j;
1097: PetscMalloc1((ninterval+1)*ndegree,&Bv);
1098: /* Point evaluation of L_p on all the source vertices */
1099: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1100: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1101: for (i=0; i<ninterval; i++) {
1102: for (j=0; j<ndegree; j++) {
1103: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1104: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1105: }
1106: }
1107: PetscFree(Bv);
1108: return(0);
1109: }
1113: /*@
1114: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
1116: Not Collective
1118: Input Arguments:
1119: + degree - degree of reconstruction polynomial
1120: . nsource - number of source intervals
1121: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1122: . ntarget - number of target intervals
1123: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
1125: Output Arguments:
1126: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
1128: Level: advanced
1130: .seealso: PetscDTLegendreEval()
1131: @*/
1132: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1133: {
1135: PetscInt i,j,k,*bdegrees,worksize;
1136: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1137: PetscScalar *tau,*work;
1143: 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);
1144: #if defined(PETSC_USE_DEBUG)
1145: for (i=0; i<nsource; i++) {
1146: 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]);
1147: }
1148: for (i=0; i<ntarget; i++) {
1149: 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]);
1150: }
1151: #endif
1152: xmin = PetscMin(sourcex[0],targetx[0]);
1153: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1154: center = (xmin + xmax)/2;
1155: hscale = (xmax - xmin)/2;
1156: worksize = nsource;
1157: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1158: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1159: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1160: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1161: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1162: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1163: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1164: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1165: for (i=0; i<ntarget; i++) {
1166: PetscReal rowsum = 0;
1167: for (j=0; j<nsource; j++) {
1168: PetscReal sum = 0;
1169: for (k=0; k<degree+1; k++) {
1170: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1171: }
1172: R[i*nsource+j] = sum;
1173: rowsum += sum;
1174: }
1175: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1176: }
1177: PetscFree4(bdegrees,sourcey,Bsource,work);
1178: PetscFree4(tau,Bsinv,targety,Btarget);
1179: return(0);
1180: }