Actual source code: dt.c
petsc-3.5.4 2015-05-23
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
8: #include <petscdt.h> /*I "petscdt.h" I*/
9: #include <petscblaslapack.h>
10: #include <petsc-private/petscimpl.h>
11: #include <petsc-private/dtimpl.h>
12: #include <petscviewer.h>
13: #include <petscdmplex.h>
14: #include <petscdmshell.h>
16: static PetscBool GaussCite = PETSC_FALSE;
17: const char GaussCitation[] = "@article{GolubWelsch1969,\n"
18: " author = {Golub and Welsch},\n"
19: " title = {Calculation of Quadrature Rules},\n"
20: " journal = {Math. Comp.},\n"
21: " volume = {23},\n"
22: " number = {106},\n"
23: " pages = {221--230},\n"
24: " year = {1969}\n}\n";
28: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
29: {
34: DMInitializePackage();
35: PetscHeaderCreate(*q,_p_PetscQuadrature,int,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
36: (*q)->dim = -1;
37: (*q)->numPoints = 0;
38: (*q)->points = NULL;
39: (*q)->weights = NULL;
40: return(0);
41: }
45: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
46: {
50: if (!*q) return(0);
52: if (--((PetscObject)(*q))->refct > 0) {
53: *q = NULL;
54: return(0);
55: }
56: PetscFree((*q)->points);
57: PetscFree((*q)->weights);
58: PetscHeaderDestroy(q);
59: return(0);
60: }
64: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
65: {
68: if (dim) {
70: *dim = q->dim;
71: }
72: if (npoints) {
74: *npoints = q->numPoints;
75: }
76: if (points) {
78: *points = q->points;
79: }
80: if (weights) {
82: *weights = q->weights;
83: }
84: return(0);
85: }
89: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
90: {
93: if (dim >= 0) q->dim = dim;
94: if (npoints >= 0) q->numPoints = npoints;
95: if (points) {
97: q->points = points;
98: }
99: if (weights) {
101: q->weights = weights;
102: }
103: return(0);
104: }
108: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
109: {
110: PetscInt q, d;
114: PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
115: PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n (", quad->numPoints);
116: for (q = 0; q < quad->numPoints; ++q) {
117: for (d = 0; d < quad->dim; ++d) {
118: if (d) PetscViewerASCIIPrintf(viewer, ", ");
119: PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
120: }
121: PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
122: }
123: return(0);
124: }
128: /*@
129: PetscDTLegendreEval - evaluate Legendre polynomial at points
131: Not Collective
133: Input Arguments:
134: + npoints - number of spatial points to evaluate at
135: . points - array of locations to evaluate at
136: . ndegree - number of basis degrees to evaluate
137: - degrees - sorted array of degrees to evaluate
139: Output Arguments:
140: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
141: . D - row-oriented derivative evaluation matrix (or NULL)
142: - D2 - row-oriented second derivative evaluation matrix (or NULL)
144: Level: intermediate
146: .seealso: PetscDTGaussQuadrature()
147: @*/
148: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
149: {
150: PetscInt i,maxdegree;
153: if (!npoints || !ndegree) return(0);
154: maxdegree = degrees[ndegree-1];
155: for (i=0; i<npoints; i++) {
156: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
157: PetscInt j,k;
158: x = points[i];
159: pm2 = 0;
160: pm1 = 1;
161: pd2 = 0;
162: pd1 = 0;
163: pdd2 = 0;
164: pdd1 = 0;
165: k = 0;
166: if (degrees[k] == 0) {
167: if (B) B[i*ndegree+k] = pm1;
168: if (D) D[i*ndegree+k] = pd1;
169: if (D2) D2[i*ndegree+k] = pdd1;
170: k++;
171: }
172: for (j=1; j<=maxdegree; j++,k++) {
173: PetscReal p,d,dd;
174: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
175: d = pd2 + (2*j-1)*pm1;
176: dd = pdd2 + (2*j-1)*pd1;
177: pm2 = pm1;
178: pm1 = p;
179: pd2 = pd1;
180: pd1 = d;
181: pdd2 = pdd1;
182: pdd1 = dd;
183: if (degrees[k] == j) {
184: if (B) B[i*ndegree+k] = p;
185: if (D) D[i*ndegree+k] = d;
186: if (D2) D2[i*ndegree+k] = dd;
187: }
188: }
189: }
190: return(0);
191: }
195: /*@
196: PetscDTGaussQuadrature - create Gauss quadrature
198: Not Collective
200: Input Arguments:
201: + npoints - number of points
202: . a - left end of interval (often-1)
203: - b - right end of interval (often +1)
205: Output Arguments:
206: + x - quadrature points
207: - w - quadrature weights
209: Level: intermediate
211: References:
212: Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969.
214: .seealso: PetscDTLegendreEval()
215: @*/
216: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
217: {
219: PetscInt i;
220: PetscReal *work;
221: PetscScalar *Z;
222: PetscBLASInt N,LDZ,info;
225: PetscCitationsRegister(GaussCitation, &GaussCite);
226: /* Set up the Golub-Welsch system */
227: for (i=0; i<npoints; i++) {
228: x[i] = 0; /* diagonal is 0 */
229: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
230: }
231: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
232: PetscBLASIntCast(npoints,&N);
233: LDZ = N;
234: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
235: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
236: PetscFPTrapPop();
237: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
239: for (i=0; i<(npoints+1)/2; i++) {
240: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
241: x[i] = (a+b)/2 - y*(b-a)/2;
242: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
244: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
245: }
246: PetscFree2(Z,work);
247: return(0);
248: }
252: /*@
253: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
255: Not Collective
257: Input Arguments:
258: + dim - The spatial dimension
259: . npoints - number of points in one dimension
260: . a - left end of interval (often-1)
261: - b - right end of interval (often +1)
263: Output Argument:
264: . q - A PetscQuadrature object
266: Level: intermediate
268: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
269: @*/
270: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
271: {
272: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
273: PetscReal *x, *w, *xw, *ww;
277: PetscMalloc1(totpoints*dim,&x);
278: PetscMalloc1(totpoints,&w);
279: /* Set up the Golub-Welsch system */
280: switch (dim) {
281: case 0:
282: PetscFree(x);
283: PetscFree(w);
284: PetscMalloc1(1, &x);
285: PetscMalloc1(1, &w);
286: x[0] = 0.0;
287: w[0] = 1.0;
288: break;
289: case 1:
290: PetscDTGaussQuadrature(npoints, a, b, x, w);
291: break;
292: case 2:
293: PetscMalloc2(npoints,&xw,npoints,&ww);
294: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
295: for (i = 0; i < npoints; ++i) {
296: for (j = 0; j < npoints; ++j) {
297: x[(i*npoints+j)*dim+0] = xw[i];
298: x[(i*npoints+j)*dim+1] = xw[j];
299: w[i*npoints+j] = ww[i] * ww[j];
300: }
301: }
302: PetscFree2(xw,ww);
303: break;
304: case 3:
305: PetscMalloc2(npoints,&xw,npoints,&ww);
306: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
307: for (i = 0; i < npoints; ++i) {
308: for (j = 0; j < npoints; ++j) {
309: for (k = 0; k < npoints; ++k) {
310: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
311: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
312: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
313: w[(i*npoints+j)*npoints+k] = ww[i] * ww[j] * ww[k];
314: }
315: }
316: }
317: PetscFree2(xw,ww);
318: break;
319: default:
320: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
321: }
322: PetscQuadratureCreate(PETSC_COMM_SELF, q);
323: PetscQuadratureSetData(*q, dim, totpoints, x, w);
324: return(0);
325: }
329: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
330: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
331: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
332: {
333: PetscReal f = 1.0;
334: PetscInt i;
337: for (i = 1; i < n+1; ++i) f *= i;
338: *factorial = f;
339: return(0);
340: }
344: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
345: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
346: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
347: {
348: PetscReal apb, pn1, pn2;
349: PetscInt k;
352: if (!n) {*P = 1.0; return(0);}
353: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
354: apb = a + b;
355: pn2 = 1.0;
356: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
357: *P = 0.0;
358: for (k = 2; k < n+1; ++k) {
359: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
360: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
361: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
362: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
364: a2 = a2 / a1;
365: a3 = a3 / a1;
366: a4 = a4 / a1;
367: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
368: pn2 = pn1;
369: pn1 = *P;
370: }
371: return(0);
372: }
376: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
377: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
378: {
379: PetscReal nP;
383: if (!n) {*P = 0.0; return(0);}
384: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
385: *P = 0.5 * (a + b + n + 1) * nP;
386: return(0);
387: }
391: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
392: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
393: {
395: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
396: *eta = y;
397: return(0);
398: }
402: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
403: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
404: {
406: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
407: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
408: *zeta = z;
409: return(0);
410: }
414: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
415: {
416: PetscInt maxIter = 100;
417: PetscReal eps = 1.0e-8;
418: PetscReal a1, a2, a3, a4, a5, a6;
419: PetscInt k;
424: a1 = PetscPowReal(2.0, a+b+1);
425: #if defined(PETSC_HAVE_TGAMMA)
426: a2 = PetscTGamma(a + npoints + 1);
427: a3 = PetscTGamma(b + npoints + 1);
428: a4 = PetscTGamma(a + b + npoints + 1);
429: #else
430: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
431: #endif
433: PetscDTFactorial_Internal(npoints, &a5);
434: a6 = a1 * a2 * a3 / a4 / a5;
435: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
436: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
437: for (k = 0; k < npoints; ++k) {
438: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
439: PetscInt j;
441: if (k > 0) r = 0.5 * (r + x[k-1]);
442: for (j = 0; j < maxIter; ++j) {
443: PetscReal s = 0.0, delta, f, fp;
444: PetscInt i;
446: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
447: PetscDTComputeJacobi(a, b, npoints, r, &f);
448: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
449: delta = f / (fp - f * s);
450: r = r - delta;
451: if (PetscAbsReal(delta) < eps) break;
452: }
453: x[k] = r;
454: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
455: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
456: }
457: return(0);
458: }
462: /*@C
463: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
465: Not Collective
467: Input Arguments:
468: + dim - The simplex dimension
469: . order - The number of points in one dimension
470: . a - left end of interval (often-1)
471: - b - right end of interval (often +1)
473: Output Argument:
474: . q - A PetscQuadrature object
476: Level: intermediate
478: References:
479: Karniadakis and Sherwin.
480: FIAT
482: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
483: @*/
484: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
485: {
486: PetscInt npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
487: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
488: PetscInt i, j, k;
492: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
493: PetscMalloc1(npoints*dim, &x);
494: PetscMalloc1(npoints, &w);
495: switch (dim) {
496: case 0:
497: PetscFree(x);
498: PetscFree(w);
499: PetscMalloc1(1, &x);
500: PetscMalloc1(1, &w);
501: x[0] = 0.0;
502: w[0] = 1.0;
503: break;
504: case 1:
505: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
506: break;
507: case 2:
508: PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
509: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
510: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
511: for (i = 0; i < order; ++i) {
512: for (j = 0; j < order; ++j) {
513: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
514: w[i*order+j] = 0.5 * wx[i] * wy[j];
515: }
516: }
517: PetscFree4(px,wx,py,wy);
518: break;
519: case 3:
520: PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
521: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
522: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
523: PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
524: for (i = 0; i < order; ++i) {
525: for (j = 0; j < order; ++j) {
526: for (k = 0; k < order; ++k) {
527: 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]);
528: w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
529: }
530: }
531: }
532: PetscFree6(px,wx,py,wy,pz,wz);
533: break;
534: default:
535: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
536: }
537: PetscQuadratureCreate(PETSC_COMM_SELF, q);
538: PetscQuadratureSetData(*q, dim, npoints, x, w);
539: return(0);
540: }
544: /* Overwrites A. Can only handle full-rank problems with m>=n
545: * A in column-major format
546: * Ainv in row-major format
547: * tau has length m
548: * worksize must be >= max(1,n)
549: */
550: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
551: {
553: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
554: PetscScalar *A,*Ainv,*R,*Q,Alpha;
557: #if defined(PETSC_USE_COMPLEX)
558: {
559: PetscInt i,j;
560: PetscMalloc2(m*n,&A,m*n,&Ainv);
561: for (j=0; j<n; j++) {
562: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
563: }
564: mstride = m;
565: }
566: #else
567: A = A_in;
568: Ainv = Ainv_out;
569: #endif
571: PetscBLASIntCast(m,&M);
572: PetscBLASIntCast(n,&N);
573: PetscBLASIntCast(mstride,&lda);
574: PetscBLASIntCast(worksize,&ldwork);
575: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
576: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
577: PetscFPTrapPop();
578: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
579: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
581: /* Extract an explicit representation of Q */
582: Q = Ainv;
583: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
584: K = N; /* full rank */
585: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
586: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
588: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
589: Alpha = 1.0;
590: ldb = lda;
591: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
592: /* Ainv is Q, overwritten with inverse */
594: #if defined(PETSC_USE_COMPLEX)
595: {
596: PetscInt i;
597: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
598: PetscFree2(A,Ainv);
599: }
600: #endif
601: return(0);
602: }
606: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
607: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
608: {
610: PetscReal *Bv;
611: PetscInt i,j;
614: PetscMalloc1((ninterval+1)*ndegree,&Bv);
615: /* Point evaluation of L_p on all the source vertices */
616: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
617: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
618: for (i=0; i<ninterval; i++) {
619: for (j=0; j<ndegree; j++) {
620: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
621: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
622: }
623: }
624: PetscFree(Bv);
625: return(0);
626: }
630: /*@
631: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
633: Not Collective
635: Input Arguments:
636: + degree - degree of reconstruction polynomial
637: . nsource - number of source intervals
638: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
639: . ntarget - number of target intervals
640: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
642: Output Arguments:
643: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
645: Level: advanced
647: .seealso: PetscDTLegendreEval()
648: @*/
649: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
650: {
652: PetscInt i,j,k,*bdegrees,worksize;
653: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
654: PetscScalar *tau,*work;
660: 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);
661: #if defined(PETSC_USE_DEBUG)
662: for (i=0; i<nsource; i++) {
663: 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]);
664: }
665: for (i=0; i<ntarget; i++) {
666: 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]);
667: }
668: #endif
669: xmin = PetscMin(sourcex[0],targetx[0]);
670: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
671: center = (xmin + xmax)/2;
672: hscale = (xmax - xmin)/2;
673: worksize = nsource;
674: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
675: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
676: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
677: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
678: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
679: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
680: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
681: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
682: for (i=0; i<ntarget; i++) {
683: PetscReal rowsum = 0;
684: for (j=0; j<nsource; j++) {
685: PetscReal sum = 0;
686: for (k=0; k<degree+1; k++) {
687: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
688: }
689: R[i*nsource+j] = sum;
690: rowsum += sum;
691: }
692: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
693: }
694: PetscFree4(bdegrees,sourcey,Bsource,work);
695: PetscFree4(tau,Bsinv,targety,Btarget);
696: return(0);
697: }