137045ce4SJed Brown /* Discretization tools */ 237045ce4SJed Brown 337045ce4SJed Brown #include <petscdt.h> /*I "petscdt.h" I*/ 437045ce4SJed Brown #include <petscblaslapack.h> 537045ce4SJed Brown 637045ce4SJed Brown #undef __FUNCT__ 737045ce4SJed Brown #define __FUNCT__ "PetscDTLegendreEval" 837045ce4SJed Brown /*@ 937045ce4SJed Brown PetscDTLegendreEval - evaluate Legendre polynomial at points 1037045ce4SJed Brown 1137045ce4SJed Brown Not Collective 1237045ce4SJed Brown 1337045ce4SJed Brown Input Arguments: 1437045ce4SJed Brown + npoints - number of spatial points to evaluate at 1537045ce4SJed Brown . points - array of locations to evaluate at 1637045ce4SJed Brown . ndegree - number of basis degrees to evaluate 1737045ce4SJed Brown - degrees - sorted array of degrees to evaluate 1837045ce4SJed Brown 1937045ce4SJed Brown Output Arguments: 2037045ce4SJed Brown + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or PETSC_NULL) 2137045ce4SJed Brown . D - row-oriented derivative evaluation matrix (or PETSC_NULL) 2237045ce4SJed Brown - D2 - row-oriented second derivative evaluation matrix (or PETSC_NULL) 2337045ce4SJed Brown 2437045ce4SJed Brown Level: intermediate 2537045ce4SJed Brown 2637045ce4SJed Brown .seealso: PetscDTGaussQuadrature() 2737045ce4SJed Brown @*/ 2837045ce4SJed Brown PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 2937045ce4SJed Brown { 3037045ce4SJed Brown PetscInt i,maxdegree; 3137045ce4SJed Brown 3237045ce4SJed Brown PetscFunctionBegin; 3337045ce4SJed Brown if (!npoints || !ndegree) PetscFunctionReturn(0); 3437045ce4SJed Brown maxdegree = degrees[ndegree-1]; 3537045ce4SJed Brown for (i=0; i<npoints; i++) { 3637045ce4SJed Brown PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x; 3737045ce4SJed Brown PetscInt j,k; 3837045ce4SJed Brown x = points[i]; 3937045ce4SJed Brown pm2 = 0; 4037045ce4SJed Brown pm1 = 1; 4137045ce4SJed Brown pd2 = 0; 4237045ce4SJed Brown pd1 = 0; 4337045ce4SJed Brown pdd2 = 0; 4437045ce4SJed Brown pdd1 = 0; 4537045ce4SJed Brown k = 0; 4637045ce4SJed Brown if (degrees[k] == 0) { 4737045ce4SJed Brown if (B) B[i*ndegree+k] = pm1; 4837045ce4SJed Brown if (D) D[i*ndegree+k] = pd1; 4937045ce4SJed Brown if (D2) D2[i*ndegree+k] = pdd1; 5037045ce4SJed Brown k++; 5137045ce4SJed Brown } 5237045ce4SJed Brown for (j=1; j<=maxdegree; j++,k++) { 5337045ce4SJed Brown PetscReal p,d,dd; 5437045ce4SJed Brown p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j; 5537045ce4SJed Brown d = pd2 + (2*j-1)*pm1; 5637045ce4SJed Brown dd = pdd2 + (2*j-1)*pd1; 5737045ce4SJed Brown pm2 = pm1; 5837045ce4SJed Brown pm1 = p; 5937045ce4SJed Brown pd2 = pd1; 6037045ce4SJed Brown pd1 = d; 6137045ce4SJed Brown pdd2 = pdd1; 6237045ce4SJed Brown pdd1 = dd; 6337045ce4SJed Brown if (degrees[k] == j) { 6437045ce4SJed Brown if (B) B[i*ndegree+k] = p; 6537045ce4SJed Brown if (D) D[i*ndegree+k] = d; 6637045ce4SJed Brown if (D2) D2[i*ndegree+k] = dd; 6737045ce4SJed Brown } 6837045ce4SJed Brown } 6937045ce4SJed Brown } 7037045ce4SJed Brown PetscFunctionReturn(0); 7137045ce4SJed Brown } 7237045ce4SJed Brown 7337045ce4SJed Brown #undef __FUNCT__ 7437045ce4SJed Brown #define __FUNCT__ "PetscDTGaussQuadrature" 7537045ce4SJed Brown /*@ 7637045ce4SJed Brown PetscDTGaussQuadrature - create Gauss quadrature 7737045ce4SJed Brown 7837045ce4SJed Brown Not Collective 7937045ce4SJed Brown 8037045ce4SJed Brown Input Arguments: 8137045ce4SJed Brown + npoints - number of points 8237045ce4SJed Brown . a - left end of interval (often-1) 8337045ce4SJed Brown - b - right end of interval (often +1) 8437045ce4SJed Brown 8537045ce4SJed Brown Output Arguments: 8637045ce4SJed Brown + x - quadrature points 8737045ce4SJed Brown - w - quadrature weights 8837045ce4SJed Brown 8937045ce4SJed Brown Level: intermediate 9037045ce4SJed Brown 9137045ce4SJed Brown References: 9237045ce4SJed Brown Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969. 9337045ce4SJed Brown 9437045ce4SJed Brown .seealso: PetscDTLegendreEval() 9537045ce4SJed Brown @*/ 9637045ce4SJed Brown PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w) 9737045ce4SJed Brown { 9837045ce4SJed Brown PetscErrorCode ierr; 9937045ce4SJed Brown PetscInt i; 10037045ce4SJed Brown PetscReal *work; 10137045ce4SJed Brown PetscScalar *Z; 10237045ce4SJed Brown PetscBLASInt N,LDZ,info; 10337045ce4SJed Brown 10437045ce4SJed Brown PetscFunctionBegin; 10537045ce4SJed Brown /* Set up the Golub-Welsch system */ 10637045ce4SJed Brown for (i=0; i<npoints; i++) { 10737045ce4SJed Brown x[i] = 0; /* diagonal is 0 */ 10837045ce4SJed Brown if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i)); 10937045ce4SJed Brown } 11037045ce4SJed Brown ierr = PetscRealView(npoints-1,w,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 11137045ce4SJed Brown ierr = PetscMalloc2(npoints*npoints,PetscScalar,&Z,PetscMax(1,2*npoints-2),PetscReal,&work);CHKERRQ(ierr); 112*c5df96a5SBarry Smith ierr = PetscBLASIntCast(npoints,&N);CHKERRQ(ierr); 11337045ce4SJed Brown LDZ = N; 11437045ce4SJed Brown ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1151c3d6f74SJed Brown LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info); 11637045ce4SJed Brown ierr = PetscFPTrapPop();CHKERRQ(ierr); 1171c3d6f74SJed Brown if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 11837045ce4SJed Brown 11937045ce4SJed Brown for (i=0; i<(npoints+1)/2; i++) { 12037045ce4SJed Brown PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */ 12137045ce4SJed Brown x[i] = (a+b)/2 - y*(b-a)/2; 12237045ce4SJed Brown x[npoints-i-1] = (a+b)/2 + y*(b-a)/2; 12337045ce4SJed Brown w[i] = w[npoints-1-i] = (b-a)*PetscSqr(0.5*PetscAbsScalar(Z[i*npoints] + Z[(npoints-i-1)*npoints])); 12437045ce4SJed Brown } 12537045ce4SJed Brown ierr = PetscFree2(Z,work);CHKERRQ(ierr); 12637045ce4SJed Brown PetscFunctionReturn(0); 12737045ce4SJed Brown } 128