Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_raw_bernstein.cc
Go to the documentation of this file.
1
21#include "basis_raw_bernstein.h"
22#include "basis_ordering.icc"
23#include "bernstein.icc"
24
25namespace rheolef {
26using namespace std;
27
28// =========================================================================
29// basis raw bernstein members
30// =========================================================================
31template<class T>
35template<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}
53template<class T>
56{
57 return reference_element::n_node (hat_K.variant(), base::_degree);
58}
59template<class T>
60void
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:
82template<class T>
83void
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}
101template<class T>
102void
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) \
123template class basis_raw_bernstein<T>;
124
126
127}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
const point_basic< T > & grad() const
Definition ad3.h:160
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value) const
size_type ndof(reference_element hat_K) const
void _initialize(reference_element hat_K) const
std::string name() const
Definition basis_raw.h:49
see the reference_element page for the full documentation
variant_type variant() const
static size_type n_node(variant_type variant, size_type order)
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
STL namespace.