Loading...
Searching...
No Matches
quadrature_Pr.cc
Go to the documentation of this file.
#include "rheolef/quadrature.h"
#include "rheolef/gauss_jacobi.h"
namespace rheolef {
using namespace std;
template<class T>
void
quadrature_on_geo<T>::init_prism (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(1/T(3),1/T(3),0), 1);
} else {
size_type n = ((r+1)*(r+2)/2)*(r+1);
T w = 1/T(int(n));
for (size_type i = 0; i <= r; i++) {
for (size_type j = 0; i+j <= r; j++) {
for (size_type k = 0; k <= r; k++) {
wx (x(T(int(i))/r,T(int(j))/r,2*T(int(k))/r-1), w);
}
}
}
}
return;
}
// -------------------------------------------------------------------------
// general case: Gauss
// -------------------------------------------------------------------------
check_macro (f == quadrature_option::gauss,
"unsupported quadrature family \"" << opt.get_family_name() << "\"");
switch (opt.get_order()) {
case 0:
case 1: {
// central point:
// better than the general case
// r=0 : 2 nodes
// r=1 : 8 nodes
// here: 1 node
T a = T(1)/3;
wx(x(a,a,0), 1);
break;
}
default: {
// Gauss-Legendre quadrature formulae
// where Legendre = Jacobi(alpha=0,beta=0) polynoms
size_t r = opt.get_order();
size_type n0 = n_node_gauss(r);
size_type n1 = n_node_gauss(r+1);
size_type n2 = n_node_gauss(r);
vector<T> zeta0(n0), omega0(n0);
vector<T> zeta1(n1), omega1(n1);
vector<T> zeta2(n2), omega2(n2);
gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
for (size_type i = 0; i < n0; i++) {
for (size_type j = 0; j < n1; j++) {
for (size_type k = 0; k < n2; k++) {
// 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 eta_2 = zeta2[k];
T J = (1-zeta1[j])/8;
wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
}
}
}
}
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template void quadrature_on_geo<T>::init_prism (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::init_prism
void init_prism(quadrature_option opt)
Definition quadrature_Pr.cc:28
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::wx
void wx(const point_basic< T > &x, const T &w)
Definition quadrature.h:116
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
f
Definition cavity_dg.h:29