NiHu  2.0
gaussian_quadrature.hpp
Go to the documentation of this file.
1 // This file is a part of NiHu, a C++ BEM template library.
2 //
3 // Copyright (C) 2012-2014 Peter Fiala <fiala@hit.bme.hu>
4 // Copyright (C) 2012-2014 Peter Rucz <rucz@hit.bme.hu>
5 //
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 //
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with this program. If not, see <http://www.gnu.org/licenses/>.
18 
25 #ifndef GAUSSIAN_QUADRATURE_HPP_INCLUDED
26 #define GAUSSIAN_QUADRATURE_HPP_INCLUDED
27 
28 #include "quadrature.hpp"
29 #include "../library/lib_domain.hpp"
30 
31 #include <stdexcept>
32 
33 namespace NiHu
34 {
35 
61 template <class scalar_t>
62 Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> gauss_impl(size_t N)
63 {
64  typedef Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic> mat_t;
65  mat_t J(N, N);
66  J.setZero();
67 
68  // Fill the diagonals of the matrix
69  for (size_t n = 0; n < N-1; ++n)
70  {
71  scalar_t v = scalar_t(n+1);
72  scalar_t u = v / sqrt(4.0*v*v - 1);
73  J(n,n+1) = u;
74  J(n+1,n) = u;
75  }
76 
77  // Get the eigenvalues and eigenvectors
78  Eigen::SelfAdjointEigenSolver<mat_t> solver(J);
79  Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> V;
80  V.resize(N,2);
81  V.col(0) = solver.eigenvalues();
82  V.col(1) = 2.0 * solver.eigenvectors().row(0).cwiseAbs2();
83 
84  return V;
85 }
86 
87 // forward declaration
88 template <class Domain>
89 class gaussian_quadrature;
90 
91 
96 template <class Domain>
97 struct quadrature_traits<gaussian_quadrature<Domain> >
98 {
100  typedef Domain domain_t;
101 };
102 
106 template <>
107 class gaussian_quadrature<line_domain> :
108  public quadrature_base<gaussian_quadrature<line_domain> >
109 {
110 public:
112  typedef quadrature_base<gaussian_quadrature<line_domain> > base_t;
114  typedef base_t::xi_t xi_t;
116  typedef base_t::scalar_t scalar_t;
117 
119  typedef gaussian_quadrature type;
120 
124  gaussian_quadrature() :
125  base_t(0)
126  {
127  }
128 
133  gaussian_quadrature(size_t degree) :
134  base_t(degree/2+1)
135  {
136  size_t N = degree/2+1;
137  // compute 1D Gaussian locations and weights
138  auto V = gauss_impl<scalar_t>(N);
139 
140  // Fill the points and weights
141  for(size_t i = 0; i < N; ++i)
142  {
143  xi_t xi;
144  xi << V(i,0);
145  push_back(quadrature_elem_t(xi, V(i,1)));
146  }
147  }
148 };
149 
150 
154 template <>
155 class gaussian_quadrature<quad_domain> :
156  public quadrature_base<gaussian_quadrature<quad_domain> >
157 {
158 public:
160  typedef quadrature_base<gaussian_quadrature<quad_domain> > base_t;
162  typedef base_t::domain_t domain_t;
164  typedef base_t::xi_t xi_t;
166  typedef base_t::scalar_t scalar_t;
167 
169  typedef gaussian_quadrature type;
170 
174  gaussian_quadrature() :
175  base_t(0)
176  {
177  }
178 
183  gaussian_quadrature(size_t degree) :
184  base_t((degree/2+1) * (degree/2+1))
185  {
186  size_t N = degree/2+1;
187  auto V = gauss_impl<scalar_t>(N);
188 
189  // Fill the points and weights
190  for(size_t i = 0; i < N; ++i)
191  {
192  for(size_t j = 0; j < N; ++j)
193  {
194  xi_t xi;
195  xi << V(i,0), V(j,0);
196  push_back(quadrature_elem_t(xi, V(i,1)*V(j,1)));
197  }
198  }
199  }
200 }; // end of class gauss_quad
201 
202 
206 static std::array<size_t, 15> const dunavant_num = {
207  1 /*0*/,
208  1 /*1*/,
209  3 /*2*/,
210  4 /*3*/,
211  6 /*4*/,
212  7 /*5*/,
213  12 /*6*/,
214  13 /*7*/,
215  16 /*8*/,
216  19 /*9*/,
217  25 /*10*/,
218  33 /*11*/,
219  33 /*12*/,
220  37 /*13*/,
221  42 /*14*/
222  };
223 
227 template <>
228 class gaussian_quadrature<tria_domain> :
229  public quadrature_base<gaussian_quadrature<tria_domain> >
230 {
231 public:
233  typedef quadrature_base<gaussian_quadrature<tria_domain> > base_t;
235  typedef base_t::quadrature_elem_t quadrature_elem_t;
237  typedef base_t::xi_t xi_t;
238 
240  typedef gaussian_quadrature type;
241 
242 
246  gaussian_quadrature() :
247  base_t(0)
248  {
249  }
250 
255  gaussian_quadrature(size_t degree) :
256  base_t(degree < dunavant_num.size() ? dunavant_num[degree] : 0)
257  {
258  switch(degree)
259  {
260  case 0:
261  case 1:
262  push_back(quadrature_elem_t(xi_t(1./3., 1./3.0), 1./2.0));
263  break;
264  case 2:
265  push_back(quadrature_elem_t(xi_t(1./6., 1./6.0), 1./6.0));
266  push_back(quadrature_elem_t(xi_t(4./6., 1./6.0), 1./6.0));
267  push_back(quadrature_elem_t(xi_t(1./6., 4./6.0), 1./6.0));
268  break;
269  case 3:
270  push_back(quadrature_elem_t(xi_t(1./3., 1./3.0), -0.281250000000000));
271  push_back(quadrature_elem_t(xi_t(1./5., 3./5.0), 0.260416666666667));
272  push_back(quadrature_elem_t(xi_t(1./5., 1./5.0), 0.260416666666667));
273  push_back(quadrature_elem_t(xi_t(3./5., 1./5.0), 0.260416666666667));
274  break;
275  case 4:
276  push_back(quadrature_elem_t(xi_t(0.445948490915965, 0.108103018168070), 0.111690794839006));
277  push_back(quadrature_elem_t(xi_t(0.445948490915965, 0.445948490915965), 0.111690794839006));
278  push_back(quadrature_elem_t(xi_t(0.108103018168070, 0.445948490915965), 0.111690794839006));
279  push_back(quadrature_elem_t(xi_t(0.091576213509771, 0.816847572980459), 0.054975871827661));
280  push_back(quadrature_elem_t(xi_t(0.091576213509771, 0.091576213509771), 0.054975871827661));
281  push_back(quadrature_elem_t(xi_t(0.816847572980459, 0.091576213509771), 0.054975871827661));
282  break;
283  case 5:
284  push_back(quadrature_elem_t(xi_t(0.333333333333333, 0.333333333333333), 0.112500000000000));
285  push_back(quadrature_elem_t(xi_t(0.470142064105115, 0.059715871789770), 0.066197076394253));
286  push_back(quadrature_elem_t(xi_t(0.470142064105115, 0.470142064105115), 0.066197076394253));
287  push_back(quadrature_elem_t(xi_t(0.059715871789770, 0.470142064105115), 0.066197076394253));
288  push_back(quadrature_elem_t(xi_t(0.101286507323456, 0.797426985353087), 0.062969590272414));
289  push_back(quadrature_elem_t(xi_t(0.101286507323456, 0.101286507323456), 0.062969590272414));
290  push_back(quadrature_elem_t(xi_t(0.797426985353087, 0.101286507323456), 0.062969590272414));
291  break;
292  case 6:
293  push_back(quadrature_elem_t(xi_t(0.249286745170910, 0.501426509658179), 0.058393137863189));
294  push_back(quadrature_elem_t(xi_t(0.249286745170910, 0.249286745170910), 0.058393137863189));
295  push_back(quadrature_elem_t(xi_t(0.501426509658179, 0.249286745170910), 0.058393137863189));
296  push_back(quadrature_elem_t(xi_t(0.063089014491502, 0.873821971016996), 0.025422453185104));
297  push_back(quadrature_elem_t(xi_t(0.063089014491502, 0.063089014491502), 0.025422453185104));
298  push_back(quadrature_elem_t(xi_t(0.873821971016996, 0.063089014491502), 0.025422453185104));
299  push_back(quadrature_elem_t(xi_t(0.310352451033784, 0.053145049844817), 0.041425537809187));
300  push_back(quadrature_elem_t(xi_t(0.636502499121399, 0.310352451033784), 0.041425537809187));
301  push_back(quadrature_elem_t(xi_t(0.053145049844817, 0.636502499121399), 0.041425537809187));
302  push_back(quadrature_elem_t(xi_t(0.053145049844817, 0.310352451033784), 0.041425537809187));
303  push_back(quadrature_elem_t(xi_t(0.310352451033784, 0.636502499121399), 0.041425537809187));
304  push_back(quadrature_elem_t(xi_t(0.636502499121399, 0.053145049844817), 0.041425537809187));
305  break;
306  case 7:
307  push_back(quadrature_elem_t(xi_t(0.333333333333333, 0.333333333333333), -0.074785022233841));
308  push_back(quadrature_elem_t(xi_t(0.260345966079040, 0.479308067841920), 0.087807628716604));
309  push_back(quadrature_elem_t(xi_t(0.260345966079040, 0.260345966079040), 0.087807628716604));
310  push_back(quadrature_elem_t(xi_t(0.479308067841920, 0.260345966079040), 0.087807628716604));
311  push_back(quadrature_elem_t(xi_t(0.065130102902216, 0.869739794195568), 0.026673617804419));
312  push_back(quadrature_elem_t(xi_t(0.065130102902216, 0.065130102902216), 0.026673617804419));
313  push_back(quadrature_elem_t(xi_t(0.869739794195568, 0.065130102902216), 0.026673617804419));
314  push_back(quadrature_elem_t(xi_t(0.312865496004874, 0.048690315425316), 0.038556880445128));
315  push_back(quadrature_elem_t(xi_t(0.638444188569810, 0.312865496004874), 0.038556880445128));
316  push_back(quadrature_elem_t(xi_t(0.048690315425316, 0.638444188569810), 0.038556880445128));
317  push_back(quadrature_elem_t(xi_t(0.048690315425316, 0.312865496004874), 0.038556880445128));
318  push_back(quadrature_elem_t(xi_t(0.312865496004874, 0.638444188569810), 0.038556880445128));
319  push_back(quadrature_elem_t(xi_t(0.638444188569810, 0.048690315425316), 0.038556880445128));
320  break;
321  case 8:
322  push_back(quadrature_elem_t(xi_t(0.333333333333333 , 0.333333333333333), 0.072157803838894));
323  push_back(quadrature_elem_t(xi_t(0.459292588292723 , 0.081414823414554), 0.047545817133642));
324  push_back(quadrature_elem_t(xi_t(0.459292588292723 , 0.459292588292723), 0.047545817133642));
325  push_back(quadrature_elem_t(xi_t(0.081414823414554 , 0.459292588292723), 0.047545817133642));
326  push_back(quadrature_elem_t(xi_t(0.170569307751760 , 0.658861384496480), 0.051608685267359));
327  push_back(quadrature_elem_t(xi_t(0.170569307751760 , 0.170569307751760), 0.051608685267359));
328  push_back(quadrature_elem_t(xi_t(0.658861384496480 , 0.170569307751760), 0.051608685267359));
329  push_back(quadrature_elem_t(xi_t(0.050547228317031 , 0.898905543365938), 0.016229248811599));
330  push_back(quadrature_elem_t(xi_t(0.050547228317031 , 0.050547228317031), 0.016229248811599));
331  push_back(quadrature_elem_t(xi_t(0.898905543365938 , 0.050547228317031), 0.016229248811599));
332  push_back(quadrature_elem_t(xi_t(0.263112829634638 , 0.008394777409958), 0.013615157087218));
333  push_back(quadrature_elem_t(xi_t(0.728492392955404 , 0.263112829634638), 0.013615157087218));
334  push_back(quadrature_elem_t(xi_t(0.008394777409958 , 0.728492392955404), 0.013615157087218));
335  push_back(quadrature_elem_t(xi_t(0.008394777409958 , 0.263112829634638), 0.013615157087218));
336  push_back(quadrature_elem_t(xi_t(0.263112829634638 , 0.728492392955404), 0.013615157087218));
337  push_back(quadrature_elem_t(xi_t(0.728492392955404 , 0.008394777409958), 0.013615157087218));
338  break;
339  case 9:
340  push_back(quadrature_elem_t(xi_t(0.333333333333333, 0.333333333333333), 0.048567898141400));
341  push_back(quadrature_elem_t(xi_t(0.489682519198738, 0.020634961602525), 0.015667350113570));
342  push_back(quadrature_elem_t(xi_t(0.489682519198738, 0.489682519198738), 0.015667350113570));
343  push_back(quadrature_elem_t(xi_t(0.020634961602525, 0.489682519198738), 0.015667350113570));
344  push_back(quadrature_elem_t(xi_t(0.437089591492937, 0.125820817014127), 0.038913770502387));
345  push_back(quadrature_elem_t(xi_t(0.437089591492937, 0.437089591492937), 0.038913770502387));
346  push_back(quadrature_elem_t(xi_t(0.125820817014127, 0.437089591492937), 0.038913770502387));
347  push_back(quadrature_elem_t(xi_t(0.188203535619033, 0.623592928761935), 0.039823869463605));
348  push_back(quadrature_elem_t(xi_t(0.188203535619033, 0.188203535619033), 0.039823869463605));
349  push_back(quadrature_elem_t(xi_t(0.623592928761935, 0.188203535619033), 0.039823869463605));
350  push_back(quadrature_elem_t(xi_t(0.044729513394453, 0.910540973211095), 0.012788837829349));
351  push_back(quadrature_elem_t(xi_t(0.044729513394453, 0.044729513394453), 0.012788837829349));
352  push_back(quadrature_elem_t(xi_t(0.910540973211095, 0.044729513394453), 0.012788837829349));
353  push_back(quadrature_elem_t(xi_t(0.221962989160766, 0.036838412054736), 0.021641769688645));
354  push_back(quadrature_elem_t(xi_t(0.741198598784498, 0.221962989160766), 0.021641769688645));
355  push_back(quadrature_elem_t(xi_t(0.036838412054736, 0.741198598784498), 0.021641769688645));
356  push_back(quadrature_elem_t(xi_t(0.036838412054736, 0.221962989160766), 0.021641769688645));
357  push_back(quadrature_elem_t(xi_t(0.221962989160766, 0.741198598784498), 0.021641769688645));
358  push_back(quadrature_elem_t(xi_t(0.741198598784498, 0.036838412054736), 0.021641769688645));
359  break;
360  case 10:
361  push_back(quadrature_elem_t(xi_t(0.333333333333333, 0.333333333333333), 0.045408995191377));
362  push_back(quadrature_elem_t(xi_t(0.028844733232685, 0.485577633383657), 0.018362978878233));
363  push_back(quadrature_elem_t(xi_t(0.485577633383657, 0.485577633383657), 0.018362978878233));
364  push_back(quadrature_elem_t(xi_t(0.485577633383657, 0.028844733232685), 0.018362978878233));
365  push_back(quadrature_elem_t(xi_t(0.781036849029926, 0.109481575485037), 0.022660529717764));
366  push_back(quadrature_elem_t(xi_t(0.109481575485037, 0.109481575485037), 0.022660529717764));
367  push_back(quadrature_elem_t(xi_t(0.109481575485037, 0.781036849029926), 0.022660529717764));
368  push_back(quadrature_elem_t(xi_t(0.141707219414880, 0.307939838764121), 0.036378958422710));
369  push_back(quadrature_elem_t(xi_t(0.307939838764121, 0.550352941820999), 0.036378958422710));
370  push_back(quadrature_elem_t(xi_t(0.550352941820999, 0.141707219414880), 0.036378958422710));
371  push_back(quadrature_elem_t(xi_t(0.307939838764121, 0.141707219414880), 0.036378958422710));
372  push_back(quadrature_elem_t(xi_t(0.550352941820999, 0.307939838764121), 0.036378958422710));
373  push_back(quadrature_elem_t(xi_t(0.141707219414880, 0.550352941820999), 0.036378958422710));
374  push_back(quadrature_elem_t(xi_t(0.025003534762686, 0.246672560639903), 0.014163621265528));
375  push_back(quadrature_elem_t(xi_t(0.246672560639903, 0.728323904597411), 0.014163621265528));
376  push_back(quadrature_elem_t(xi_t(0.728323904597411, 0.025003534762686), 0.014163621265528));
377  push_back(quadrature_elem_t(xi_t(0.246672560639903, 0.025003534762686), 0.014163621265528));
378  push_back(quadrature_elem_t(xi_t(0.728323904597411, 0.246672560639903), 0.014163621265528));
379  push_back(quadrature_elem_t(xi_t(0.025003534762686, 0.728323904597411), 0.014163621265528));
380  push_back(quadrature_elem_t(xi_t(0.009540815400299, 0.066803251012200), 0.004710833481867));
381  push_back(quadrature_elem_t(xi_t(0.066803251012200, 0.923655933587500), 0.004710833481867));
382  push_back(quadrature_elem_t(xi_t(0.923655933587500, 0.009540815400299), 0.004710833481867));
383  push_back(quadrature_elem_t(xi_t(0.066803251012200, 0.009540815400299), 0.004710833481867));
384  push_back(quadrature_elem_t(xi_t(0.923655933587500, 0.066803251012200), 0.004710833481867));
385  push_back(quadrature_elem_t(xi_t(0.009540815400299, 0.923655933587500), 0.004710833481867));
386  break;
387  case 11: // NOTE: 11 same as 12 because of outliers at order 11
388  case 12:
389  push_back(quadrature_elem_t(xi_t(0.023565220452390, 0.488217389773805), 0.012865533220227));
390  push_back(quadrature_elem_t(xi_t(0.488217389773805, 0.488217389773805), 0.012865533220227));
391  push_back(quadrature_elem_t(xi_t(0.488217389773805, 0.023565220452390), 0.012865533220227));
392  push_back(quadrature_elem_t(xi_t(0.120551215411079, 0.439724392294460), 0.021846272269019));
393  push_back(quadrature_elem_t(xi_t(0.439724392294460, 0.439724392294460), 0.021846272269019));
394  push_back(quadrature_elem_t(xi_t(0.439724392294460, 0.120551215411079), 0.021846272269019));
395  push_back(quadrature_elem_t(xi_t(0.457579229975768, 0.271210385012116), 0.031429112108943));
396  push_back(quadrature_elem_t(xi_t(0.271210385012116, 0.271210385012116), 0.031429112108943));
397  push_back(quadrature_elem_t(xi_t(0.271210385012116, 0.457579229975768), 0.031429112108943));
398  push_back(quadrature_elem_t(xi_t(0.744847708916828, 0.127576145541586), 0.017398056465355));
399  push_back(quadrature_elem_t(xi_t(0.127576145541586, 0.127576145541586), 0.017398056465355));
400  push_back(quadrature_elem_t(xi_t(0.127576145541586, 0.744847708916828), 0.017398056465355));
401  push_back(quadrature_elem_t(xi_t(0.957365299093579, 0.021317350453210), 0.003083130525780));
402  push_back(quadrature_elem_t(xi_t(0.021317350453210, 0.021317350453210), 0.003083130525780));
403  push_back(quadrature_elem_t(xi_t(0.021317350453210, 0.957365299093579), 0.003083130525780));
404  push_back(quadrature_elem_t(xi_t(0.115343494534698, 0.275713269685514), 0.020185778883191));
405  push_back(quadrature_elem_t(xi_t(0.275713269685514, 0.608943235779788), 0.020185778883191));
406  push_back(quadrature_elem_t(xi_t(0.608943235779788, 0.115343494534698), 0.020185778883191));
407  push_back(quadrature_elem_t(xi_t(0.275713269685514, 0.115343494534698), 0.020185778883191));
408  push_back(quadrature_elem_t(xi_t(0.608943235779788, 0.275713269685514), 0.020185778883191));
409  push_back(quadrature_elem_t(xi_t(0.115343494534698, 0.608943235779788), 0.020185778883191));
410  push_back(quadrature_elem_t(xi_t(0.022838332222257, 0.281325580989940), 0.011178386601152));
411  push_back(quadrature_elem_t(xi_t(0.281325580989940, 0.695836086787803), 0.011178386601152));
412  push_back(quadrature_elem_t(xi_t(0.695836086787803, 0.022838332222257), 0.011178386601152));
413  push_back(quadrature_elem_t(xi_t(0.281325580989940, 0.022838332222257), 0.011178386601152));
414  push_back(quadrature_elem_t(xi_t(0.695836086787803, 0.281325580989940), 0.011178386601152));
415  push_back(quadrature_elem_t(xi_t(0.022838332222257, 0.695836086787803), 0.011178386601152));
416  push_back(quadrature_elem_t(xi_t(0.025734050548330, 0.116251915907597), 0.008658115554329));
417  push_back(quadrature_elem_t(xi_t(0.116251915907597, 0.858014033544073), 0.008658115554329));
418  push_back(quadrature_elem_t(xi_t(0.858014033544073, 0.025734050548330), 0.008658115554329));
419  push_back(quadrature_elem_t(xi_t(0.116251915907597, 0.025734050548330), 0.008658115554329));
420  push_back(quadrature_elem_t(xi_t(0.858014033544073, 0.116251915907597), 0.008658115554329));
421  push_back(quadrature_elem_t(xi_t(0.025734050548330, 0.858014033544073), 0.008658115554329));
422  break;
423  case 13:
424  push_back(quadrature_elem_t(xi_t(0.333333333333333, 0.333333333333333), 0.026260461700401));
425  push_back(quadrature_elem_t(xi_t(0.009903630120591, 0.495048184939705), 0.005640072604665));
426  push_back(quadrature_elem_t(xi_t(0.495048184939705, 0.495048184939705), 0.005640072604665));
427  push_back(quadrature_elem_t(xi_t(0.495048184939705, 0.009903630120591), 0.005640072604665));
428  push_back(quadrature_elem_t(xi_t(0.062566729780852, 0.468716635109574), 0.015711759181227));
429  push_back(quadrature_elem_t(xi_t(0.468716635109574, 0.468716635109574), 0.015711759181227));
430  push_back(quadrature_elem_t(xi_t(0.468716635109574, 0.062566729780852), 0.015711759181227));
431  push_back(quadrature_elem_t(xi_t(0.170957326397447, 0.414521336801277), 0.023536251252097));
432  push_back(quadrature_elem_t(xi_t(0.414521336801277, 0.414521336801277), 0.023536251252097));
433  push_back(quadrature_elem_t(xi_t(0.414521336801277, 0.170957326397447), 0.023536251252097));
434  push_back(quadrature_elem_t(xi_t(0.541200855914337, 0.229399572042831), 0.023681793268178));
435  push_back(quadrature_elem_t(xi_t(0.229399572042831, 0.229399572042831), 0.023681793268178));
436  push_back(quadrature_elem_t(xi_t(0.229399572042831, 0.541200855914337), 0.023681793268178));
437  push_back(quadrature_elem_t(xi_t(0.771151009607340, 0.114424495196330), 0.015583764522897));
438  push_back(quadrature_elem_t(xi_t(0.114424495196330, 0.114424495196330), 0.015583764522897));
439  push_back(quadrature_elem_t(xi_t(0.114424495196330, 0.771151009607340), 0.015583764522897));
440  push_back(quadrature_elem_t(xi_t(0.950377217273082, 0.024811391363459), 0.003987885732537));
441  push_back(quadrature_elem_t(xi_t(0.024811391363459, 0.024811391363459), 0.003987885732537));
442  push_back(quadrature_elem_t(xi_t(0.024811391363459, 0.950377217273082), 0.003987885732537));
443  push_back(quadrature_elem_t(xi_t(0.094853828379579, 0.268794997058761), 0.018424201364366));
444  push_back(quadrature_elem_t(xi_t(0.268794997058761, 0.636351174561660), 0.018424201364366));
445  push_back(quadrature_elem_t(xi_t(0.636351174561660, 0.094853828379579), 0.018424201364366));
446  push_back(quadrature_elem_t(xi_t(0.268794997058761, 0.094853828379579), 0.018424201364366));
447  push_back(quadrature_elem_t(xi_t(0.636351174561660, 0.268794997058761), 0.018424201364366));
448  push_back(quadrature_elem_t(xi_t(0.094853828379579, 0.636351174561660), 0.018424201364366));
449  push_back(quadrature_elem_t(xi_t(0.018100773278807, 0.291730066734288), 0.008700731651911));
450  push_back(quadrature_elem_t(xi_t(0.291730066734288, 0.690169159986905), 0.008700731651911));
451  push_back(quadrature_elem_t(xi_t(0.690169159986905, 0.018100773278807), 0.008700731651911));
452  push_back(quadrature_elem_t(xi_t(0.291730066734288, 0.018100773278807), 0.008700731651911));
453  push_back(quadrature_elem_t(xi_t(0.690169159986905, 0.291730066734288), 0.008700731651911));
454  push_back(quadrature_elem_t(xi_t(0.018100773278807, 0.690169159986905), 0.008700731651911));
455  push_back(quadrature_elem_t(xi_t(0.022233076674090, 0.126357385491669), 0.007760893419522));
456  push_back(quadrature_elem_t(xi_t(0.126357385491669, 0.851409537834241), 0.007760893419522));
457  push_back(quadrature_elem_t(xi_t(0.851409537834241, 0.022233076674090), 0.007760893419522));
458  push_back(quadrature_elem_t(xi_t(0.126357385491669, 0.022233076674090), 0.007760893419522));
459  push_back(quadrature_elem_t(xi_t(0.851409537834241, 0.126357385491669), 0.007760893419522));
460  push_back(quadrature_elem_t(xi_t(0.022233076674090, 0.851409537834241), 0.007760893419522));
461  break;
462  case 14:
463  push_back(quadrature_elem_t(xi_t(0.022072179275643, 0.488963910362179), 0.010941790684715));
464  push_back(quadrature_elem_t(xi_t(0.488963910362179, 0.488963910362179), 0.010941790684715));
465  push_back(quadrature_elem_t(xi_t(0.488963910362179, 0.022072179275643), 0.010941790684715));
466  push_back(quadrature_elem_t(xi_t(0.164710561319092, 0.417644719340454), 0.016394176772063));
467  push_back(quadrature_elem_t(xi_t(0.417644719340454, 0.417644719340454), 0.016394176772063));
468  push_back(quadrature_elem_t(xi_t(0.417644719340454, 0.164710561319092), 0.016394176772063));
469  push_back(quadrature_elem_t(xi_t(0.453044943382323, 0.273477528308839), 0.025887052253646));
470  push_back(quadrature_elem_t(xi_t(0.273477528308839, 0.273477528308839), 0.025887052253646));
471  push_back(quadrature_elem_t(xi_t(0.273477528308839, 0.453044943382323), 0.025887052253646));
472  push_back(quadrature_elem_t(xi_t(0.645588935174913, 0.177205532412543), 0.021081294368497));
473  push_back(quadrature_elem_t(xi_t(0.177205532412543, 0.177205532412543), 0.021081294368497));
474  push_back(quadrature_elem_t(xi_t(0.177205532412543, 0.645588935174913), 0.021081294368497));
475  push_back(quadrature_elem_t(xi_t(0.876400233818255, 0.061799883090873), 0.007216849834889));
476  push_back(quadrature_elem_t(xi_t(0.061799883090873, 0.061799883090873), 0.007216849834889));
477  push_back(quadrature_elem_t(xi_t(0.061799883090873, 0.876400233818255), 0.007216849834889));
478  push_back(quadrature_elem_t(xi_t(0.961218077502598, 0.019390961248701), 0.002461701801200));
479  push_back(quadrature_elem_t(xi_t(0.019390961248701, 0.019390961248701), 0.002461701801200));
480  push_back(quadrature_elem_t(xi_t(0.019390961248701, 0.961218077502598), 0.002461701801200));
481  push_back(quadrature_elem_t(xi_t(0.057124757403648, 0.172266687821356), 0.012332876606282));
482  push_back(quadrature_elem_t(xi_t(0.172266687821356, 0.770608554774996), 0.012332876606282));
483  push_back(quadrature_elem_t(xi_t(0.770608554774996, 0.057124757403648), 0.012332876606282));
484  push_back(quadrature_elem_t(xi_t(0.172266687821356, 0.057124757403648), 0.012332876606282));
485  push_back(quadrature_elem_t(xi_t(0.770608554774996, 0.172266687821356), 0.012332876606282));
486  push_back(quadrature_elem_t(xi_t(0.057124757403648, 0.770608554774996), 0.012332876606282));
487  push_back(quadrature_elem_t(xi_t(0.092916249356972, 0.336861459796345), 0.019285755393531));
488  push_back(quadrature_elem_t(xi_t(0.336861459796345, 0.570222290846683), 0.019285755393531));
489  push_back(quadrature_elem_t(xi_t(0.570222290846683, 0.092916249356972), 0.019285755393531));
490  push_back(quadrature_elem_t(xi_t(0.336861459796345, 0.092916249356972), 0.019285755393531));
491  push_back(quadrature_elem_t(xi_t(0.570222290846683, 0.336861459796345), 0.019285755393531));
492  push_back(quadrature_elem_t(xi_t(0.092916249356972, 0.570222290846683), 0.019285755393531));
493  push_back(quadrature_elem_t(xi_t(0.014646950055654, 0.298372882136258), 0.007218154056767));
494  push_back(quadrature_elem_t(xi_t(0.298372882136258, 0.686980167808088), 0.007218154056767));
495  push_back(quadrature_elem_t(xi_t(0.686980167808088, 0.014646950055654), 0.007218154056767));
496  push_back(quadrature_elem_t(xi_t(0.298372882136258, 0.014646950055654), 0.007218154056767));
497  push_back(quadrature_elem_t(xi_t(0.686980167808088, 0.298372882136258), 0.007218154056767));
498  push_back(quadrature_elem_t(xi_t(0.014646950055654, 0.686980167808088), 0.007218154056767));
499  push_back(quadrature_elem_t(xi_t(0.001268330932872, 0.118974497696957), 0.002505114419250));
500  push_back(quadrature_elem_t(xi_t(0.118974497696957, 0.879757171370171), 0.002505114419250));
501  push_back(quadrature_elem_t(xi_t(0.879757171370171, 0.001268330932872), 0.002505114419250));
502  push_back(quadrature_elem_t(xi_t(0.118974497696957, 0.001268330932872), 0.002505114419250));
503  push_back(quadrature_elem_t(xi_t(0.879757171370171, 0.118974497696957), 0.002505114419250));
504  push_back(quadrature_elem_t(xi_t(0.001268330932872, 0.879757171370171), 0.002505114419250));
505  break;
506  default:
507  throw std::out_of_range("unsupported dunavant degree: " + std::to_string(degree) +
508  " (max degree: " + std::to_string(dunavant_num.size()-1) + ")");
509  break;
510  }
511  }
512 };
513 
517 struct gauss_family_tag { typedef gauss_family_tag type; };
518 
520 template<class Domain>
521 struct quadrature_type<gauss_family_tag, Domain> :
522  gaussian_quadrature<Domain>
523 {
524 };
525 
526 
528 class log_gaussian_quadrature;
529 
531 template <>
532 struct quadrature_traits<log_gaussian_quadrature>
533 {
536  typedef NiHu::line_domain domain_t;
537 };
538 
540 class log_gaussian_quadrature :
541  public quadrature_base<log_gaussian_quadrature>
542 {
543 public:
545  typedef quadrature_base<log_gaussian_quadrature> base_t;
547  typedef base_t::xi_t xi_t;
549  typedef base_t::scalar_t scalar_t;
550 
552  typedef log_gaussian_quadrature type;
553 
558  static Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> log_gauss_impl(size_t N)
559  {
560  Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> V(N, 2);
561 
562  switch (N)
563  {
564  case 1:
565  V <<
566  2.5000000000000000e-01, 1.0000000000000000e+00;
567  break;
568  case 2:
569  V <<
570  1.1200880616697619e-01, 7.1853931903038448e-01,
571  6.0227690811873813e-01, 2.8146068096961557e-01;
572  break;
573  case 3:
574  V <<
575  6.3890793087325398e-02, 5.1340455223236336e-01,
576  3.6899706371561874e-01, 3.9198004120148755e-01,
577  7.6688030393894147e-01, 9.4615406566149127e-02;
578  break;
579  case 4:
580  V <<
581  4.1448480199383221e-02, 3.8346406814513512e-01,
582  2.4527491432060225e-01, 3.8687531777476264e-01,
583  5.5616545356027580e-01, 1.9043512695014242e-01,
584  8.4898239453298519e-01, 3.9225487129959831e-02;
585  break;
586  case 5:
587  V <<
588  2.9134472151972100e-02, 2.9789347178289399e-01,
589  1.7397721332089799e-01, 3.4977622651322399e-01,
590  4.1170252028490201e-01, 2.3448829004405200e-01,
591  6.7731417458281995e-01, 9.8930459516633096e-02,
592  8.9477136103100796e-01, 1.8911552143195801e-02;
593  break;
594  case 6:
595  case 7:
596  case 8:
597  case 9:
598  V.resize(10, 2);
599  case 10:
600  V <<
601  9.0426309621996510e-03, 1.2095513195457051e-01,
602  5.3971266222500633e-02, 1.8636354256407187e-01,
603  1.3531182463925079e-01, 1.9566087327775999e-01,
604  2.4705241628715982e-01, 1.7357714218290693e-01,
605  3.8021253960933232e-01, 1.3569567299548421e-01,
606  5.2379231797184322e-01, 9.3646758538110525e-02,
607  6.6577520551642455e-01, 5.5787727351415871e-02,
608  7.9419041601196627e-01, 2.7159810899233330e-02,
609  8.9816109121900356e-01, 9.5151826028485147e-03,
610  9.6884798871863353e-01, 1.6381576335982632e-03;
611  break;
612  case 11:
613  case 12:
614  case 13:
615  case 14:
616  V.resize(15, 2);
617  case 15:
618  V <<
619  4.3831101754754041e-03, 6.7009978916493712e-02,
620  2.5935898105330615e-02, 1.1226415028670574e-01,
621  6.5596095412316244e-02, 1.3176017703967990e-01,
622  1.2210193407333160e-01, 1.3521764906193473e-01,
623  1.9339526237400712e-01, 1.2788179864568036e-01,
624  2.7677283870610203e-01, 1.1353290749021942e-01,
625  3.6901512713974294e-01, 9.5205239784358658e-02,
626  4.6652432896470658e-01, 7.5389314167395957e-02,
627  5.6547347379181734e-01, 5.6078424492653718e-02,
628  6.6196291901245641e-01, 3.8768295375018233e-02,
629  7.5217888337878580e-01, 2.4451483268750077e-02,
630  8.3254803386618959e-01, 1.3624630138238846e-02,
631  8.9988205012089806e-01, 6.3164475985907614e-03,
632  9.5150618874340986e-01, 2.1388899159444715e-03,
633  9.8536446812213196e-01, 3.6061381833540662e-04;
634  break;
635  case 16:
636  case 17:
637  case 18:
638  case 19:
639  V.resize(20, 2);
640  case 20:
641  V <<
642  2.5883279559219554e-03, 4.3142752133208076e-02,
643  1.5209662349560232e-02, 7.5383709908589364e-02,
644  3.8536550372165329e-02, 9.3053267451663049e-02,
645  7.2181613815873902e-02, 1.0145671184982975e-01,
646  1.1546052648763315e-01, 1.0320176205607207e-01,
647  1.6744285627532968e-01, 1.0002254980527317e-01,
648  2.2698378726020249e-01, 9.3259799300297680e-02,
649  2.9275496094154585e-01, 8.4028952871941051e-02,
650  3.6327742985785888e-01, 7.3285589130030734e-02,
651  4.3695714009076830e-01, 6.1850336913730292e-02,
652  5.1212259467896737e-01, 5.0416604438374681e-02,
653  5.8706404491440989e-01, 3.9551370005298382e-02,
654  6.6007341331490943e-01, 2.9694077895812843e-02,
655  7.2948408392968744e-01, 2.1156315355427099e-02,
656  7.9370967198708586e-01, 1.4123732938964021e-02,
657  8.5128089278912578e-01, 8.6609745043354988e-03,
658  9.0087968085441761e-01, 4.7199401462036054e-03,
659  9.4136974912909166e-01, 2.1513974039652061e-03,
660  9.7182274107526323e-01, 7.1972821465320265e-04,
661  9.9153808143871203e-01, 1.2042767633021674e-04;
662  break;
663  default:
664  throw std::out_of_range("unsupported log Gauss degree");
665  }
666  return V;
667  }
668 
672  log_gaussian_quadrature() :
673  base_t(0)
674  {
675  }
676 
681  log_gaussian_quadrature(size_t degree) :
682  base_t(degree / 2 + 1)
683  {
684  size_t N = degree / 2 + 1;
685  // compute 1D Gaussian locations and weights
686  auto V = log_gauss_impl(N);
687 
688  // Fill the points and weights
689  for (size_t i = 0; i < V.rows(); ++i)
690  push_back(quadrature_elem_t(xi_t(V(i, 0)), V(i, 1)));
691  }
692 }; // class log_gaussian_quadrature
693 
694 } // end of namespace NiHu
695 
696 #endif // GAUSSIAN_QUADRATURE_HPP_INCLUDED
NiHu::bessel::J
T J(T const &z)
Bessel function J_nu(z)
Definition: math_functions.hpp:223
NiHu::tria_domain
a 2D triangle domain
Definition: lib_domain.hpp:101
NiHu::gaussian_quadrature
Definition: gaussian_quadrature.hpp:89
NiHu::quadrature_traits
Definition: quadrature.hpp:145
NiHu::gaussian_quadrature< line_domain >::xi_t
base_t::xi_t xi_t
the location type
Definition: gaussian_quadrature.hpp:114
NiHu::gaussian_quadrature< quad_domain >::type
gaussian_quadrature type
self-returning
Definition: gaussian_quadrature.hpp:169
NiHu::log_gaussian_quadrature::scalar_t
base_t::scalar_t scalar_t
the scalar type
Definition: gaussian_quadrature.hpp:549
quadrature.hpp
implementation of class NiHu::quadrature_elem, NiHu::quadrature_base
NiHu::gaussian_quadrature< line_domain >::scalar_t
base_t::scalar_t scalar_t
the scalar type
Definition: gaussian_quadrature.hpp:116
NiHu::gaussian_quadrature< quad_domain >::gaussian_quadrature
gaussian_quadrature()
default constructor creating an empty quadrature
Definition: gaussian_quadrature.hpp:174
NiHu::quadrature_base< gaussian_quadrature< line_domain > >::scalar_t
domain_t::scalar_t scalar_t
local scalar type
Definition: quadrature.hpp:179
NiHu::gaussian_quadrature< quad_domain >::domain_t
base_t::domain_t domain_t
the domain type
Definition: gaussian_quadrature.hpp:162
NiHu::gaussian_quadrature< line_domain >::type
gaussian_quadrature type
self-returning
Definition: gaussian_quadrature.hpp:119
NiHu::gaussian_quadrature< quad_domain >::xi_t
base_t::xi_t xi_t
the location type
Definition: gaussian_quadrature.hpp:164
NiHu::quadrature_base< gaussian_quadrature< line_domain > >::xi_t
domain_t::xi_t xi_t
local coordinate type
Definition: quadrature.hpp:177
NiHu::quadrature_base< gaussian_quadrature< quad_domain > >::domain_t
traits_t::domain_t domain_t
domain type
Definition: quadrature.hpp:175
NiHu::quadrature_traits< log_gaussian_quadrature >::domain_t
NiHu::line_domain domain_t
type of the domain
Definition: gaussian_quadrature.hpp:536
NiHu::gaussian_quadrature< quad_domain >::scalar_t
base_t::scalar_t scalar_t
the scalar type
Definition: gaussian_quadrature.hpp:166
NiHu::gaussian_quadrature< tria_domain >::xi_t
base_t::xi_t xi_t
the quadrature location type
Definition: gaussian_quadrature.hpp:237
NiHu::quadrature_base
CRTP base class of all quadratures.
Definition: quadrature.hpp:162
NiHu::gaussian_quadrature< quad_domain >::base_t
quadrature_base< gaussian_quadrature< quad_domain > > base_t
the base class
Definition: gaussian_quadrature.hpp:160
NiHu::log_gaussian_quadrature::xi_t
base_t::xi_t xi_t
the location type
Definition: gaussian_quadrature.hpp:547
NiHu::fmm::fmm_matrix
Matrix representation of the FMM method.
Definition: fmm_matrix.hpp:69
NiHu::line_domain
a 1D line domain
Definition: lib_domain.hpp:53
NiHu::quadrature_elem
a quadrature element is a base point and a weight
Definition: quadrature.hpp:39
NiHu::log_gaussian_quadrature::log_gauss_impl
static Eigen::Matrix< scalar_t, Eigen::Dynamic, 2 > log_gauss_impl(size_t N)
return 1D N-point log Gaussian quadrature
Definition: gaussian_quadrature.hpp:558
NiHu::gauss_impl
Eigen::Matrix< scalar_t, Eigen::Dynamic, 2 > gauss_impl(size_t N)
return 1D N-point Guassian quadrature
Definition: gaussian_quadrature.hpp:62
NiHu::quadrature_type
metafunction to assign a quadrature type to a quadrature family and a domain
Definition: quadrature.hpp:326
NiHu::quadrature_traits< gaussian_quadrature< Domain > >::domain_t
Domain domain_t
type of the domain
Definition: gaussian_quadrature.hpp:100
NiHu::quadrature_base< log_gaussian_quadrature >::quadrature_elem_t
quadrature_elem< xi_t, scalar_t > quadrature_elem_t
quadrature elem type
Definition: quadrature.hpp:181
NiHu::gaussian_quadrature< line_domain >::gaussian_quadrature
gaussian_quadrature(size_t degree)
constructor for a given polynomial degree
Definition: gaussian_quadrature.hpp:133
NiHu::gaussian_quadrature< tria_domain >::type
gaussian_quadrature type
self-returning
Definition: gaussian_quadrature.hpp:240
NiHu::log_gaussian_quadrature::log_gaussian_quadrature
log_gaussian_quadrature(size_t degree)
constructor for a given polynomial degree
Definition: gaussian_quadrature.hpp:681
NiHu::log_gaussian_quadrature::type
log_gaussian_quadrature type
self-returning
Definition: gaussian_quadrature.hpp:552
NiHu::gaussian_quadrature< tria_domain >::quadrature_elem_t
base_t::quadrature_elem_t quadrature_elem_t
the quadrature elem
Definition: gaussian_quadrature.hpp:235
NiHu::gaussian_quadrature< tria_domain >::gaussian_quadrature
gaussian_quadrature()
default constructor creating an empty quadrature
Definition: gaussian_quadrature.hpp:246
NiHu::gaussian_quadrature< tria_domain >::gaussian_quadrature
gaussian_quadrature(size_t degree)
constructor for a given polynomial order
Definition: gaussian_quadrature.hpp:255
NiHu::gaussian_quadrature< tria_domain >::base_t
quadrature_base< gaussian_quadrature< tria_domain > > base_t
base class
Definition: gaussian_quadrature.hpp:233
NiHu::gaussian_quadrature< line_domain >::base_t
quadrature_base< gaussian_quadrature< line_domain > > base_t
the base class
Definition: gaussian_quadrature.hpp:112
NiHu::log_gaussian_quadrature::base_t
quadrature_base< log_gaussian_quadrature > base_t
the base class
Definition: gaussian_quadrature.hpp:545
NiHu::log_gaussian_quadrature::log_gaussian_quadrature
log_gaussian_quadrature()
default constructor creating an empty quadrature
Definition: gaussian_quadrature.hpp:672
NiHu::gaussian_quadrature< line_domain >::gaussian_quadrature
gaussian_quadrature()
default constructor creating an empty quadrature
Definition: gaussian_quadrature.hpp:124
NiHu::gaussian_quadrature< quad_domain >::gaussian_quadrature
gaussian_quadrature(size_t degree)
constructor for a given polynomial order
Definition: gaussian_quadrature.hpp:183
NiHu::quad_domain
a 2D quad domain
Definition: lib_domain.hpp:157
NiHu::log_gaussian_quadrature
Log-Gaussian quadrature over a line domain.
Definition: gaussian_quadrature.hpp:540
NiHu::gauss_family_tag
tag for the family of Gaussian quadratures
Definition: gaussian_quadrature.hpp:517