Rheolef  7.1
an efficient C++ finite element environment
quadrature_Te.cc
Go to the documentation of this file.
1 #include "rheolef/quadrature.h"
22 #include "rheolef/gauss_jacobi.h"
23 namespace rheolef {
24 using namespace std;
25 
26 template<class T>
27 void
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) \
182 template void quadrature_on_geo<T>::init_tetrahedron (quadrature_option);
183 
185 
186 }// namespace rheolef
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
rheolef::gauss_jacobi
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)
Definition: gauss_jacobi.icc:29
rheolef::point_basic
Definition: point.h:87
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_contraction.c2
c2
Definition: mkgeo_contraction.sh:199
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
a
Definition: diffusion_isotropic.h:25
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition: quadrature.h:88
Float
see the Float page for the full documentation
rheolef::quadrature_on_geo::init_tetrahedron
void init_tetrahedron(quadrature_option opt)
Definition: quadrature_Te.cc:28
f
Definition: cavity_dg.h:29
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
T
Expr1::float_type T
Definition: field_expr.h:261