Loading...
Searching...
No Matches
gauss_radau_jacobi.icc
Go to the documentation of this file.
#include "rheolef/compiler.h"
#include "rheolef/gamma.h"
#include "rheolef/jacobi.h"
#include "rheolef/jacobi_roots.h"
#include <iterator>
namespace rheolef {
template <class Size, class OutputIterator1, class OutputIterator2>
void gauss_radau_jacobi (Size R,
typename std::iterator_traits<OutputIterator1>::value_type alpha,
typename std::iterator_traits<OutputIterator1>::value_type beta,
OutputIterator1 zeta, OutputIterator2 omega)
{
typedef typename std::iterator_traits<OutputIterator1>::value_type T;
check_macro (R >= 1, "gauss_radau_jacobi: node number " << R << " may be >= 1");
T num = pow(T(2), alpha+beta)/(beta+T(1.*R));
T w0 = pow(T(2), alpha+beta+1)*(beta+1);
if (alpha == floor(alpha) && beta == floor(beta)) {
num *= T(1.*R)/((alpha+T(1.*R))*(beta+T(1.*R)));
w0 *= 1/(T(1.*R)*(alpha+T(1.*R)));
for (Size k = 1; k <= size_t(static_cast<int>(beta)); k++) {
num *= T(1.*R+k)/(alpha+T(1.*R+k));
w0 *= T(sqr(T(int(k))))/(T(1.*R+k)*(alpha+R+k));
}
} else {
num *= (my_gamma(alpha+R)/my_gamma(alpha+beta+R+1))
*(my_gamma(beta+R)/my_gamma(T(1.*R)));
w0 *= sqr(my_gamma(beta+1))
*(my_gamma(T(1.*R))/my_gamma(beta+R+1))
*(my_gamma(alpha+R)/my_gamma(alpha+beta+R+1));
}
zeta [0] = -1;
omega[0] = w0;
jacobi_roots (R-1, alpha, beta+1, zeta+1);
jacobi<T> P1 (R-1, alpha, beta);
for (Size r = 1; r < R; ++r)
omega[r] = num*(1-zeta[r])/sqr(P1(zeta[r]));
}
} // namespace rheolef
rheolef::jacobi
Definition jacobi.icc:24
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::gauss_radau_jacobi
void gauss_radau_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)
Definition gauss_radau_jacobi.icc:28
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
rheolef::jacobi_roots
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
Definition jacobi_roots.icc:26
rheolef::my_gamma
T my_gamma(const T &x)
Definition gamma.icc:25