Rheolef  7.1
an efficient C++ finite element environment
basis_raw_bernstein.cc
Go to the documentation of this file.
1 #include "basis_raw_bernstein.h"
22 #include "basis_ordering.icc"
23 #include "bernstein.icc"
24 
25 namespace rheolef {
26 using namespace std;
27 
28 // =========================================================================
29 // basis raw bernstein members
30 // =========================================================================
31 template<class T>
33 {
34 }
35 template<class T>
37  : basis_raw_rep<T> (name),
38  _factorial(),
39  _power_index(),
40  _lambda_pow(),
41  _lambda_ad_pow()
42 {
43  if ((name.length()) > 0 && (name[0] == 'B')) {
44  // TODO: check also that name fits "Bk" where is an k integer
45  base::_degree = atoi(name.c_str()+1);
46  } else if (name.length() > 0) { // missing 'B' !
47  error_macro ("invalid polynomial name `"<<name<<"' for the Bk raw polynomial set");
48  } else {
49  // empty name : default cstor
50  base::_degree = 0;
51  }
52 }
53 template<class T>
56 {
57  return reference_element::n_node (hat_K.variant(), base::_degree);
58 }
59 template<class T>
60 void
62 {
63  if (_factorial.size() == 0) {
64  _factorial.resize (base::_degree+1);
65  precompute_factorial (base::_degree, _factorial);
66  }
67  // for each bernstein with index loc_idof:
68  // p(x,y) = x^a*y^b*(1-x-y)^c
69  // power_index[loc_idof] contains the set of power indexes (a,b) with c=degree-a-b
70  // note: power_index is computed here one time for all
71  // as it does not depend upon hat_x but only upon (hat_K,degree)
72  // note: Bernstein polynoms are not hierarchized by degree
73  // => cannot be used easily for building RTk basis
74  make_power_indexes_sorted_by_inodes (hat_K, base::_degree, _power_index[hat_K.variant()]);
75  // note: lambda_pow depends upon hat_x, but is here allocated
76  // one time for all, as a working array,
77  // as its size does not depend upon hat_x but only upon (hat_K,degree)
78  _lambda_pow [hat_K.variant()].resize(base::_degree+1);
79  _lambda_ad_pow[hat_K.variant()].resize(base::_degree+1);
80 }
81 // evaluation of all basis functions at hat_x:
82 template<class T>
83 void
85  reference_element hat_K,
86  const point_basic<T>& hat_x,
87  Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
88 {
89  base::_initialize_guard (hat_K);
90  size_t d = hat_K.dimension();
91  // each x^a and y^b and (1-x-y)^c are computed first
92  // for all a,b,c in the lambda_pow array
93  // => this avoid imbricated loops
94  precompute_power_bernstein (hat_K, d, hat_x, base::_degree, _lambda_pow[hat_K.variant()]);
95  size_t loc_ndof = _power_index[hat_K.variant()].size();
96  value.resize(loc_ndof);
97  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
98  value[loc_idof] = eval_bernstein_internal (hat_K, d, _lambda_pow[hat_K.variant()], _factorial, _power_index[hat_K.variant()][loc_idof], base::_degree);
99  }
100 }
101 template<class T>
102 void
104  reference_element hat_K,
105  const point_basic<T>& hat_x,
106  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
107 {
108  base::_initialize_guard (hat_K);
109  point_basic<ad3_basic<T> > hat_x_ad = ad3::point (hat_x);
110  size_t d = hat_K.dimension();
111  precompute_power_bernstein (hat_K, d, hat_x_ad, base::_degree, _lambda_ad_pow[hat_K.variant()]);
112  size_t loc_ndof = _power_index[hat_K.variant()].size();
113  value.resize(loc_ndof);
114  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
115  ad3_basic<T> bx = eval_bernstein_internal (hat_K, d, _lambda_ad_pow[hat_K.variant()], _factorial, _power_index[hat_K.variant()][loc_idof], base::_degree);
116  value[loc_idof] = bx.grad();
117  }
118 }
119 // ----------------------------------------------------------------------------
120 // instanciation in library
121 // ----------------------------------------------------------------------------
122 #define _RHEOLEF_instanciation(T) \
123 template class basis_raw_bernstein<T>;
124 
126 
127 }// namespace rheolef
rheolef::basis_raw_bernstein::size_type
base::size_type size_type
Definition: basis_raw_bernstein.h:41
rheolef::basis_raw_bernstein::_initialize
void _initialize(reference_element hat_K) const
Definition: basis_raw_bernstein.cc:61
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
rheolef::point_basic
Definition: point.h:87
rheolef::basis_raw_bernstein::~basis_raw_bernstein
~basis_raw_bernstein()
Definition: basis_raw_bernstein.cc:32
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::value
rheolef::std value
rheolef::ad3_basic::grad
const point_basic< T > & grad() const
Definition: ad3.h:160
rheolef::ad3_basic
Definition: ad3.h:50
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef::basis_raw_bernstein::evaluate
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
Definition: basis_raw_bernstein.cc:84
rheolef::basis_raw_bernstein::grad_evaluate
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value) const
Definition: basis_raw_bernstein.cc:103
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::basis_raw_rep
Definition: basis_raw.h:36
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::basis_raw_bernstein::basis_raw_bernstein
basis_raw_bernstein(std::string name)
Definition: basis_raw_bernstein.cc:36
rheolef::basis_raw_rep::_degree
size_type _degree
Definition: basis_raw.h:80
rheolef::basis_raw_bernstein::ndof
size_type ndof(reference_element hat_K) const
Definition: basis_raw_bernstein.cc:55
bernstein.icc
basis_raw_bernstein.h
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
basis_ordering.icc
rheolef::basis_raw_rep::name
std::string name() const
Definition: basis_raw.h:49
T
Expr1::float_type T
Definition: field_expr.h:218