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