Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_H.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
26template<class T>
27void
29{
30 quadrature_option::family_type f = opt.get_family();
31 // -------------------------------------------------------------------------
32 // special case : equispaced, for irregular (e.g. Heaviside) functions
33 // -------------------------------------------------------------------------
34 if (f == quadrature_option::equispaced) {
35 size_type r = opt.get_order();
36 if (r == 0) {
37 wx (x(0,0,0), 8);
38 } else {
39 size_type n = (r+1)*(r+1)*(r+1);
40 T w = 8/T(int(n));
41 for (size_type i = 0; i <= r; i++) {
42 for (size_type j = 0; j <= r; j++) {
43 for (size_type k = 0; k <= r; k++) {
44 wx (x(2*T(int(i))/r-1,2*T(int(j))/r-1,2*T(int(k))/r-1), w);
45 }
46 }
47 }
48 }
49 return;
50 }
51 // -------------------------------------------------------------------------
52 // TODO: special case : superconvergent patch recovery nodes & weights
53 // -------------------------------------------------------------------------
54
55 // -------------------------------------------------------------------------
56 // default: gauss
57 // -------------------------------------------------------------------------
58 check_macro (f == quadrature_option::gauss,
59 "unsupported quadrature family \"" << opt.get_family_name() << "\"");
60
61 // Gauss-Legendre quadrature formulae
62 // where Legendre = Jacobi(alpha=0,beta=0) polynoms
63 // when using n nodes : quadrature formulae order is 2*r-1
64 size_type n = n_node_gauss(opt.get_order());
65 vector<T> zeta(n), omega(n);
66 gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
67 for (size_type i = 0; i < n; i++)
68 for (size_type j = 0; j < n; j++)
69 for (size_type k = 0; k < n; k++)
70 wx (x(zeta[i], zeta[j], zeta[k]), omega[i]*omega[j]*omega[k]);
71}
72// ----------------------------------------------------------------------------
73// instanciation in library
74// ----------------------------------------------------------------------------
75#define _RHEOLEF_instanciation(T) \
76template void quadrature_on_geo<T>::init_hexahedron (quadrature_option);
77
79
80}// 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_hexahedron(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