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
Float
see the Float page for the full documentation
rheolef::point_basic
Definition: point.h:87
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition: quadrature.h:88
rheolef::quadrature_on_geo::init_edge
void init_edge(quadrature_option opt)
Definition: quadrature_e.cc:38
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::gauss_jacobi
void gauss_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_jacobi.icc:29
f
Definition: cavity_dg.h:29