Rheolef  7.2
an efficient C++ finite element environment
gauss_lobatto_jacobi.icc
Go to the documentation of this file.
1 #include "rheolef/compiler.h"
22 #include "rheolef/gamma.h"
23 #include "rheolef/jacobi.h"
24 #include "rheolef/jacobi_roots.h"
25 #include <iterator>
26 namespace rheolef {
27 template <class Size, class OutputIterator1, class OutputIterator2>
28 void gauss_lobatto_jacobi (Size R,
29  typename std::iterator_traits<OutputIterator1>::value_type alpha,
30  typename std::iterator_traits<OutputIterator1>::value_type beta,
31  OutputIterator1 zeta, OutputIterator2 omega)
32 {
33  typedef typename std::iterator_traits<OutputIterator1>::value_type T;
34  check_macro (R >= 2, "gauss_lobatto_jacobi: node number " << R << " may be >= 2");
35  T num = pow(T(2), alpha+beta+1)/T(1.*R-1);
36  T w0 = pow(T(2), alpha+beta+1)*(beta+1);
37  T wf = pow(T(2), alpha+beta+1)*(alpha+1);
38  if (alpha == floor(alpha) && beta == floor(beta)) {
39  num *= T(1.*R)/((alpha+R)*(beta+R));
40  w0 *= 1/(T(R-1.)*(alpha+R));
41  wf *= 1/(T(R-1.)*(beta+R));
42  for (Size k = 1; k <= size_t(static_cast<int>(beta)); k++) {
43  num *= T(1.*R+k)/(alpha+T(1.*R+k));
44  w0 *= T(sqr(int(k)))/(T(1.*R-1+k)*(alpha+R+k));
45  }
46  for (Size k = 1; k <= size_t(static_cast<int>(alpha)); k++) {
47  wf *= T(sqr(int(k)))/(T(1.*R-1+k)*(beta+T(1.*R+k)));
48  }
49  } else {
50  num *= (my_gamma(alpha+R)/my_gamma(alpha+beta+R+1))
51  *(my_gamma(beta+R)/my_gamma(T(1.*R)));
52  w0 *= (my_gamma(alpha+R)/my_gamma(alpha+beta+R+1))
53  *(sqr(my_gamma(beta+1))*my_gamma(T(1.*R-1))/my_gamma(beta+R));
54  wf *= (my_gamma(beta+R)/my_gamma(alpha+beta+R+1))
55  *(sqr(my_gamma(alpha+1))*my_gamma(T(1.*R-1))/my_gamma(alpha+R));
56  }
57  zeta [0] = -1;
58  omega[0] = w0;
59  jacobi_roots (R-2, alpha+1, beta+1, zeta+1);
60  jacobi<T> P1 (R-1, alpha, beta);
61  for (Size r = 1; r < R-1; ++r)
62  omega[r] = num/sqr(P1(zeta[r]));
63  zeta [R-1] = 1;
64  omega[R-1] = wf;
65 }
66 } // namespace rheolef
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
This file is part of Rheolef.
void gauss_lobatto_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
T my_gamma(const T &x)
Definition: gamma.icc:25
Float beta[][pmax+1]