basis_fem_tensor.cc

Go to the documentation of this file.

#include "basis_fem_tensor.h"
#include "rheolef/rheostream.h"
 
namespace rheolef {
 
using namespace std;
 
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_tensor<T>::~basis_fem_tensor()
{
}
template<class T>
basis_fem_tensor<T>::basis_fem_tensor (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::tensor);
  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(),
    "tensor(basis): basis.option.map_dimension should be initialized for component number");
  size_type d = base::option().dimension();
  _n_comp = d*(d+1)/2;
  if (d == 2 &&
      (base::option().coordinate_system() == space_constant::axisymmetric_rz ||
       base::option().coordinate_system() == space_constant::axisymmetric_zr)) {
    _n_comp++;
  }
  _initialize_cstor_sizes();
}
#ifdef TO_CLEAN
template<class T>
typename basis_fem_tensor<T>::size_type
basis_fem_tensor<T>::n_component(size_type map_d) const
{
  const size_type unset = std::numeric_limits<basis_option::size_type>::max();
  size_type d = (base::option().dimension() == unset) ? map_d : base::option().dimension();
  size_type n_comp = d*(d+1)/2;
  if (d == 2 &&
      (base::option().coordinate_system() == space_constant::axisymmetric_rz ||
       base::option().coordinate_system() == space_constant::axisymmetric_zr)) {
    n_comp++;
  }
  return n_comp;
}
#endif // TO_CLEAN
template<class T>
void
basis_fem_tensor<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_tensor<T>::_initialize_data (reference_element hat_K) const
{
}
template<class T>
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
basis_fem_tensor<T>::hat_node (reference_element hat_K) const
{
  return _scalar_basis.hat_node (hat_K);
}
template<class T>
void
basis_fem_tensor<T>::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);
  space_constant::coordinate_type sys_coord = base::option().coordinate_system();
  _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 (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 ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
      size_type loc_idof = _n_comp*loc_comp_idof + ij_comp;
      std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
      size_type i_comp = ij.first;
      size_type j_comp = ij.second;
      value[loc_idof] (i_comp,j_comp) = _scalar_value[loc_comp_idof];
      if (ij.first == ij.second) continue;
      value[loc_idof] (j_comp,i_comp) = _scalar_value[loc_comp_idof];
    }
  }
}
template<class T>
void
basis_fem_tensor<T>::grad_evaluate (
    reference_element                                 hat_K,
    const point_basic<T>&                             hat_x,
    Eigen::Matrix<tensor3_basic<T>,Eigen::Dynamic,1>& value) const
{
  base::_initialize_data_guard (hat_K);
  space_constant::coordinate_type sys_coord = base::option().coordinate_system();
  _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 (tensor3_basic<T>()); // do not remove !
  for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
    for (size_type ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
      size_type loc_idof = _n_comp*loc_comp_idof + ij_comp;
      std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
      size_type i_comp = ij.first;
      size_type j_comp = ij.second;
      for (size_type k_comp = 0; k_comp < _n_comp; ++k_comp) {
        value[loc_idof](i_comp,j_comp,k_comp) = _vector_value[loc_comp_idof][k_comp];
        if (i_comp == j_comp) continue;
        value[loc_idof](j_comp,i_comp,k_comp) = _vector_value[loc_comp_idof][k_comp];
      }
    }
  }
}
// dofs for a scalar-valued function
template<class T>
void
basis_fem_tensor<T>::_compute_dofs (
  reference_element                                      hat_K,
  const Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& f_xnod, 
        Eigen::Matrix<T,Eigen::Dynamic,1>&               dof) const
{
  base::_initialize_data_guard (hat_K);
  space_constant::coordinate_type sys_coord = base::option().coordinate_system();
  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 ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
    std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
    size_type i_comp = ij.first;
    size_type j_comp = ij.second;
    for (size_type loc_comp_inod = 0; loc_comp_inod < loc_comp_nnod; ++loc_comp_inod) {
      // interpolate the symmetric value of the tensor-valued function
      f_comp_xnod [loc_comp_inod]
        = 0.5*( f_xnod [loc_comp_inod] (i_comp,j_comp)
              + f_xnod [loc_comp_inod] (j_comp,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 + ij_comp;
      dof [loc_idof] = comp_dof [loc_comp_idof];
    }
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_fem_tensor<T>;
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef