piola_on_pointset.cc

Go to the documentation of this file.

//
#include "rheolef/piola_on_pointset.h"
#include "rheolef/piola_util.h"
 
namespace rheolef {
 
// ----------------------------------------------------------------------------
// initializers
// ----------------------------------------------------------------------------
template<class T>
void
piola_on_pointset_rep<T>::initialize (
    const basis_basic<T>&   piola_basis,
    const quadrature<T>&    quad,
    const integrate_option& iopt)
{
  _bops.set (quad, piola_basis);
  _ignore_sys_coord = iopt.ignore_sys_coord;
  _last_visited_geo.fill ("");
  _last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
}
template<class T>
void
piola_on_pointset_rep<T>::initialize (
    const basis_basic<T>&   piola_basis,
    const basis_basic<T>&   nodal_basis,
    const integrate_option& iopt)
{
  _bops.set (nodal_basis, piola_basis);
  _ignore_sys_coord = iopt.ignore_sys_coord;
  _last_visited_geo.fill ("");
  _last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
}
// ----------------------------------------------------------------------------
// piola transformation on a pointset
// ----------------------------------------------------------------------------
template<class T>
template<class M>
void
piola_on_pointset_rep<T>::_update (const geo_basic<T,M>& omega_K, const geo_element& K) const
{
  // memorize the last K call: avoid multiple _piola evaluation
  reference_element hat_K = K;
  if (_last_visited_geo    [hat_K.variant()] == omega_K.name() &&
      _last_visited_dis_ie [hat_K.variant()] == K.dis_ie()) {
    return;
  }
  _last_visited_geo    [hat_K.variant()] = omega_K.name();
  _last_visited_dis_ie [hat_K.variant()] = K.dis_ie();
 
  // here, compute piola and weight on K:
  Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola  = _piola  [hat_K.variant()];
  Eigen::Matrix<T,Eigen::Dynamic,1>&        weight = _weight [hat_K.variant()];
 
  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
      s_phij_xi = _bops.template evaluate<T> (hat_K);
  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
      grad_phij_xi = _bops.template grad_evaluate<point_basic<T>> (hat_K);
  size_type loc_nnod = s_phij_xi.rows();
  size_type loc_ndof = s_phij_xi.cols();
  size_type     d = omega_K.dimension();
  size_type map_d = hat_K.dimension();
  space_constant::coordinate_type sys_coord = omega_K.coordinate_system();
 
  // pass 0: reset
  piola .resize (loc_nnod);
  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
    piola[loc_inod].clear();
    piola[loc_inod].d         = d;
    piola[loc_inod].map_d     = map_d;
    piola[loc_inod].sys_coord = sys_coord;
    piola[loc_inod].ignore_sys_coord = _ignore_sys_coord;
  }
  // pass 1: F and DF
  omega_K.dis_inod (K, _dis_inod_K);
  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
    // dis_node: in outer loop: could require more time with external node
    const point_basic<T>& xjnod = omega_K.dis_node (_dis_inod_K[loc_jdof]);
    for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
      // step 1.1: F(hat_x)
      for (size_type alpha = 0; alpha < d; alpha++) {
        piola[loc_inod].F[alpha] += s_phij_xi (loc_inod,loc_jdof)*xjnod[alpha];
      }
      // step 1.2: DF(hat_x)
      cumul_otimes (piola[loc_inod].DF, xjnod, grad_phij_xi(loc_inod,loc_jdof), d, map_d);
 
      // step 1.3: surface projector
      if (map_d > 0 && map_d + 1 == d) {
        map_projector (piola[loc_inod].DF, d, map_d, piola[loc_inod].P);
      }
    }
  }
  // pass 2: invDF and detDF
  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
    piola[loc_inod].invDF = pseudo_inverse_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
    piola[loc_inod].detDF = det_jacobian_piola_transformation            (piola[loc_inod].DF, d, map_d);
  }
  // pass 3: quadrature weight, when pointset is a quadrature
  if (!_bops.has_quadrature()) return;
  weight.resize (loc_nnod);
  const quadrature<T>& quad = _bops.get_quadrature();
  size_type i_comp_axi = (omega_K.coordinate_system() == space_constant::axisymmetric_rz) ? 0 : 1;
  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
    T weight_q = quad (hat_K, loc_inod).w;
    weight[loc_inod] = piola[loc_inod].detDF*weight_q;
    if (omega_K.coordinate_system() != space_constant::cartesian && !_ignore_sys_coord) {
      weight[loc_inod] *= piola[loc_inod].F [i_comp_axi];
    }
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
 
#define _RHEOLEF_instanciation(T)                                               \
template class piola_on_pointset_rep<T>;                                        \
 
#define _RHEOLEF_instanciation_update(T,M)                                      \
template void piola_on_pointset_rep<T>::_update (                               \
  const geo_basic<T,M>& omega_K,                                                \
  const geo_element&    K) const;                                               \
 
_RHEOLEF_instanciation(Float)
_RHEOLEF_instanciation_update(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation_update(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef