Loading...
Searching...
No Matches
quadrature_H.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_hexahedron (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,0,0), 8);
} else {
size_type n = (r+1)*(r+1)*(r+1);
T w = 8/T(int(n));
for (size_type i = 0; i <= r; i++) {
for (size_type j = 0; j <= r; j++) {
for (size_type k = 0; k <= r; k++) {
wx (x(2*T(int(i))/r-1,2*T(int(j))/r-1,2*T(int(k))/r-1), w);
}
}
}
}
return;
}
// -------------------------------------------------------------------------
// TODO: special case : superconvergent patch recovery nodes & weights
// -------------------------------------------------------------------------
// -------------------------------------------------------------------------
// default: gauss
// -------------------------------------------------------------------------
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
// when using n nodes : quadrature formulae order is 2*r-1
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++)
for (size_type j = 0; j < n; j++)
for (size_type k = 0; k < n; k++)
wx (x(zeta[i], zeta[j], zeta[k]), omega[i]*omega[j]*omega[k]);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template void quadrature_on_geo<T>::init_hexahedron (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::point_basic
Definition point.h:86
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition quadrature.h:88
rheolef::quadrature_on_geo::init_hexahedron
void init_hexahedron(quadrature_option opt)
Definition quadrature_H.cc:28
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
std
STL namespace.
f
Definition cavity_dg.h:29