Rheolef  7.1
an efficient C++ finite element environment
quadrature_t.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 297
37  * 2) M. Crouzeix et A. L. Mignot
38  * Exercices d'analyse des equations differentielles,
39  * Masson,
40  * 1986,
41  * p. 61
42  * PROBLEM: direct optimized formulae have hard-coded
43  * coefficients in double precision.
44  * TODO: compute numeric value at the fly why full
45  * precision.
46  * ------------------------------------------
47  */
48 template<class T>
49 void
51 {
52  // -------------------------------------------------------------------------
53  // special case : superconvergent patch recovery nodes & weights
54  // -------------------------------------------------------------------------
55  quadrature_option::family_type f = opt.get_family();
56  if (f == quadrature_option::superconvergent) {
57  switch (opt.get_order()) {
58  case 0:
59  case 1:
60  f = quadrature_option::gauss;
61  break;
62  case 2:
63  f = quadrature_option::middle_edge;
64  break;
65  default:
66  error_macro ("unsupported superconvergent("<<opt.get_order()<<")");
67  }
68  }
69  // -------------------------------------------------------------------------
70  // special case : Middle Edge formulae
71  // e.g. when using special option in riesz_representer
72  // -------------------------------------------------------------------------
73  if (f == quadrature_option::middle_edge) {
74  switch (opt.get_order()) {
75  case 0 :
76  case 1 :
77  case 2 :
78  // middle edge formulae:
79  wx(x(T(0.5), T(0.0)), T(1)/6);
80  wx(x(T(0.0), T(0.5)), T(1)/6);
81  wx(x(T(0.5), T(0.5)), T(1)/6);
82  break;
83  default:
84  error_macro ("unsupported Middle-Edge("<<opt.get_order()<<")");
85  }
86  return;
87  }
88  // -------------------------------------------------------------------------
89  // special case : Gauss-Lobatto points
90  // e.g. when using special option in riesz() function
91  // -------------------------------------------------------------------------
92  if (f == quadrature_option::gauss_lobatto) {
93  switch (opt.get_order()) {
94  case 0 :
95  case 1 :
96  // trapezes:
97  wx(x(T(0), T(0)), T(1)/6);
98  wx(x(T(1), T(0)), T(1)/6);
99  wx(x(T(0), T(1)), T(1)/6);
100  break;
101  case 2 :
102  // simpson:
103  wx(x(T(0), T(0)), 3/T(120));
104  wx(x(T(1), T(0)), 3/T(120));
105  wx(x(T(0), T(1)), 3/T(120));
106  wx(x(T(0.5), T(0.0)), 8/T(120));
107  wx(x(T(0.0), T(0.5)), 8/T(120));
108  wx(x(T(0.5), T(0.5)), 8/T(120));
109  wx(x(1/T(3), 1/T(3)), 27/T(120));
110  break;
111  default:
112  error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
113  }
114  return;
115  }
116  // -------------------------------------------------------------------------
117  // special case : equispaced, for irregular (e.g. Heaviside) functions
118  // -------------------------------------------------------------------------
119  if (f == quadrature_option::equispaced) {
120  size_type r = opt.get_order();
121  if (r == 0) {
122  wx (x(1/T(3),1/T(3)), 0.5);
123  } else {
124  size_type n = (r+1)*(r+2)/2;
125  T w = 0.5/T(int(n));
126  for (size_type i = 0; i <= r; i++) {
127  for (size_type j = 0; i+j <= r; j++) {
128  wx (x(T(int(i))/r,T(int(j))/r), w);
129  }
130  }
131  }
132  return;
133  }
134  // -------------------------------------------------------------------------
135  // default & general case : Gauss points
136  // -------------------------------------------------------------------------
137  check_macro (f == quadrature_option::gauss,
138  "unsupported quadrature family \"" << opt.get_family_name() << "\"");
139 
140  switch (opt.get_order()) {
141  case 0 :
142  case 1 :
143  // better than the general case:
144  // r=0 => 2 nodes
145  // r=1 => 4 nodes
146  // here : 1 node
147  wx(x(1/T(3), 1/T(3)), 0.5);
148  break;
149  case 2 :
150  // better than the general case:
151  // r=2 => 6 nodes
152  // here : 3 node
153  wx(x(1/T(6), 1/T(6)), 1/T(6));
154  wx(x(4/T(6), 1/T(6)), 1/T(6));
155  wx(x(1/T(6), 4/T(6)), 1/T(6));
156  break;
157 #if !(defined(_RHEOLEF_HAVE_CLN) || defined(_RHEOLEF_HAVE_LONG_DOUBLE) || defined(_RHEOLEF_HAVE_FLOAT128))
158  // numerical values are inlined in double precision : 15 digits !
159  case 3 :
160  case 4 : {
161  // better than the general case:
162  // r=3 => 9 nodes
163  // r=4 => 12 nodes
164  // here : 6 node
165  T a = 0.445948490915965;
166  T b = 0.091576213509771;
167  T wa = 0.111690794839005;
168  T wb = 0.054975871827661;
169  wx(x(a, a), wa);
170  wx(x(1-2*a, a), wa);
171  wx(x(a, 1-2*a), wa);
172  wx(x(b, b), wb);
173  wx(x(1-2*b, b), wb);
174  wx(x(b, 1-2*b), wb);
175  break;
176  }
177  case 5 :
178  case 6 : {
179  // better than the general case:
180  // r=5 => 16 nodes
181  // r=6 => 20 nodes
182  // here : 12 node
183  T a = 0.063089014491502;
184  T b = 0.249286745170910;
185  T c = 0.310352451033785;
186  T d = 0.053145049844816;
187  T wa = 0.025422453185103;
188  T wb = 0.058393137863189;
189  T wc = 0.041425537809187;
190  wx(x(a, a), wa);
191  wx(x(1-2*a, a), wa);
192  wx(x(a, 1-2*a), wa);
193  wx(x(b, b), wb);
194  wx(x(1-2*b, b), wb);
195  wx(x(b, 1-2*b), wb);
196  wx(x(c, d), wc);
197  wx(x(d, c), wc);
198  wx(x(1-(c+d), c), wc);
199  wx(x(1-(c+d), d), wc);
200  wx(x(c, 1-(c+d)), wc);
201  wx(x(d, 1-(c+d)), wc);
202  break;
203  }
204 #endif // BIGFLOAT
205  default: {
206  // Gauss-Legendre quadrature formulae
207  // where Legendre = Jacobi(alpha=0,beta=0) polynoms
208  size_type r = opt.get_order();
209  size_type n0 = n_node_gauss(r);
210  size_type n1 = n_node_gauss(r+1);
211  vector<T> zeta0(n0), omega0(n0);
212  vector<T> zeta1(n1), omega1(n1);
213  gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
214  gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
215 
216  for (size_type i = 0; i < n0; i++) {
217  for (size_type j = 0; j < n1; j++) {
218  // we transform the square into the triangle
219  T eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
220  T eta_1 = (1+zeta1[j])/2;
221  T J = (1-zeta1[j])/8;
222  wx (x(eta_0,eta_1), J*omega0[i]*omega1[j]);
223  }
224  }
225  }
226  }
227 }
228 // ----------------------------------------------------------------------------
229 // instanciation in library
230 // ----------------------------------------------------------------------------
231 #define _RHEOLEF_instanciation(T) \
232 template void quadrature_on_geo<T>::init_triangle (quadrature_option);
233 
235 
236 }// namespace rheolef
rheolef::quadrature_on_geo::init_triangle
void init_triangle(quadrature_option opt)
Definition: quadrature_t.cc:50
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
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
a
Definition: diffusion_isotropic.h:25
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
f
Definition: cavity_dg.h:29
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
T
Expr1::float_type T
Definition: field_expr.h:261