1 #ifndef DUNE_GEOMETRY_QUADRATURERULES_JACOBI_N_0_H 2 #define DUNE_GEOMETRY_QUADRATURERULES_JACOBI_N_0_H 3 4 #include <vector> 5 #include <type_traits> 6 7 #include <dune/common/fvector.hh> 8 #include <dune/common/fmatrix.hh> 9 #include <dune/common/exceptions.hh> 10 11 #include <dune/common/dynmatrix.hh> 12 #include <dune/common/dynvector.hh> 13 #include <dune/common/dynmatrixev.hh> 14 15 16 namespace Dune { 17 18 /* 19 * This class calculates 1D quadrature rules of dynamic order, that respect the 20 * weightfunctions resulting from the conical product. 21 * For quadrature rules with order <= highest order (61) and dimension <= 3 exact quadrature rules are used. 22 * Else the quadrature rules are computed by Lapack and thus a floating point type is required. 23 * For more information see the comment in `tensorproductquadrature.hh` 24 * 25 */ 26 27 template<typename ct, int dim> 28 class JacobiNQuadratureRule; 29 30 template<typename ct> 31 using JacobiNQuadratureRule1D = JacobiNQuadratureRule<ct,1>; 32 33 template<typename ct> 34 class JacobiNQuadratureRule<ct,1> : public QuadratureRule<ct,1> 35 { 36 public: 37 // compile time parameters 38 enum { dim=1 }; 39 40 private: 41 42 typedef QuadratureRule<ct, dim> Rule; 43 44 friend class QuadratureRuleFactory<ct,dim>; 45 46 template< class ctype, int dimension> 47 friend class TensorProductQuadratureRule; 48 JacobiNQuadratureRule(int const order,int const alpha=0)49 explicit JacobiNQuadratureRule (int const order, int const alpha=0) 50 : Rule( GeometryTypes::line ) 51 { 52 if (unsigned(order) > maxOrder()) 53 DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule " << order << " not supported!"); 54 auto&& rule = decideRule(order,alpha); 55 for( auto qpoint : rule ) 56 this->push_back(qpoint); 57 this->delivered_order = 2*rule.size()-1; 58 59 } 60 maxOrder()61 static unsigned maxOrder() 62 { 63 return 127; // can be changed 64 } 65 decideRule(int const degree,int const alpha)66 QuadratureRule<ct,1> decideRule(int const degree, int const alpha) 67 { 68 const auto maxOrder = std::min( unsigned(GaussQuadratureRule1D<ct>::highest_order), 69 std::min( unsigned(Jacobi1QuadratureRule1D<ct>::highest_order), 70 unsigned(Jacobi2QuadratureRule1D<ct>::highest_order)) 71 ); 72 return unsigned(degree) < maxOrder ? decideRuleExponent(degree,alpha) : UseLapackOrError<ct>(degree, alpha); 73 } 74 75 #if HAVE_LAPACK 76 template<typename type, std::enable_if_t<!(std::is_floating_point<type>::value)>* = nullptr> UseLapackOrError(int const degree,int const alpha)77 QuadratureRule<ct,1> UseLapackOrError( int const degree, int const alpha) 78 { 79 DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule with degree: "<< degree << " and jacobi exponent: "<< alpha<< " is not supported for this type!"); 80 } 81 82 template<typename type, std::enable_if_t<std::is_floating_point<type>::value>* = nullptr> UseLapackOrError(int const degree,int const alpha)83 QuadratureRule<ct,1> UseLapackOrError( int const degree, int const alpha) 84 { 85 return jacobiNodesWeights<ct>(degree, alpha); 86 } 87 #else 88 template<typename type> UseLapackOrError(int const,int const)89 QuadratureRule<ct,1> UseLapackOrError( int const, int const) 90 { 91 DUNE_THROW(NotImplemented, "LAPACK must be enable to use JacobiN quadrature rules."); 92 } 93 #endif 94 decideRuleExponent(int const degree,int const alpha)95 QuadratureRule<ct,1> decideRuleExponent(int const degree, int const alpha) 96 { 97 switch(alpha) 98 { 99 case 0 : return QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussLegendre); 100 case 1 : { 101 // scale already existing weights by 0.5 102 auto&& rule = QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussJacobi_1_0); 103 QuadratureRule<ct,1> quadratureRule; 104 quadratureRule.reserve(rule.size()); 105 for( auto qpoint : rule ) 106 quadratureRule.push_back(QuadraturePoint<ct,1>(qpoint.position(),ct(0.5)*qpoint.weight())); 107 return quadratureRule; 108 } 109 case 2 : { 110 // scale already existing weights by 0.25 111 auto&& rule = QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussJacobi_2_0); 112 QuadratureRule<ct,1> quadratureRule; 113 quadratureRule.reserve(rule.size()); 114 for( auto qpoint : rule ) 115 quadratureRule.push_back(QuadraturePoint<ct,1>(qpoint.position(),ct(0.25)*qpoint.weight())); 116 return quadratureRule; 117 } 118 default : return UseLapackOrError<ct>(degree,alpha); 119 } 120 } 121 122 #if HAVE_LAPACK 123 // computes the nodes and weights for the weight function (1-x)^alpha, which are exact for given polynomials with degree "degree" 124 template<typename ctype, std::enable_if_t<std::is_floating_point<ctype>::value>* = nullptr> jacobiNodesWeights(int const degree,int const alpha)125 QuadratureRule<ct,1> jacobiNodesWeights(int const degree, int const alpha) 126 { 127 using std::sqrt; 128 129 // compute the degree of the needed jacobi polynomial 130 const int n = degree/2 +1; 131 132 DynamicMatrix<double> J(n,n,0); 133 134 J[0][0] = -double(alpha)/(2 + alpha); 135 for(int i=1; i<n; ++i) 136 { 137 ct a_i = double(2*i*(i + alpha)) /( (2*i + alpha - 1)*(2*i + alpha) ); 138 ct c_i = double(2*i*(i + alpha)) /( (2*i + alpha + 1)*(2*i + alpha) ); 139 140 J[i][i] = -double(alpha*alpha) /( (2*i + alpha + 2)*(2*i + alpha) ); 141 J[i][i-1] = sqrt(a_i*c_i); 142 J[i-1][i] = J[i][i-1]; 143 } 144 145 DynamicVector<std::complex<double> > eigenValues(n,0); 146 std::vector<DynamicVector<double> > eigenVectors(n, DynamicVector<double>(n,0)); 147 148 DynamicMatrixHelp::eigenValuesNonSym(J, eigenValues, &eigenVectors); 149 150 double mu = 1.0/(alpha + 1); 151 152 QuadratureRule<ct,1> quadratureRule; 153 quadratureRule.reserve(n); 154 for (int i=0; i<n; ++i) 155 { 156 auto&& eV0 = eigenVectors[i][0]; 157 ct weight = mu * eV0*eV0; 158 DynamicVector<ct> node(1,0.5*eigenValues[i].real() + 0.5); 159 160 // bundle the nodes and the weights 161 QuadraturePoint<ct,1> temp(node, weight); 162 quadratureRule.push_back(temp); 163 } 164 165 return quadratureRule; 166 167 } 168 #endif 169 170 }; 171 172 } 173 174 #endif 175