Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_Te.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(1/T(3),1/T(3),1/(3)), T(1)/6);
38 } else {
39 size_type n = (r+1)*(r+2)*(r+3)/6;
40 T w = (1/T(6))/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; i+j+k <= r; k++) {
44 wx (x(T(int(i))/r,T(int(j))/r,T(int(k))/r), w);
45 }
46 }
47 }
48 }
49 return;
50 }
51 // -------------------------------------------------------------------------
52 // special case : superconvergent patch recovery nodes & weights
53 // -------------------------------------------------------------------------
54 if (f == quadrature_option::superconvergent) {
55 switch (opt.get_order()) {
56 case 0:
57 case 1:
58 case 2:
59 f = quadrature_option::gauss;
60 break;
61 default:
62 error_macro ("unsupported superconvergent("<<opt.get_order()<<")");
63 }
64 }
65 // -------------------------------------------------------------------------
66 // special case : Gauss-Lobatto points
67 // e.g. when using special option in riesz() function
68 // -------------------------------------------------------------------------
69 if (f == quadrature_option::gauss_lobatto) {
70 switch (opt.get_order()) {
71 case 0 :
72 case 1 :
73 // trapezes:
74 wx(x(T(0), T(0), T(0)), T(1)/24);
75 wx(x(T(1), T(0), T(0)), T(1)/24);
76 wx(x(T(0), T(1), T(0)), T(1)/24);
77 wx(x(T(0), T(0), T(1)), T(1)/24);
78 break;
79 case 2 :
80 // TODO: simpson
81 // break;
82 default:
83 error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
84 }
85 return;
86 }
87 // -------------------------------------------------------------------------
88 // general case: Gauss
89 // -------------------------------------------------------------------------
90 check_macro (f == quadrature_option::gauss,
91 "unsupported quadrature family \"" << opt.get_family_name() << "\"");
92
93 switch (opt.get_order()) {
94 case 0:
95 case 1:
96 // central point:
97 // better than the general case
98 // r=0 : 4 nodes
99 // r=1 : 12 nodes
100 // here: 1 node
101 wx (x(0.25, 0.25, 0.25), T(1)/6);
102 break;
103 case 2: {
104 // better than the general case
105 // r=2 : 18 nodes
106 // here: 4 nodes
107 T a = (5 - sqrt(T(5)))/20;
108 T b = (5 + 3*sqrt(T(5)))/20;
109 wx (x(a, a, a), T(1)/24);
110 wx (x(a, a, b), T(1)/24);
111 wx (x(a, b, a), T(1)/24);
112 wx (x(b, a, a), T(1)/24);
113 break;
114 }
115 case 3:
116 case 4:
117 case 5: {
118 // better than the general case
119 // r=3 : 36 nodes
120 // r=4 : 48 nodes
121 // r=5 : 80 nodes
122 // here: 15 nodes
123 T a = 0.25;
124 T b1 = (7 + sqrt(T(15)))/34;
125 T b2 = (7 - sqrt(T(15)))/34;
126 T c1 = (13 - 3*sqrt(T(15)))/34;
127 T c2 = (13 + 3*sqrt(T(15)))/34;
128 T w1 = (2665 - 14*sqrt(T(15)))/226800;
129 T w2 = (2665 + 14*sqrt(T(15)))/226800;
130 T d = (5 - sqrt(T(15)))/20;
131 T e = (5 + sqrt(T(15)))/20;
132 wx (x(a, a, a), T(8)/405);
133 wx (x(b1, b1, b1), w1);
134 wx (x(b1, b1, c1), w1);
135 wx (x(b1, c1, b1), w1);
136 wx (x(c1, b1, b1), w1);
137 wx (x(b2, b2, b2), w2);
138 wx (x(b2, b2, c2), w2);
139 wx (x(b2, c2, b2), w2);
140 wx (x(c2, b2, b2), w2);
141 wx (x(d, d, e), T(5)/567);
142 wx (x(d, e, d), T(5)/567);
143 wx (x(e, d, d), T(5)/567);
144 wx (x(d, e, e), T(5)/567);
145 wx (x(e, d, e), T(5)/567);
146 wx (x(e, e, d), T(5)/567);
147 break;
148 }
149 default: {
150 // Gauss-Legendre quadrature formulae
151 // where Legendre = Jacobi(alpha=0,beta=0) polynoms
152 size_t r = opt.get_order();
153 size_type n0 = n_node_gauss(r);
154 size_type n1 = n_node_gauss(r+1);
155 size_type n2 = n_node_gauss(r+2);
156 vector<T> zeta0(n0), omega0(n0);
157 vector<T> zeta1(n1), omega1(n1);
158 vector<T> zeta2(n2), omega2(n2);
159 gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
160 gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
161 gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
162
163 for (size_type i = 0; i < n0; i++) {
164 for (size_type j = 0; j < n1; j++) {
165 for (size_type k = 0; k < n2; k++) {
166 // we transform the cube into the tetra
167 T eta_0 = (1+zeta0[i])*(1-zeta1[j])*(1-zeta2[k])/8;
168 T eta_1 = (1+zeta1[j])*(1-zeta2[k])/4;
169 T eta_2 = (1+zeta2[k])/2;
170 T J = (1-zeta1[j])*sqr(1-zeta2[k])/64;
171 wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
172 }
173 }
174 }
175 }
176 }
177}
178// ----------------------------------------------------------------------------
179// instanciation in library
180// ----------------------------------------------------------------------------
181#define _RHEOLEF_instanciation(T) \
182template void quadrature_on_geo<T>::init_tetrahedron (quadrature_option);
183
185
186}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
void init_tetrahedron(quadrature_option opt)
base::size_type size_type
Definition quadrature.h:88
#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