Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_Pr.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
28quadrature_on_geo<T>::init_prism (quadrature_option opt)
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(1/T(3),1/T(3),0), 1);
38 } else {
39 size_type n = ((r+1)*(r+2)/2)*(r+1);
40 T w = 1/T(int(n));
41 for (size_type i = 0; i <= r; i++) {
42 for (size_type j = 0; i+j <= r; j++) {
43 for (size_type k = 0; k <= r; k++) {
44 wx (x(T(int(i))/r,T(int(j))/r,2*T(int(k))/r-1), w);
45 }
46 }
47 }
48 }
49 return;
50 }
51 // -------------------------------------------------------------------------
52 // general case: Gauss
53 // -------------------------------------------------------------------------
54 check_macro (f == quadrature_option::gauss,
55 "unsupported quadrature family \"" << opt.get_family_name() << "\"");
56
57 switch (opt.get_order()) {
58 case 0:
59 case 1: {
60 // central point:
61 // better than the general case
62 // r=0 : 2 nodes
63 // r=1 : 8 nodes
64 // here: 1 node
65 T a = T(1)/3;
66 wx(x(a,a,0), 1);
67 break;
68 }
69 default: {
70 // Gauss-Legendre quadrature formulae
71 // where Legendre = Jacobi(alpha=0,beta=0) polynoms
72 size_t r = opt.get_order();
73 size_type n0 = n_node_gauss(r);
74 size_type n1 = n_node_gauss(r+1);
75 size_type n2 = n_node_gauss(r);
76 vector<T> zeta0(n0), omega0(n0);
77 vector<T> zeta1(n1), omega1(n1);
78 vector<T> zeta2(n2), omega2(n2);
79 gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
80 gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
81 gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
82
83 for (size_type i = 0; i < n0; i++) {
84 for (size_type j = 0; j < n1; j++) {
85 for (size_type k = 0; k < n2; k++) {
86 // we transform the square into the triangle
87 T eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
88 T eta_1 = (1+zeta1[j])/2;
89 T eta_2 = zeta2[k];
90 T J = (1-zeta1[j])/8;
91 wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
92 }
93 }
94 }
95 }
96 }
97}
98// ----------------------------------------------------------------------------
99// instanciation in library
100// ----------------------------------------------------------------------------
101#define _RHEOLEF_instanciation(T) \
102template void quadrature_on_geo<T>::init_prism (quadrature_option);
103
105
106}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
void init_prism(quadrature_option opt)
base::size_type size_type
Definition quadrature.h:88
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