basis_fem_Pk_lagrange.cc

Go to the documentation of this file.

#include "basis_fem_Pk_lagrange.h"
#include "rheolef/rheostream.h"
 
#include "piola_fem_lagrange.h"
#include "equispaced.icc"
#include "warburton.icc"
#include "fekete.icc"
#include "Pk_get_local_idof_on_side.icc"
#include "basis_on_pointset_evaluate.icc"
#include "eigen_util.h"
 
namespace rheolef {
 
using namespace std;
 
// =========================================================================
// utilities
// =========================================================================
template <class T>
void
basis_fem_Pk_lagrange<T>::initialize_local_first (
  size_type k,
  bool is_continuous,
  std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo_internal,
  std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo,
  std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo_internal,
  std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo,
  std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension_internal,
  std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension,
  std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension_internal,
  std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension)
{
  // 1) ndof_on_subgeo
  if (k == size_t(-1)) {
    // P_{-1} == empty polynomial space
    for (size_type map_d = 0; map_d < 4; ++map_d) {
      ndof_on_subgeo_internal [map_d].fill (0);
    }
  } else if (k == 0) {
    // P0 initialization
    for (size_type map_d = 0; map_d < 4; ++map_d) {
      ndof_on_subgeo_internal [map_d].fill (0);
      for (size_type variant = reference_element::first_variant_by_dimension (map_d);
                     variant < reference_element:: last_variant_by_dimension (map_d);
                   ++variant) {
        ndof_on_subgeo_internal [map_d] [variant] = 1;
      }
    }
  } else { // k > 0
    for (size_type map_d = 0; map_d < 4; ++map_d) {
      reference_element::init_local_nnode_by_variant (k, ndof_on_subgeo_internal [map_d]);
      // clean upper-dimensional subgeos, since init_local_nnode_by_variant do not consider dimension
      for (size_type variant = reference_element::first_variant_by_dimension (map_d+1);
                     variant < reference_element::max_variant; ++variant) {
        ndof_on_subgeo_internal [map_d] [variant] = 0;
      }
    }
  }
  // when discontinuous, fix, but conserve the subgeo structure in _internal:
  base::_helper_make_discontinuous_ndof_on_subgeo (is_continuous, ndof_on_subgeo_internal, ndof_on_subgeo);
 
  // 2) deduce automatically first_idof_by_dimension:
  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo_internal, first_idof_by_dimension_internal);
  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo,          first_idof_by_dimension);
  // 3) Lagrange nodes follow dofs:
  nnod_on_subgeo_internal          = ndof_on_subgeo_internal;
  nnod_on_subgeo                   = ndof_on_subgeo;
  first_inod_by_dimension_internal = first_idof_by_dimension_internal;
  first_inod_by_dimension          = first_idof_by_dimension;
}
// =========================================================================
// basis members
// =========================================================================
template<class T>
basis_fem_Pk_lagrange<T>::~basis_fem_Pk_lagrange()
{
}
template<class T>
basis_fem_Pk_lagrange<T>::basis_fem_Pk_lagrange (size_type degree, const basis_option& sopt) 
  : basis_rep<T> (sopt),
    _raw_basis(),
    _hat_node(),
    _vdm(),
    _inv_vdm()
{
  base::_name = base::standard_naming (family_name(), degree, base::_sopt);
  string R = "?";
  switch (base::_sopt.get_raw_polynomial()) {
    case basis_option::monomial:  R = "M"; break;
    case basis_option::bernstein: R = "B"; break;
    case basis_option::dubiner:   R = "D"; break;
    default: error_macro ("unsupported polynomial: "<<sopt.get_raw_polynomial_name());
  }
  _raw_basis = basis_raw_basic<T> (R+std::to_string(degree));
  _initialize_cstor_sizes();
 
  // piola FEM transformation:
  typedef piola_fem_lagrange<T> piola_fem_type;
  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
}
template<class T>
bool
basis_fem_Pk_lagrange<T>::is_nodal() const
{
  return true;
}
template<class T>
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
basis_fem_Pk_lagrange<T>::hat_node (reference_element hat_K) const
{
  base::_initialize_data_guard (hat_K);
  return _hat_node [hat_K.variant()];
}
template<class T>
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
basis_fem_Pk_lagrange<T>::vdm (reference_element hat_K) const
{
  base::_initialize_data_guard (hat_K);
  return _vdm [hat_K.variant()];
}
template<class T>
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
basis_fem_Pk_lagrange<T>::inv_vdm (reference_element hat_K) const
{
  base::_initialize_data_guard (hat_K);
  return _inv_vdm [hat_K.variant()];
}
template<class T>
void
basis_fem_Pk_lagrange<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_lagrange<T>::_initialize_data (reference_element hat_K) const
{
  size_type k = degree();
  size_type variant = hat_K.variant();
 
  // nodes:
  switch (base::_sopt.get_node()) {
    case basis_option::equispaced:
          pointset_lagrange_equispaced (hat_K, k, _hat_node[variant]); break;
    case basis_option::warburton:
          pointset_lagrange_warburton  (hat_K, k, _hat_node[variant]); break;
    case basis_option::fekete:
          pointset_lagrange_fekete     (hat_K, k, _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, _hat_node[variant], _vdm[variant]);
  check_macro (invert(_vdm[variant], _inv_vdm[variant]),
        "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
}
// evaluation of all basis functions at hat_x:
template<class T>
void
basis_fem_Pk_lagrange<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> raw_value;
  _raw_basis.evaluate (hat_K, hat_x, raw_value);
  value = _inv_vdm[hat_K.variant()].transpose()*raw_value;
}
// evaluate the gradient:
template<class T>
void
basis_fem_Pk_lagrange<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);
  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm [hat_K.variant()];
  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> raw_v_value;
  _raw_basis.grad_evaluate (hat_K, hat_x, raw_v_value);
  size_type loc_ndof = raw_v_value.size(); 
  value.resize (loc_ndof);
  // TODO: point-valued blas2 : value := trans(inv_vdm)*raw_value
  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
    value [loc_idof] = point_basic<T>(0,0,0);
    for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
      value [loc_idof] += inv_vdm (loc_jdof,loc_idof)*raw_v_value[loc_jdof];
    }
  }
}
// extract local dof-indexes on a side
template<class T>
typename basis_fem_Pk_lagrange<T>::size_type
basis_fem_Pk_lagrange<T>::local_ndof_on_side (
  reference_element            hat_K,
  const side_information_type& sid) const
{
  return base::ndof (sid.hat);
}
template<class T>
void
basis_fem_Pk_lagrange<T>::local_idof_on_side (
  reference_element            hat_K,
  const side_information_type& sid,
  Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
{
  details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);
}
// dofs for a scalar-valued function
template<class T>
void
basis_fem_Pk_lagrange<T>::_compute_dofs (
  reference_element     hat_K,
  const Eigen::Matrix<T,Eigen::Dynamic,1>&   f_xnod, 
        Eigen::Matrix<T,Eigen::Dynamic,1>&   dof) const
{
  dof = f_xnod; // TODO: could avoid a physical copy: check it in compute_dof(f) with is_nodal() 
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_fem_Pk_lagrange<T>;
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef