Rheolef  7.2
an efficient C++ finite element environment
bernstein.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_BERNSTEIN_ICC
2 #define _RHEOLEF_BERNSTEIN_ICC
23 //
24 // Bernstein basis, as initial local FEM basis
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 2 september 2017
29 //
30 #include "rheolef/reference_element.h"
31 
32 namespace rheolef {
33 using namespace std;
34 
35 // pre-compute factorial[i] = i! for i=0..degree
36 // thus, avoid some inner loops
37 template<class T>
38 static
39 void
40 precompute_factorial (
41  size_t degree,
42  std::vector<T>& factorial)
43 {
45  factorial.resize (degree+1);
46  factorial[0] = 1;
47  for (size_type i = 1; i <= degree; i++) {
48  factorial[i] = factorial[i-1]*i;
49  }
50 }
51 // lambda[1][mu] = i-th barycentric coordinate, mu=0..d-1
52 // lambda[i][mu] = (lambda[1][mu])^i, i=0..degree, mu=0..d-1
53 // manage also variants for tensor-product elements (qPH)
54 template<class T>
55 static
56 void
57 precompute_power_bernstein (
58  reference_element hat_K,
59  size_t d,
60  const point_basic<T>& x,
61  size_t degree,
62  std::vector<std::array<T,6> >& lambda)
63 {
64  lambda.resize (degree+1);
65  size_t n_comp = 0;
66  switch (hat_K.variant()) {
67  case reference_element::p: n_comp = 0; break;
70  case reference_element::T: n_comp = d+1; break;
72  case reference_element::H: n_comp = 2*d; break;
73  case reference_element::P: n_comp = 5; break;
74  default: error_macro ("unsupported element: "<<hat_K.name());
75  }
76  for (size_t mu = 0; mu < n_comp; mu++) {
77  lambda[0][mu] = 1;
78  }
79  if (degree == 0) return;
80  switch (hat_K.variant()) {
82  break;
84  lambda[1][0] = x[0];
85  lambda[1][1] = 1-x[0];
86  break;
88  lambda[1][0] = x[0];
89  lambda[1][1] = x[1];
90  lambda[1][2] = 1-x[0]-x[1];
91  break;
93  lambda[1][0] = x[0];
94  lambda[1][1] = x[1];
95  lambda[1][2] = x[2];
96  lambda[1][3] = 1-x[0]-x[1]-x[2];
97  break;
98  case reference_element::q: {
99  T x0 = (1+x[0])/2,
100  x1 = (1+x[1])/2;
101  lambda[1][0] = x0;
102  lambda[1][1] = 1-x0;
103  lambda[1][2] = x1;
104  lambda[1][3] = 1-x1;
105  break;
106  }
107  case reference_element::P: {
108  T x_2 = (1+x[2])/2;
109  lambda[1][0] = x[0];
110  lambda[1][1] = x[1];
111  lambda[1][2] = 1-x[0]-x[1];
112  lambda[1][3] = x_2;
113  lambda[1][4] = 1-x_2;
114  break;
115  }
116  case reference_element::H: {
117  T x0 = (1+x[0])/2,
118  x1 = (1+x[1])/2,
119  x2 = (1+x[2])/2;
120  lambda[1][0] = x0;
121  lambda[1][1] = 1-x0;
122  lambda[1][2] = x1;
123  lambda[1][3] = 1-x1;
124  lambda[1][4] = x2;
125  lambda[1][5] = 1-x2;
126  break;
127  }
128  }
129  for (size_t i = 2; i <= degree; i++) {
130  for (size_t mu = 0; mu < n_comp; mu++) {
131  lambda[i][mu] = lambda[i-1][mu]*lambda[1][mu];
132  }
133  }
134 }
135 // ----------------------------------------------
136 // Bernstein basis evaluation
137 // result = c * prod_{i=0}^{d} lambda[i]^m[i]
138 // with C = s!/(prod m[i]!)
139 // s = sum_{i=0}^d m[i]!
140 // m[0] = sum_{i=1}^d m[i]!
141 // use precomputed i! and x[i]^m for efficiency
142 // Bernstein basis functions on quadrilaterals
143 // or hexahedra are defined using
144 // a tensor product construction [ChaWar-2016]
145 // ----------------------------------------------
146 template<class T, class U>
147 static
148 U
149 eval_bernstein_internal (
150  reference_element hat_K,
151  size_t d,
152  const std::vector<std::array<U,6> >& lambda_m,
153  const std::vector<T>& factorial,
154  const point_basic<size_t>& m,
155  size_t degree)
156 {
157  switch (hat_K.variant()) {
158  case reference_element::p: {
159  return 1;
160  }
161  case reference_element::e: {
162  T deno = factorial[ m[0]]
163  *factorial[degree-m[0]];
164  T c = factorial[degree]/deno;
165  return c*lambda_m[ m[0]][0]
166  *lambda_m[degree-m[0]][1];
167  }
168  case reference_element::t: {
169  T deno = factorial[ m[0] ]
170  *factorial[ m[1]]
171  *factorial[degree-m[0]-m[1]];
172  T c = factorial[degree]/deno;
173  return c*lambda_m[ m[0] ][0]
174  *lambda_m[ m[1]][1]
175  *lambda_m[degree-m[0]-m[1]][2];
176  }
177  case reference_element::T: {
178  T deno = factorial[ m[0] ]
179  *factorial[ m[1] ]
180  *factorial[ m[2]]
181  *factorial[degree-m[0]-m[1]-m[2]];
182  T c = factorial[degree]/deno;
183  return c*lambda_m[ m[0] ][0]
184  *lambda_m[ m[1] ][1]
185  *lambda_m[ m[2]][2]
186  *lambda_m[degree-m[0]-m[1]-m[2]][3];
187  }
188  case reference_element::q: {
189  T deno = factorial[ m[0]]
190  *factorial[degree-m[0]]
191  *factorial[ m[1]]
192  *factorial[degree-m[1]];
193  T c = sqr(factorial[degree])/deno;
194  return c*lambda_m[ m[0]][0]
195  *lambda_m[degree-m[0]][1]
196  *lambda_m[ m[1]][2]
197  *lambda_m[degree-m[1]][3];
198  }
199  case reference_element::H: {
200  T deno = factorial[ m[0]]
201  *factorial[degree-m[0]]
202  *factorial[ m[1]]
203  *factorial[degree-m[1]]
204  *factorial[ m[2]]
205  *factorial[degree-m[2]];
206  T c = pow(factorial[degree],3)/deno;
207  return c*lambda_m[ m[0]][0]
208  *lambda_m[degree-m[0]][1]
209  *lambda_m[ m[1]][2]
210  *lambda_m[degree-m[1]][3]
211  *lambda_m[ m[2]][4]
212  *lambda_m[degree-m[2]][5];
213  }
214  case reference_element::P: {
215  T deno = factorial[ m[0] ]
216  *factorial[ m[1]]
217  *factorial[degree-m[0]-m[1]]
218  *factorial[ m[2]]
219  *factorial[degree-m[2]];
220  T c = sqr(factorial[degree])/deno;
221  return c*lambda_m[ m[0] ][0]
222  *lambda_m[ m[1]][1]
223  *lambda_m[degree-m[0]-m[1]][2]
224  *lambda_m[ m[2] ][3]
225  *lambda_m[degree-m[2] ][4];
226  }
227  default:
228  error_macro ("unsupported element: "<<hat_K.name());
229  return 0;
230  }
231 }
232 
233 }// namespace rheolef
234 #endif // _RHEOLEF_BERNSTEIN_ICC
field::size_type size_type
Definition: branch.cc:430
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
std::vector< int >::size_type size_type
static const variant_type T
static const variant_type P
static const variant_type t
size_t size_type
Definition: basis_get.cc:76
point_basic< T >
Definition: piola_fem.h:135
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
This file is part of Rheolef.
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120