Rheolef  7.1
an efficient C++ finite element environment
quadrature_q.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  * Gauss quadrature formulae on the triangle
28  * order r : exact for all polynoms of P_r
29  *
30  * References for optimal formulae :
31  * 1) G. Dhatt and G. Touzot,
32  * Une presentation de la methode des elements finis,
33  * Maloine editeur, Paris
34  * Deuxieme edition,
35  * 1984,
36  * page 293
37  * 2) M. Crouzeix et A. L. Mignot
38  * Exercices d'analyse des equations differentielles,
39  * Masson,
40  * 1986,
41  * p. 63
42  * ------------------------------------------
43  */
44 template<class T>
45 void
47 {
48  quadrature_option::family_type f = opt.get_family();
49  // -------------------------------------------------------------------------
50  // special case : equispaced, for irregular (e.g. Heaviside) functions
51  // -------------------------------------------------------------------------
52  if (f == quadrature_option::equispaced) {
53  size_type r = opt.get_order();
54  if (r == 0) {
55  wx (x(0,0), 4);
56  } else {
57  size_type n = (r+1)*(r+1);
58  T w = 4/T(int(n));
59  for (size_type i = 0; i <= r; i++) {
60  for (size_type j = 0; j <= r; j++) {
61  wx (x(2*T(int(i))/r-1,2*T(int(j))/r-1), w);
62  }
63  }
64  }
65  return;
66  }
67  // -------------------------------------------------------------------------
68  // TODO: special case : superconvergent patch recovery nodes & weights
69  // -------------------------------------------------------------------------
70 
71  // -------------------------------------------------------------------------
72  // default: gauss
73  // -------------------------------------------------------------------------
74  check_macro (f == quadrature_option::gauss,
75  "unsupported quadrature family \"" << opt.get_family_name() << "\"");
76 
77  // Gauss-Legendre quadrature formulae
78  // where Legendre = Jacobi(alpha=0,beta=0) polynoms
79  // when using n nodes : quadrature formulae order is 2*r-1
80  size_type n = n_node_gauss(opt.get_order());
81  vector<T> zeta(n), omega(n);
82  gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
83  for (size_type i = 0; i < n; i++)
84  for (size_type j = 0; j < n; j++)
85  wx (x(zeta[i], zeta[j]), omega[i]*omega[j]);
86 }
87 // ----------------------------------------------------------------------------
88 // instanciation in library
89 // ----------------------------------------------------------------------------
90 #define _RHEOLEF_instanciation(T) \
91 template void quadrature_on_geo<T>::init_square (quadrature_option);
92 
94 
95 }// 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
rheolef::quadrature_on_geo::init_square
void init_square(quadrature_option opt)
Definition: quadrature_q.cc:46
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