Rheolef  7.1
an efficient C++ finite element environment
sherwin.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_SHERWIN_ICC
2 #define _RHEOLEF_SHERWIN_ICC
3 //
24 // evaluate the Dubiner-Sherwin basis on a point of the reference element
25 // see SheKar-2005, pages 52 (1d) & 114 (2d) & app. D
26 //
27 // author: Pierre.Saramito@imag.fr
28 //
29 // date: 5 september 2017
30 //
31 #include "rheolef/quadrature.h"
32 
33 #include "basis_option.h"
35 
36 #ifdef SHERWIN_USE_GINAC
37 #include <ginac/ginac.h>
38 #endif // SHERWIN_USE_GINAC
39 
40 namespace rheolef {
41 
42 // --------------------------------------------------------------------------
43 // compute all [0:order] Jacobi polynomials
44 // --------------------------------------------------------------------------
45 namespace details {
46 
47 template<class T0, class T, class Container>
48 static
49 void
50 eval_jacobi (
51  const T0& tilde_x, // in [-1,1]
52  size_t order,
53  const T& alpha,
54  const T& beta,
55  Container& value)
56 {
58  value.resize (order+1);
59  value[0] = 1;
60  if (order == 0) return;
61  value[1] = ((alpha+beta+2)*tilde_x + alpha - beta)/2;
62  for (size_t r = 1; r < order; r++) {
63  T a = (beta*beta - alpha*alpha)/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+2)));
64  T b = 2*(alpha+T(1.*r))*(beta+T(1.*r))/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+1)));
65  T c = 2*T(r+1.)*(alpha+beta+T(r+1.))/((alpha+beta+T(2.*r+1))*(alpha+beta+T(2.*r+2)));
66  value[r+1] = ((tilde_x-a)*value[r] - b*value[r-1])/c;
67  }
68 }
69 
70 } // namespace details
71 // --------------------------------------------------------------------------
72 // Note: SheKar-2005 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
73 // --------------------------------------------------------------------------
74 template<class T0, class T, class Container>
75 static
76 void
77 eval_sherwin_basis_e (
78  const T0& tilde_x, // in [-1,1]
79  size_t degree,
80  const T& alpha,
81  const T& beta,
82  Container& jacobi, // working array
83  Container& value)
84 {
86  typedef typename Container::value_type value_type;
88  value.resize (degree+1);
89  size_t loc_idof;
90  if (degree == 0) {
91  loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
92  value [loc_idof] = 1;
93  return;
94  }
95  loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
96  value [loc_idof] = (1-tilde_x)/2;
97  loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(degree));
98  value [loc_idof] = (1+tilde_x)/2;
99  if (degree == 1) return;
100  details::eval_jacobi (tilde_x, degree-2, alpha, beta, jacobi);
101  for (size_type i = 1; i+1 <= degree; i++) {
102  loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
103  value [loc_idof] = 0.25*(1-tilde_x)*(1+tilde_x)*jacobi [i-1];
104  }
105 }
106 template<class T, class Container>
107 inline
108 bool
110  const point_basic<T>& tilde_x, // in [-1,1]^2
111  reference_element hat_K,
112  size_t degree,
113  Container& value)
114 {
115  typedef typename Container::value_type value_type;
118  size_t d = hat_K.dimension();
119  if (d == 0) return false;
120  T distance_to_corner = fabs(tilde_x[d-1]-1);
121  if (distance_to_corner >= eps) {
122  return false;
123  }
124  // we are at top singular vertex:
125  for (size_t i = 0; i < size_t(value.size()); ++i) {
126  value[i] = value_type(0);
127  }
128  point_basic<size_t> ilat (0,0,0);
129  ilat [d] = degree;
130  size_t loc_idof = ilat2loc_inod (hat_K, degree, ilat);
131  value[loc_idof] = 1;
132  return true;
133 }
134 #ifdef SHERWIN_USE_GINAC
135 // is_singular is not relevant by ginac (symbolic computations):
136 template<>
137 inline
138 bool
139 eval_sherwin_basis_is_singular_point<GiNaC::ex,std::vector<GiNaC::ex> > (
140  const point_basic<GiNaC::ex>& tilde_x, // in [-1,1]^2
141  reference_element hat_K,
142  size_t degree,
143  std::vector<GiNaC::ex>& value)
144 {
145  return false;
146 }
147 #endif // SHERWIN_USE_GINAC
148 // Dubiner polynoms on a triangle with C0 junctions
149 // See KarShe-2005 p. 585
150 template<class T0, class T, class Container>
151 static
152 void
153 eval_sherwin_basis_t (
154  const point_basic<T0>& tilde_x, // in [-1,1]^2
155  size_t degree,
156  const T& alpha,
157  const T& beta,
158  Container& jacobi_x, // working array
159  Container& jacobi_y, // working array
160  Container& value)
161 {
163  typedef typename Container::value_type value_type;
165  value.resize ((degree+1)*(degree+2)/2);
166  if (degree == 0) {
167  value[0] = 1;
168  return;
169  }
170  if (eval_sherwin_basis_is_singular_point (tilde_x, reference_element(reference_element::t), degree, value)) return;
171  value_type eta0 = 2*(1+tilde_x[0])/(1-tilde_x[1])-1,
172  eta1 = tilde_x[1];
173  size_t loc_idof;
174 
175  // 1. verticies polynomials
176  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
177  value [loc_idof] = 0.25*(1-eta0)*(1-eta1);
178  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree,0));
179  value [loc_idof] = 0.25*(1+eta0)*(1-eta1);
180  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,degree));
181  value [loc_idof] = 0.5*(1+eta1);
182 
183  // 2. edge polynomials
184  if (degree <= 1) return;
185  // 2.1. 1st edge: (-1,-1) --> (1,-1)
186  details::eval_jacobi (eta0, degree-2, T(1), T(1), jacobi_x);
187  value_type coef = 0.5*(1-eta1);
188  for (size_type i = 1; i <= degree-1; ++i) {
189  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
190  coef *= 0.5*(1-eta1);
191  value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*coef*jacobi_x[i-1];
192  }
193  // 2.2. 2nd edge: (1,-1) --> (-1,1)
194  details::eval_jacobi (eta1, degree-2, T(1), T(1), jacobi_y);
195  coef = 0.5*(1+eta1);
196  for (size_type j = 1; j <= degree-1; ++j) {
197  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree-j,j));
198  value [loc_idof] = 0.25*(1+eta0)*(1-eta1)*coef*jacobi_y[j-1];
199  }
200  // 2.3. 3rd edge: (-1,1) --> (-1,-1)
201  for (size_type j = 1; j <= degree-1; ++j) {
202  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,j));
203  value [loc_idof] = 0.25*(1-eta0)*(1-eta1)*coef*jacobi_y[j-1];
204  }
205  // 3. interior polynomials
206  coef = 0.5*(1-eta1);
207  for (size_type i = 1; i <= degree-1; ++i) {
208  coef *= 0.5*(1-eta1);
209  details::eval_jacobi (eta1, degree-2, T(2*int(i)+1), T(1), jacobi_y);
210  for (size_type j = 1; i+j <= degree-1; ++j) {
211  loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
212  value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*jacobi_x[i-1]*coef
213  *0.5*(1+eta1)*jacobi_y[j-1];
214  }
215  }
216 }
217 #ifdef TODO
218 template<class T>
219 static
220 void
221 eval_dubiner_basis_T (
222  const point_basic<T>& tilde_x, // in [-1,1]^3
223  size_t degree,
224  const std::vector<size_t>& perm,
225  Eigen::Matrix<T,Eigen::Dynamic,1>& value,
226  std::vector<point_basic<size_t> >& power_index)
227 {
228 }
229 #endif // TODO
230 // --------------------------------------------------
231 // main call: evaluate all the basis at a given point
232 // --------------------------------------------------
233 template<class T0, class T, class Container>
234 static
235 void
236 eval_sherwin_basis (
237  const point_basic<T0>& hat_x,
238  reference_element hat_K,
239  size_t degree,
240  const T& alpha,
241  const T& beta,
242  Container& jacobi_x, // working array
243  Container& jacobi_y, // working array
244  Container& jacobi_z, // working array
245  Container& value,
246  bool do_tilde = true)
247 {
249  typedef typename Container::value_type value_type;
251  value.resize (reference_element::n_node(hat_K.variant(), degree));
252  switch (hat_K.variant()) {
253  case reference_element::p: {
254  value [0] = 1;
255  break;
256  }
257  case reference_element::e: {
258  // HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
259  value_type tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
260  eval_sherwin_basis_e (tilde_x, degree, alpha, beta, jacobi_x, value);
261  break;
262  }
263  case reference_element::t: {
264  // Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
265  point_basic<value_type> tilde_x
266  = do_tilde ? point_basic<value_type>(2*hat_x[0]-1, 2*hat_x[1]-1)
267  : point_basic<value_type>( hat_x[0], hat_x[1]);
268  eval_sherwin_basis_t (tilde_x, degree, alpha, beta, jacobi_x, jacobi_y, value);
269  break;
270  }
271 #ifdef TODO
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  Eigen::Matrix<T,Eigen::Dynamic,1> 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  Eigen::Matrix<T,Eigen::Dynamic,1> 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  Eigen::Matrix<T,Eigen::Dynamic,1> 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 #endif // TODO
340  default:
341  error_macro ("unsupported element: "<<hat_K.name());
342  }
343 }
344 
345 }// namespace rheolef
346 #endif // _RHEOLEF_SHERWIN_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
Definition: point.h:87
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::value
rheolef::std value
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
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
reference_element_aux.icc
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
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
rheolef::value_type
result_type value_type
Definition: field_expr_recursive.h:499
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_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
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
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
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::eval_sherwin_basis_is_singular_point
bool eval_sherwin_basis_is_singular_point(const point_basic< T > &tilde_x, reference_element hat_K, size_t degree, Container &value)
Definition: sherwin.icc:109
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