Loading...
Searching...
No Matches
quadrature_t.cc
Go to the documentation of this file.
#include "rheolef/quadrature.h"
#include "rheolef/gauss_jacobi.h"
namespace rheolef{
using namespace std;
/* ------------------------------------------
* Gauss quadrature formulae on the triangle
* order r : exact for all polynoms of P_r
*
* References for optimal formulae :
* 1) G. Dhatt and G. Touzot,
* Une presentation de la methode des elements finis,
* Maloine editeur, Paris
* Deuxieme edition,
* 1984,
* page 297
* 2) M. Crouzeix et A. L. Mignot
* Exercices d'analyse des equations differentielles,
* Masson,
* 1986,
* p. 61
* PROBLEM: direct optimized formulae have hard-coded
* coefficients in double precision.
* TODO: compute numeric value at the fly why full
* precision.
* ------------------------------------------
*/
template<class T>
void
quadrature_on_geo<T>::init_triangle (quadrature_option opt)
{
// -------------------------------------------------------------------------
// special case : superconvergent patch recovery nodes & weights
// -------------------------------------------------------------------------
quadrature_option::family_type f = opt.get_family();
if (f == quadrature_option::superconvergent) {
switch (opt.get_order()) {
case 0:
case 1:
f = quadrature_option::gauss;
break;
case 2:
f = quadrature_option::middle_edge;
break;
default:
error_macro ("unsupported superconvergent("<<opt.get_order()<<")");
}
}
// -------------------------------------------------------------------------
// special case : Middle Edge formulae
// e.g. when using special option in riesz_representer
// -------------------------------------------------------------------------
if (f == quadrature_option::middle_edge) {
switch (opt.get_order()) {
case 0 :
case 1 :
case 2 :
// middle edge formulae:
wx(x(T(0.5), T(0.0)), T(1)/6);
wx(x(T(0.0), T(0.5)), T(1)/6);
wx(x(T(0.5), T(0.5)), T(1)/6);
break;
default:
error_macro ("unsupported Middle-Edge("<<opt.get_order()<<")");
}
return;
}
// -------------------------------------------------------------------------
// special case : Gauss-Lobatto points
// e.g. when using special option in riesz() function
// -------------------------------------------------------------------------
if (f == quadrature_option::gauss_lobatto) {
switch (opt.get_order()) {
case 0 :
case 1 :
// trapezes:
wx(x(T(0), T(0)), T(1)/6);
wx(x(T(1), T(0)), T(1)/6);
wx(x(T(0), T(1)), T(1)/6);
break;
case 2 :
// simpson:
wx(x(T(0), T(0)), 3/T(120));
wx(x(T(1), T(0)), 3/T(120));
wx(x(T(0), T(1)), 3/T(120));
wx(x(T(0.5), T(0.0)), 8/T(120));
wx(x(T(0.0), T(0.5)), 8/T(120));
wx(x(T(0.5), T(0.5)), 8/T(120));
wx(x(1/T(3), 1/T(3)), 27/T(120));
break;
default:
error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
}
return;
}
// -------------------------------------------------------------------------
// 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(1/T(3),1/T(3)), 0.5);
} else {
size_type n = (r+1)*(r+2)/2;
T w = 0.5/T(int(n));
for (size_type i = 0; i <= r; i++) {
for (size_type j = 0; i+j <= r; j++) {
wx (x(T(int(i))/r,T(int(j))/r), w);
}
}
}
return;
}
// -------------------------------------------------------------------------
// default & general case : Gauss points
// -------------------------------------------------------------------------
check_macro (f == quadrature_option::gauss,
"unsupported quadrature family \"" << opt.get_family_name() << "\"");
switch (opt.get_order()) {
case 0 :
case 1 :
// better than the general case:
// r=0 => 2 nodes
// r=1 => 4 nodes
// here : 1 node
wx(x(1/T(3), 1/T(3)), 0.5);
break;
case 2 :
// better than the general case:
// r=2 => 6 nodes
// here : 3 node
wx(x(1/T(6), 1/T(6)), 1/T(6));
wx(x(4/T(6), 1/T(6)), 1/T(6));
wx(x(1/T(6), 4/T(6)), 1/T(6));
break;
#if !(defined(_RHEOLEF_HAVE_CLN) || defined(_RHEOLEF_HAVE_LONG_DOUBLE) || defined(_RHEOLEF_HAVE_FLOAT128))
// numerical values are inlined in double precision : 15 digits !
case 3 :
case 4 : {
// better than the general case:
// r=3 => 9 nodes
// r=4 => 12 nodes
// here : 6 node
T a = 0.445948490915965;
T b = 0.091576213509771;
T wa = 0.111690794839005;
T wb = 0.054975871827661;
wx(x(a, a), wa);
wx(x(1-2*a, a), wa);
wx(x(a, 1-2*a), wa);
wx(x(b, b), wb);
wx(x(1-2*b, b), wb);
wx(x(b, 1-2*b), wb);
break;
}
case 5 :
case 6 : {
// better than the general case:
// r=5 => 16 nodes
// r=6 => 20 nodes
// here : 12 node
T a = 0.063089014491502;
T b = 0.249286745170910;
T c = 0.310352451033785;
T d = 0.053145049844816;
T wa = 0.025422453185103;
T wb = 0.058393137863189;
T wc = 0.041425537809187;
wx(x(a, a), wa);
wx(x(1-2*a, a), wa);
wx(x(a, 1-2*a), wa);
wx(x(b, b), wb);
wx(x(1-2*b, b), wb);
wx(x(b, 1-2*b), wb);
wx(x(c, d), wc);
wx(x(d, c), wc);
wx(x(1-(c+d), c), wc);
wx(x(1-(c+d), d), wc);
wx(x(c, 1-(c+d)), wc);
wx(x(d, 1-(c+d)), wc);
break;
}
#endif // BIGFLOAT
default: {
// Gauss-Legendre quadrature formulae
// where Legendre = Jacobi(alpha=0,beta=0) polynoms
size_type r = opt.get_order();
size_type n0 = n_node_gauss(r);
size_type n1 = n_node_gauss(r+1);
vector<T> zeta0(n0), omega0(n0);
vector<T> zeta1(n1), omega1(n1);
gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
for (size_type i = 0; i < n0; i++) {
for (size_type j = 0; j < n1; j++) {
// we transform the square into the triangle
T eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
T eta_1 = (1+zeta1[j])/2;
T J = (1-zeta1[j])/8;
wx (x(eta_0,eta_1), J*omega0[i]*omega1[j]);
}
}
}
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template void quadrature_on_geo<T>::init_triangle (quadrature_option);
_RHEOLEF_instanciation(Float)
}// namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
Float
see the Float page for the full documentation
rheolef::quadrature_on_geo::x
point_basic< T > x
Definition quadrature.h:89
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition quadrature.h:88
rheolef::quadrature_on_geo::n_node_gauss
static size_type n_node_gauss(size_type r)
Definition quadrature.h:245
rheolef::quadrature_on_geo::init_triangle
void init_triangle(quadrature_option opt)
Definition quadrature_t.cc:50
rheolef::quadrature_on_geo::wx
void wx(const point_basic< T > &x, const T &w)
Definition quadrature.h:116
error_macro
#define error_macro(message)
Definition dis_macros.h:49
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_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
d
size_t d
Definition neumann-laplace-lambda.cc:28
f
Definition cavity_dg.h:29