basis_raw_monomial.cc

Go to the documentation of this file.

#include "basis_raw_monomial.h"
#include "monomial.icc"
 
namespace rheolef {
using namespace std;
 
// =========================================================================
// basis raw monomial members
// =========================================================================
template<class T>
basis_raw_monomial<T>::~basis_raw_monomial()
{
}
template<class T>
basis_raw_monomial<T>::basis_raw_monomial (std::string name)
  : basis_raw_rep<T> (name),
    _power_index(),
    _hat_x_pow(),
    _hat_x_ad_pow()
{
trace_macro ("basis_raw_monomial: name="<<name);
  if ((name.length()) > 0 && (name[0] == 'M')) {
    // TODO: check also that name fits "Mk" where is an k integer
    base::_degree = atoi(name.c_str()+1);
  } else if (name.length() > 0) { // missing 'M' !
    error_macro ("invalid polynomial name `"<<name<<"' for the Mk polynomial set");
  } else {
    // empty name : default cstor
    base::_degree = 0;
  }
trace_macro ("basis_raw_monomial: name()="<<base::name());
}
template<class T>
typename basis_raw_monomial<T>::size_type
basis_raw_monomial<T>::ndof (reference_element hat_K) const
{
  return reference_element::n_node (hat_K.variant(), base::_degree);
}
template<class T>
void
basis_raw_monomial<T>::_initialize (reference_element hat_K) const
{
  // for each monomial with index loc_idof:
  //    p(x,y) = x^a*y^b
  // power_index[loc_idof] contains the set of power indexes (a,b)
  // note: power_index is computed here one time for all 
  // as it does not depend upon hat_x but onlhy upon (hat_K,degree)
  make_power_indexes_sorted_by_degrees (hat_K, base::_degree, _power_index[hat_K.variant()]);
  // note: hat_x_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)
  _hat_x_pow   [hat_K.variant()].resize (base::_degree+1);
  _hat_x_ad_pow[hat_K.variant()].resize (base::_degree+1);
}
// evaluation of all basis functions at hat_x:
template<class T>
void
basis_raw_monomial<T>::evaluate (
    reference_element                 hat_K,
    const point_basic<T>&             hat_x,
    Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
{
trace_macro ("basis_raw_monomial: evaluate (name="<<base::name()<<")...");
  base::_initialize_guard (hat_K);
  size_t d = hat_K.dimension();
  // each x^a and y^b are computed first
  // for all a,b=0..degree in the _hat_x_pow array
  // => this avoid imbricated loops
  precompute_power_monomial (hat_K, d, hat_x, base::_degree, _hat_x_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_monomial_internal (hat_K, d, _hat_x_pow[hat_K.variant()], _power_index[hat_K.variant()][loc_idof]);
  }
trace_macro ("basis_raw_monomial: evaluate (name="<<base::name()<<") done");
}
template<class T>
void
basis_raw_monomial<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_monomial (hat_K, d, hat_x_ad, base::_degree, _hat_x_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_monomial_internal (hat_K, d, _hat_x_ad_pow[hat_K.variant()], _power_index[hat_K.variant()][loc_idof]);
    value[loc_idof] = bx.grad();
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_raw_monomial<T>;
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef