space_numbering.cc

Go to the documentation of this file.

#include "rheolef/space_numbering.h"
#include "rheolef/geo.h"
#include "rheolef/rheostream.h"
namespace rheolef { namespace space_numbering {
 
template <class T>
size_type
ndof (
  const basis_basic<T>&  b,
  const geo_size&        gs, 
  size_type              map_dim)
{
  size_type ndof = 0;
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
    ndof += gs.ownership_by_variant [variant].size() * b.ndof_on_subgeo (map_dim, variant);
  }
  return ndof;
}
template <class T>
size_type
dis_ndof (
  const basis_basic<T>&  b,
  const geo_size&        gs, 
  size_type              map_dim)
{
  size_type dis_ndof = 0;
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
    dis_ndof += gs.ownership_by_variant [variant].dis_size() * b.ndof_on_subgeo (map_dim, variant);
  }
  return dis_ndof;
}
template <class T>
size_type
nnod (
  const basis_basic<T>&  b,
  const geo_size&        gs, 
  size_type              map_dim)
{
  size_type nnod = 0;
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
    nnod += gs.ownership_by_variant [variant].size() * b.nnod_on_subgeo (map_dim, variant);
  }
  return nnod;
}
template <class T>
size_type
dis_nnod (
  const basis_basic<T>&  b,
  const geo_size&        gs, 
  size_type              map_dim)
{
  size_type dis_nnod = 0;
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
    dis_nnod += gs.ownership_by_variant [variant].dis_size() * b.nnod_on_subgeo (map_dim, variant);
  }
  return dis_nnod;
}
// --------------------------------------------------------------------------
// dis_inod/dis_idof shared internal function
// --------------------------------------------------------------------------
template <class Callback>
static
void
dis_ixxx (
  const Callback&       nxxx_on_subgeo,
  const geo_size&       gs, 
  const geo_element&    K, 
  typename std::vector<size_t>::iterator dis_ixxx_tab)
{
  typedef size_t size_type;
  reference_element hat_K = K;
  size_type K_map_dim = hat_K.dimension();
  size_type n_comp    = nxxx_on_subgeo.n_component();
  size_type first_loc_ixxx = 0;
  for (size_type subgeo_dim = 0; subgeo_dim <= K_map_dim; ++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) {
      if (subgeo_dim == K_map_dim && subgeo_variant != hat_K.variant()) continue;
      reference_element hat_S (subgeo_variant);
      size_type loc_ngev = hat_K.n_subgeo_by_variant (hat_S.variant());
      if (loc_ngev == 0) continue;
      size_type loc_nxxx_on_subgeo = nxxx_on_subgeo (hat_S.variant());
      size_type loc_comp_nxxx_on_subgeo = loc_nxxx_on_subgeo / n_comp;
      // => reduce the basis interface
      if (loc_nxxx_on_subgeo == 0) continue;
      // loop on all subgeos S like hat_S and that belongs to K
      for (size_type loc_ige = 0, loc_nge = hat_K.n_subgeo (hat_S.dimension()); loc_ige < loc_nge; ++ loc_ige) {
        reference_element hat_S2 = hat_K.subgeo (hat_S.dimension(), loc_ige);
        if (hat_S2.variant() != hat_S.variant()) continue;
        // 0) compute loc_igev
        size_type loc_igev = hat_K.local_subgeo_index2index_by_variant (hat_S.variant(), loc_ige);
        // 1) compute dis_igev
        size_type dis_ige = K.subgeo_dis_index (hat_S.dimension(), loc_ige);
        if (hat_S.dimension() == 0) { dis_ige = gs.dis_inod2dis_iv (dis_ige); }
        size_type dis_igev = gs.dis_ige2dis_igev_by_variant (hat_S.variant(), dis_ige);
        // 2) compute dis_ixxx_loc_start
        size_type dis_ixxx_loc_start = 0;
        size_type iproc = gs.ownership_by_variant [hat_S.variant()].find_owner (dis_igev);
        for (size_type prev_subgeo_variant = 0;
                       prev_subgeo_variant < hat_S.variant(); 
                       prev_subgeo_variant++) {
          dis_ixxx_loc_start
          += gs.ownership_by_variant [prev_subgeo_variant]. last_index(iproc)
           * nxxx_on_subgeo (prev_subgeo_variant);
        }
        dis_ixxx_loc_start += dis_igev*nxxx_on_subgeo (hat_S.variant());
        for (size_type next_subgeo_variant = hat_S.variant()+1;
                       next_subgeo_variant < reference_element::max_variant;
                       next_subgeo_variant++) {
          dis_ixxx_loc_start
          += gs.ownership_by_variant [next_subgeo_variant].first_index(iproc)
           * nxxx_on_subgeo (next_subgeo_variant);
        }
        // 3) set all ixxxs
        size_type loc_ixxx_loc_start = first_loc_ixxx + loc_igev*loc_nxxx_on_subgeo;
        for (size_type loc_comp_ixxx_on_subgeo = 0; loc_comp_ixxx_on_subgeo < loc_comp_nxxx_on_subgeo; ++loc_comp_ixxx_on_subgeo) {
          size_type loc_comp_ixxx_on_subgeo_fixed = geo_element::fix_indirect (K, hat_S.variant(), loc_ige, loc_comp_ixxx_on_subgeo, nxxx_on_subgeo.degree()); 
          for (size_type i_comp = 0; i_comp < n_comp; ++i_comp) {
            size_type loc_ixxx      = loc_ixxx_loc_start + loc_comp_ixxx_on_subgeo      *n_comp + i_comp;
            dis_ixxx_tab [loc_ixxx] = dis_ixxx_loc_start + loc_comp_ixxx_on_subgeo_fixed*n_comp + i_comp;
          }
        }
      }
      first_loc_ixxx += loc_ngev*loc_nxxx_on_subgeo;
    }
  }
}
// --------------------------------------------------------------------------
// dis_inod/dis_idof call-back
// --------------------------------------------------------------------------
// codes for these functions are very similar
// thus create a call-back that allows one to share one main code between them
template<class T>
struct nnod_callback {
  nnod_callback (const basis_basic<T>& b, size_t map_d) 
  : _b(b), _map_d(map_d) {}
  size_t operator() (size_t variant) const {
    return _b.nnod_on_subgeo (_map_d, variant); }
  size_t n_component() const { return 1; }
  size_t degree() const { return _b.degree(); }
  basis_basic<T> _b;
  size_t         _map_d;
};
template<class T>
struct ndof_callback {
  ndof_callback (const basis_basic<T>& b, size_t map_d) 
  : _b(b), _map_d(map_d) {}
  size_t operator() (size_t variant) const {
    return _b.ndof_on_subgeo (_map_d, variant); }
  size_t n_component() const { return _b.size(); }
  size_t degree() const { return _b.degree(); }
  basis_basic<T> _b;
  size_t         _map_d;
};
template <class T>
void
dis_inod (
  const basis_basic<T>& b,
  const geo_size&       gs, 
  const geo_element&    K, 
  typename std::vector<size_type>::iterator dis_inod_tab)
{
  dis_ixxx (nnod_callback<T> (b, gs.map_dimension()), gs, K, dis_inod_tab);
}
template <class T>
void
dis_idof (
  const basis_basic<T>& b,
  const geo_size&       gs, 
  const geo_element&    K, 
  typename std::vector<size_type>::iterator dis_idof_tab)
{
  dis_ixxx (ndof_callback<T> (b, gs.map_dimension()), gs, K, dis_idof_tab);
}
template <class T>
void
dis_inod (
  const basis_basic<T>& b,
  const geo_size&       gs, 
  const geo_element&    K, 
  std::vector<size_type>& dis_inod_tab)
{
  dis_inod_tab.resize (b.nnod (K.variant()));
#ifdef _RHEOLEF_PARANO
  std::fill (dis_inod_tab.begin(), dis_inod_tab.end(), std::numeric_limits<size_type>::max());
#endif // _RHEOLEF_PARANO
  dis_inod (b, gs, K, dis_inod_tab.begin());
}
template <class T>
void
dis_idof (
  const basis_basic<T>& b,
  const geo_size&       gs, 
  const geo_element&    K, 
  std::vector<size_type>& dis_idof_tab)
{
  dis_idof_tab.resize (b.ndof (K.variant()));
#ifdef _RHEOLEF_PARANO
  std::fill (dis_idof_tab.begin(), dis_idof_tab.end(), std::numeric_limits<size_type>::max());
#endif // _RHEOLEF_PARANO
  dis_idof (b, gs, K, dis_idof_tab.begin());
}
#ifdef _RHEOLEF_HAVE_MPI
// -------------------------------------------------------------
// permut to/from ios dof numbering, for distributed environment
// --------------------------------------------------------------------------
// set permutation for idofs: ios_idof2dis_idof & idof2ios_dis_idof
// and redistribute ios_dof[] into _dof[]
// applied for initialiing distributed mesh:
//       posibly curved with Pk piola basis.
// --------------------------------------------------------------------------
template <class T>
void
generic_set_ios_permutation (
  const basis_basic<T>&         b,
  size_t                        map_d,
  const geo_size&               gs,
  const std::array<disarray<size_t,distributed>,reference_element::max_variant>&
                                igev2ios_dis_igev,
  disarray<size_t,distributed>& idof2ios_dis_idof)
{
  typedef geo_element::size_type size_type;
  // 1) count & resize
  size_type ndof = 0;
  for (size_type variant = 0,
                 variant_last = reference_element::last_variant_by_dimension(map_d);
                 variant < variant_last;
                 variant++) {
    ndof += gs.ownership_by_variant [variant].size() * b.ndof_on_subgeo (map_d, variant);
  }
  communicator comm = gs.ownership_by_variant[0].comm();
  distributor dof_ownership (distributor::decide, comm, ndof);
  idof2ios_dis_idof.resize (dof_ownership);
  // 2) set idof2ios_dis_idof[]
  size_type first_dis_idof= dof_ownership.first_index();
  size_type idof = 0;
  size_type first_ios_dis_v = 0;
  for (size_type variant = 0,
                 variant_last = reference_element::last_variant_by_dimension(map_d);
                 variant < variant_last;
                 variant++) {
    size_type loc_ndof = b.ndof_on_subgeo (map_d, variant);
    for (size_type igev = 0, ngev = gs.ownership_by_variant [variant].size(); igev < ngev; igev++) {
      size_type ios_dis_igev = igev2ios_dis_igev [variant] [igev];
      for (size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++, idof++) {
        size_type ios_dis_idof = first_ios_dis_v + ios_dis_igev*loc_ndof + loc_idof;
        idof2ios_dis_idof [idof] = ios_dis_idof;
      }
    }
    first_ios_dis_v += gs.ownership_by_variant [variant].dis_size() * loc_ndof;
  }
}
template <class T>
void
set_ios_permutations (
    const basis_basic<T>&              b,
    const geo_basic<T,distributed>&    omega,
    disarray<size_type,distributed>&   idof2ios_dis_idof,
    disarray<size_type,distributed>&   ios_idof2dis_idof)
{
  generic_set_ios_permutation (b, omega.map_dimension(), omega.sizes(), omega.get_igev2ios_dis_igev(), idof2ios_dis_idof);
  size_type dis_ndof = idof2ios_dis_idof.ownership().dis_size();
  distributor ios_dof_ownership (dis_ndof, idof2ios_dis_idof.comm(), distributor::decide);
  ios_idof2dis_idof.resize (ios_dof_ownership, std::numeric_limits<size_type>::max());
  idof2ios_dis_idof.reverse_permutation (ios_idof2dis_idof);
}
#endif // _RHEOLEF_HAVE_MPI
 
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciate(T)                                                 \
template size_type ndof (                                                       \
  const basis_basic<T>&  b,                                                     \
  const geo_size&        gs,                                                    \
  size_type              map_dim);                                              \
template size_type nnod (                                                       \
  const basis_basic<T>&  b,                                                     \
  const geo_size&        gs,                                                    \
  size_type              map_dim);                                              \
template size_type dis_ndof (                                                   \
  const basis_basic<T>&  b,                                                     \
  const geo_size&        gs,                                                    \
  size_type              map_dim);                                              \
template size_type dis_nnod (                                                   \
  const basis_basic<T>&  b,                                                     \
  const geo_size&        gs,                                                    \
  size_type              map_dim);                                              \
template void dis_idof (                                                        \
  const basis_basic<T>& b,                                                      \
  const geo_size&       gs,                                                     \
  const geo_element&    K,                                                      \
  typename std::vector<size_type>::iterator dis_inod_tab);                      \
template void dis_inod (                                                        \
  const basis_basic<T>& b,                                                      \
  const geo_size&       gs,                                                     \
  const geo_element&    K,                                                      \
  typename std::vector<size_type>::iterator dis_inod_tab);                      \
template void dis_idof (                                                        \
  const basis_basic<T>& b,                                                      \
  const geo_size&       gs,                                                     \
  const geo_element&    K,                                                      \
  std::vector<size_type>& dis_inod_tab);                                        \
template void dis_inod (                                                        \
  const basis_basic<T>& b,                                                      \
  const geo_size&       gs,                                                     \
  const geo_element&    K,                                                      \
  std::vector<size_type>& dis_inod_tab);                                        \
 
_RHEOLEF_instanciate(Float)
 
#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciate_distributed(T)                                     \
template void generic_set_ios_permutation (                                     \
  const basis_basic<T>&         b,                                              \
  size_t                        map_d,                                          \
  const geo_size&               gs,                                             \
  const std::array<disarray<size_t,distributed>,reference_element::max_variant>& \
                                igev2ios_dis_igev,                              \
  disarray<size_t,distributed>& idof2ios_dis_idof);                             \
template void set_ios_permutations (                                            \
    const basis_basic<T>&              b,                                       \
    const geo_basic<T,distributed>&    omega,                                   \
    disarray<size_type,distributed>&   idof2ios_dis_idof,                       \
    disarray<size_type,distributed>&   ios_idof2dis_idof);                      \
 
_RHEOLEF_instanciate_distributed(Float)
#endif // _RHEOLEF_HAVE_MPI
 
}} // namespace rheolef::space_numbering