2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
6
/// Rheolef 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 2 of the License, or
9
/// (at your option) any later version.
11
/// Rheolef 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.
16
/// You should have received a copy of the GNU General Public License
17
/// along with Rheolef; if not, write to the Free Software
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
/// =========================================================================
21
#include "rheolef/quadrature.h"
22
#include "rheolef/gauss_jacobi.h"
23
using namespace rheolef;
26
/* ------------------------------------------
27
* Gauss quadrature formulae on the triangle
28
* order r : exact for all polynoms of P_r
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
37
* 2) M. Crouzeix et A. L. Mignot
38
* Exercices d'analyse des equations differentielles,
42
* PROBLEM: direct optimized formulae have hard-coded
43
* coefficients in double precision.
44
* TODO: compute numeric value at the fly why full
46
* ------------------------------------------
50
quadrature_on_geo::init_triangle (quadrature_option_type opt)
52
// -------------------------------------------------------------------------
53
// special case : superconvergent patch recovery nodes & weights
54
// -------------------------------------------------------------------------
55
quadrature_option_type::family_type f = opt.get_family();
56
if (f == quadrature_option_type::superconvergent) {
57
switch (opt.get_order()) {
60
f = quadrature_option_type::gauss;
63
f = quadrature_option_type::middle_edge;
66
error_macro ("unsupported superconvergent("<<opt.get_order()<<")");
69
// -------------------------------------------------------------------------
70
// special case : Middle Edge formulae
71
// e.g. when using special option in riesz_representer
72
// -------------------------------------------------------------------------
73
if (f == quadrature_option_type::middle_edge) {
74
switch (opt.get_order()) {
78
// middle edge formulae:
79
wx(x(Float(0.5), Float(0.0)), 1./6);
80
wx(x(Float(0.0), Float(0.5)), 1./6);
81
wx(x(Float(0.5), Float(0.5)), 1./6);
84
error_macro ("unsupported Middle-Edge("<<opt.get_order()<<")");
88
// -------------------------------------------------------------------------
89
// special case : Gauss-Lobatto points
90
// e.g. when using special option in riesz_representer
91
// -------------------------------------------------------------------------
92
if (f == quadrature_option_type::gauss_lobatto) {
93
switch (opt.get_order()) {
97
wx(x(Float(0), Float(0)), 1./6);
98
wx(x(Float(1), Float(0)), 1./6);
99
wx(x(Float(0), Float(1)), 1./6);
103
wx(x(Float(0), Float(0)), 3/Float(120));
104
wx(x(Float(1), Float(0)), 3/Float(120));
105
wx(x(Float(0), Float(1)), 3/Float(120));
106
wx(x(Float(0.5), Float(0.0)), 8/Float(120));
107
wx(x(Float(0.0), Float(0.5)), 8/Float(120));
108
wx(x(Float(0.5), Float(0.5)), 8/Float(120));
109
wx(x(1/Float(3), 1/Float(3)), 27/Float(120));
112
error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
116
// -------------------------------------------------------------------------
117
// default & general case : Gauss points
118
// -------------------------------------------------------------------------
119
assert_macro (f == quadrature_option_type::gauss,
120
"unsupported quadrature family \"" << opt.get_family_name() << "\"");
122
switch (opt.get_order()) {
125
// better than the general case:
129
wx(x(1/Float(3), 1/Float(3)), 1./2);
132
// better than the general case:
135
wx(x(1/Float(6), 1/Float(6)), 1/Float(6));
136
wx(x(4/Float(6), 1/Float(6)), 1/Float(6));
137
wx(x(1/Float(6), 4/Float(6)), 1/Float(6));
139
#if !(defined(_RHEOLEF_HAVE_CLN) || defined(_RHEOLEF_HAVE_DOUBLEDOUBLE) || defined(_RHEOLEF_HAVE_LONG_DOUBLE))
140
// numerical values are inlined in double precision : 15 digits !
143
// better than the general case:
147
Float a = 0.445948490915965;
148
Float b = 0.091576213509771;
149
Float wa = 0.111690794839005;
150
Float wb = 0.054975871827661;
161
// better than the general case:
165
Float a = 0.063089014491502;
166
Float b = 0.249286745170910;
167
Float c = 0.310352451033785;
168
Float d = 0.053145049844816;
169
Float wa = 0.025422453185103;
170
Float wb = 0.058393137863189;
171
Float wc = 0.041425537809187;
180
wx(x(1-(c+d), c), wc);
181
wx(x(1-(c+d), d), wc);
182
wx(x(c, 1-(c+d)), wc);
183
wx(x(d, 1-(c+d)), wc);
188
// Gauss-Legendre quadrature formulae
189
// where Legendre = Jacobi(alpha=0,beta=0) polynoms
190
size_t r = opt.get_order();
191
size_type n0 = n_node_gauss(r);
192
size_type n1 = n_node_gauss(r+1);
193
vector<Float> zeta0(n0), omega0(n0);
194
vector<Float> zeta1(n1), omega1(n1);
195
gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
196
gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
198
for (size_type i = 0; i < n0; i++) {
199
for (size_type j = 0; j < n1; j++) {
200
// we transform the square into the triangle
201
Float eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
202
Float eta_1 = (1+zeta1[j])/2;
203
Float J = (1-zeta1[j])/8;
204
wx (x(eta_0,eta_1), J*omega0[i]*omega1[j]);
added
removed
nfem/basis/quadrature_t.cc
