Actual source code: dt.c
petsc-3.12.5 2020-03-29
1: /* Discretization tools */
3: #include <petscdt.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/dtimpl.h>
7: #include <petscviewer.h>
8: #include <petscdmplex.h>
9: #include <petscdmshell.h>
11: #if defined(PETSC_HAVE_MPFR)
12: #include <mpfr.h>
13: #endif
15: static PetscBool GaussCite = PETSC_FALSE;
16: const char GaussCitation[] = "@article{GolubWelsch1969,\n"
17: " author = {Golub and Welsch},\n"
18: " title = {Calculation of Quadrature Rules},\n"
19: " journal = {Math. Comp.},\n"
20: " volume = {23},\n"
21: " number = {106},\n"
22: " pages = {221--230},\n"
23: " year = {1969}\n}\n";
25: /*@
26: PetscQuadratureCreate - Create a PetscQuadrature object
28: Collective
30: Input Parameter:
31: . comm - The communicator for the PetscQuadrature object
33: Output Parameter:
34: . q - The PetscQuadrature object
36: Level: beginner
38: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
39: @*/
40: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
41: {
46: PetscSysInitializePackage();
47: PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
48: (*q)->dim = -1;
49: (*q)->Nc = 1;
50: (*q)->order = -1;
51: (*q)->numPoints = 0;
52: (*q)->points = NULL;
53: (*q)->weights = NULL;
54: return(0);
55: }
57: /*@
58: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
60: Collective on q
62: Input Parameter:
63: . q - The PetscQuadrature object
65: Output Parameter:
66: . r - The new PetscQuadrature object
68: Level: beginner
70: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
71: @*/
72: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
73: {
74: PetscInt order, dim, Nc, Nq;
75: const PetscReal *points, *weights;
76: PetscReal *p, *w;
77: PetscErrorCode ierr;
81: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
82: PetscQuadratureGetOrder(q, &order);
83: PetscQuadratureSetOrder(*r, order);
84: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
85: PetscMalloc1(Nq*dim, &p);
86: PetscMalloc1(Nq*Nc, &w);
87: PetscArraycpy(p, points, Nq*dim);
88: PetscArraycpy(w, weights, Nc * Nq);
89: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
90: return(0);
91: }
93: /*@
94: PetscQuadratureDestroy - Destroys a PetscQuadrature object
96: Collective on q
98: Input Parameter:
99: . q - The PetscQuadrature object
101: Level: beginner
103: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
104: @*/
105: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
106: {
110: if (!*q) return(0);
112: if (--((PetscObject)(*q))->refct > 0) {
113: *q = NULL;
114: return(0);
115: }
116: PetscFree((*q)->points);
117: PetscFree((*q)->weights);
118: PetscHeaderDestroy(q);
119: return(0);
120: }
122: /*@
123: PetscQuadratureGetOrder - Return the order of the method
125: Not collective
127: Input Parameter:
128: . q - The PetscQuadrature object
130: Output Parameter:
131: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
133: Level: intermediate
135: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
136: @*/
137: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
138: {
142: *order = q->order;
143: return(0);
144: }
146: /*@
147: PetscQuadratureSetOrder - Return the order of the method
149: Not collective
151: Input Parameters:
152: + q - The PetscQuadrature object
153: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
155: Level: intermediate
157: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
158: @*/
159: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
160: {
163: q->order = order;
164: return(0);
165: }
167: /*@
168: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
170: Not collective
172: Input Parameter:
173: . q - The PetscQuadrature object
175: Output Parameter:
176: . Nc - The number of components
178: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
180: Level: intermediate
182: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
183: @*/
184: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
185: {
189: *Nc = q->Nc;
190: return(0);
191: }
193: /*@
194: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
196: Not collective
198: Input Parameters:
199: + q - The PetscQuadrature object
200: - Nc - The number of components
202: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
204: Level: intermediate
206: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
207: @*/
208: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
209: {
212: q->Nc = Nc;
213: return(0);
214: }
216: /*@C
217: PetscQuadratureGetData - Returns the data defining the quadrature
219: Not collective
221: Input Parameter:
222: . q - The PetscQuadrature object
224: Output Parameters:
225: + dim - The spatial dimension
226: . Nc - The number of components
227: . npoints - The number of quadrature points
228: . points - The coordinates of each quadrature point
229: - weights - The weight of each quadrature point
231: Level: intermediate
233: Fortran Notes:
234: From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
236: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
237: @*/
238: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
239: {
242: if (dim) {
244: *dim = q->dim;
245: }
246: if (Nc) {
248: *Nc = q->Nc;
249: }
250: if (npoints) {
252: *npoints = q->numPoints;
253: }
254: if (points) {
256: *points = q->points;
257: }
258: if (weights) {
260: *weights = q->weights;
261: }
262: return(0);
263: }
265: /*@C
266: PetscQuadratureSetData - Sets the data defining the quadrature
268: Not collective
270: Input Parameters:
271: + q - The PetscQuadrature object
272: . dim - The spatial dimension
273: . Nc - The number of components
274: . npoints - The number of quadrature points
275: . points - The coordinates of each quadrature point
276: - weights - The weight of each quadrature point
278: Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
280: Level: intermediate
282: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
283: @*/
284: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
285: {
288: if (dim >= 0) q->dim = dim;
289: if (Nc >= 0) q->Nc = Nc;
290: if (npoints >= 0) q->numPoints = npoints;
291: if (points) {
293: q->points = points;
294: }
295: if (weights) {
297: q->weights = weights;
298: }
299: return(0);
300: }
302: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
303: {
304: PetscInt q, d, c;
305: PetscViewerFormat format;
306: PetscErrorCode ierr;
309: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);}
310: else {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);}
311: PetscViewerGetFormat(v, &format);
312: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
313: for (q = 0; q < quad->numPoints; ++q) {
314: PetscViewerASCIIPrintf(v, "p%D (", q);
315: PetscViewerASCIIUseTabs(v, PETSC_FALSE);
316: for (d = 0; d < quad->dim; ++d) {
317: if (d) {PetscViewerASCIIPrintf(v, ", ");}
318: PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
319: }
320: PetscViewerASCIIPrintf(v, ") ");
321: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "w%D (", q);}
322: for (c = 0; c < quad->Nc; ++c) {
323: if (c) {PetscViewerASCIIPrintf(v, ", ");}
324: PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
325: }
326: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
327: PetscViewerASCIIPrintf(v, "\n");
328: PetscViewerASCIIUseTabs(v, PETSC_TRUE);
329: }
330: return(0);
331: }
333: /*@C
334: PetscQuadratureView - Views a PetscQuadrature object
336: Collective on quad
338: Input Parameters:
339: + quad - The PetscQuadrature object
340: - viewer - The PetscViewer object
342: Level: beginner
344: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
345: @*/
346: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
347: {
348: PetscBool iascii;
354: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
355: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
356: PetscViewerASCIIPushTab(viewer);
357: if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
358: PetscViewerASCIIPopTab(viewer);
359: return(0);
360: }
362: /*@C
363: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
365: Not collective
367: Input Parameter:
368: + q - The original PetscQuadrature
369: . numSubelements - The number of subelements the original element is divided into
370: . v0 - An array of the initial points for each subelement
371: - jac - An array of the Jacobian mappings from the reference to each subelement
373: Output Parameters:
374: . dim - The dimension
376: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
378: Not available from Fortran
380: Level: intermediate
382: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
383: @*/
384: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
385: {
386: const PetscReal *points, *weights;
387: PetscReal *pointsRef, *weightsRef;
388: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
389: PetscErrorCode ierr;
396: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
397: PetscQuadratureGetOrder(q, &order);
398: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
399: npointsRef = npoints*numSubelements;
400: PetscMalloc1(npointsRef*dim,&pointsRef);
401: PetscMalloc1(npointsRef*Nc, &weightsRef);
402: for (c = 0; c < numSubelements; ++c) {
403: for (p = 0; p < npoints; ++p) {
404: for (d = 0; d < dim; ++d) {
405: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
406: for (e = 0; e < dim; ++e) {
407: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
408: }
409: }
410: /* Could also use detJ here */
411: for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
412: }
413: }
414: PetscQuadratureSetOrder(*qref, order);
415: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
416: return(0);
417: }
419: /*@
420: PetscDTLegendreEval - evaluate Legendre polynomial at points
422: Not Collective
424: Input Arguments:
425: + npoints - number of spatial points to evaluate at
426: . points - array of locations to evaluate at
427: . ndegree - number of basis degrees to evaluate
428: - degrees - sorted array of degrees to evaluate
430: Output Arguments:
431: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
432: . D - row-oriented derivative evaluation matrix (or NULL)
433: - D2 - row-oriented second derivative evaluation matrix (or NULL)
435: Level: intermediate
437: .seealso: PetscDTGaussQuadrature()
438: @*/
439: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
440: {
441: PetscInt i,maxdegree;
444: if (!npoints || !ndegree) return(0);
445: maxdegree = degrees[ndegree-1];
446: for (i=0; i<npoints; i++) {
447: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
448: PetscInt j,k;
449: x = points[i];
450: pm2 = 0;
451: pm1 = 1;
452: pd2 = 0;
453: pd1 = 0;
454: pdd2 = 0;
455: pdd1 = 0;
456: k = 0;
457: if (degrees[k] == 0) {
458: if (B) B[i*ndegree+k] = pm1;
459: if (D) D[i*ndegree+k] = pd1;
460: if (D2) D2[i*ndegree+k] = pdd1;
461: k++;
462: }
463: for (j=1; j<=maxdegree; j++,k++) {
464: PetscReal p,d,dd;
465: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
466: d = pd2 + (2*j-1)*pm1;
467: dd = pdd2 + (2*j-1)*pd1;
468: pm2 = pm1;
469: pm1 = p;
470: pd2 = pd1;
471: pd1 = d;
472: pdd2 = pdd1;
473: pdd1 = dd;
474: if (degrees[k] == j) {
475: if (B) B[i*ndegree+k] = p;
476: if (D) D[i*ndegree+k] = d;
477: if (D2) D2[i*ndegree+k] = dd;
478: }
479: }
480: }
481: return(0);
482: }
484: /*@
485: PetscDTGaussQuadrature - create Gauss quadrature
487: Not Collective
489: Input Arguments:
490: + npoints - number of points
491: . a - left end of interval (often-1)
492: - b - right end of interval (often +1)
494: Output Arguments:
495: + x - quadrature points
496: - w - quadrature weights
498: Level: intermediate
500: References:
501: . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
503: .seealso: PetscDTLegendreEval()
504: @*/
505: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
506: {
508: PetscInt i;
509: PetscReal *work;
510: PetscScalar *Z;
511: PetscBLASInt N,LDZ,info;
514: PetscCitationsRegister(GaussCitation, &GaussCite);
515: /* Set up the Golub-Welsch system */
516: for (i=0; i<npoints; i++) {
517: x[i] = 0; /* diagonal is 0 */
518: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
519: }
520: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
521: PetscBLASIntCast(npoints,&N);
522: LDZ = N;
523: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
524: #if defined(PETSC_MISSING_LAPACK_STEQR)
525: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEQR - Lapack routine is unavailable.");
526: #else
527: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
528: #endif
529: PetscFPTrapPop();
530: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
532: for (i=0; i<(npoints+1)/2; i++) {
533: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
534: x[i] = (a+b)/2 - y*(b-a)/2;
535: if (x[i] == -0.0) x[i] = 0.0;
536: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
538: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
539: }
540: PetscFree2(Z,work);
541: return(0);
542: }
544: static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln)
545: /*
546: Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in
547: addition to L_N(x) as these are needed for computing the GLL points via Newton's method.
548: Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms
549: for Scientists and Engineers" by David A. Kopriva.
550: */
551: {
552: PetscInt k;
554: PetscReal Lnp;
555: PetscReal Lnp1, Lnp1p;
556: PetscReal Lnm1, Lnm1p;
557: PetscReal Lnm2, Lnm2p;
559: Lnm1 = 1.0;
560: *Ln = x;
561: Lnm1p = 0.0;
562: Lnp = 1.0;
564: for (k=2; k<=n; ++k) {
565: Lnm2 = Lnm1;
566: Lnm1 = *Ln;
567: Lnm2p = Lnm1p;
568: Lnm1p = Lnp;
569: *Ln = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2;
570: Lnp = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1;
571: }
572: k = n+1;
573: Lnp1 = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1;
574: Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln);
575: *q = Lnp1 - Lnm1;
576: *qp = Lnp1p - Lnm1p;
577: }
579: /*@C
580: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
581: nodes of a given size on the domain [-1,1]
583: Not Collective
585: Input Parameter:
586: + n - number of grid nodes
587: - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON
589: Output Arguments:
590: + x - quadrature points
591: - w - quadrature weights
593: Notes:
594: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
595: close enough to the desired solution
597: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
599: See https://epubs.siam.org/doi/abs/10.1137/110855442 https://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes
601: Level: intermediate
603: .seealso: PetscDTGaussQuadrature()
605: @*/
606: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
607: {
611: if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
613: if (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA) {
614: PetscReal *M,si;
615: PetscBLASInt bn,lierr;
616: PetscReal x0,z0,z1,z2;
617: PetscInt i,p = npoints - 1,nn;
619: x[0] =-1.0;
620: x[npoints-1] = 1.0;
621: if (npoints-2 > 0){
622: PetscMalloc1(npoints-1,&M);
623: for (i=0; i<npoints-2; i++) {
624: si = ((PetscReal)i)+1.0;
625: M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
626: }
627: PetscBLASIntCast(npoints-2,&bn);
628: PetscArrayzero(&x[1],bn);
629: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
630: x0=0;
631: #if defined(PETSC_MISSING_LAPACK_STEQR)
632: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEQR - Lapack routine is unavailable.");
633: #else
634: PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&x[1],M,&x0,&bn,M,&lierr));
635: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
636: #endif
637: PetscFPTrapPop();
638: PetscFree(M);
639: }
640: if ((npoints-1)%2==0) {
641: x[(npoints-1)/2] = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
642: }
644: w[0] = w[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
645: z2 = -1.; /* Dummy value to avoid -Wmaybe-initialized */
646: for (i=1; i<p; i++) {
647: x0 = x[i];
648: z0 = 1.0;
649: z1 = x0;
650: for (nn=1; nn<p; nn++) {
651: z2 = x0*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
652: z0 = z1;
653: z1 = z2;
654: }
655: w[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
656: }
657: } else {
658: PetscInt j,m;
659: PetscReal z1,z,q,qp,Ln;
660: PetscReal *pt;
661: PetscMalloc1(npoints,&pt);
663: if (npoints > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON produces incorrect answers for n > 30");
664: x[0] = -1.0;
665: x[npoints-1] = 1.0;
666: w[0] = w[npoints-1] = 2./(((PetscReal)npoints)*(((PetscReal)npoints)-1.0));
667: m = (npoints-1)/2; /* The roots are symmetric, so we only find half of them. */
668: for (j=1; j<=m; j++) { /* Loop over the desired roots. */
669: z = -1.0*PetscCosReal((PETSC_PI*((PetscReal)j)+0.25)/(((PetscReal)npoints)-1.0))-(3.0/(8.0*(((PetscReal)npoints)-1.0)*PETSC_PI))*(1.0/(((PetscReal)j)+0.25));
670: /* Starting with the above approximation to the ith root, we enter */
671: /* the main loop of refinement by Newton's method. */
672: do {
673: qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);
674: z1 = z;
675: z = z1-q/qp; /* Newton's method. */
676: } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
677: qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);
679: x[j] = z;
680: x[npoints-1-j] = -z; /* and put in its symmetric counterpart. */
681: w[j] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln); /* Compute the weight */
682: w[npoints-1-j] = w[j]; /* and its symmetric counterpart. */
683: pt[j]=qp;
684: }
686: if ((npoints-1)%2==0) {
687: qAndLEvaluation(npoints-1,0.0,&q,&qp,&Ln);
688: x[(npoints-1)/2] = 0.0;
689: w[(npoints-1)/2] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln);
690: }
691: PetscFree(pt);
692: }
693: return(0);
694: }
696: /*@
697: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
699: Not Collective
701: Input Arguments:
702: + dim - The spatial dimension
703: . Nc - The number of components
704: . npoints - number of points in one dimension
705: . a - left end of interval (often-1)
706: - b - right end of interval (often +1)
708: Output Argument:
709: . q - A PetscQuadrature object
711: Level: intermediate
713: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
714: @*/
715: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
716: {
717: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
718: PetscReal *x, *w, *xw, *ww;
722: PetscMalloc1(totpoints*dim,&x);
723: PetscMalloc1(totpoints*Nc,&w);
724: /* Set up the Golub-Welsch system */
725: switch (dim) {
726: case 0:
727: PetscFree(x);
728: PetscFree(w);
729: PetscMalloc1(1, &x);
730: PetscMalloc1(Nc, &w);
731: x[0] = 0.0;
732: for (c = 0; c < Nc; ++c) w[c] = 1.0;
733: break;
734: case 1:
735: PetscMalloc1(npoints,&ww);
736: PetscDTGaussQuadrature(npoints, a, b, x, ww);
737: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
738: PetscFree(ww);
739: break;
740: case 2:
741: PetscMalloc2(npoints,&xw,npoints,&ww);
742: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
743: for (i = 0; i < npoints; ++i) {
744: for (j = 0; j < npoints; ++j) {
745: x[(i*npoints+j)*dim+0] = xw[i];
746: x[(i*npoints+j)*dim+1] = xw[j];
747: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
748: }
749: }
750: PetscFree2(xw,ww);
751: break;
752: case 3:
753: PetscMalloc2(npoints,&xw,npoints,&ww);
754: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
755: for (i = 0; i < npoints; ++i) {
756: for (j = 0; j < npoints; ++j) {
757: for (k = 0; k < npoints; ++k) {
758: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
759: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
760: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
761: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
762: }
763: }
764: }
765: PetscFree2(xw,ww);
766: break;
767: default:
768: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
769: }
770: PetscQuadratureCreate(PETSC_COMM_SELF, q);
771: PetscQuadratureSetOrder(*q, 2*npoints-1);
772: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
773: PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
774: return(0);
775: }
777: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
778: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
779: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
780: {
781: PetscReal f = 1.0;
782: PetscInt i;
785: for (i = 1; i < n+1; ++i) f *= i;
786: *factorial = f;
787: return(0);
788: }
790: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
791: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
792: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
793: {
794: PetscReal apb, pn1, pn2;
795: PetscInt k;
798: if (!n) {*P = 1.0; return(0);}
799: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
800: apb = a + b;
801: pn2 = 1.0;
802: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
803: *P = 0.0;
804: for (k = 2; k < n+1; ++k) {
805: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
806: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
807: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
808: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
810: a2 = a2 / a1;
811: a3 = a3 / a1;
812: a4 = a4 / a1;
813: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
814: pn2 = pn1;
815: pn1 = *P;
816: }
817: return(0);
818: }
820: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
821: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
822: {
823: PetscReal nP;
827: if (!n) {*P = 0.0; return(0);}
828: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
829: *P = 0.5 * (a + b + n + 1) * nP;
830: return(0);
831: }
833: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
834: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
835: {
837: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
838: *eta = y;
839: return(0);
840: }
842: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
843: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
844: {
846: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
847: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
848: *zeta = z;
849: return(0);
850: }
852: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
853: {
854: PetscInt maxIter = 100;
855: PetscReal eps = 1.0e-8;
856: PetscReal a1, a2, a3, a4, a5, a6;
857: PetscInt k;
862: a1 = PetscPowReal(2.0, a+b+1);
863: #if defined(PETSC_HAVE_TGAMMA)
864: a2 = PetscTGamma(a + npoints + 1);
865: a3 = PetscTGamma(b + npoints + 1);
866: a4 = PetscTGamma(a + b + npoints + 1);
867: #else
868: {
869: PetscInt ia, ib;
871: ia = (PetscInt) a;
872: ib = (PetscInt) b;
873: 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 */
874: PetscDTFactorial_Internal(ia + npoints, &a2);
875: PetscDTFactorial_Internal(ib + npoints, &a3);
876: PetscDTFactorial_Internal(ia + ib + npoints, &a4);
877: } else {
878: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
879: }
880: }
881: #endif
883: PetscDTFactorial_Internal(npoints, &a5);
884: a6 = a1 * a2 * a3 / a4 / a5;
885: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
886: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
887: for (k = 0; k < npoints; ++k) {
888: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
889: PetscInt j;
891: if (k > 0) r = 0.5 * (r + x[k-1]);
892: for (j = 0; j < maxIter; ++j) {
893: PetscReal s = 0.0, delta, f, fp;
894: PetscInt i;
896: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
897: PetscDTComputeJacobi(a, b, npoints, r, &f);
898: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
899: delta = f / (fp - f * s);
900: r = r - delta;
901: if (PetscAbsReal(delta) < eps) break;
902: }
903: x[k] = r;
904: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
905: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
906: }
907: return(0);
908: }
910: /*@
911: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
913: Not Collective
915: Input Arguments:
916: + dim - The simplex dimension
917: . Nc - The number of components
918: . npoints - The number of points in one dimension
919: . a - left end of interval (often-1)
920: - b - right end of interval (often +1)
922: Output Argument:
923: . q - A PetscQuadrature object
925: Level: intermediate
927: References:
928: . 1. - Karniadakis and Sherwin. FIAT
930: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
931: @*/
932: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
933: {
934: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
935: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
936: PetscInt i, j, k, c;
940: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
941: PetscMalloc1(totpoints*dim, &x);
942: PetscMalloc1(totpoints*Nc, &w);
943: switch (dim) {
944: case 0:
945: PetscFree(x);
946: PetscFree(w);
947: PetscMalloc1(1, &x);
948: PetscMalloc1(Nc, &w);
949: x[0] = 0.0;
950: for (c = 0; c < Nc; ++c) w[c] = 1.0;
951: break;
952: case 1:
953: PetscMalloc1(npoints,&wx);
954: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
955: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
956: PetscFree(wx);
957: break;
958: case 2:
959: PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
960: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
961: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
962: for (i = 0; i < npoints; ++i) {
963: for (j = 0; j < npoints; ++j) {
964: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
965: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
966: }
967: }
968: PetscFree4(px,wx,py,wy);
969: break;
970: case 3:
971: PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
972: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
973: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
974: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
975: for (i = 0; i < npoints; ++i) {
976: for (j = 0; j < npoints; ++j) {
977: for (k = 0; k < npoints; ++k) {
978: 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]);
979: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
980: }
981: }
982: }
983: PetscFree6(px,wx,py,wy,pz,wz);
984: break;
985: default:
986: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
987: }
988: PetscQuadratureCreate(PETSC_COMM_SELF, q);
989: PetscQuadratureSetOrder(*q, 2*npoints-1);
990: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
991: PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
992: return(0);
993: }
995: /*@
996: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
998: Not Collective
1000: Input Arguments:
1001: + dim - The cell dimension
1002: . level - The number of points in one dimension, 2^l
1003: . a - left end of interval (often-1)
1004: - b - right end of interval (often +1)
1006: Output Argument:
1007: . q - A PetscQuadrature object
1009: Level: intermediate
1011: .seealso: PetscDTGaussTensorQuadrature()
1012: @*/
1013: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1014: {
1015: const PetscInt p = 16; /* Digits of precision in the evaluation */
1016: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1017: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1018: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1019: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
1020: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
1021: PetscReal *x, *w;
1022: PetscInt K, k, npoints;
1023: PetscErrorCode ierr;
1026: if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1027: if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1028: /* Find K such that the weights are < 32 digits of precision */
1029: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1030: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1031: }
1032: PetscQuadratureCreate(PETSC_COMM_SELF, q);
1033: PetscQuadratureSetOrder(*q, 2*K+1);
1034: npoints = 2*K-1;
1035: PetscMalloc1(npoints*dim, &x);
1036: PetscMalloc1(npoints, &w);
1037: /* Center term */
1038: x[0] = beta;
1039: w[0] = 0.5*alpha*PETSC_PI;
1040: for (k = 1; k < K; ++k) {
1041: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1042: xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1043: x[2*k-1] = -alpha*xk+beta;
1044: w[2*k-1] = wk;
1045: x[2*k+0] = alpha*xk+beta;
1046: w[2*k+0] = wk;
1047: }
1048: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1049: return(0);
1050: }
1052: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1053: {
1054: const PetscInt p = 16; /* Digits of precision in the evaluation */
1055: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1056: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1057: PetscReal h = 1.0; /* Step size, length between x_k */
1058: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1059: PetscReal osum = 0.0; /* Integral on last level */
1060: PetscReal psum = 0.0; /* Integral on the level before the last level */
1061: PetscReal sum; /* Integral on current level */
1062: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1063: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1064: PetscReal wk; /* Quadrature weight at x_k */
1065: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1066: PetscInt d; /* Digits of precision in the integral */
1069: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1070: /* Center term */
1071: func(beta, &lval);
1072: sum = 0.5*alpha*PETSC_PI*lval;
1073: /* */
1074: do {
1075: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1076: PetscInt k = 1;
1078: ++l;
1079: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1080: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1081: psum = osum;
1082: osum = sum;
1083: h *= 0.5;
1084: sum *= 0.5;
1085: do {
1086: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1087: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1088: lx = -alpha*(1.0 - yk)+beta;
1089: rx = alpha*(1.0 - yk)+beta;
1090: func(lx, &lval);
1091: func(rx, &rval);
1092: lterm = alpha*wk*lval;
1093: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1094: sum += lterm;
1095: rterm = alpha*wk*rval;
1096: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1097: sum += rterm;
1098: ++k;
1099: /* Only need to evaluate every other point on refined levels */
1100: if (l != 1) ++k;
1101: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
1103: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1104: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1105: d3 = PetscLog10Real(maxTerm) - p;
1106: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1107: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1108: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1109: } while (d < digits && l < 12);
1110: *sol = sum;
1112: return(0);
1113: }
1115: #if defined(PETSC_HAVE_MPFR)
1116: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1117: {
1118: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
1119: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1120: mpfr_t alpha; /* Half-width of the integration interval */
1121: mpfr_t beta; /* Center of the integration interval */
1122: mpfr_t h; /* Step size, length between x_k */
1123: mpfr_t osum; /* Integral on last level */
1124: mpfr_t psum; /* Integral on the level before the last level */
1125: mpfr_t sum; /* Integral on current level */
1126: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1127: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1128: mpfr_t wk; /* Quadrature weight at x_k */
1129: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1130: PetscInt d; /* Digits of precision in the integral */
1131: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
1134: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1135: /* Create high precision storage */
1136: 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);
1137: /* Initialization */
1138: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1139: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
1140: mpfr_set_d(osum, 0.0, MPFR_RNDN);
1141: mpfr_set_d(psum, 0.0, MPFR_RNDN);
1142: mpfr_set_d(h, 1.0, MPFR_RNDN);
1143: mpfr_const_pi(pi2, MPFR_RNDN);
1144: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1145: /* Center term */
1146: func(0.5*(b+a), &lval);
1147: mpfr_set(sum, pi2, MPFR_RNDN);
1148: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1149: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1150: /* */
1151: do {
1152: PetscReal d1, d2, d3, d4;
1153: PetscInt k = 1;
1155: ++l;
1156: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1157: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1158: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1159: mpfr_set(psum, osum, MPFR_RNDN);
1160: mpfr_set(osum, sum, MPFR_RNDN);
1161: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
1162: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1163: do {
1164: mpfr_set_si(kh, k, MPFR_RNDN);
1165: mpfr_mul(kh, kh, h, MPFR_RNDN);
1166: /* Weight */
1167: mpfr_set(wk, h, MPFR_RNDN);
1168: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1169: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1170: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1171: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1172: mpfr_sqr(tmp, tmp, MPFR_RNDN);
1173: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1174: mpfr_div(wk, wk, tmp, MPFR_RNDN);
1175: /* Abscissa */
1176: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1177: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1178: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1179: mpfr_exp(tmp, msinh, MPFR_RNDN);
1180: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1181: /* Quadrature points */
1182: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1183: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1184: mpfr_add(lx, lx, beta, MPFR_RNDU);
1185: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1186: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1187: mpfr_add(rx, rx, beta, MPFR_RNDD);
1188: /* Evaluation */
1189: func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1190: func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1191: /* Update */
1192: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1193: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1194: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1195: mpfr_abs(tmp, tmp, MPFR_RNDN);
1196: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1197: mpfr_set(curTerm, tmp, MPFR_RNDN);
1198: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1199: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1200: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1201: mpfr_abs(tmp, tmp, MPFR_RNDN);
1202: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1203: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1204: ++k;
1205: /* Only need to evaluate every other point on refined levels */
1206: if (l != 1) ++k;
1207: mpfr_log10(tmp, wk, MPFR_RNDN);
1208: mpfr_abs(tmp, tmp, MPFR_RNDN);
1209: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1210: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1211: mpfr_abs(tmp, tmp, MPFR_RNDN);
1212: mpfr_log10(tmp, tmp, MPFR_RNDN);
1213: d1 = mpfr_get_d(tmp, MPFR_RNDN);
1214: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1215: mpfr_abs(tmp, tmp, MPFR_RNDN);
1216: mpfr_log10(tmp, tmp, MPFR_RNDN);
1217: d2 = mpfr_get_d(tmp, MPFR_RNDN);
1218: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1219: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1220: mpfr_log10(tmp, curTerm, MPFR_RNDN);
1221: d4 = mpfr_get_d(tmp, MPFR_RNDN);
1222: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1223: } while (d < digits && l < 8);
1224: *sol = mpfr_get_d(sum, MPFR_RNDN);
1225: /* Cleanup */
1226: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1227: return(0);
1228: }
1229: #else
1231: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1232: {
1233: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1234: }
1235: #endif
1237: /* Overwrites A. Can only handle full-rank problems with m>=n
1238: * A in column-major format
1239: * Ainv in row-major format
1240: * tau has length m
1241: * worksize must be >= max(1,n)
1242: */
1243: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1244: {
1246: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1247: PetscScalar *A,*Ainv,*R,*Q,Alpha;
1250: #if defined(PETSC_USE_COMPLEX)
1251: {
1252: PetscInt i,j;
1253: PetscMalloc2(m*n,&A,m*n,&Ainv);
1254: for (j=0; j<n; j++) {
1255: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1256: }
1257: mstride = m;
1258: }
1259: #else
1260: A = A_in;
1261: Ainv = Ainv_out;
1262: #endif
1264: PetscBLASIntCast(m,&M);
1265: PetscBLASIntCast(n,&N);
1266: PetscBLASIntCast(mstride,&lda);
1267: PetscBLASIntCast(worksize,&ldwork);
1268: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1269: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1270: PetscFPTrapPop();
1271: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1272: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1274: /* Extract an explicit representation of Q */
1275: Q = Ainv;
1276: PetscArraycpy(Q,A,mstride*n);
1277: K = N; /* full rank */
1278: #if defined(PETSC_MISSING_LAPACK_ORGQR)
1279: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ORGQR - Lapack routine is unavailable.");
1280: #else
1281: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1282: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1283: #endif
1285: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1286: Alpha = 1.0;
1287: ldb = lda;
1288: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1289: /* Ainv is Q, overwritten with inverse */
1291: #if defined(PETSC_USE_COMPLEX)
1292: {
1293: PetscInt i;
1294: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1295: PetscFree2(A,Ainv);
1296: }
1297: #endif
1298: return(0);
1299: }
1301: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1302: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1303: {
1305: PetscReal *Bv;
1306: PetscInt i,j;
1309: PetscMalloc1((ninterval+1)*ndegree,&Bv);
1310: /* Point evaluation of L_p on all the source vertices */
1311: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1312: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1313: for (i=0; i<ninterval; i++) {
1314: for (j=0; j<ndegree; j++) {
1315: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1316: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1317: }
1318: }
1319: PetscFree(Bv);
1320: return(0);
1321: }
1323: /*@
1324: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
1326: Not Collective
1328: Input Arguments:
1329: + degree - degree of reconstruction polynomial
1330: . nsource - number of source intervals
1331: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1332: . ntarget - number of target intervals
1333: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
1335: Output Arguments:
1336: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
1338: Level: advanced
1340: .seealso: PetscDTLegendreEval()
1341: @*/
1342: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1343: {
1345: PetscInt i,j,k,*bdegrees,worksize;
1346: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1347: PetscScalar *tau,*work;
1353: 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);
1354: #if defined(PETSC_USE_DEBUG)
1355: for (i=0; i<nsource; i++) {
1356: 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]);
1357: }
1358: for (i=0; i<ntarget; i++) {
1359: 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]);
1360: }
1361: #endif
1362: xmin = PetscMin(sourcex[0],targetx[0]);
1363: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1364: center = (xmin + xmax)/2;
1365: hscale = (xmax - xmin)/2;
1366: worksize = nsource;
1367: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1368: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1369: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1370: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1371: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1372: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1373: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1374: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1375: for (i=0; i<ntarget; i++) {
1376: PetscReal rowsum = 0;
1377: for (j=0; j<nsource; j++) {
1378: PetscReal sum = 0;
1379: for (k=0; k<degree+1; k++) {
1380: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1381: }
1382: R[i*nsource+j] = sum;
1383: rowsum += sum;
1384: }
1385: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1386: }
1387: PetscFree4(bdegrees,sourcey,Bsource,work);
1388: PetscFree4(tau,Bsinv,targety,Btarget);
1389: return(0);
1390: }
1392: /*@C
1393: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
1395: Not Collective
1397: Input Parameter:
1398: + n - the number of GLL nodes
1399: . nodes - the GLL nodes
1400: . weights - the GLL weights
1401: . f - the function values at the nodes
1403: Output Parameter:
1404: . in - the value of the integral
1406: Level: beginner
1408: .seealso: PetscDTGaussLobattoLegendreQuadrature()
1410: @*/
1411: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
1412: {
1413: PetscInt i;
1416: *in = 0.;
1417: for (i=0; i<n; i++) {
1418: *in += f[i]*f[i]*weights[i];
1419: }
1420: return(0);
1421: }
1423: /*@C
1424: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
1426: Not Collective
1428: Input Parameter:
1429: + n - the number of GLL nodes
1430: . nodes - the GLL nodes
1431: . weights - the GLL weights
1433: Output Parameter:
1434: . A - the stiffness element
1436: Level: beginner
1438: Notes:
1439: Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()
1441: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)
1443: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
1445: @*/
1446: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1447: {
1448: PetscReal **A;
1449: PetscErrorCode ierr;
1450: const PetscReal *gllnodes = nodes;
1451: const PetscInt p = n-1;
1452: PetscReal z0,z1,z2 = -1,x,Lpj,Lpr;
1453: PetscInt i,j,nn,r;
1456: PetscMalloc1(n,&A);
1457: PetscMalloc1(n*n,&A[0]);
1458: for (i=1; i<n; i++) A[i] = A[i-1]+n;
1460: for (j=1; j<p; j++) {
1461: x = gllnodes[j];
1462: z0 = 1.;
1463: z1 = x;
1464: for (nn=1; nn<p; nn++) {
1465: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1466: z0 = z1;
1467: z1 = z2;
1468: }
1469: Lpj=z2;
1470: for (r=1; r<p; r++) {
1471: if (r == j) {
1472: A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
1473: } else {
1474: x = gllnodes[r];
1475: z0 = 1.;
1476: z1 = x;
1477: for (nn=1; nn<p; nn++) {
1478: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1479: z0 = z1;
1480: z1 = z2;
1481: }
1482: Lpr = z2;
1483: A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
1484: }
1485: }
1486: }
1487: for (j=1; j<p+1; j++) {
1488: x = gllnodes[j];
1489: z0 = 1.;
1490: z1 = x;
1491: for (nn=1; nn<p; nn++) {
1492: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1493: z0 = z1;
1494: z1 = z2;
1495: }
1496: Lpj = z2;
1497: A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
1498: A[0][j] = A[j][0];
1499: }
1500: for (j=0; j<p; j++) {
1501: x = gllnodes[j];
1502: z0 = 1.;
1503: z1 = x;
1504: for (nn=1; nn<p; nn++) {
1505: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1506: z0 = z1;
1507: z1 = z2;
1508: }
1509: Lpj=z2;
1511: A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
1512: A[j][p] = A[p][j];
1513: }
1514: A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
1515: A[p][p]=A[0][0];
1516: *AA = A;
1517: return(0);
1518: }
1520: /*@C
1521: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
1523: Not Collective
1525: Input Parameter:
1526: + n - the number of GLL nodes
1527: . nodes - the GLL nodes
1528: . weights - the GLL weightss
1529: - A - the stiffness element
1531: Level: beginner
1533: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()
1535: @*/
1536: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1537: {
1541: PetscFree((*AA)[0]);
1542: PetscFree(*AA);
1543: *AA = NULL;
1544: return(0);
1545: }
1547: /*@C
1548: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
1550: Not Collective
1552: Input Parameter:
1553: + n - the number of GLL nodes
1554: . nodes - the GLL nodes
1555: . weights - the GLL weights
1557: Output Parameter:
1558: . AA - the stiffness element
1559: - AAT - the transpose of AA (pass in NULL if you do not need this array)
1561: Level: beginner
1563: Notes:
1564: Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()
1566: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
1568: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
1570: @*/
1571: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1572: {
1573: PetscReal **A, **AT = NULL;
1574: PetscErrorCode ierr;
1575: const PetscReal *gllnodes = nodes;
1576: const PetscInt p = n-1;
1577: PetscReal q,qp,Li, Lj,d0;
1578: PetscInt i,j;
1581: PetscMalloc1(n,&A);
1582: PetscMalloc1(n*n,&A[0]);
1583: for (i=1; i<n; i++) A[i] = A[i-1]+n;
1585: if (AAT) {
1586: PetscMalloc1(n,&AT);
1587: PetscMalloc1(n*n,&AT[0]);
1588: for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
1589: }
1591: if (n==1) {A[0][0] = 0.;}
1592: d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
1593: for (i=0; i<n; i++) {
1594: for (j=0; j<n; j++) {
1595: A[i][j] = 0.;
1596: qAndLEvaluation(p,gllnodes[i],&q,&qp,&Li);
1597: qAndLEvaluation(p,gllnodes[j],&q,&qp,&Lj);
1598: if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
1599: if ((j==i) && (i==0)) A[i][j] = -d0;
1600: if (j==i && i==p) A[i][j] = d0;
1601: if (AT) AT[j][i] = A[i][j];
1602: }
1603: }
1604: if (AAT) *AAT = AT;
1605: *AA = A;
1606: return(0);
1607: }
1609: /*@C
1610: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()
1612: Not Collective
1614: Input Parameter:
1615: + n - the number of GLL nodes
1616: . nodes - the GLL nodes
1617: . weights - the GLL weights
1618: . AA - the stiffness element
1619: - AAT - the transpose of the element
1621: Level: beginner
1623: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()
1625: @*/
1626: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1627: {
1631: PetscFree((*AA)[0]);
1632: PetscFree(*AA);
1633: *AA = NULL;
1634: if (*AAT) {
1635: PetscFree((*AAT)[0]);
1636: PetscFree(*AAT);
1637: *AAT = NULL;
1638: }
1639: return(0);
1640: }
1642: /*@C
1643: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
1645: Not Collective
1647: Input Parameter:
1648: + n - the number of GLL nodes
1649: . nodes - the GLL nodes
1650: . weights - the GLL weightss
1652: Output Parameter:
1653: . AA - the stiffness element
1655: Level: beginner
1657: Notes:
1658: Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()
1660: This is the same as the Gradient operator multiplied by the diagonal mass matrix
1662: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
1664: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()
1666: @*/
1667: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1668: {
1669: PetscReal **D;
1670: PetscErrorCode ierr;
1671: const PetscReal *gllweights = weights;
1672: const PetscInt glln = n;
1673: PetscInt i,j;
1676: PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
1677: for (i=0; i<glln; i++){
1678: for (j=0; j<glln; j++) {
1679: D[i][j] = gllweights[i]*D[i][j];
1680: }
1681: }
1682: *AA = D;
1683: return(0);
1684: }
1686: /*@C
1687: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
1689: Not Collective
1691: Input Parameter:
1692: + n - the number of GLL nodes
1693: . nodes - the GLL nodes
1694: . weights - the GLL weights
1695: - A - advection
1697: Level: beginner
1699: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()
1701: @*/
1702: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1703: {
1707: PetscFree((*AA)[0]);
1708: PetscFree(*AA);
1709: *AA = NULL;
1710: return(0);
1711: }
1713: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1714: {
1715: PetscReal **A;
1716: PetscErrorCode ierr;
1717: const PetscReal *gllweights = weights;
1718: const PetscInt glln = n;
1719: PetscInt i,j;
1722: PetscMalloc1(glln,&A);
1723: PetscMalloc1(glln*glln,&A[0]);
1724: for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
1725: if (glln==1) {A[0][0] = 0.;}
1726: for (i=0; i<glln; i++) {
1727: for (j=0; j<glln; j++) {
1728: A[i][j] = 0.;
1729: if (j==i) A[i][j] = gllweights[i];
1730: }
1731: }
1732: *AA = A;
1733: return(0);
1734: }
1736: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1737: {
1741: PetscFree((*AA)[0]);
1742: PetscFree(*AA);
1743: *AA = NULL;
1744: return(0);
1745: }