xref: /petsc/src/dm/dt/interface/dt.c (revision 37045ce4b26a202cc7a954b115583b9ceffe94d5) 
1*37045ce4SJed Brown /* Discretization tools */
2*37045ce4SJed Brown 
3*37045ce4SJed Brown #include <petscdt.h>            /*I "petscdt.h" I*/
4*37045ce4SJed Brown #include <petscblaslapack.h>
5*37045ce4SJed Brown 
6*37045ce4SJed Brown #undef __FUNCT__
7*37045ce4SJed Brown #define __FUNCT__ "PetscDTLegendreEval"
8*37045ce4SJed Brown /*@
9*37045ce4SJed Brown    PetscDTLegendreEval - evaluate Legendre polynomial at points
10*37045ce4SJed Brown 
11*37045ce4SJed Brown    Not Collective
12*37045ce4SJed Brown 
13*37045ce4SJed Brown    Input Arguments:
14*37045ce4SJed Brown +  npoints - number of spatial points to evaluate at
15*37045ce4SJed Brown .  points - array of locations to evaluate at
16*37045ce4SJed Brown .  ndegree - number of basis degrees to evaluate
17*37045ce4SJed Brown -  degrees - sorted array of degrees to evaluate
18*37045ce4SJed Brown 
19*37045ce4SJed Brown    Output Arguments:
20*37045ce4SJed Brown +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or PETSC_NULL)
21*37045ce4SJed Brown .  D - row-oriented derivative evaluation matrix (or PETSC_NULL)
22*37045ce4SJed Brown -  D2 - row-oriented second derivative evaluation matrix (or PETSC_NULL)
23*37045ce4SJed Brown 
24*37045ce4SJed Brown    Level: intermediate
25*37045ce4SJed Brown 
26*37045ce4SJed Brown .seealso: PetscDTGaussQuadrature()
27*37045ce4SJed Brown @*/
28*37045ce4SJed Brown PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
29*37045ce4SJed Brown {
30*37045ce4SJed Brown   PetscInt i,maxdegree;
31*37045ce4SJed Brown 
32*37045ce4SJed Brown   PetscFunctionBegin;
33*37045ce4SJed Brown   if (!npoints || !ndegree) PetscFunctionReturn(0);
34*37045ce4SJed Brown   maxdegree = degrees[ndegree-1];
35*37045ce4SJed Brown   for (i=0; i<npoints; i++) {
36*37045ce4SJed Brown     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
37*37045ce4SJed Brown     PetscInt j,k;
38*37045ce4SJed Brown     x = points[i];
39*37045ce4SJed Brown     pm2 = 0;
40*37045ce4SJed Brown     pm1 = 1;
41*37045ce4SJed Brown     pd2 = 0;
42*37045ce4SJed Brown     pd1 = 0;
43*37045ce4SJed Brown     pdd2 = 0;
44*37045ce4SJed Brown     pdd1 = 0;
45*37045ce4SJed Brown     k = 0;
46*37045ce4SJed Brown     if (degrees[k] == 0) {
47*37045ce4SJed Brown       if (B) B[i*ndegree+k] = pm1;
48*37045ce4SJed Brown       if (D) D[i*ndegree+k] = pd1;
49*37045ce4SJed Brown       if (D2) D2[i*ndegree+k] = pdd1;
50*37045ce4SJed Brown       k++;
51*37045ce4SJed Brown     }
52*37045ce4SJed Brown     for (j=1; j<=maxdegree; j++,k++) {
53*37045ce4SJed Brown       PetscReal p,d,dd;
54*37045ce4SJed Brown       p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
55*37045ce4SJed Brown       d = pd2 + (2*j-1)*pm1;
56*37045ce4SJed Brown       dd = pdd2 + (2*j-1)*pd1;
57*37045ce4SJed Brown       pm2 = pm1;
58*37045ce4SJed Brown       pm1 = p;
59*37045ce4SJed Brown       pd2 = pd1;
60*37045ce4SJed Brown       pd1 = d;
61*37045ce4SJed Brown       pdd2 = pdd1;
62*37045ce4SJed Brown       pdd1 = dd;
63*37045ce4SJed Brown       if (degrees[k] == j) {
64*37045ce4SJed Brown         if (B) B[i*ndegree+k] = p;
65*37045ce4SJed Brown         if (D) D[i*ndegree+k] = d;
66*37045ce4SJed Brown         if (D2) D2[i*ndegree+k] = dd;
67*37045ce4SJed Brown       }
68*37045ce4SJed Brown     }
69*37045ce4SJed Brown   }
70*37045ce4SJed Brown   PetscFunctionReturn(0);
71*37045ce4SJed Brown }
72*37045ce4SJed Brown 
73*37045ce4SJed Brown #undef __FUNCT__
74*37045ce4SJed Brown #define __FUNCT__ "PetscDTGaussQuadrature"
75*37045ce4SJed Brown /*@
76*37045ce4SJed Brown    PetscDTGaussQuadrature - create Gauss quadrature
77*37045ce4SJed Brown 
78*37045ce4SJed Brown    Not Collective
79*37045ce4SJed Brown 
80*37045ce4SJed Brown    Input Arguments:
81*37045ce4SJed Brown +  npoints - number of points
82*37045ce4SJed Brown .  a - left end of interval (often-1)
83*37045ce4SJed Brown -  b - right end of interval (often +1)
84*37045ce4SJed Brown 
85*37045ce4SJed Brown    Output Arguments:
86*37045ce4SJed Brown +  x - quadrature points
87*37045ce4SJed Brown -  w - quadrature weights
88*37045ce4SJed Brown 
89*37045ce4SJed Brown    Level: intermediate
90*37045ce4SJed Brown 
91*37045ce4SJed Brown    References:
92*37045ce4SJed Brown    Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969.
93*37045ce4SJed Brown 
94*37045ce4SJed Brown .seealso: PetscDTLegendreEval()
95*37045ce4SJed Brown @*/
96*37045ce4SJed Brown PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
97*37045ce4SJed Brown {
98*37045ce4SJed Brown   PetscErrorCode ierr;
99*37045ce4SJed Brown   PetscInt i;
100*37045ce4SJed Brown   PetscReal *work;
101*37045ce4SJed Brown   PetscScalar *Z;
102*37045ce4SJed Brown #if defined(PETSC_HAVE_COMPLEX)
103*37045ce4SJed Brown   char compz[] = "I";           /* compute the eigenvectors of the tridiagonal matrix directly */
104*37045ce4SJed Brown #else
105*37045ce4SJed Brown   char compz[] = "V";
106*37045ce4SJed Brown #endif
107*37045ce4SJed Brown   PetscBLASInt N,LDZ,info;
108*37045ce4SJed Brown 
109*37045ce4SJed Brown   PetscFunctionBegin;
110*37045ce4SJed Brown   /* Set up the Golub-Welsch system */
111*37045ce4SJed Brown   for (i=0; i<npoints; i++) {
112*37045ce4SJed Brown     x[i] = 0;                   /* diagonal is 0 */
113*37045ce4SJed Brown     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
114*37045ce4SJed Brown   }
115*37045ce4SJed Brown   ierr = PetscRealView(npoints-1,w,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
116*37045ce4SJed Brown   ierr = PetscMalloc2(npoints*npoints,PetscScalar,&Z,PetscMax(1,2*npoints-2),PetscReal,&work);CHKERRQ(ierr);
117*37045ce4SJed Brown   N = PetscBLASIntCast(npoints);
118*37045ce4SJed Brown   LDZ = N;
119*37045ce4SJed Brown   ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
120*37045ce4SJed Brown   LAPACKstev_(compz,&N,x,w,Z,&LDZ,work,&info);
121*37045ce4SJed Brown   ierr = PetscFPTrapPop();CHKERRQ(ierr);
122*37045ce4SJed Brown   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEV/xSTEQR error");
123*37045ce4SJed Brown 
124*37045ce4SJed Brown   for (i=0; i<(npoints+1)/2; i++) {
125*37045ce4SJed Brown     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
126*37045ce4SJed Brown     x[i] = (a+b)/2 - y*(b-a)/2;
127*37045ce4SJed Brown     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
128*37045ce4SJed Brown     w[i] = w[npoints-1-i] = (b-a)*PetscSqr(0.5*PetscAbsScalar(Z[i*npoints] + Z[(npoints-i-1)*npoints]));
129*37045ce4SJed Brown   }
130*37045ce4SJed Brown   ierr = PetscFree2(Z,work);CHKERRQ(ierr);
131*37045ce4SJed Brown   PetscFunctionReturn(0);
132*37045ce4SJed Brown }
133