basis_fem_Pk_sherwin.cc

Go to the documentation of this file.

#include "basis_fem_Pk_sherwin.h"
#include "basis_fem_Pk_lagrange.h"
#include "piola_fem_lagrange.h"
#include "rheolef/rheostream.h"
#include "equispaced.icc"
#include "warburton.icc"
#include "sherwin.icc"
#include "eigen_util.h"
#include "basis_on_pointset_evaluate.icc"
 
namespace rheolef {
using namespace std;
 
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_Pk_sherwin<T>::~basis_fem_Pk_sherwin()
{
}
template<class T>
basis_fem_Pk_sherwin<T>::basis_fem_Pk_sherwin (size_type degree, const basis_option& sopt) 
  : base (sopt),
    _degree(degree),
    _alpha(1.0),
    _beta(1.0)
#ifdef TODO
    ,
    _value_ad(),
    _work0_ad(),
    _work1_ad(),
    _work2_ad()
    _work0(),
    _work1(),
    _work2()
#endif // TODO
{
  base::_name = base::standard_naming (family_name(), degree, base::_sopt);
warning_macro ("cstor...");
  _alpha = _beta = 1; // 1 or 2, see SheKar-2005 p. 53
  _initialize_cstor_sizes();
 
  // piola FEM transformatiion:
  typedef piola_fem_lagrange<T> piola_fem_type;
  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
warning_macro ("cstor done");
}
template<class T>
void
basis_fem_Pk_sherwin<T>::_initialize_cstor_sizes() const
{
  basis_fem_Pk_lagrange<T>::initialize_local_first (
    degree(),
    base::is_continuous(),
    base::_ndof_on_subgeo_internal,
    base::_ndof_on_subgeo,
    base::_nnod_on_subgeo_internal,
    base::_nnod_on_subgeo,
    base::_first_idof_by_dimension_internal,
    base::_first_idof_by_dimension,
    base::_first_inod_by_dimension_internal,
    base::_first_inod_by_dimension);
}
template<class T>
void
basis_fem_Pk_sherwin<T>::_initialize_data (reference_element hat_K) const
{
  // initialization is similar to Pk-Lagrange
  size_type k = degree();
  size_type variant = hat_K.variant();
 
  fatal_macro ("sherwin: not yet");
#ifdef TODO // requires a _raw_basis here
  // nodes: TODO: volume-internal nodes should be quadrature, because of integral dofs
  switch (base::_sopt.get_node()) {
    case basis_option::equispaced:
          pointset_lagrange_equispaced (hat_K, k, base::_hat_node[variant]);
          break;
    case basis_option::warburton:
          pointset_lagrange_warburton  (hat_K, k, base::_hat_node[variant]); break;
    default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
  }
  // vdm:
  details::basis_on_pointset_evaluate (_raw_basis, hat_K, base::_hat_node[variant], base::_vdm[variant]);
  check_macro (invert(base::_vdm[variant], base::_inv_vdm[variant]),
        "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
#endif // TODO
}
// evaluation of all basis functions at hat_x:
template<class T>
void
basis_fem_Pk_sherwin<T>::evaluate (
    reference_element                 hat_K,
    const point_basic<T>&             hat_x,
    Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
{
  base::_initialize_data_guard (hat_K);
 
  Eigen::Matrix<T,Eigen::Dynamic,1> _work0, _work1, _work2; // TODO: allocate one time for all in class ?
 
  eval_sherwin_basis (hat_x, hat_K, degree(), _alpha, _beta, _work0, _work1, _work2, value);
}
// evaluate the gradient:
template<class T>
void
basis_fem_Pk_sherwin<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_data_guard (hat_K);
  point_basic<ad3_basic<T> > hat_x_ad = ad3::point (hat_x);
 
  std::vector<ad3_basic<T> > _work0_ad, _work1_ad, _work2_ad; // TODO: allocate one time for all in class ?
  std::array<std::vector<ad3_basic<T> >,
             reference_element::max_variant>       _value_ad; // TODO: idem
 
  std::vector<ad3_basic<T> >& value_ad = _value_ad [hat_K.variant()];
 
  eval_sherwin_basis (hat_x_ad, hat_K, degree(), _alpha, _beta, _work0_ad, _work1_ad, _work2_ad, value_ad);
  size_t loc_ndof = value_ad.size();
  value.resize(loc_ndof);
  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
    value[loc_idof] = value_ad[loc_idof].grad();
  }
}
// dofs for a scalar-valued function
template<class T>
void
basis_fem_Pk_sherwin<T>::_compute_dofs (
  reference_element     hat_K,
  const Eigen::Matrix<T,Eigen::Dynamic,1>&   f_xnod, 
        Eigen::Matrix<T,Eigen::Dynamic,1>&   dof) const
{
#ifdef TODO
  // TODO: nodal values + integrals inside, instead of Lagrange node
  base::_initialize_data_guard (hat_K);
  dof = base::_inv_vdm[hat_K.variant()]*f_xnod;
#endif // TODO
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_fem_Pk_sherwin<T>;
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef