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