Actual source code: dt.c
petsc-3.6.1 2015-08-06
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: /*@
29: PetscQuadratureCreate - Create a PetscQuadrature object
31: Collective on MPI_Comm
33: Input Parameter:
34: . comm - The communicator for the PetscQuadrature object
36: Output Parameter:
37: . q - The PetscQuadrature object
39: Level: beginner
41: .keywords: PetscQuadrature, quadrature, create
42: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
43: @*/
44: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
45: {
50: DMInitializePackage();
51: PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
52: (*q)->dim = -1;
53: (*q)->order = -1;
54: (*q)->numPoints = 0;
55: (*q)->points = NULL;
56: (*q)->weights = NULL;
57: return(0);
58: }
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, 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, &Nq, &points, &weights);
91: PetscMalloc1(Nq*dim, &p);
92: PetscMalloc1(Nq, &w);
93: PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
94: PetscMemcpy(w, weights, Nq * sizeof(PetscReal));
95: PetscQuadratureSetData(*r, dim, Nq, p, w);
96: return(0);
97: }
101: /*@
102: PetscQuadratureDestroy - Destroys a PetscQuadrature object
104: Collective on PetscQuadrature
106: Input Parameter:
107: . q - The PetscQuadrature object
109: Level: beginner
111: .keywords: PetscQuadrature, quadrature, destroy
112: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
113: @*/
114: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
115: {
119: if (!*q) return(0);
121: if (--((PetscObject)(*q))->refct > 0) {
122: *q = NULL;
123: return(0);
124: }
125: PetscFree((*q)->points);
126: PetscFree((*q)->weights);
127: PetscHeaderDestroy(q);
128: return(0);
129: }
133: /*@
134: PetscQuadratureGetOrder - Return the quadrature information
136: Not collective
138: Input Parameter:
139: . q - The PetscQuadrature object
141: Output Parameter:
142: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
144: Output Parameter:
146: Level: intermediate
148: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
149: @*/
150: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
151: {
155: *order = q->order;
156: return(0);
157: }
161: /*@
162: PetscQuadratureSetOrder - Return the quadrature information
164: Not collective
166: Input Parameters:
167: + q - The PetscQuadrature object
168: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
170: Level: intermediate
172: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
173: @*/
174: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
175: {
178: q->order = order;
179: return(0);
180: }
184: /*@C
185: PetscQuadratureGetData - Returns the data defining the quadrature
187: Not collective
189: Input Parameter:
190: . q - The PetscQuadrature object
192: Output Parameters:
193: + dim - The spatial dimension
194: . npoints - The number of quadrature points
195: . points - The coordinates of each quadrature point
196: - weights - The weight of each quadrature point
198: Level: intermediate
200: .keywords: PetscQuadrature, quadrature
201: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
202: @*/
203: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
204: {
207: if (dim) {
209: *dim = q->dim;
210: }
211: if (npoints) {
213: *npoints = q->numPoints;
214: }
215: if (points) {
217: *points = q->points;
218: }
219: if (weights) {
221: *weights = q->weights;
222: }
223: return(0);
224: }
228: /*@C
229: PetscQuadratureSetData - Sets the data defining the quadrature
231: Not collective
233: Input Parameters:
234: + q - The PetscQuadrature object
235: . dim - The spatial dimension
236: . npoints - The number of quadrature points
237: . points - The coordinates of each quadrature point
238: - weights - The weight of each quadrature point
240: Level: intermediate
242: .keywords: PetscQuadrature, quadrature
243: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
244: @*/
245: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
246: {
249: if (dim >= 0) q->dim = dim;
250: if (npoints >= 0) q->numPoints = npoints;
251: if (points) {
253: q->points = points;
254: }
255: if (weights) {
257: q->weights = weights;
258: }
259: return(0);
260: }
264: /*@C
265: PetscQuadratureView - Views a PetscQuadrature object
267: Collective on PetscQuadrature
269: Input Parameters:
270: + q - The PetscQuadrature object
271: - viewer - The PetscViewer object
273: Level: beginner
275: .keywords: PetscQuadrature, quadrature, view
276: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
277: @*/
278: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
279: {
280: PetscInt q, d;
284: PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
285: PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n (", quad->numPoints);
286: for (q = 0; q < quad->numPoints; ++q) {
287: for (d = 0; d < quad->dim; ++d) {
288: if (d) PetscViewerASCIIPrintf(viewer, ", ");
289: PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
290: }
291: PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
292: }
293: return(0);
294: }
298: /*@C
299: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
301: Not collective
303: Input Parameter:
304: + q - The original PetscQuadrature
305: . numSubelements - The number of subelements the original element is divided into
306: . v0 - An array of the initial points for each subelement
307: - jac - An array of the Jacobian mappings from the reference to each subelement
309: Output Parameters:
310: . dim - The dimension
312: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
314: Level: intermediate
316: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
317: @*/
318: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
319: {
320: const PetscReal *points, *weights;
321: PetscReal *pointsRef, *weightsRef;
322: PetscInt dim, order, npoints, npointsRef, c, p, d, e;
323: PetscErrorCode ierr;
330: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
331: PetscQuadratureGetOrder(q, &order);
332: PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
333: npointsRef = npoints*numSubelements;
334: PetscMalloc1(npointsRef*dim,&pointsRef);
335: PetscMalloc1(npointsRef,&weightsRef);
336: for (c = 0; c < numSubelements; ++c) {
337: for (p = 0; p < npoints; ++p) {
338: for (d = 0; d < dim; ++d) {
339: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
340: for (e = 0; e < dim; ++e) {
341: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
342: }
343: }
344: /* Could also use detJ here */
345: weightsRef[c*npoints+p] = weights[p]/numSubelements;
346: }
347: }
348: PetscQuadratureSetOrder(*qref, order);
349: PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
350: return(0);
351: }
355: /*@
356: PetscDTLegendreEval - evaluate Legendre polynomial at points
358: Not Collective
360: Input Arguments:
361: + npoints - number of spatial points to evaluate at
362: . points - array of locations to evaluate at
363: . ndegree - number of basis degrees to evaluate
364: - degrees - sorted array of degrees to evaluate
366: Output Arguments:
367: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
368: . D - row-oriented derivative evaluation matrix (or NULL)
369: - D2 - row-oriented second derivative evaluation matrix (or NULL)
371: Level: intermediate
373: .seealso: PetscDTGaussQuadrature()
374: @*/
375: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
376: {
377: PetscInt i,maxdegree;
380: if (!npoints || !ndegree) return(0);
381: maxdegree = degrees[ndegree-1];
382: for (i=0; i<npoints; i++) {
383: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
384: PetscInt j,k;
385: x = points[i];
386: pm2 = 0;
387: pm1 = 1;
388: pd2 = 0;
389: pd1 = 0;
390: pdd2 = 0;
391: pdd1 = 0;
392: k = 0;
393: if (degrees[k] == 0) {
394: if (B) B[i*ndegree+k] = pm1;
395: if (D) D[i*ndegree+k] = pd1;
396: if (D2) D2[i*ndegree+k] = pdd1;
397: k++;
398: }
399: for (j=1; j<=maxdegree; j++,k++) {
400: PetscReal p,d,dd;
401: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
402: d = pd2 + (2*j-1)*pm1;
403: dd = pdd2 + (2*j-1)*pd1;
404: pm2 = pm1;
405: pm1 = p;
406: pd2 = pd1;
407: pd1 = d;
408: pdd2 = pdd1;
409: pdd1 = dd;
410: if (degrees[k] == j) {
411: if (B) B[i*ndegree+k] = p;
412: if (D) D[i*ndegree+k] = d;
413: if (D2) D2[i*ndegree+k] = dd;
414: }
415: }
416: }
417: return(0);
418: }
422: /*@
423: PetscDTGaussQuadrature - create Gauss quadrature
425: Not Collective
427: Input Arguments:
428: + npoints - number of points
429: . a - left end of interval (often-1)
430: - b - right end of interval (often +1)
432: Output Arguments:
433: + x - quadrature points
434: - w - quadrature weights
436: Level: intermediate
438: References:
439: Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969.
441: .seealso: PetscDTLegendreEval()
442: @*/
443: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
444: {
446: PetscInt i;
447: PetscReal *work;
448: PetscScalar *Z;
449: PetscBLASInt N,LDZ,info;
452: PetscCitationsRegister(GaussCitation, &GaussCite);
453: /* Set up the Golub-Welsch system */
454: for (i=0; i<npoints; i++) {
455: x[i] = 0; /* diagonal is 0 */
456: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
457: }
458: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
459: PetscBLASIntCast(npoints,&N);
460: LDZ = N;
461: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
462: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
463: PetscFPTrapPop();
464: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
466: for (i=0; i<(npoints+1)/2; i++) {
467: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
468: x[i] = (a+b)/2 - y*(b-a)/2;
469: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
471: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
472: }
473: PetscFree2(Z,work);
474: return(0);
475: }
479: /*@
480: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
482: Not Collective
484: Input Arguments:
485: + dim - The spatial dimension
486: . npoints - number of points in one dimension
487: . a - left end of interval (often-1)
488: - b - right end of interval (often +1)
490: Output Argument:
491: . q - A PetscQuadrature object
493: Level: intermediate
495: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
496: @*/
497: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
498: {
499: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
500: PetscReal *x, *w, *xw, *ww;
504: PetscMalloc1(totpoints*dim,&x);
505: PetscMalloc1(totpoints,&w);
506: /* Set up the Golub-Welsch system */
507: switch (dim) {
508: case 0:
509: PetscFree(x);
510: PetscFree(w);
511: PetscMalloc1(1, &x);
512: PetscMalloc1(1, &w);
513: x[0] = 0.0;
514: w[0] = 1.0;
515: break;
516: case 1:
517: PetscDTGaussQuadrature(npoints, a, b, x, w);
518: break;
519: case 2:
520: PetscMalloc2(npoints,&xw,npoints,&ww);
521: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
522: for (i = 0; i < npoints; ++i) {
523: for (j = 0; j < npoints; ++j) {
524: x[(i*npoints+j)*dim+0] = xw[i];
525: x[(i*npoints+j)*dim+1] = xw[j];
526: w[i*npoints+j] = ww[i] * ww[j];
527: }
528: }
529: PetscFree2(xw,ww);
530: break;
531: case 3:
532: PetscMalloc2(npoints,&xw,npoints,&ww);
533: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
534: for (i = 0; i < npoints; ++i) {
535: for (j = 0; j < npoints; ++j) {
536: for (k = 0; k < npoints; ++k) {
537: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
538: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
539: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
540: w[(i*npoints+j)*npoints+k] = ww[i] * ww[j] * ww[k];
541: }
542: }
543: }
544: PetscFree2(xw,ww);
545: break;
546: default:
547: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
548: }
549: PetscQuadratureCreate(PETSC_COMM_SELF, q);
550: PetscQuadratureSetOrder(*q, npoints);
551: PetscQuadratureSetData(*q, dim, totpoints, x, w);
552: return(0);
553: }
557: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
558: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
559: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
560: {
561: PetscReal f = 1.0;
562: PetscInt i;
565: for (i = 1; i < n+1; ++i) f *= i;
566: *factorial = f;
567: return(0);
568: }
572: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
573: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
574: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
575: {
576: PetscReal apb, pn1, pn2;
577: PetscInt k;
580: if (!n) {*P = 1.0; return(0);}
581: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
582: apb = a + b;
583: pn2 = 1.0;
584: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
585: *P = 0.0;
586: for (k = 2; k < n+1; ++k) {
587: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
588: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
589: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
590: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
592: a2 = a2 / a1;
593: a3 = a3 / a1;
594: a4 = a4 / a1;
595: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
596: pn2 = pn1;
597: pn1 = *P;
598: }
599: return(0);
600: }
604: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
605: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
606: {
607: PetscReal nP;
611: if (!n) {*P = 0.0; return(0);}
612: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
613: *P = 0.5 * (a + b + n + 1) * nP;
614: return(0);
615: }
619: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
620: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
621: {
623: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
624: *eta = y;
625: return(0);
626: }
630: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
631: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
632: {
634: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
635: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
636: *zeta = z;
637: return(0);
638: }
642: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
643: {
644: PetscInt maxIter = 100;
645: PetscReal eps = 1.0e-8;
646: PetscReal a1, a2, a3, a4, a5, a6;
647: PetscInt k;
652: a1 = PetscPowReal(2.0, a+b+1);
653: #if defined(PETSC_HAVE_TGAMMA)
654: a2 = PetscTGamma(a + npoints + 1);
655: a3 = PetscTGamma(b + npoints + 1);
656: a4 = PetscTGamma(a + b + npoints + 1);
657: #else
658: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
659: #endif
661: PetscDTFactorial_Internal(npoints, &a5);
662: a6 = a1 * a2 * a3 / a4 / a5;
663: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
664: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
665: for (k = 0; k < npoints; ++k) {
666: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
667: PetscInt j;
669: if (k > 0) r = 0.5 * (r + x[k-1]);
670: for (j = 0; j < maxIter; ++j) {
671: PetscReal s = 0.0, delta, f, fp;
672: PetscInt i;
674: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
675: PetscDTComputeJacobi(a, b, npoints, r, &f);
676: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
677: delta = f / (fp - f * s);
678: r = r - delta;
679: if (PetscAbsReal(delta) < eps) break;
680: }
681: x[k] = r;
682: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
683: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
684: }
685: return(0);
686: }
690: /*@C
691: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
693: Not Collective
695: Input Arguments:
696: + dim - The simplex dimension
697: . order - The number of points in one dimension
698: . a - left end of interval (often-1)
699: - b - right end of interval (often +1)
701: Output Argument:
702: . q - A PetscQuadrature object
704: Level: intermediate
706: References:
707: Karniadakis and Sherwin.
708: FIAT
710: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
711: @*/
712: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
713: {
714: PetscInt npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
715: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
716: PetscInt i, j, k;
720: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
721: PetscMalloc1(npoints*dim, &x);
722: PetscMalloc1(npoints, &w);
723: switch (dim) {
724: case 0:
725: PetscFree(x);
726: PetscFree(w);
727: PetscMalloc1(1, &x);
728: PetscMalloc1(1, &w);
729: x[0] = 0.0;
730: w[0] = 1.0;
731: break;
732: case 1:
733: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
734: break;
735: case 2:
736: PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
737: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
738: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
739: for (i = 0; i < order; ++i) {
740: for (j = 0; j < order; ++j) {
741: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
742: w[i*order+j] = 0.5 * wx[i] * wy[j];
743: }
744: }
745: PetscFree4(px,wx,py,wy);
746: break;
747: case 3:
748: PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
749: PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
750: PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
751: PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
752: for (i = 0; i < order; ++i) {
753: for (j = 0; j < order; ++j) {
754: for (k = 0; k < order; ++k) {
755: 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]);
756: w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
757: }
758: }
759: }
760: PetscFree6(px,wx,py,wy,pz,wz);
761: break;
762: default:
763: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
764: }
765: PetscQuadratureCreate(PETSC_COMM_SELF, q);
766: PetscQuadratureSetOrder(*q, order);
767: PetscQuadratureSetData(*q, dim, npoints, x, w);
768: return(0);
769: }
773: /* Overwrites A. Can only handle full-rank problems with m>=n
774: * A in column-major format
775: * Ainv in row-major format
776: * tau has length m
777: * worksize must be >= max(1,n)
778: */
779: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
780: {
782: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
783: PetscScalar *A,*Ainv,*R,*Q,Alpha;
786: #if defined(PETSC_USE_COMPLEX)
787: {
788: PetscInt i,j;
789: PetscMalloc2(m*n,&A,m*n,&Ainv);
790: for (j=0; j<n; j++) {
791: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
792: }
793: mstride = m;
794: }
795: #else
796: A = A_in;
797: Ainv = Ainv_out;
798: #endif
800: PetscBLASIntCast(m,&M);
801: PetscBLASIntCast(n,&N);
802: PetscBLASIntCast(mstride,&lda);
803: PetscBLASIntCast(worksize,&ldwork);
804: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
805: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
806: PetscFPTrapPop();
807: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
808: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
810: /* Extract an explicit representation of Q */
811: Q = Ainv;
812: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
813: K = N; /* full rank */
814: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
815: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
817: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
818: Alpha = 1.0;
819: ldb = lda;
820: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
821: /* Ainv is Q, overwritten with inverse */
823: #if defined(PETSC_USE_COMPLEX)
824: {
825: PetscInt i;
826: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
827: PetscFree2(A,Ainv);
828: }
829: #endif
830: return(0);
831: }
835: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
836: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
837: {
839: PetscReal *Bv;
840: PetscInt i,j;
843: PetscMalloc1((ninterval+1)*ndegree,&Bv);
844: /* Point evaluation of L_p on all the source vertices */
845: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
846: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
847: for (i=0; i<ninterval; i++) {
848: for (j=0; j<ndegree; j++) {
849: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
850: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
851: }
852: }
853: PetscFree(Bv);
854: return(0);
855: }
859: /*@
860: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
862: Not Collective
864: Input Arguments:
865: + degree - degree of reconstruction polynomial
866: . nsource - number of source intervals
867: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
868: . ntarget - number of target intervals
869: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
871: Output Arguments:
872: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
874: Level: advanced
876: .seealso: PetscDTLegendreEval()
877: @*/
878: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
879: {
881: PetscInt i,j,k,*bdegrees,worksize;
882: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
883: PetscScalar *tau,*work;
889: 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);
890: #if defined(PETSC_USE_DEBUG)
891: for (i=0; i<nsource; i++) {
892: 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]);
893: }
894: for (i=0; i<ntarget; i++) {
895: 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]);
896: }
897: #endif
898: xmin = PetscMin(sourcex[0],targetx[0]);
899: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
900: center = (xmin + xmax)/2;
901: hscale = (xmax - xmin)/2;
902: worksize = nsource;
903: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
904: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
905: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
906: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
907: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
908: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
909: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
910: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
911: for (i=0; i<ntarget; i++) {
912: PetscReal rowsum = 0;
913: for (j=0; j<nsource; j++) {
914: PetscReal sum = 0;
915: for (k=0; k<degree+1; k++) {
916: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
917: }
918: R[i*nsource+j] = sum;
919: rowsum += sum;
920: }
921: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
922: }
923: PetscFree4(bdegrees,sourcey,Bsource,work);
924: PetscFree4(tau,Bsinv,targety,Btarget);
925: return(0);
926: }