Rheolef  7.1
an efficient C++ finite element environment
quadrature_e.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 /*
27  * quadrature formulae
28  * refs: G. Dhatt and G. Touzot,
29  * Une presentation de la methode des elements finis,
30  * Maloine editeur, Paris
31  * Deuxieme edition,
32  * 1984,
33  * page 288
34  */
35 
36 template<class T>
37 void
39 {
40  quadrature_option::family_type f = opt.get_family();
41  // -------------------------------------------------------------------------
42  // special case : equispaced, for irregular (e.g. Heaviside) functions
43  // -------------------------------------------------------------------------
44  if (f == quadrature_option::equispaced) {
45  size_type r = opt.get_order();
46  if (r == 0) {
47  wx (x(0.5), 1);
48  } else {
49  size_type n = r+1;
50  T w = T(1)/T(int(n));
51  for (size_type i = 0; i <= r; i++) {
52  wx (x(T(int(i))/r), w);
53  }
54  }
55  return;
56  }
57  // -------------------------------------------------------------------------
58  // TODO: special case : superconvergent patch recovery nodes & weights
59  // -------------------------------------------------------------------------
60 
61  // -------------------------------------------------------------------------
62  // special case : Gauss-Lobatto points
63  // e.g. when using special option in riesz_representer
64  // -------------------------------------------------------------------------
65  if (f == quadrature_option::gauss_lobatto) {
66  switch (opt.get_order()) {
67  case 0 :
68  case 1 :
69  // trapezes:
70  wx(x(T(0)), 0.5);
71  wx(x(T(1)), 0.5);
72  break;
73  case 2 :
74  case 3 :
75  // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
76  wx(x(T(0)), T(1)/T(6));
77  wx(x(T(0.5)), T(4)/T(6));
78  wx(x(T(1)), T(1)/T(6));
79  break;
80  case 4 :
81  case 5 :
82  // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
83  wx(x(T(0)), T(1)/T(12));
84  wx(x(0.5-0.5/sqrt(T(5))), T(5)/T(12));
85  wx(x(0.5+0.5/sqrt(T(5))), T(5)/T(12));
86  wx(x(T(1)), T(1)/T(12));
87  break;
88  case 6 :
89  case 7 :
90  // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
91  wx(x(T(0)), T( 9)/T(180));
92  wx(x(0.5-0.5/sqrt(T(3)/T(7))), T(49)/T(180));
93  wx(x(T(0.5)), T(64)/T(180));
94  wx(x(0.5+0.5/sqrt(T(3)/T(7))), T(49)/T(180));
95  wx(x(T(1)), T( 9)/T(180));
96  break;
97  default:
98  error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
99  }
100  return;
101  }
102  // -------------------------------------------------------------------------
103  // default & general case : Gauss points
104  // -------------------------------------------------------------------------
105  check_macro (f == quadrature_option::gauss,
106  "unsupported quadrature family \"" << opt.get_family_name() << "\"");
107 
108  // Gauss-Legendre quadrature formulae
109  // where Legendre = Jacobi(alpha=0,beta=0) polynoms
110  size_type n = n_node_gauss(opt.get_order());
111  vector<T> zeta(n), omega(n);
112  gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
113  for (size_type i = 0; i < n; i++)
114  wx (x((1+zeta[i])/2), omega[i]/2);
115 }
116 // ----------------------------------------------------------------------------
117 // instanciation in library
118 // ----------------------------------------------------------------------------
119 #define _RHEOLEF_instanciation(T) \
120 template void quadrature_on_geo<T>::init_edge (quadrature_option);
121 
123 
124 }// namespace rheolef
rheolef::quadrature_on_geo::init_edge
void init_edge(quadrature_option opt)
Definition: quadrature_e.cc:38
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
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
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