Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_e.cc
Go to the documentation of this file.
1
21#include "rheolef/quadrature.h"
22#include "rheolef/gauss_jacobi.h"
23namespace rheolef {
24using 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
36template<class T>
37void
38quadrature_on_geo<T>::init_edge (quadrature_option opt)
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) \
120template void quadrature_on_geo<T>::init_edge (quadrature_option);
121
123
124}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
base::size_type size_type
Definition quadrature.h:88
void init_edge(quadrature_option opt)
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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)
STL namespace.
Definition cavity_dg.h:29