Rheolef  7.1
an efficient C++ finite element environment
quadrature_H.cc
Go to the documentation of this file.
1 #include "rheolef/quadrature.h"
22 #include "rheolef/gauss_jacobi.h"
23 namespace rheolef {
24 using namespace std;
25 
26 template<class T>
27 void
29 {
30  quadrature_option::family_type f = opt.get_family();
31  // -------------------------------------------------------------------------
32  // special case : equispaced, for irregular (e.g. Heaviside) functions
33  // -------------------------------------------------------------------------
34  if (f == quadrature_option::equispaced) {
35  size_type r = opt.get_order();
36  if (r == 0) {
37  wx (x(0,0,0), 8);
38  } else {
39  size_type n = (r+1)*(r+1)*(r+1);
40  T w = 8/T(int(n));
41  for (size_type i = 0; i <= r; i++) {
42  for (size_type j = 0; j <= r; j++) {
43  for (size_type k = 0; k <= r; k++) {
44  wx (x(2*T(int(i))/r-1,2*T(int(j))/r-1,2*T(int(k))/r-1), w);
45  }
46  }
47  }
48  }
49  return;
50  }
51  // -------------------------------------------------------------------------
52  // TODO: special case : superconvergent patch recovery nodes & weights
53  // -------------------------------------------------------------------------
54 
55  // -------------------------------------------------------------------------
56  // default: gauss
57  // -------------------------------------------------------------------------
58  check_macro (f == quadrature_option::gauss,
59  "unsupported quadrature family \"" << opt.get_family_name() << "\"");
60 
61  // Gauss-Legendre quadrature formulae
62  // where Legendre = Jacobi(alpha=0,beta=0) polynoms
63  // when using n nodes : quadrature formulae order is 2*r-1
64  size_type n = n_node_gauss(opt.get_order());
65  vector<T> zeta(n), omega(n);
66  gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
67  for (size_type i = 0; i < n; i++)
68  for (size_type j = 0; j < n; j++)
69  for (size_type k = 0; k < n; k++)
70  wx (x(zeta[i], zeta[j], zeta[k]), omega[i]*omega[j]*omega[k]);
71 }
72 // ----------------------------------------------------------------------------
73 // instanciation in library
74 // ----------------------------------------------------------------------------
75 #define _RHEOLEF_instanciation(T) \
76 template void quadrature_on_geo<T>::init_hexahedron (quadrature_option);
77 
79 
80 }// namespace rheolef
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
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
rheolef::point_basic
Definition: point.h:87
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition: quadrature.h:88
Float
see the Float page for the full documentation
rheolef::quadrature_on_geo::init_hexahedron
void init_hexahedron(quadrature_option opt)
Definition: quadrature_H.cc:28
f
Definition: cavity_dg.h:29
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
T
Expr1::float_type T
Definition: field_expr.h:261