Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_q.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 293
37 * 2) M. Crouzeix et A. L. Mignot
38 * Exercices d'analyse des equations differentielles,
39 * Masson,
40 * 1986,
41 * p. 63
42 * ------------------------------------------
43 */
44template<class T>
45void
46quadrature_on_geo<T>::init_square (quadrature_option opt)
47{
48 quadrature_option::family_type f = opt.get_family();
49 // -------------------------------------------------------------------------
50 // special case : equispaced, for irregular (e.g. Heaviside) functions
51 // -------------------------------------------------------------------------
52 if (f == quadrature_option::equispaced) {
53 size_type r = opt.get_order();
54 if (r == 0) {
55 wx (x(0,0), 4);
56 } else {
57 size_type n = (r+1)*(r+1);
58 T w = 4/T(int(n));
59 for (size_type i = 0; i <= r; i++) {
60 for (size_type j = 0; j <= r; j++) {
61 wx (x(2*T(int(i))/r-1,2*T(int(j))/r-1), w);
62 }
63 }
64 }
65 return;
66 }
67 // -------------------------------------------------------------------------
68 // TODO: special case : superconvergent patch recovery nodes & weights
69 // -------------------------------------------------------------------------
70
71 // -------------------------------------------------------------------------
72 // default: gauss
73 // -------------------------------------------------------------------------
74 check_macro (f == quadrature_option::gauss,
75 "unsupported quadrature family \"" << opt.get_family_name() << "\"");
76
77 // Gauss-Legendre quadrature formulae
78 // where Legendre = Jacobi(alpha=0,beta=0) polynoms
79 // when using n nodes : quadrature formulae order is 2*r-1
80 size_type n = n_node_gauss(opt.get_order());
81 vector<T> zeta(n), omega(n);
82 gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
83 for (size_type i = 0; i < n; i++)
84 for (size_type j = 0; j < n; j++)
85 wx (x(zeta[i], zeta[j]), omega[i]*omega[j]);
86}
87// ----------------------------------------------------------------------------
88// instanciation in library
89// ----------------------------------------------------------------------------
90#define _RHEOLEF_instanciation(T) \
91template void quadrature_on_geo<T>::init_square (quadrature_option);
92
94
95}// 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_square(quadrature_option opt)
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