basis_fem_vector.cc

Go to the documentation of this file.

//
// vectors
//
#include "basis_fem_vector.h"
#include "rheolef/rheostream.h"
 
namespace rheolef {
 
using namespace std;
 
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_vector<T>::~basis_fem_vector()
{
}
template<class T>
basis_fem_vector<T>::basis_fem_vector (const basis_basic<T>& scalar_basis, const basis_option& sopt)
  : basis_rep<T>(sopt),
    _n_comp(0),
    _scalar_basis(scalar_basis),
    _scalar_value(),
    _vector_value()
{
  base::_sopt.set_valued_tag (space_constant::vector);
  base::_name      = base::standard_naming (family_name(), family_index(), base::_sopt);
  base::_piola_fem = _scalar_basis.get_piola_fem();
  check_macro (base::option().dimension() != std::numeric_limits<basis_option::size_type>::max(),
    "vector(basis): basis.option.map_dimension should be initialized for component number");
  _n_comp = base::option().dimension();
  _initialize_cstor_sizes();
}
#ifdef TO_CLEAN
template<class T>
typename basis_fem_vector<T>::size_type
basis_fem_vector<T>::n_component(size_type map_d) const
{
  const size_type unset = std::numeric_limits<basis_option::size_type>::max();
  return (base::option().dimension() == unset) ? map_d : base::option().dimension();
}
#endif // TO_CLEAN
template<class T>
void
basis_fem_vector<T>::_initialize_cstor_sizes () const
{
  for (size_type map_d = 0; map_d < 4; ++map_d) {
    for (size_type subgeo_variant = 0; subgeo_variant < reference_element::max_variant; ++subgeo_variant) {
      base::_ndof_on_subgeo [map_d][subgeo_variant] = _n_comp*_scalar_basis.ndof_on_subgeo (map_d, subgeo_variant);
      base::_nnod_on_subgeo [map_d][subgeo_variant] =         _scalar_basis.nnod_on_subgeo (map_d, subgeo_variant);
    }
  }
  for (size_type variant = 0; variant < reference_element::max_variant; ++variant) {
    reference_element hat_K (variant);
    for (size_type subgeo_d = 0; subgeo_d < 5; ++subgeo_d) {
      base::_first_idof_by_dimension [variant][subgeo_d] = _n_comp*_scalar_basis.first_idof_by_dimension (hat_K, subgeo_d);
      base::_first_inod_by_dimension [variant][subgeo_d] =         _scalar_basis.first_inod_by_dimension (hat_K, subgeo_d);
    }
  }
}
template<class T>
void
basis_fem_vector<T>::_initialize_data (reference_element hat_K) const
{
}
template<class T>
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
basis_fem_vector<T>::hat_node (reference_element hat_K) const
{
  return _scalar_basis.hat_node (hat_K);
}
template<class T>
void
basis_fem_vector<T>::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);
  _scalar_basis.evaluate (hat_K, hat_x, _scalar_value); 
  size_type loc_comp_ndof = _scalar_value.size();
  size_type loc_ndof      = _n_comp*loc_comp_ndof;
  value.resize (loc_ndof);
  value.fill (point_basic<T>()); // do not remove !
  for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
    for (size_type i_comp = 0; i_comp < _n_comp; ++i_comp) {
      size_type loc_idof = _n_comp*loc_comp_idof + i_comp;
      value[loc_idof][i_comp] = _scalar_value[loc_comp_idof];
    }
  }
}
template<class T>
void
basis_fem_vector<T>::grad_evaluate (
    reference_element                                hat_K,
    const point_basic<T>&                            hat_x,
    Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value) const
{
  base::_initialize_data_guard (hat_K);
  _scalar_basis.grad_evaluate (hat_K, hat_x, _vector_value); 
  size_type loc_comp_ndof = _vector_value.size();
  size_type loc_ndof      = _n_comp*loc_comp_ndof;
  value.resize (loc_ndof);
  value.fill (tensor_basic<T>()); // do not remove !
  for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
    for (size_type i_comp = 0; i_comp < _n_comp; ++i_comp) {
      size_type loc_idof = _n_comp*loc_comp_idof + i_comp;
      for (size_type j_comp = 0; j_comp < _n_comp; ++j_comp) {
        value[loc_idof] (i_comp,j_comp) = _vector_value[loc_comp_idof][j_comp];
      }
    }
  }
}
template<class T>
void
basis_fem_vector<T>::_compute_dofs (
  reference_element                                      hat_K,
  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&  f_xnod, 
        Eigen::Matrix<T,Eigen::Dynamic,1>&               dof) const
{
  base::_initialize_data_guard (hat_K);
  size_type loc_comp_ndof = _scalar_basis.ndof (hat_K);
  size_type loc_comp_nnod = _scalar_basis.nnod (hat_K);
  size_type loc_ndof      = _n_comp*loc_comp_ndof;
  Eigen::Matrix<T,Eigen::Dynamic,1> f_comp_xnod (loc_comp_nnod); // TODO: class working array
  Eigen::Matrix<T,Eigen::Dynamic,1> comp_dof    (loc_comp_ndof); // TODO: class working array
  dof.resize (loc_ndof);
  for (size_type i_comp = 0; i_comp < _n_comp; ++i_comp) {
    for (size_type loc_comp_inod = 0; loc_comp_inod < loc_comp_nnod; ++loc_comp_inod) {
      f_comp_xnod [loc_comp_inod] = f_xnod [loc_comp_inod] [i_comp];
    }
    _scalar_basis.compute_dofs (hat_K, f_comp_xnod, comp_dof);
    for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
      size_type loc_idof = _n_comp*loc_comp_idof + i_comp;
      dof [loc_idof] = comp_dof [loc_comp_idof];
    }
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_fem_vector<T>;
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef