Rheolef  7.1
an efficient C++ finite element environment
dubiner.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_DUBINER_ICC
2 #define _RHEOLEF_DUBINER_ICC
3 //
24 // evaluate the Proriol-Dubiner basis on a point of the reference element
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 5 september 2017
29 //
30 #include "rheolef/quadrature.h"
31 
32 #include "basis_option.h"
33 
34 namespace rheolef {
35 
36 // ----------------------------------------------
37 // coefficients in the recurrence formula of the Jacobi polynomials:
38 // see War-2010, page 5:4
39 // ----------------------------------------------
40 namespace details {
41 
42 namespace legendre {
43 template<class T>
44 static inline T a (size_t i) {
45  return T(int(i))/sqrt(T(2*int(i)+1)*T(2*int(i)-1));
46 }
47 } // namespace legendre
48 
49 namespace jacobi {
50 template<class T>
51 static inline T a (T alpha, T beta, size_t degree) {
52  return (2*degree + 1 + alpha + beta)*(2*degree + 2 + alpha + beta)
53  /(2*(degree + 1)*(degree + 1 + alpha + beta));
54 }
55 template<class T>
56 static inline T b (T alpha, T beta, size_t degree) {
57  return (sqr(alpha) - sqr(beta))*(2*degree + 1 + alpha + beta)
58  /(2*(degree + 1)*(2*degree + alpha + beta)*(degree + 1 + alpha + beta));
59 }
60 template<class T>
61 static inline T c (T alpha, T beta, size_t degree) {
62  return (degree + alpha)*(degree + beta)*(2*degree + 2 + alpha + beta)
63  /((degree + 1)*(degree + 1 + alpha + beta)*(2*degree + alpha + beta));
64 }
65 } // namespace jacobi
66 } // namespace details
67 
68 // ----------------------------------------------
69 // Legendre/Dubiner basis evaluation:
70 // - Dubiner basis [Kir-2010] when K=t,T
71 // - Legendre basis when K=e,q,H
72 // - product basis when K=P
73 // ----------------------------------------------
74 // ortho-normalized Legendre polynoms on a segment: see e.g. HesWar-2008, p. 45
75 // Note: HesWar-2008 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
76 template<class T, class Container>
77 static
78 void
79 eval_dubiner_basis_e (
80  const T& tilde_x, // in [-1,1]
81  size_t degree,
82  const std::vector<size_t>& perm,
83  Container& value,
84  std::vector<point_basic<size_t> >& power_index)
85 {
87  typedef typename float_traits<T>::type float_type;
88  value.resize (degree+1);
89  power_index.resize (value.size());
90  size_t loc_idof;
92  power_index[perm[loc_idof]] = point_basic<size_t>(0);
93  value [perm[loc_idof]] = 1/sqrt(float_type(2));
94  if (degree == 0) return;
96  power_index[perm[loc_idof]] = point_basic<size_t>(1);
97  value [perm[loc_idof]] = sqrt(float_type(3)/float_type(2))*tilde_x;
98  for (size_type i = 1; i+1 <= degree; i++) {
99  namespace L = rheolef::details::legendre;
103  power_index[perm[loc_idof]] = point_basic<size_t>(i+1);
104  value [perm[loc_idof]] = (tilde_x*value[perm[loc_idof_i]] - L::a<float_type>(i)*value[perm[loc_idof_im1]])/L::a<float_type>(i+1);
105  }
106 }
107 // orthogonal Dubiner polynoms on a triangle: see Kir-2010, p. 5:7
108 // Notes: polynoms are not normalized => mass matrix is diag but not identity
109 // Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
110 template<class T, class Container>
111 static
112 void
113 eval_dubiner_basis_t (
114  const point_basic<T>& tilde_x, // in [-1,1]^2
115  size_t degree,
116  const std::vector<size_t>& perm,
117  Container& value,
118  std::vector<point_basic<size_t> >& power_index)
119 {
121  typedef typename float_traits<T>::type float_type;
122  value.resize ((degree+1)*(degree+2)/2);
123  power_index.resize (value.size());
124  // Dubiner, by algorithm Kir-2010, page 5:7
125  T xi0 = tilde_x[0],
126  xi1 = tilde_x[1];
127  size_t loc_idof;
128  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
129  power_index[perm[loc_idof]] = point_basic<size_t>(0,0);
130  value [perm[loc_idof]] = 1;
131  if (degree == 0) return;
132  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(1,0));
133  power_index[perm[loc_idof]] = point_basic<size_t>(1,0);
134  value [perm[loc_idof]] = (1 + 2*xi0 + xi1)/float_type(2);
135  for (size_type i = 1; i < degree; i++) {
136  size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i ,0));
137  size_t loc_idof_im1_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i-1,0));
138  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i+1,0));
139  power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0);
140  value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*(1 + 2*xi0 + xi1)/float_type(2)*value[perm[loc_idof_i_0]]
141  - float_type( int(i) )/float_type(int(i)+1)*sqr((1 - xi1)/float_type(2)) *value[perm[loc_idof_im1_0]];
142  }
143  for (size_type i = 0; i < degree; i++) {
144  size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
145  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,1));
146  power_index[perm[loc_idof]] = point_basic<size_t>(i,1);
147  value [perm[loc_idof]] = value[perm[loc_idof_i_0]]*(1 + 2*i + (3+2*i)*xi1)/float_type(2);
148  }
149  for (size_type i = 0; i < degree; i++) {
150  for (size_type j = 1; i+j < degree; j++) {
151  namespace J = rheolef::details::jacobi;
152  size_t loc_idof_i_j = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j ));
153  size_t loc_idof_i_jm1 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j-1));
154  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j+1));
155  power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1);
156  value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*xi1 + J::b<float_type>(2*int(i)+1,0,j))*value[perm[loc_idof_i_j]]
157  - J::c<float_type>(2*int(i)+1,0,j) *value[perm[loc_idof_i_jm1]];
158  }
159  }
160 }
161 template<class T, class Container>
162 static
163 void
164 eval_dubiner_basis_T (
165  const point_basic<T>& tilde_x, // in [-1,1]^3
166  size_t degree,
167  const std::vector<size_t>& perm,
168  Container& value,
169  std::vector<point_basic<size_t> >& power_index)
170 {
172  typedef typename float_traits<T>::type float_type;
174  power_index.resize (value.size());
175  // Dubiner, by algorithm Kir-2010, page 5:7
176  T xi0 = tilde_x[0],
177  xi1 = tilde_x[1],
178  xi2 = tilde_x[2];
179  T F1 = (2 + 2*xi0 + xi1 + xi2)/2,
180  F2 = sqr((xi1 + xi2)/2),
181  F3 = (2 + 3*xi1 + xi2)/2,
182  F4 = (1 + 2*xi1 + xi2)/2,
183  F5 = (1 - xi2)/2;
184  size_t loc_idof;
185  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(0,0,0));
186  power_index[perm[loc_idof]] = point_basic<size_t>(0,0,0);
187  value [perm[loc_idof]] = 1;
188  if (degree == 0) return;
189  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(1,0,0));
190  power_index[perm[loc_idof]] = point_basic<size_t>(1,0,0);
191  value [perm[loc_idof]] = F1;
192  for (size_type i = 1; i+1 <= degree; i++) {
193  size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i ,0,0));
194  size_t loc_idof_im1_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i-1,0,0));
195  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i+1,0,0));
196  power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0,0);
197  value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*F1*value[perm[loc_idof_i_0_0]]
198  - float_type( int(i) )/float_type(int(i)+1)*F2*value[perm[loc_idof_im1_0_0]];
199  }
200  for (size_type i = 0; i+1 <= degree; i++) {
201  size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,0,0));
202  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,1,0));
203  power_index[perm[loc_idof]] = point_basic<size_t>(i,1,0);
204  value [perm[loc_idof]] = (i*(1 + xi1) + F3)*value[perm[loc_idof_i_0_0]];
205  }
206  for (size_type i = 0; i+2 <= degree; i++) {
207  for (size_type j = 1; i+j+1 <= degree; j++) {
208  namespace J = rheolef::details::jacobi;
209  size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j ,0));
210  size_t loc_idof_i_jm1_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j-1,0));
211  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j+1,0));
212  power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1,0);
213  value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*F4 + J::b<float_type>(2*int(i)+1,0,j)*F5)*value[perm[loc_idof_i_j_0]]
214  - J::c<float_type>(2*int(i)+1,0,j)*sqr(F5) *value[perm[loc_idof_i_jm1_0]];
215  }}
216  for (size_type i = 0; i+1 <= degree; i++) {
217  for (size_type j = 0; i+j+1 <= degree; j++) {
218  size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,0));
219  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,1));
220  power_index[perm[loc_idof]] = point_basic<size_t>(i,j,1);
221  value [perm[loc_idof]] = (1+i+j + (2+i+j)*xi2)*value[perm[loc_idof_i_j_0]];
222  }}
223  for (size_type i = 0; i+2 <= degree; i++) {
224  for (size_type j = 0; i+j+2 <= degree; j++) {
225  for (size_type k = 1; i+j+k+1 <= degree; k++) {
226  namespace J = rheolef::details::jacobi;
227  size_t loc_idof_i_j_k = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
228  size_t loc_idof_i_j_km1 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k-1));
229  loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k+1));
230  power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k+1);
231  value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+2*int(j)+2,0,k)*xi2 + J::b<float_type>(2*int(i)+2*int(j)+2,0,k))*value[perm[loc_idof_i_j_k]]
232  - J::c<float_type>(2*int(i)+2*int(j)+2,0,k) *value[perm[loc_idof_i_j_km1]];
233  }}}
234 }
235 // --------------------------------------------------
236 // main call: evaluate all the basis at a given point
237 // --------------------------------------------------
238 template<class T, class Container>
239 static
240 void
241 eval_dubiner_basis (
242  const point_basic<T>& hat_x,
243  reference_element hat_K,
244  size_t degree,
245  const std::vector<size_t>& perm,
246  Container& value,
247  std::vector<point_basic<size_t> >& power_index,
248  bool do_tilde = true)
249 {
251  value.resize (reference_element::n_node(hat_K.variant(), degree));
252  power_index.resize (value.size());
253  switch (hat_K.variant()) {
254  case reference_element::p: {
255  power_index[perm[0]] = point_basic<size_t>();
256  value [perm[0]] = 1;
257  break;
258  }
259  case reference_element::e: {
260  // HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
261  T tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
262  eval_dubiner_basis_e (tilde_x, degree, perm, value, power_index);
263  break;
264  }
265  case reference_element::t: {
266  // Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
267  point_basic<T> tilde_x
268  = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1) : hat_x;
269  eval_dubiner_basis_t (tilde_x, degree, perm, value, power_index);
270  break;
271  }
272  case reference_element::T: {
273  // Dubiner, by algorithm Kir-2010, page 5:7
274  // Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
275  point_basic<T> tilde_x
276  = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1, 2*hat_x[2]-1) : hat_x;
277  eval_dubiner_basis_T (tilde_x, degree, perm, value, power_index);
278  break;
279  }
280  case reference_element::q: {
281  Container value0, value1;
282  std::vector<size_t> id_e (degree+1);
283  for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
284  std::vector<point_basic<size_t> > power_index_e;
285  eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
286  eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
287  for (size_type i = 0; i <= degree; i++) {
288  for (size_type j = 0; j <= degree; j++) {
289  size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
290  size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
291  size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
292  power_index[perm[loc_idof]] = point_basic<size_t>(i,j);
293  value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e];
294  }}
295  break;
296  }
297  case reference_element::H: {
298  Container value0, value1, value2;
299  std::vector<size_t> id_e (degree+1);
300  for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
301  std::vector<point_basic<size_t> > power_index_e;
302  eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
303  eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
304  eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
305  for (size_type i = 0; i <= degree; i++) {
306  for (size_type j = 0; j <= degree; j++) {
307  for (size_type k = 0; k <= degree; k++) {
308  size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
309  size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
310  size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
311  size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
312  power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
313  value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e]*value2[loc_kdof_e];
314  }}}
315  break;
316  }
317  case reference_element::P: {
318  Container value01, value2;
319  std::vector<size_t> id_e (degree+1), id_t ((degree+1)*(degree+2)/2);
320  for (size_type iloc = 0; iloc < id_e.size(); iloc ++) { id_e [iloc ] = iloc ; }
321  for (size_type iloc = 0; iloc < id_t.size(); iloc ++) { id_t [iloc ] = iloc ; }
322  point_basic<T> tilde_x
323  = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1)
324  : point_basic<T>( hat_x[0], hat_x[1]);
325  std::vector<point_basic<size_t> > power_index_e, power_index_t;
326  eval_dubiner_basis_t (tilde_x, degree, id_t, value01, power_index_t);
327  eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
328  for (size_type i = 0; i <= degree; i++) {
329  for (size_type j = 0; i+j <= degree; j++) {
330  for (size_type k = 0; k <= degree; k++) {
331  size_t loc_ijdof_t = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
332  size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
333  size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
334  power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
335  value [perm[loc_idof]] = value01[loc_ijdof_t]*value2[loc_kdof_e];
336  }}}
337  break;
338  }
339  default:
340  error_macro ("unsupported element: "<<hat_K.name());
341  }
342 }
343 // omit the power_index last arg
344 template<class T, class Container>
345 static
346 void
347 eval_dubiner_basis (
348  const point_basic<T>& hat_x,
349  reference_element hat_K,
350  size_t degree,
351  const std::vector<size_t>& perm,
352  Container& value,
353  bool do_tilde = true)
354 {
355  std::vector<point_basic<size_t> > power_index;
356  eval_dubiner_basis (hat_x, hat_K, degree, perm, value, power_index, do_tilde);
357 }
358 
359 
360 }// namespace rheolef
361 #endif // _RHEOLEF_DUBINER_ICC
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::point_basic< size_t >
rheolef::reference_element_P::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:695
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
rheolef::value
rheolef::std value
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:481
a
Definition: diffusion_isotropic.h:25
basis_option.h
basis_option - finite element method options
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::float_traits::type
T type
Definition: Float.h:94
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::reference_element_H::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:790
mkgeo_obstacle.L
L
Definition: mkgeo_obstacle.sh:133
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::jacobi
Definition: jacobi.icc:24
rheolef::details::legendre
Definition: dubiner.icc:42
rheolef::details::jacobi
Definition: dubiner.icc:49
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
T
Expr1::float_type T
Definition: field_expr.h:218