Actual source code: dt.c
petsc-3.8.4 2018-03-24
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>
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";
29: /*@
30: PetscQuadratureCreate - Create a PetscQuadrature object
32: Collective on MPI_Comm
34: Input Parameter:
35: . comm - The communicator for the PetscQuadrature object
37: Output Parameter:
38: . q - The PetscQuadrature object
40: Level: beginner
42: .keywords: PetscQuadrature, quadrature, create
43: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
44: @*/
45: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
46: {
51: PetscSysInitializePackage();
52: PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
53: (*q)->dim = -1;
54: (*q)->Nc = 1;
55: (*q)->order = -1;
56: (*q)->numPoints = 0;
57: (*q)->points = NULL;
58: (*q)->weights = NULL;
59: return(0);
60: }
62: /*@
63: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
65: Collective on PetscQuadrature
67: Input Parameter:
68: . q - The PetscQuadrature object
70: Output Parameter:
71: . r - The new PetscQuadrature object
73: Level: beginner
75: .keywords: PetscQuadrature, quadrature, clone
76: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
77: @*/
78: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
79: {
80: PetscInt order, dim, Nc, Nq;
81: const PetscReal *points, *weights;
82: PetscReal *p, *w;
83: PetscErrorCode ierr;
87: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
88: PetscQuadratureGetOrder(q, &order);
89: PetscQuadratureSetOrder(*r, order);
90: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
91: PetscMalloc1(Nq*dim, &p);
92: PetscMalloc1(Nq*Nc, &w);
93: PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
94: PetscMemcpy(w, weights, Nc * Nq * sizeof(PetscReal));
95: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
96: return(0);
97: }
99: /*@
100: PetscQuadratureDestroy - Destroys a PetscQuadrature object
102: Collective on PetscQuadrature
104: Input Parameter:
105: . q - The PetscQuadrature object
107: Level: beginner
109: .keywords: PetscQuadrature, quadrature, destroy
110: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
111: @*/
112: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
113: {
117: if (!*q) return(0);
119: if (--((PetscObject)(*q))->refct > 0) {
120: *q = NULL;
121: return(0);
122: }
123: PetscFree((*q)->points);
124: PetscFree((*q)->weights);
125: PetscHeaderDestroy(q);
126: return(0);
127: }
129: /*@
130: PetscQuadratureGetOrder - Return the order of the method
132: Not collective
134: Input Parameter:
135: . q - The PetscQuadrature object
137: Output Parameter:
138: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
140: Level: intermediate
142: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
143: @*/
144: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
145: {
149: *order = q->order;
150: return(0);
151: }
153: /*@
154: PetscQuadratureSetOrder - Return the order of the method
156: Not collective
158: Input Parameters:
159: + q - The PetscQuadrature object
160: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
162: Level: intermediate
164: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
165: @*/
166: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
167: {
170: q->order = order;
171: return(0);
172: }
174: /*@
175: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
177: Not collective
179: Input Parameter:
180: . q - The PetscQuadrature object
182: Output Parameter:
183: . Nc - The number of components
185: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
187: Level: intermediate
189: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
190: @*/
191: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
192: {
196: *Nc = q->Nc;
197: return(0);
198: }
200: /*@
201: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
203: Not collective
205: Input Parameters:
206: + q - The PetscQuadrature object
207: - Nc - The number of components
209: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
211: Level: intermediate
213: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
214: @*/
215: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
216: {
219: q->Nc = Nc;
220: return(0);
221: }
223: /*@C
224: PetscQuadratureGetData - Returns the data defining the quadrature
226: Not collective
228: Input Parameter:
229: . q - The PetscQuadrature object
231: Output Parameters:
232: + dim - The spatial dimension
233: , Nc - The number of components
234: . npoints - The number of quadrature points
235: . points - The coordinates of each quadrature point
236: - weights - The weight of each quadrature point
238: Level: intermediate
240: .keywords: PetscQuadrature, quadrature
241: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
242: @*/
243: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
244: {
247: if (dim) {
249: *dim = q->dim;
250: }
251: if (Nc) {
253: *Nc = q->Nc;
254: }
255: if (npoints) {
257: *npoints = q->numPoints;
258: }
259: if (points) {
261: *points = q->points;
262: }
263: if (weights) {
265: *weights = q->weights;
266: }
267: return(0);
268: }
270: /*@C
271: PetscQuadratureSetData - Sets the data defining the quadrature
273: Not collective
275: Input Parameters:
276: + q - The PetscQuadrature object
277: . dim - The spatial dimension
278: , Nc - The number of components
279: . npoints - The number of quadrature points
280: . points - The coordinates of each quadrature point
281: - weights - The weight of each quadrature point
283: Level: intermediate
285: .keywords: PetscQuadrature, quadrature
286: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
287: @*/
288: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
289: {
292: if (dim >= 0) q->dim = dim;
293: if (Nc >= 0) q->Nc = Nc;
294: if (npoints >= 0) q->numPoints = npoints;
295: if (points) {
297: q->points = points;
298: }
299: if (weights) {
301: q->weights = weights;
302: }
303: return(0);
304: }
306: /*@C
307: PetscQuadratureView - Views a PetscQuadrature object
309: Collective on PetscQuadrature
311: Input Parameters:
312: + q - The PetscQuadrature object
313: - viewer - The PetscViewer object
315: Level: beginner
317: .keywords: PetscQuadrature, quadrature, view
318: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
319: @*/
320: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
321: {
322: PetscInt q, d, c;
326: PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
327: if (quad->Nc > 1) {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points with %D components\n (", quad->numPoints, quad->Nc);}
328: else {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points\n (", quad->numPoints);}
329: for (q = 0; q < quad->numPoints; ++q) {
330: for (d = 0; d < quad->dim; ++d) {
331: if (d) PetscViewerASCIIPrintf(viewer, ", ");
332: PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
333: }
334: if (quad->Nc > 1) {
335: PetscViewerASCIIPrintf(viewer, ") (");
336: for (c = 0; c < quad->Nc; ++c) {
337: if (c) PetscViewerASCIIPrintf(viewer, ", ");
338: PetscViewerASCIIPrintf(viewer, "%g", (double)quad->weights[q*quad->Nc+c]);
339: }
340: PetscViewerASCIIPrintf(viewer, ")\n");
341: } else {
342: PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
343: }
344: }
345: return(0);
346: }
348: /*@C
349: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
351: Not collective
353: Input Parameter:
354: + q - The original PetscQuadrature
355: . numSubelements - The number of subelements the original element is divided into
356: . v0 - An array of the initial points for each subelement
357: - jac - An array of the Jacobian mappings from the reference to each subelement
359: Output Parameters:
360: . dim - The dimension
362: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
364: Level: intermediate
366: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
367: @*/
368: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
369: {
370: const PetscReal *points, *weights;
371: PetscReal *pointsRef, *weightsRef;
372: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
373: PetscErrorCode ierr;
380: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
381: PetscQuadratureGetOrder(q, &order);
382: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
383: npointsRef = npoints*numSubelements;
384: PetscMalloc1(npointsRef*dim,&pointsRef);
385: PetscMalloc1(npointsRef*Nc, &weightsRef);
386: for (c = 0; c < numSubelements; ++c) {
387: for (p = 0; p < npoints; ++p) {
388: for (d = 0; d < dim; ++d) {
389: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
390: for (e = 0; e < dim; ++e) {
391: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
392: }
393: }
394: /* Could also use detJ here */
395: for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
396: }
397: }
398: PetscQuadratureSetOrder(*qref, order);
399: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
400: return(0);
401: }
403: /*@
404: PetscDTLegendreEval - evaluate Legendre polynomial at points
406: Not Collective
408: Input Arguments:
409: + npoints - number of spatial points to evaluate at
410: . points - array of locations to evaluate at
411: . ndegree - number of basis degrees to evaluate
412: - degrees - sorted array of degrees to evaluate
414: Output Arguments:
415: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
416: . D - row-oriented derivative evaluation matrix (or NULL)
417: - D2 - row-oriented second derivative evaluation matrix (or NULL)
419: Level: intermediate
421: .seealso: PetscDTGaussQuadrature()
422: @*/
423: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
424: {
425: PetscInt i,maxdegree;
428: if (!npoints || !ndegree) return(0);
429: maxdegree = degrees[ndegree-1];
430: for (i=0; i<npoints; i++) {
431: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
432: PetscInt j,k;
433: x = points[i];
434: pm2 = 0;
435: pm1 = 1;
436: pd2 = 0;
437: pd1 = 0;
438: pdd2 = 0;
439: pdd1 = 0;
440: k = 0;
441: if (degrees[k] == 0) {
442: if (B) B[i*ndegree+k] = pm1;
443: if (D) D[i*ndegree+k] = pd1;
444: if (D2) D2[i*ndegree+k] = pdd1;
445: k++;
446: }
447: for (j=1; j<=maxdegree; j++,k++) {
448: PetscReal p,d,dd;
449: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
450: d = pd2 + (2*j-1)*pm1;
451: dd = pdd2 + (2*j-1)*pd1;
452: pm2 = pm1;
453: pm1 = p;
454: pd2 = pd1;
455: pd1 = d;
456: pdd2 = pdd1;
457: pdd1 = dd;
458: if (degrees[k] == j) {
459: if (B) B[i*ndegree+k] = p;
460: if (D) D[i*ndegree+k] = d;
461: if (D2) D2[i*ndegree+k] = dd;
462: }
463: }
464: }
465: return(0);
466: }
468: /*@
469: PetscDTGaussQuadrature - create Gauss quadrature
471: Not Collective
473: Input Arguments:
474: + npoints - number of points
475: . a - left end of interval (often-1)
476: - b - right end of interval (often +1)
478: Output Arguments:
479: + x - quadrature points
480: - w - quadrature weights
482: Level: intermediate
484: References:
485: . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
487: .seealso: PetscDTLegendreEval()
488: @*/
489: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
490: {
492: PetscInt i;
493: PetscReal *work;
494: PetscScalar *Z;
495: PetscBLASInt N,LDZ,info;
498: PetscCitationsRegister(GaussCitation, &GaussCite);
499: /* Set up the Golub-Welsch system */
500: for (i=0; i<npoints; i++) {
501: x[i] = 0; /* diagonal is 0 */
502: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
503: }
504: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
505: PetscBLASIntCast(npoints,&N);
506: LDZ = N;
507: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
508: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
509: PetscFPTrapPop();
510: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
512: for (i=0; i<(npoints+1)/2; i++) {
513: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
514: x[i] = (a+b)/2 - y*(b-a)/2;
515: if (x[i] == -0.0) x[i] = 0.0;
516: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
518: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
519: }
520: PetscFree2(Z,work);
521: return(0);
522: }
524: /*@
525: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
527: Not Collective
529: Input Arguments:
530: + dim - The spatial dimension
531: . Nc - The number of components
532: . npoints - number of points in one dimension
533: . a - left end of interval (often-1)
534: - b - right end of interval (often +1)
536: Output Argument:
537: . q - A PetscQuadrature object
539: Level: intermediate
541: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
542: @*/
543: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
544: {
545: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
546: PetscReal *x, *w, *xw, *ww;
550: PetscMalloc1(totpoints*dim,&x);
551: PetscMalloc1(totpoints*Nc,&w);
552: /* Set up the Golub-Welsch system */
553: switch (dim) {
554: case 0:
555: PetscFree(x);
556: PetscFree(w);
557: PetscMalloc1(1, &x);
558: PetscMalloc1(Nc, &w);
559: x[0] = 0.0;
560: for (c = 0; c < Nc; ++c) w[c] = 1.0;
561: break;
562: case 1:
563: PetscMalloc1(npoints,&ww);
564: PetscDTGaussQuadrature(npoints, a, b, x, ww);
565: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
566: PetscFree(ww);
567: break;
568: case 2:
569: PetscMalloc2(npoints,&xw,npoints,&ww);
570: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
571: for (i = 0; i < npoints; ++i) {
572: for (j = 0; j < npoints; ++j) {
573: x[(i*npoints+j)*dim+0] = xw[i];
574: x[(i*npoints+j)*dim+1] = xw[j];
575: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
576: }
577: }
578: PetscFree2(xw,ww);
579: break;
580: case 3:
581: PetscMalloc2(npoints,&xw,npoints,&ww);
582: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
583: for (i = 0; i < npoints; ++i) {
584: for (j = 0; j < npoints; ++j) {
585: for (k = 0; k < npoints; ++k) {
586: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
587: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
588: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
589: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
590: }
591: }
592: }
593: PetscFree2(xw,ww);
594: break;
595: default:
596: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
597: }
598: PetscQuadratureCreate(PETSC_COMM_SELF, q);
599: PetscQuadratureSetOrder(*q, npoints-1);
600: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
601: return(0);
602: }
604: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
605: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
606: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
607: {
608: PetscReal f = 1.0;
609: PetscInt i;
612: for (i = 1; i < n+1; ++i) f *= i;
613: *factorial = f;
614: return(0);
615: }
617: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
618: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
619: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
620: {
621: PetscReal apb, pn1, pn2;
622: PetscInt k;
625: if (!n) {*P = 1.0; return(0);}
626: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
627: apb = a + b;
628: pn2 = 1.0;
629: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
630: *P = 0.0;
631: for (k = 2; k < n+1; ++k) {
632: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
633: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
634: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
635: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
637: a2 = a2 / a1;
638: a3 = a3 / a1;
639: a4 = a4 / a1;
640: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
641: pn2 = pn1;
642: pn1 = *P;
643: }
644: return(0);
645: }
647: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
648: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
649: {
650: PetscReal nP;
654: if (!n) {*P = 0.0; return(0);}
655: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
656: *P = 0.5 * (a + b + n + 1) * nP;
657: return(0);
658: }
660: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
661: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
662: {
664: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
665: *eta = y;
666: return(0);
667: }
669: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
670: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
671: {
673: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
674: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
675: *zeta = z;
676: return(0);
677: }
679: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
680: {
681: PetscInt maxIter = 100;
682: PetscReal eps = 1.0e-8;
683: PetscReal a1, a2, a3, a4, a5, a6;
684: PetscInt k;
689: a1 = PetscPowReal(2.0, a+b+1);
690: #if defined(PETSC_HAVE_TGAMMA)
691: a2 = PetscTGamma(a + npoints + 1);
692: a3 = PetscTGamma(b + npoints + 1);
693: a4 = PetscTGamma(a + b + npoints + 1);
694: #else
695: {
696: PetscInt ia, ib;
698: ia = (PetscInt) a;
699: ib = (PetscInt) b;
700: if (ia == a && ib == b && ia + npoints + 1 > 0 && ib + npoints + 1 > 0 && ia + ib + npoints + 1 > 0) { /* All gamma(x) terms are (x-1)! terms */
701: PetscDTFactorial_Internal(ia + npoints, &a2);
702: PetscDTFactorial_Internal(ib + npoints, &a3);
703: PetscDTFactorial_Internal(ia + ib + npoints, &a4);
704: } else {
705: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
706: }
707: }
708: #endif
710: PetscDTFactorial_Internal(npoints, &a5);
711: a6 = a1 * a2 * a3 / a4 / a5;
712: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
713: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
714: for (k = 0; k < npoints; ++k) {
715: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
716: PetscInt j;
718: if (k > 0) r = 0.5 * (r + x[k-1]);
719: for (j = 0; j < maxIter; ++j) {
720: PetscReal s = 0.0, delta, f, fp;
721: PetscInt i;
723: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
724: PetscDTComputeJacobi(a, b, npoints, r, &f);
725: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
726: delta = f / (fp - f * s);
727: r = r - delta;
728: if (PetscAbsReal(delta) < eps) break;
729: }
730: x[k] = r;
731: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
732: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
733: }
734: return(0);
735: }
737: /*@C
738: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
740: Not Collective
742: Input Arguments:
743: + dim - The simplex dimension
744: . Nc - The number of components
745: . npoints - The number of points in one dimension
746: . a - left end of interval (often-1)
747: - b - right end of interval (often +1)
749: Output Argument:
750: . q - A PetscQuadrature object
752: Level: intermediate
754: References:
755: . 1. - Karniadakis and Sherwin. FIAT
757: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
758: @*/
759: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
760: {
761: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
762: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
763: PetscInt i, j, k, c;
767: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
768: PetscMalloc1(totpoints*dim, &x);
769: PetscMalloc1(totpoints*Nc, &w);
770: switch (dim) {
771: case 0:
772: PetscFree(x);
773: PetscFree(w);
774: PetscMalloc1(1, &x);
775: PetscMalloc1(Nc, &w);
776: x[0] = 0.0;
777: for (c = 0; c < Nc; ++c) w[c] = 1.0;
778: break;
779: case 1:
780: PetscMalloc1(npoints,&wx);
781: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
782: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
783: PetscFree(wx);
784: break;
785: case 2:
786: PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
787: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
788: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
789: for (i = 0; i < npoints; ++i) {
790: for (j = 0; j < npoints; ++j) {
791: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
792: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
793: }
794: }
795: PetscFree4(px,wx,py,wy);
796: break;
797: case 3:
798: PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
799: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
800: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
801: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
802: for (i = 0; i < npoints; ++i) {
803: for (j = 0; j < npoints; ++j) {
804: for (k = 0; k < npoints; ++k) {
805: PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);
806: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
807: }
808: }
809: }
810: PetscFree6(px,wx,py,wy,pz,wz);
811: break;
812: default:
813: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
814: }
815: PetscQuadratureCreate(PETSC_COMM_SELF, q);
816: PetscQuadratureSetOrder(*q, npoints-1);
817: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
818: return(0);
819: }
821: /*@C
822: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
824: Not Collective
826: Input Arguments:
827: + dim - The cell dimension
828: . level - The number of points in one dimension, 2^l
829: . a - left end of interval (often-1)
830: - b - right end of interval (often +1)
832: Output Argument:
833: . q - A PetscQuadrature object
835: Level: intermediate
837: .seealso: PetscDTGaussTensorQuadrature()
838: @*/
839: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
840: {
841: const PetscInt p = 16; /* Digits of precision in the evaluation */
842: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
843: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
844: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
845: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
846: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
847: PetscReal *x, *w;
848: PetscInt K, k, npoints;
849: PetscErrorCode ierr;
852: if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
853: if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
854: /* Find K such that the weights are < 32 digits of precision */
855: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
856: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
857: }
858: PetscQuadratureCreate(PETSC_COMM_SELF, q);
859: PetscQuadratureSetOrder(*q, 2*K+1);
860: npoints = 2*K-1;
861: PetscMalloc1(npoints*dim, &x);
862: PetscMalloc1(npoints, &w);
863: /* Center term */
864: x[0] = beta;
865: w[0] = 0.5*alpha*PETSC_PI;
866: for (k = 1; k < K; ++k) {
867: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
868: xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
869: x[2*k-1] = -alpha*xk+beta;
870: w[2*k-1] = wk;
871: x[2*k+0] = alpha*xk+beta;
872: w[2*k+0] = wk;
873: }
874: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
875: return(0);
876: }
878: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
879: {
880: const PetscInt p = 16; /* Digits of precision in the evaluation */
881: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
882: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
883: PetscReal h = 1.0; /* Step size, length between x_k */
884: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
885: PetscReal osum = 0.0; /* Integral on last level */
886: PetscReal psum = 0.0; /* Integral on the level before the last level */
887: PetscReal sum; /* Integral on current level */
888: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
889: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
890: PetscReal wk; /* Quadrature weight at x_k */
891: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
892: PetscInt d; /* Digits of precision in the integral */
895: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
896: /* Center term */
897: func(beta, &lval);
898: sum = 0.5*alpha*PETSC_PI*lval;
899: /* */
900: do {
901: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
902: PetscInt k = 1;
904: ++l;
905: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
906: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
907: psum = osum;
908: osum = sum;
909: h *= 0.5;
910: sum *= 0.5;
911: do {
912: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
913: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
914: lx = -alpha*(1.0 - yk)+beta;
915: rx = alpha*(1.0 - yk)+beta;
916: func(lx, &lval);
917: func(rx, &rval);
918: lterm = alpha*wk*lval;
919: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
920: sum += lterm;
921: rterm = alpha*wk*rval;
922: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
923: sum += rterm;
924: ++k;
925: /* Only need to evaluate every other point on refined levels */
926: if (l != 1) ++k;
927: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
929: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
930: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
931: d3 = PetscLog10Real(maxTerm) - p;
932: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
933: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
934: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
935: } while (d < digits && l < 12);
936: *sol = sum;
938: return(0);
939: }
941: #ifdef PETSC_HAVE_MPFR
942: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
943: {
944: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
945: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
946: mpfr_t alpha; /* Half-width of the integration interval */
947: mpfr_t beta; /* Center of the integration interval */
948: mpfr_t h; /* Step size, length between x_k */
949: mpfr_t osum; /* Integral on last level */
950: mpfr_t psum; /* Integral on the level before the last level */
951: mpfr_t sum; /* Integral on current level */
952: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
953: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
954: mpfr_t wk; /* Quadrature weight at x_k */
955: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
956: PetscInt d; /* Digits of precision in the integral */
957: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
960: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
961: /* Create high precision storage */
962: 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);
963: /* Initialization */
964: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
965: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
966: mpfr_set_d(osum, 0.0, MPFR_RNDN);
967: mpfr_set_d(psum, 0.0, MPFR_RNDN);
968: mpfr_set_d(h, 1.0, MPFR_RNDN);
969: mpfr_const_pi(pi2, MPFR_RNDN);
970: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
971: /* Center term */
972: func(0.5*(b+a), &lval);
973: mpfr_set(sum, pi2, MPFR_RNDN);
974: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
975: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
976: /* */
977: do {
978: PetscReal d1, d2, d3, d4;
979: PetscInt k = 1;
981: ++l;
982: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
983: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
984: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
985: mpfr_set(psum, osum, MPFR_RNDN);
986: mpfr_set(osum, sum, MPFR_RNDN);
987: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
988: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
989: do {
990: mpfr_set_si(kh, k, MPFR_RNDN);
991: mpfr_mul(kh, kh, h, MPFR_RNDN);
992: /* Weight */
993: mpfr_set(wk, h, MPFR_RNDN);
994: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
995: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
996: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
997: mpfr_cosh(tmp, msinh, MPFR_RNDN);
998: mpfr_sqr(tmp, tmp, MPFR_RNDN);
999: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1000: mpfr_div(wk, wk, tmp, MPFR_RNDN);
1001: /* Abscissa */
1002: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1003: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1004: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1005: mpfr_exp(tmp, msinh, MPFR_RNDN);
1006: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1007: /* Quadrature points */
1008: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1009: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1010: mpfr_add(lx, lx, beta, MPFR_RNDU);
1011: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1012: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1013: mpfr_add(rx, rx, beta, MPFR_RNDD);
1014: /* Evaluation */
1015: func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1016: func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1017: /* Update */
1018: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1019: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1020: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1021: mpfr_abs(tmp, tmp, MPFR_RNDN);
1022: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1023: mpfr_set(curTerm, tmp, MPFR_RNDN);
1024: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1025: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1026: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1027: mpfr_abs(tmp, tmp, MPFR_RNDN);
1028: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1029: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1030: ++k;
1031: /* Only need to evaluate every other point on refined levels */
1032: if (l != 1) ++k;
1033: mpfr_log10(tmp, wk, MPFR_RNDN);
1034: mpfr_abs(tmp, tmp, MPFR_RNDN);
1035: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1036: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1037: mpfr_abs(tmp, tmp, MPFR_RNDN);
1038: mpfr_log10(tmp, tmp, MPFR_RNDN);
1039: d1 = mpfr_get_d(tmp, MPFR_RNDN);
1040: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1041: mpfr_abs(tmp, tmp, MPFR_RNDN);
1042: mpfr_log10(tmp, tmp, MPFR_RNDN);
1043: d2 = mpfr_get_d(tmp, MPFR_RNDN);
1044: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1045: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1046: mpfr_log10(tmp, curTerm, MPFR_RNDN);
1047: d4 = mpfr_get_d(tmp, MPFR_RNDN);
1048: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1049: } while (d < digits && l < 8);
1050: *sol = mpfr_get_d(sum, MPFR_RNDN);
1051: /* Cleanup */
1052: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1053: return(0);
1054: }
1055: #else
1057: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1058: {
1059: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1060: }
1061: #endif
1063: /* Overwrites A. Can only handle full-rank problems with m>=n
1064: * A in column-major format
1065: * Ainv in row-major format
1066: * tau has length m
1067: * worksize must be >= max(1,n)
1068: */
1069: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1070: {
1072: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1073: PetscScalar *A,*Ainv,*R,*Q,Alpha;
1076: #if defined(PETSC_USE_COMPLEX)
1077: {
1078: PetscInt i,j;
1079: PetscMalloc2(m*n,&A,m*n,&Ainv);
1080: for (j=0; j<n; j++) {
1081: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1082: }
1083: mstride = m;
1084: }
1085: #else
1086: A = A_in;
1087: Ainv = Ainv_out;
1088: #endif
1090: PetscBLASIntCast(m,&M);
1091: PetscBLASIntCast(n,&N);
1092: PetscBLASIntCast(mstride,&lda);
1093: PetscBLASIntCast(worksize,&ldwork);
1094: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1095: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1096: PetscFPTrapPop();
1097: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1098: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1100: /* Extract an explicit representation of Q */
1101: Q = Ainv;
1102: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1103: K = N; /* full rank */
1104: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1105: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1107: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1108: Alpha = 1.0;
1109: ldb = lda;
1110: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1111: /* Ainv is Q, overwritten with inverse */
1113: #if defined(PETSC_USE_COMPLEX)
1114: {
1115: PetscInt i;
1116: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1117: PetscFree2(A,Ainv);
1118: }
1119: #endif
1120: return(0);
1121: }
1123: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1124: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1125: {
1127: PetscReal *Bv;
1128: PetscInt i,j;
1131: PetscMalloc1((ninterval+1)*ndegree,&Bv);
1132: /* Point evaluation of L_p on all the source vertices */
1133: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1134: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1135: for (i=0; i<ninterval; i++) {
1136: for (j=0; j<ndegree; j++) {
1137: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1138: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1139: }
1140: }
1141: PetscFree(Bv);
1142: return(0);
1143: }
1145: /*@
1146: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
1148: Not Collective
1150: Input Arguments:
1151: + degree - degree of reconstruction polynomial
1152: . nsource - number of source intervals
1153: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1154: . ntarget - number of target intervals
1155: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
1157: Output Arguments:
1158: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
1160: Level: advanced
1162: .seealso: PetscDTLegendreEval()
1163: @*/
1164: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1165: {
1167: PetscInt i,j,k,*bdegrees,worksize;
1168: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1169: PetscScalar *tau,*work;
1175: 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);
1176: #if defined(PETSC_USE_DEBUG)
1177: for (i=0; i<nsource; i++) {
1178: 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]);
1179: }
1180: for (i=0; i<ntarget; i++) {
1181: 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]);
1182: }
1183: #endif
1184: xmin = PetscMin(sourcex[0],targetx[0]);
1185: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1186: center = (xmin + xmax)/2;
1187: hscale = (xmax - xmin)/2;
1188: worksize = nsource;
1189: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1190: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1191: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1192: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1193: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1194: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1195: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1196: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1197: for (i=0; i<ntarget; i++) {
1198: PetscReal rowsum = 0;
1199: for (j=0; j<nsource; j++) {
1200: PetscReal sum = 0;
1201: for (k=0; k<degree+1; k++) {
1202: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1203: }
1204: R[i*nsource+j] = sum;
1205: rowsum += sum;
1206: }
1207: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1208: }
1209: PetscFree4(bdegrees,sourcey,Bsource,work);
1210: PetscFree4(tau,Bsinv,targety,Btarget);
1211: return(0);
1212: }