Loading...
Searching...
No Matches
basis_raw_bernstein.cc
Go to the documentation of this file.
#include "basis_raw_bernstein.h"
#include "basis_ordering.icc"
#include "bernstein.icc"
namespace rheolef {
using namespace std;
// =========================================================================
// basis raw bernstein members
// =========================================================================
template<class T>
basis_raw_bernstein<T>::~basis_raw_bernstein()
{
}
template<class T>
basis_raw_bernstein<T>::basis_raw_bernstein (std::string name)
: basis_raw_rep<T> (name),
_factorial(),
_power_index(),
_lambda_pow(),
_lambda_ad_pow()
{
if ((name.length()) > 0 && (name[0] == 'B')) {
// TODO: check also that name fits "Bk" where is an k integer
base::_degree = atoi(name.c_str()+1);
} else if (name.length() > 0) { // missing 'B' !
error_macro ("invalid polynomial name `"<<name<<"' for the Bk raw polynomial set");
} else {
// empty name : default cstor
base::_degree = 0;
}
}
template<class T>
typename basis_raw_bernstein<T>::size_type
basis_raw_bernstein<T>::ndof (reference_element hat_K) const
{
return reference_element::n_node (hat_K.variant(), base::_degree);
}
template<class T>
void
basis_raw_bernstein<T>::_initialize (reference_element hat_K) const
{
if (_factorial.size() == 0) {
_factorial.resize (base::_degree+1);
precompute_factorial (base::_degree, _factorial);
}
// for each bernstein with index loc_idof:
// p(x,y) = x^a*y^b*(1-x-y)^c
// power_index[loc_idof] contains the set of power indexes (a,b) with c=degree-a-b
// note: power_index is computed here one time for all
// as it does not depend upon hat_x but only upon (hat_K,degree)
// note: Bernstein polynoms are not hierarchized by degree
// => cannot be used easily for building RTk basis
make_power_indexes_sorted_by_inodes (hat_K, base::_degree, _power_index[hat_K.variant()]);
// note: lambda_pow depends upon hat_x, but is here allocated
// one time for all, as a working array,
// as its size does not depend upon hat_x but only upon (hat_K,degree)
_lambda_pow [hat_K.variant()].resize(base::_degree+1);
_lambda_ad_pow[hat_K.variant()].resize(base::_degree+1);
}
// evaluation of all basis functions at hat_x:
template<class T>
void
basis_raw_bernstein<T>::evaluate (
reference_element hat_K,
const point_basic<T>& hat_x,
Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
{
base::_initialize_guard (hat_K);
size_t d = hat_K.dimension();
// each x^a and y^b and (1-x-y)^c are computed first
// for all a,b,c in the lambda_pow array
// => this avoid imbricated loops
precompute_power_bernstein (hat_K, d, hat_x, base::_degree, _lambda_pow[hat_K.variant()]);
size_t loc_ndof = _power_index[hat_K.variant()].size();
value.resize(loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
value[loc_idof] = eval_bernstein_internal (hat_K, d, _lambda_pow[hat_K.variant()], _factorial, _power_index[hat_K.variant()][loc_idof], base::_degree);
}
}
template<class T>
void
basis_raw_bernstein<T>::grad_evaluate (
reference_element hat_K,
const point_basic<T>& hat_x,
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
{
base::_initialize_guard (hat_K);
point_basic<ad3_basic<T> > hat_x_ad = ad3::point (hat_x);
size_t d = hat_K.dimension();
precompute_power_bernstein (hat_K, d, hat_x_ad, base::_degree, _lambda_ad_pow[hat_K.variant()]);
size_t loc_ndof = _power_index[hat_K.variant()].size();
value.resize(loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
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);
value[loc_idof] = bx.grad();
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class basis_raw_bernstein<T>;
_RHEOLEF_instanciation(Float)
}// namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
basis_ordering.icc
basis_raw_bernstein.h
bernstein.icc
Float
see the Float page for the full documentation
rheolef::ad3_basic
Definition ad3.h:50
rheolef::ad3_basic::grad
const point_basic< T > & grad() const
Definition ad3.h:160
rheolef::basis_raw_bernstein::~basis_raw_bernstein
~basis_raw_bernstein()
Definition basis_raw_bernstein.cc:32
rheolef::basis_raw_bernstein::size_type
base::size_type size_type
Definition basis_raw_bernstein.h:41
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::_lambda_ad_pow
std::array< std::vector< std::array< ad3_basic< T >, 6 > >, reference_element::max_variant > _lambda_ad_pow
Definition basis_raw_bernstein.h:81
rheolef::basis_raw_bernstein::_lambda_pow
std::array< std::vector< std::array< T, 6 > >, reference_element::max_variant > _lambda_pow
Definition basis_raw_bernstein.h:79
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::basis_raw_bernstein::basis_raw_bernstein
basis_raw_bernstein(std::string name)
Definition basis_raw_bernstein.cc:36
rheolef::basis_raw_bernstein::ndof
size_type ndof(reference_element hat_K) const
Definition basis_raw_bernstein.cc:55
rheolef::basis_raw_bernstein::_power_index
std::array< std::vector< point_basic< size_type > >, reference_element::max_variant > _power_index
Definition basis_raw_bernstein.h:77
rheolef::basis_raw_bernstein::_initialize
void _initialize(reference_element hat_K) const
Definition basis_raw_bernstein.cc:61
rheolef::basis_raw_bernstein::_factorial
std::vector< T > _factorial
Definition basis_raw_bernstein.h:75
rheolef::basis_raw_rep::name
std::string name() const
Definition basis_raw.h:49
rheolef::basis_raw_rep::_degree
size_type _degree
Definition basis_raw.h:80
rheolef::basis_raw_rep::basis_raw_rep
basis_raw_rep(std::string name="")
Definition basis_raw.cc:30
rheolef::basis_raw_rep::_initialize_guard
void _initialize_guard(reference_element hat_K) const
Definition basis_raw.h:90
rheolef::point_basic
Definition point.h:86
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::dimension
size_type dimension() const
Definition reference_element.h:101
rheolef::reference_element::variant
variant_type variant() const
Definition reference_element.h:99
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition reference_element.cc:201
T
Expr1::float_type T
Definition field_expr.h:230
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::if
if(!(yy_init))
Definition basis_get.cc:598
d
size_t d
Definition neumann-laplace-lambda.cc:28