quadrature_e.cc

Go to the documentation of this file.

#include "rheolef/quadrature.h"
#include "rheolef/gauss_jacobi.h"
namespace rheolef {
using namespace std;
 
/*
 * quadrature formulae
 * refs: G. Dhatt and G. Touzot,
 *       Une presentation de la methode des elements finis,
 *       Maloine editeur, Paris
 *       Deuxieme edition,
 *       1984,
 *       page 288
 */
 
template<class T>
void
quadrature_on_geo<T>::init_edge (quadrature_option opt)
{
  quadrature_option::family_type f = opt.get_family();
  // -------------------------------------------------------------------------
  // special case : equispaced, for irregular (e.g. Heaviside) functions
  // -------------------------------------------------------------------------
  if (f == quadrature_option::equispaced) {
    size_type r = opt.get_order();
    if (r == 0) {
      wx (x(0.5), 1);
    } else {
      size_type n = r+1;
      T w = T(1)/T(int(n));
      for (size_type i = 0; i <= r; i++) {
        wx (x(T(int(i))/r), w);
      }
    }
    return;
  }
  // -------------------------------------------------------------------------
  // TODO: special case : superconvergent patch recovery nodes & weights
  // -------------------------------------------------------------------------
 
  // -------------------------------------------------------------------------
  // special case : Gauss-Lobatto points
  // e.g. when using special option in riesz_representer
  // -------------------------------------------------------------------------
  if (f == quadrature_option::gauss_lobatto) {
    switch (opt.get_order()) {
     case 0 :
     case 1 :
            // trapezes:
            wx(x(T(0)), 0.5);
            wx(x(T(1)), 0.5);
            break;
     case 2 :
     case 3 :
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
            wx(x(T(0)),   T(1)/T(6));
            wx(x(T(0.5)), T(4)/T(6));
            wx(x(T(1)),   T(1)/T(6));
            break;
     case 4 :
     case 5 :
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
            wx(x(T(0)),               T(1)/T(12));
            wx(x(0.5-0.5/sqrt(T(5))), T(5)/T(12));
            wx(x(0.5+0.5/sqrt(T(5))), T(5)/T(12));
            wx(x(T(1)),               T(1)/T(12));
            break;
     case 6 :
     case 7 :
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
            wx(x(T(0)),                    T( 9)/T(180));
            wx(x(0.5-0.5/sqrt(T(3)/T(7))), T(49)/T(180));
            wx(x(T(0.5)),                  T(64)/T(180));
            wx(x(0.5+0.5/sqrt(T(3)/T(7))), T(49)/T(180));
            wx(x(T(1)),                    T( 9)/T(180));
            break;
     default:
            error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
    }
    return;
  }
  // -------------------------------------------------------------------------
  // default & general case : Gauss points
  // -------------------------------------------------------------------------
  check_macro (f == quadrature_option::gauss, 
        "unsupported quadrature family \"" << opt.get_family_name() << "\"");
 
  // Gauss-Legendre quadrature formulae 
  //  where Legendre = Jacobi(alpha=0,beta=0) polynoms
  size_type n = n_node_gauss(opt.get_order());
  vector<T> zeta(n), omega(n);
  gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
  for (size_type i = 0; i < n; i++)
    wx (x((1+zeta[i])/2), omega[i]/2);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template void quadrature_on_geo<T>::init_edge (quadrature_option);
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef