geo_neighbour.cc

Go to the documentation of this file.

// Banded level set routines 
//
// Author: Pierre Saramito
//
#include "rheolef/geo.h"
 
namespace rheolef {
 
// ----------------------------------------------------------------------------
// 1) utilities for neighbour initialization
// ----------------------------------------------------------------------------
// Implementation note: as member function in geo_base_rep<T,M> could be simpler
// but function member of template class cannot be specialized
// and the sequential case do not define scatter_map_type and such
// so use a specialization of a non-member function
template<class T, class M>
static
void
add_ball_externals (
  const geo_base_rep<T,M>&      omega,
  const disarray<index_set,M>&  ball)
{
}
// explicit M=sequential specialization (otherwise: linker problems)
template<class T>
static
void
add_ball_externals (
  const geo_base_rep<T,sequential>&      omega,
  const disarray<index_set,sequential>&  ball)
{
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
void
add_ball_externals (
  const geo_base_rep<T,distributed>&      omega,
  const disarray<index_set,distributed>&  ball)
{
  typedef typename geo_base_rep<T,distributed>::size_type   size_type;
  index_set ext_iv_set;
  // faudrait attraper _geo_element[0].get_dis_map_entries() : suit les vertices au lieu nodes
  const hack_array<geo_element_hack,distributed>& vertex = omega._geo_element[0];
  typedef typename hack_array<geo_element_hack,distributed>::scatter_map_type map_t;
  const map_t& map_vertex = vertex.get_dis_map_entries();
  size_type dis_nv   = omega.dis_size (0);
  for (typename map_t::const_iterator
        iter = map_vertex.begin(),
        last = map_vertex.end(); iter != last; ++iter) {
    size_type dis_iv = (*iter).first;
    check_macro (dis_iv < dis_nv, "dis_iv="<<dis_iv<<" out of range [0:"<<dis_nv<<"[");
    ext_iv_set += dis_iv;
  }
  ball.set_dis_indexes (ext_iv_set);
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// 2) (lazy) initialization of neighbour informations
// ----------------------------------------------------------------------------
template<class T, class M>
void
geo_base_rep<T,M>::init_neighbour() const
{
  if (is_broken()) return; // skip broken domain e.g. "sides" or "internal_sides": no neighbours
  size_type map_dim = map_dimension();
  if (map_dim <= 0) return; // a set of points have no neighbours
  size_type sid_dim = map_dim-1;
 
  // 0.0) clean-up neighbour information (could come from a copy of elements, from domain, e.g. band)
  for (const_iterator iter_S = begin(sid_dim), last_S = end(sid_dim); iter_S != last_S; ++iter_S) {
    const geo_element& S = *iter_S;
    S.set_master (0, std::numeric_limits<size_type>::max());
    S.set_master (1, std::numeric_limits<size_type>::max());
  }
  // 0.1) manage all external sides S that appears on some local elements K
  distributor is_ownership = _gs.ownership_by_dimension[map_dim-1];
  std::array<index_set, reference_element::max_variant> ext_isv_set; // external sides, by variants
  index_set ext_is_set; // external sides, all variants merged
  for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
    const geo_element& K = *iter_K;
    for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
      size_type dis_is = (map_dim == 1) ? K[loc_is] : ((map_dim == 2) ? K.edge(loc_is) : K.face(loc_is));
      if (is_ownership.is_owned (dis_is)) continue;
      // have an external side S: convert its dis_is to dis_isv and get its variant
      size_type variant;
      size_type dis_isv = _gs.dis_ige2dis_igev_by_dimension (sid_dim, dis_is, variant);
      ext_isv_set [variant] += dis_isv;
      ext_is_set            += dis_is;
    }
  }
  // for external side S that appears in a local element K may be available
  // => enlarge the external set of _geo_element for sid_dim
  // but ext_is_set may be classified by variants and all dis_is may be
  // converted to dis_isv
  for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
                 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
    _geo_element [variant].append_dis_indexes (ext_isv_set [variant]);
  }
  // 1) build ball(xi) := { K; xi is a vertex of K }
  index_set emptyset;
  distributor vertex_ownership  = sizes().ownership_by_dimension[0];
  disarray<index_set,M> ball (vertex_ownership, emptyset);
  for (const_iterator iter = begin(map_dim), last = end(map_dim); iter != last; ++iter) {
    const geo_element& K = *iter;
    index_set dis_ie_set;
    dis_ie_set += K.dis_ie();
    for (size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
      size_type dis_inod = K [loc_inod];
      size_type dis_iv = dis_inod2dis_iv (dis_inod);
      ball.dis_entry (dis_iv) += dis_ie_set; // union with {dis_ie}
    }
  }
  ball.dis_entry_assembly();
  // 2) set ball available at partition boundaries
  add_ball_externals (*this, ball);
  // 3) set master elements on all sides
  distributor   ie_ownership = _gs.ownership_by_dimension[map_dim];
  distributor inod_ownership = _gs.node_ownership;
  size_type first_dis_ie = ie_ownership.first_index();
  std::array<index_set, reference_element::max_variant> ext_igev_set; // external elts, by variants
  index_set ext_ie_set; // external elements, all variants merged
  index_set ext_inod_set;
  disarray<size_type, M> side_marked (is_ownership, 0); // boolean: 0,1 
  side_marked.set_dis_indexes (ext_is_set);
  for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
    const geo_element& K = *iter_K;
    for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
      size_type dis_is = (sid_dim == 0) ? K[loc_is] : ((sid_dim == 1) ? K.edge(loc_is) : K.face(loc_is));
      if (side_marked.dis_at (dis_is)) continue;
      side_marked.dis_entry (dis_is) = 1;
      const geo_element& S = dis_get_geo_element (sid_dim, dis_is);
      size_type dis_iv0 = dis_inod2dis_iv (S[0]);
      index_set ball_intersect = ball.dis_at (dis_iv0);
      for (size_type loc_iv = 1, loc_nv = S.size(); loc_iv < loc_nv; loc_iv++) {
        size_type dis_iv = dis_inod2dis_iv (S[loc_iv]);
        ball_intersect.inplace_intersection (ball.dis_at (dis_iv));
      }
      switch (ball_intersect.size()) {
        case 1: { // one master K: S is on the boundary
          index_set::const_iterator iter = ball_intersect.begin();
          size_type dis_ie = *iter;
          S.set_master (0, dis_ie);
          S.set_master (1, std::numeric_limits<size_type>::max());
          break;
        }
        case 2: { // two masters K1, K2: S is an internal side
          index_set::const_iterator iter = ball_intersect.begin();
          size_type dis_ie1 = *iter++;
          size_type dis_ie2 = *iter;
          // one of K1 and K2 is available in distributed environment
          // the other is not a priori available
          if (! ie_ownership.is_owned (dis_ie1)) std::swap (dis_ie1, dis_ie2);
          if (! ie_ownership.is_owned (dis_ie2)) {
            // convert dis_ie2 to dis_igev2 and get its variant2:
            size_type variant2;
            size_type dis_igev2 = _gs.dis_ige2dis_igev_by_dimension (map_dim, dis_ie2, variant2);
            ext_igev_set [variant2] += dis_igev2;
            ext_ie_set              += dis_ie2;
          }
          check_macro (ie_ownership.is_owned (dis_ie1), "unexpected external K1");
          // let K1 indicated by the normal on S given by S orient
          // let K2, the other one
          size_type ie1 = dis_ie1 - first_dis_ie;
          const geo_element& K1 = get_geo_element (map_dim, ie1);
          geo_element::orientation_type orient = K1.get_side_orientation (S);
          if (orient < 0) std::swap (dis_ie1, dis_ie2);
          S.set_master (0, dis_ie1);
          S.set_master (1, dis_ie2);
          break;
        }
        default: {
          error_macro ("neighbour: side #"<< ball_intersect.size()<<" connectivity problem");
        }
      }
    }
  }
  // side values have been changed:
  // propagate these modifs for sides at the partition boundaries
  for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
                 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
    _geo_element [variant].update_dis_entries();
  }
  // for all S in the partition boundary, both K1 & K2 may be available
  // => enlarge the external set of _geo_element for map_dim
  // but ext_ie_set may be classified by variants and all dis_ie may be
  // converted to dis_igev
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
    _geo_element [variant].append_dis_indexes (ext_igev_set [variant]);
  }
  // manage also external nodes of external elements K:
  for (index_set::const_iterator iter = ext_ie_set.begin(), last = ext_ie_set.end(); iter != last; ++iter) {
    size_type dis_ie = *iter;
    const geo_element& K = dis_get_geo_element (map_dim, dis_ie);
    for (size_type loc_inod = 0, loc_nnod = K.n_node(); loc_inod < loc_nnod; loc_inod++) {
      size_type dis_inod = K[loc_inod];
      if (! inod_ownership.is_owned (dis_inod)) {
        ext_inod_set += dis_inod;
      }
    }
  }
  _node.append_dis_indexes (ext_inod_set);
}
// ----------------------------------------------------------------------------
// 3) compute the neighbour K2 of an element K1 across a side S
// ----------------------------------------------------------------------------
// returns the dis_ie2 global index associated to K2, 
// the neighbour element of K1, with index ie1,
// across the side S of local index loc_isid in K1
// when S is on the boundary, returns size_t(-1)
template<class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::neighbour (size_type ie1, size_type loc_isid) const
{
  check_macro (_have_neighbour, "neighbour_guard() may be called globally before local neighbour()");
  // 1) get the i-th side S in K1
  size_type map_dim = map_dimension();
  check_macro (ie1 < size(map_dim), "neighbour: element index ie=" << ie1
        << " is out of range [0:" << size(map_dim) << "[");
  const geo_element& K1 = get_geo_element (map_dim, ie1);
  size_type dis_isid;
  switch (map_dim) {
    case 3: dis_isid = K1.face (loc_isid); break;
    case 2: dis_isid = K1.edge (loc_isid); break;
    case 1: {
        size_type dis_inod = K1[loc_isid];
        size_type dis_iv = dis_inod2dis_iv (dis_inod);
        dis_isid = dis_iv;
        break;
    }
    case 0:
    default: return std::numeric_limits<size_type>::max();
  }
  size_type sid_dim = map_dim-1;
  const geo_element& S = dis_get_geo_element (sid_dim, dis_isid);
 
  // 2) get the neighbour element K2 to K1 across S
  size_type dis_ie2 = S.master (0);
  size_type first_dis_ie = _gs.ownership_by_dimension[map_dim].first_index();
  size_type dis_ie1 = first_dis_ie + ie1;
  if (dis_ie2 != dis_ie1) return dis_ie2;
  return S.master(1);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M)                                             \
template void geo_base_rep<T,M>::init_neighbour() const;                        \
template geo_base_rep<T,M>::size_type                                           \
        geo_base_rep<T,M>::neighbour (size_type, size_type) const;
 
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace