Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_t.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 * 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 */
48template<class T>
49void
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) \
232template void quadrature_on_geo<T>::init_triangle (quadrature_option);
233
235
236}// 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_triangle(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