1 #ifndef _RHEOLEF_DUBINER_ICC
2 #define _RHEOLEF_DUBINER_ICC
30 #include "rheolef/quadrature.h"
44 static inline T a (
size_t i) {
45 return T(
int(i))/sqrt(
T(2*
int(i)+1)*
T(2*
int(i)-1));
53 /(2*(degree + 1)*(degree + 1 +
alpha +
beta));
76 template<
class T,
class Container>
79 eval_dubiner_basis_e (
82 const std::vector<size_t>& perm,
88 value.resize (degree+1);
89 power_index.resize (
value.size());
94 if (degree == 0)
return;
98 for (
size_type i = 1; i+1 <= degree; i++) {
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);
110 template<
class T,
class Container>
113 eval_dubiner_basis_t (
116 const std::vector<size_t>& perm,
122 value.resize ((degree+1)*(degree+2)/2);
123 power_index.resize (
value.size());
129 power_index[perm[loc_idof]] = point_basic<size_t>(0,0);
130 value [perm[loc_idof]] = 1;
131 if (degree == 0)
return;
133 power_index[perm[loc_idof]] = point_basic<size_t>(1,0);
139 power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0);
146 power_index[perm[loc_idof]] = point_basic<size_t>(i,1);
150 for (
size_type j = 1; i+j < degree; j++) {
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]];
161 template<
class T,
class Container>
164 eval_dubiner_basis_T (
167 const std::vector<size_t>& perm,
174 power_index.resize (
value.size());
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,
186 power_index[perm[loc_idof]] = point_basic<size_t>(0,0,0);
187 value [perm[loc_idof]] = 1;
188 if (degree == 0)
return;
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++) {
196 power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0,0);
200 for (
size_type i = 0; i+1 <= degree; i++) {
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]];
206 for (
size_type i = 0; i+2 <= degree; i++) {
207 for (
size_type j = 1; i+j+1 <= degree; j++) {
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]];
216 for (
size_type i = 0; i+1 <= degree; i++) {
217 for (
size_type j = 0; i+j+1 <= degree; j++) {
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]];
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++) {
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]];
238 template<
class T,
class Container>
243 reference_element hat_K,
245 const std::vector<size_t>& perm,
248 bool do_tilde =
true)
252 power_index.resize (
value.size());
253 switch (hat_K.variant()) {
255 power_index[perm[0]] = point_basic<size_t>();
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);
269 eval_dubiner_basis_t (tilde_x, degree, perm,
value, power_index);
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);
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++) {
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];
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++) {
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];
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 ; }
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++) {
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];
340 error_macro (
"unsupported element: "<<hat_K.name());
344 template<
class T,
class Container>
349 reference_element hat_K,
351 const std::vector<size_t>& perm,
353 bool do_tilde =
true)
355 std::vector<point_basic<size_t> > power_index;
356 eval_dubiner_basis (hat_x, hat_K, degree, perm,
value, power_index, do_tilde);
basis_option - finite element method options
field::size_type size_type
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
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 size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
typename float_traits< value_type >::type float_type
#define error_macro(message)
Float alpha[pmax+1][pmax+1]
This file is part of Rheolef.