Rheolef  7.1
an efficient C++ finite element environment
quadrature_Pr.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(1/T(3),1/T(3),0), 1);
38  } else {
39  size_type n = ((r+1)*(r+2)/2)*(r+1);
40  T w = 1/T(int(n));
41  for (size_type i = 0; i <= r; i++) {
42  for (size_type j = 0; i+j <= r; j++) {
43  for (size_type k = 0; k <= r; k++) {
44  wx (x(T(int(i))/r,T(int(j))/r,2*T(int(k))/r-1), w);
45  }
46  }
47  }
48  }
49  return;
50  }
51  // -------------------------------------------------------------------------
52  // general case: Gauss
53  // -------------------------------------------------------------------------
54  check_macro (f == quadrature_option::gauss,
55  "unsupported quadrature family \"" << opt.get_family_name() << "\"");
56 
57  switch (opt.get_order()) {
58  case 0:
59  case 1: {
60  // central point:
61  // better than the general case
62  // r=0 : 2 nodes
63  // r=1 : 8 nodes
64  // here: 1 node
65  T a = T(1)/3;
66  wx(x(a,a,0), 1);
67  break;
68  }
69  default: {
70  // Gauss-Legendre quadrature formulae
71  // where Legendre = Jacobi(alpha=0,beta=0) polynoms
72  size_t r = opt.get_order();
73  size_type n0 = n_node_gauss(r);
74  size_type n1 = n_node_gauss(r+1);
75  size_type n2 = n_node_gauss(r);
76  vector<T> zeta0(n0), omega0(n0);
77  vector<T> zeta1(n1), omega1(n1);
78  vector<T> zeta2(n2), omega2(n2);
79  gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
80  gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
81  gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
82 
83  for (size_type i = 0; i < n0; i++) {
84  for (size_type j = 0; j < n1; j++) {
85  for (size_type k = 0; k < n2; k++) {
86  // we transform the square into the triangle
87  T eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
88  T eta_1 = (1+zeta1[j])/2;
89  T eta_2 = zeta2[k];
90  T J = (1-zeta1[j])/8;
91  wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
92  }
93  }
94  }
95  }
96  }
97 }
98 // ----------------------------------------------------------------------------
99 // instanciation in library
100 // ----------------------------------------------------------------------------
101 #define _RHEOLEF_instanciation(T) \
102 template void quadrature_on_geo<T>::init_prism (quadrature_option);
103 
105 
106 }// namespace rheolef
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
a
Definition: diffusion_isotropic.h:25
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
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
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:218
rheolef::quadrature_on_geo::init_prism
void init_prism(quadrature_option opt)
Definition: quadrature_Pr.cc:28