basis_rep.cc

Go to the documentation of this file.

#include "rheolef/basis.h"
#include "rheolef/basis_get.h"
 
namespace rheolef {
using namespace std;
 
// =========================================================================
// naming scheme for standard FEM families
// =========================================================================
static
bool
is_family_no_index (const std::string& family_name)
{
  // TODO: how to merge this keywords list with basis_lex.flex ?
  // basis::have_index_parameter() givzes the response, but requires a
  //  previously builded basis
  if (family_name == "bubble") return true;
  if (family_name == "P1qd")   return true;
  if (family_name == "empty")  return true;
  return false;
}
template <class T>
std::string
basis_rep<T>::standard_naming (
  std::string         family_name,
  size_type           index,
  const basis_option& sopt)
{
  string opt_d = (sopt.is_continuous() || (family_name == "P" && index == 0)) ? "" : "d";
  string basename    = is_family_no_index (family_name) ? family_name : family_name + std::to_string(index) + opt_d;
  string name_opt = basename + sopt.stamp();
  if (sopt.valued_tag() == space_constant::scalar) {
    trace_macro("standard_naming("<<family_name<<","<<index<<",opts)=\""<<name_opt<<"\"");
    if (! sopt.is_trace_n()) {
      return name_opt;
    } else {
      return "trace_n(" + name_opt + ")";
    }
  }
  // here: non-scalar basis
  const size_type unset = std::numeric_limits<basis_option::size_type>::max();
  string opt_dim = (sopt.dimension() == unset) ? "" : "d=" + std::to_string(sopt.dimension());
  string opt_sc  = (sopt.valued_tag() == space_constant::vector || sopt.coordinate_system() == space_constant::cartesian) ?
                 "" : space_constant::coordinate_system_name(sopt.coordinate_system());
  string coma = (opt_dim == "" || opt_sc == "") ? "" : ",";
  string opt_valued = (opt_dim == "" && opt_sc == "") ? "" : "[" + opt_dim + coma + opt_sc + "]";
  string braced_name = sopt.valued() + opt_valued + "(" + name_opt + ")";
  if (! sopt.is_trace_n()) {
    return braced_name;
  } else {
    return "trace_n(" + braced_name + ")";
  }
}
// =========================================================================
// basis members (should be inlined)
// =========================================================================
template<class T>
void
basis_basic<T>::reset (std::string& name_in)
{
  if (name_in == "") {
    base::operator= (0);
  } else {
    // strip name_in as name: drop "scalar(Pk)" as "Pk" or "P0d" as "P0", to get unique name scheme
    family_index_option_type fio;
    basis_parse_from_string (name_in, fio);
    std::string name = basis_rep<T>::standard_naming (fio.family, fio.index, fio.option);
    base::operator= (persistent_table<basis_basic<T>>::load (name));
    _clear();
  }
}
template<class T>
void
basis_basic<T>::reset_family_index (size_type k)
{
  std::string fname = family_name();
  basis_option bopt = option();
  std::string name = basis_rep<T>::standard_naming (fname, k, bopt); \
  reset (name);
}
template <class T>
basis_rep<T>::~basis_rep()
{
  persistent_table<basis_basic<T>>::unload (_name);
}
template<class T>
basis_rep<T>::basis_rep (const basis_option& sopt)
: _name(""),
  _sopt(sopt),
  _piola_fem(),
  _have_initialize_data(),
  _ndof_on_subgeo_internal(),
  _ndof_on_subgeo(),
  _nnod_on_subgeo_internal(),
  _nnod_on_subgeo(),
  _first_idof_by_dimension_internal(),
  _first_idof_by_dimension(),
  _first_inod_by_dimension_internal(),
  _first_inod_by_dimension()
{
  _clear();
}
template<class T>
void
basis_rep<T>::_clear() const
{
  _have_initialize_data.fill (false);
}
template <class T>
void
basis_rep<T>::_initialize_data_guard (reference_element hat_K) const
{
  if (_have_initialize_data [hat_K.variant()]) return;
  _have_initialize_data [hat_K.variant()] = true;
  _initialize_data (hat_K);
}
// =========================================================================
// dof & nof: hierarchial organization by increasing subgeo dimension
// =========================================================================
// inplace change nxxx_on_subgeo for discontinuous elements
template <class T>
void
basis_rep<T>::_helper_make_discontinuous_ndof_on_subgeo (
  bool is_continuous,
  const
  std::array<
    std::array<
      size_type
     ,reference_element::max_variant>
    ,4>&                                    nxxx_on_subgeo_internal,
  std::array<
    std::array<
      size_type
     ,reference_element::max_variant>
    ,4>&                                    nxxx_on_subgeo)
{
  if (is_continuous) {
    nxxx_on_subgeo = nxxx_on_subgeo_internal;
    return; // no changes
  }
  for (size_type map_dim = 0; map_dim < 4; ++map_dim) {
    nxxx_on_subgeo [map_dim].fill (0);
  }
  for (size_type variant = 0;
                 variant < reference_element::max_variant;
                 variant++)
  {
    reference_element hat_K (variant);
    size_type map_dim = hat_K.dimension();
    size_type sum = 0;
    for (size_type subgeo_variant = 0;
                   subgeo_variant < reference_element:: last_variant_by_dimension(map_dim);
                   subgeo_variant++)
    {
      size_type n_subgeo = hat_K.n_subgeo_by_variant (subgeo_variant);
      sum += n_subgeo*nxxx_on_subgeo_internal [map_dim][subgeo_variant];
    }
    nxxx_on_subgeo [map_dim][variant] = sum;
  }
}
// deduce automatically first_ixxx_by_dimension from nxxx_on_subgeo
template <class T>
void
basis_rep<T>::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (
  const std::array<
    std::array<
      size_type
     ,reference_element::max_variant>
    ,4>&                                    nxxx_on_subgeo,
  std::array<
    std::array<
      size_type
      ,5>
    ,reference_element::max_variant>&       first_ixxx_by_dimension)
{
  for (size_type variant = 0;
                 variant < reference_element::max_variant;
                 variant++)
  {
    reference_element hat_K (variant);
    size_type map_dim = hat_K.dimension();
    first_ixxx_by_dimension [variant].fill(0);
    for (size_type subgeo_dim = 0; subgeo_dim <= map_dim; ++subgeo_dim) {
      size_type sum = first_ixxx_by_dimension [variant][subgeo_dim];
      for (size_type subgeo_variant = reference_element::first_variant_by_dimension(subgeo_dim);
                     subgeo_variant < reference_element:: last_variant_by_dimension(subgeo_dim);
                     subgeo_variant++)
      {
        size_type n_subgeo = hat_K.n_subgeo_by_variant (subgeo_variant);
        sum += n_subgeo*nxxx_on_subgeo [map_dim][subgeo_variant];
      }
      first_ixxx_by_dimension [variant][subgeo_dim+1] = sum;
    }
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template class basis_rep<T>;                                                    \
template void basis_basic<T>::reset (std::string&);                             \
template void basis_basic<T>::reset_family_index (size_type);                   \
 
_RHEOLEF_instanciation(Float)
 
}// namespace rheolef