geo_mpi_get.cc

Go to the documentation of this file.

#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MPI
#include "rheolef/geo.h"
#include "rheolef/geo_domain.h"
#include "rheolef/space_numbering.h"
#include "rheolef/dis_macros.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/index_set.h"
 
namespace rheolef {
 
extern
disarray<size_t>
geo_mpi_partition (
  const std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                                  ios_geo_element,
  const distributor&              ownership_by_dimension,
  size_t                          map_dim,
  size_t                          dis_nv);
 
// --------------------------------------------------------------------------
// utility for geo get
// --------------------------------------------------------------------------
static
distributor
build_true_ios_ge_ownership_by_dimension (
    const std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                                               ios_geo_element,
          size_t                               side_dim)
{
  using namespace std;
  typedef geo_element::size_type size_type;
  size_type first_variant = reference_element::first_variant_by_dimension(side_dim);
  size_type  last_variant = reference_element:: last_variant_by_dimension(side_dim);
  if (first_variant >= last_variant) return distributor();
  const communicator& comm = ios_geo_element [first_variant].ownership().comm();
 
  size_type nge = 0;
  size_type dis_nge = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
    nge     += ios_geo_element [variant].size();
    dis_nge += ios_geo_element [variant].dis_size();
  }
  return distributor (dis_nge, comm, nge);
}
static
void
build_apparent_ios_ownership(
// input:
  const std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                                  ios_geo_element,
  size_t                          side_dim,
// output:
  distributor&                    ios_ownership,
  distributor                     ios_ownership_by_variant [reference_element::max_variant])
{
  using namespace std;
  typedef geo_element::size_type size_type;
  size_type first_variant = reference_element::first_variant_by_dimension(side_dim);
  size_type  last_variant = reference_element:: last_variant_by_dimension(side_dim);
  if (first_variant >= last_variant) return;
  const communicator& comm = ios_geo_element [first_variant].ownership().comm();
 
  // 1) build ios_ownership = ios_geo_element ownership by dimension
  size_type dis_nge = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
    dis_nge += ios_geo_element [variant].dis_size();
  }
  ios_ownership = distributor (dis_nge, comm, distributor::decide);
 
  // 2) build ios_ownership_by_variant as if it was readed in sequence, by dimension
  // Note: used by P1d and P2 element numbering, for output of results in the same
  //       vertex, edges, faces & volume as these entities was readed
  size_type first_dis_v = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
 
    size_type dis_ngev = ios_geo_element [variant].dis_size();
    size_type last_dis_v = first_dis_v + dis_ngev;
    size_type first_ios_dis_nge = ios_ownership.first_index();
    size_type  last_ios_dis_nge = ios_ownership. last_index();
    size_type first_ios_dis_ngev = min (max (first_ios_dis_nge, first_dis_v),  last_ios_dis_nge);
    size_type  last_ios_dis_ngev = max (min ( last_ios_dis_nge,  last_dis_v), first_ios_dis_nge);
    size_type ios_ngev = last_ios_dis_ngev - first_ios_dis_ngev;
    ios_ownership_by_variant [variant] = distributor (dis_ngev, comm, ios_ngev);
    first_dis_v = last_dis_v;
  }
}  
// ----------------------------------------------------------------------------
// new renumbering:
// - a vertex belongs to the same proc as at least an  edge that contains this vertex
// - an edge  belongs to the same proc as at least a   face that contains this edge
// - a face   belongs to the same proc as at least a volume that contains this face
// => a side S belongs to the same proc as at least a super-side K that contains it
// ----------------------------------------------------------------------------
void
geo_element_renumbering_part1_new (
// input:
  const std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                               ios_geo_element,
  const geo_size&              ios_gs, // TODO: not used...
  size_t                       S_dim,
// output:
  std::array<std::vector<size_t>, 4>&                           massive_partition_by_dimension,
  std::array<disarray<size_t>, reference_element::max_variant>& partition_by_variant)
{
  typedef size_t size_type;
  distributor ios_vertex_ownership = ios_geo_element [reference_element::p].ownership();
  size_type dis_nv           = ios_vertex_ownership.dis_size();
  size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
  // ------------------------------------------------------------------------
  // 1) K_ball(X) := { K; X is a vertex of K }
  // ------------------------------------------------------------------------
  index_set empty_set; // TODO: add a global allocator to empty_set
  disarray<index_set,distributed> K_ball (ios_vertex_ownership, empty_set);
  size_type K_dim = S_dim + 1;
  size_type K_ios_ige = 0;
  distributor true_K_ios_ge_ownership = build_true_ios_ge_ownership_by_dimension (ios_geo_element, K_dim);
  size_type first_K_ios_dis_ige = true_K_ios_ge_ownership.first_index();
  for (size_type K_variant = reference_element::first_variant_by_dimension(K_dim);
                 K_variant < reference_element:: last_variant_by_dimension(K_dim); K_variant++) {
    distributor K_ios_gev_ownership  = ios_geo_element [K_variant].ownership();
    for (size_type K_ios_igev = 0, K_ios_ngev = K_ios_gev_ownership.size(); K_ios_igev < K_ios_ngev; K_ios_igev++, K_ios_ige++) {
      const geo_element& K = ios_geo_element [K_variant] [K_ios_igev];
      size_type K_ios_dis_ige  = first_K_ios_dis_ige  + K_ios_ige;
      for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
        size_type ios_dis_iv = K[iloc];
        assert_macro (ios_dis_iv < dis_nv, "K={"<<K<<"}: "<<iloc<<"-vertex index " << ios_dis_iv << " out of range [0:"<<dis_nv<<"[");
        if (ios_vertex_ownership.is_owned(ios_dis_iv)) {
          size_type ios_iv = ios_dis_iv - first_ios_dis_iv;
          K_ball[ios_iv] += K_ios_dis_ige;
        } else {
          index_set K_ios_dis_ige_set;
          K_ios_dis_ige_set += K_ios_dis_ige;
          K_ball.dis_entry (ios_dis_iv) += K_ios_dis_ige_set; // not so efficient: union with {dis_iegv}
        }
      }
    }
  }
  K_ball.dis_entry_assembly();
  // ------------------------------------------------------------------------
  // 2) import ball[ios_dis_iv] when dis_iv is referenced by a side S
  // ------------------------------------------------------------------------
  index_set ext_ios_dis_iv;
  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
                 S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
    distributor ios_gev_ownership  = ios_geo_element [S_variant].ownership();
    for (size_type ios_igev = 0, ios_ngev = ios_gev_ownership.size(); ios_igev < ios_ngev; ios_igev++) {
      const geo_element& S = ios_geo_element [S_variant] [ios_igev];
      for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
        size_type ios_dis_iv = S[iloc];
        if (! ios_vertex_ownership.is_owned (ios_dis_iv)) {
          ext_ios_dis_iv += ios_dis_iv;
        }
      }
    }
  }
  K_ball.set_dis_indexes (ext_ios_dis_iv);
  // ------------------------------------------------------------------------
  // 3) owner(S) := max { owner(K);  K is a super-subgeo containing S }
  // ------------------------------------------------------------------------
  distributor true_S_ios_ge_ownership = build_true_ios_ge_ownership_by_dimension (ios_geo_element, S_dim);
  size_type S_dis_nge           = true_S_ios_ge_ownership.dis_size();
  size_type first_S_ios_dis_ige = true_S_ios_ge_ownership.first_index();
  std::vector<size_type> tmp_S_massive_partition (S_dis_nge, 0); // TODO: massive memory area
  size_type S_ios_ige = 0;
  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
                 S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
    distributor ios_gev_ownership  = ios_geo_element [S_variant].ownership();
    size_type   first_ios_dis_igev = ios_gev_ownership.first_index();
    for (size_type ios_igev = 0, ios_ngev = ios_gev_ownership.size(); ios_igev < ios_ngev; ios_igev++, S_ios_ige++) {
      const geo_element& S = ios_geo_element [S_variant] [ios_igev];
      index_set K_ios_ige_set = K_ball.dis_at (S[0]);
      for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
        K_ios_ige_set.inplace_intersection (K_ball.dis_at(S[iloc]));
      }
      check_macro (K_ios_ige_set.size() >  0, "connectivity: S={"<<S<<"} not found in the side set");
      size_type S_owner = 0;
      for (index_set::const_iterator iter = K_ios_ige_set.begin(), last = K_ios_ige_set.end(); iter != last; iter++) {
        size_type K_ios_ige = *iter;
        S_owner = std::max(S_owner, massive_partition_by_dimension[K_dim][K_ios_ige]);
      }
      size_type S_ios_dis_ige = first_S_ios_dis_ige + S_ios_ige;
      tmp_S_massive_partition [S_ios_dis_ige] = S_owner;
    }
  }
  // ------------------------------------------------------------------------
  // 4) all_reduce:  TODO: massive memory area comms
  // ------------------------------------------------------------------------
  massive_partition_by_dimension[S_dim].resize (tmp_S_massive_partition.size(), 0);
  communicator comm = ios_vertex_ownership.comm();
  mpi::all_reduce (
        comm,
        tmp_S_massive_partition.begin().operator->(),
        tmp_S_massive_partition.size(),
        massive_partition_by_dimension[S_dim].begin().operator->(),
        mpi::maximum<size_type>());
  // ------------------------------------------------------------------------
  // 5) copy the massive_partition_by_dimension into a distributed disarray "partition", variant by variant
  // ------------------------------------------------------------------------
  size_type S_ios_dis_ige = first_S_ios_dis_ige;
  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
                 S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
    partition_by_variant [S_variant].resize (ios_geo_element [S_variant].ownership());
    for (size_type S_ios_igev = 0, S_ios_negv = partition_by_variant [S_variant].size(); 
                   S_ios_igev < S_ios_negv; S_ios_igev++, S_ios_dis_ige++) {
      partition_by_variant [S_variant][S_ios_igev] = massive_partition_by_dimension[S_dim] [S_ios_dis_ige];
    }
  }
}
// --------------------------------------------------------------------------
// edges & faces renumbering subroutine 
// --------------------------------------------------------------------------
void
geo_element_renumbering_part2 (
 // input:
    const std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                                               ios_geo_element,
    const geo_size&                            ios_gs,
          size_t                               dis_nv,
          size_t                               side_dim,
 // output:
          std::array<hack_array<geo_element_hack>,reference_element::max_variant>&
                                               geo_element,
          geo_size&                            gs,
          std::array<disarray<size_t>,reference_element::max_variant>&
                                               igev2ios_dis_igev,
          std::array<disarray<size_t>,reference_element::max_variant>&
                                               ios_igev2dis_igev,
          std::array<disarray<size_t>,4>&      ios_ige2dis_ige,
 // tmp:
          std::array<disarray<size_t>, reference_element::max_variant>& 
                                               partition_by_variant)
{
  typedef geo_element::size_type size_type;
  typedef hack_array<geo_element_hack>::iterator             iterator_by_variant;
  typedef hack_array<geo_element_hack>::const_iterator const_iterator_by_variant;
 
  communicator comm = ios_geo_element[0].ownership().comm();
  //
  // 3) build geo_elemen [variant]:
  //
  // MERGE: could be done from here:
  size_type nge = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
 
    ios_geo_element [variant].repartition (
          partition_by_variant [variant],
          geo_element [variant],
          igev2ios_dis_igev [variant],
          ios_igev2dis_igev [variant]);
 
    gs.ownership_by_variant [variant] = geo_element [variant].ownership();
    nge += geo_element [variant].size();
  }
  gs.ownership_by_dimension [side_dim] = distributor (distributor::decide, comm, nge);
  //
  // 4) build first_by_variant [variant]
  //
  {
    size_type first_v = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                   variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
  
        gs.first_by_variant[variant] = distributor (distributor::decide, comm, first_v);
        first_v += geo_element [variant].size();
    } 
  } 
  //
  // 5) set element number
  //   & build ios_ige2dis_ige[map_dim] from ios_igev2dis_igev[variant]
  // TODO: a way to build directly ios_ige2dis_ige without assembly (and many comms)
  //
  ios_ige2dis_ige [side_dim].resize (ios_gs.ownership_by_dimension [side_dim]); // OK
  {
    size_type first_dis_ige = gs.ownership_by_dimension [side_dim].first_index();
    size_type ige = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                   variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
      size_type first_v = gs.first_by_variant [variant].size();
      for (size_type igev = 0, ngev = geo_element [variant].size(); igev < ngev; igev++, ige++) {
        ::rheolef::geo_element& S = geo_element [variant] [igev]; // GNU C++ 4.5 BUG : add ::rheolef::
        size_type dis_ige  = first_dis_ige + ige;
        S.set_dis_ie (dis_ige);
        size_type ios_dis_ige = S.ios_dis_ie();
        ios_ige2dis_ige [side_dim].dis_entry (ios_dis_ige) = dis_ige;
      }
    }
    ios_ige2dis_ige [side_dim].dis_entry_assembly();
  }
}  
void
geo_element_renumbering_propagate (
// input:
    const std::vector<geo_element::size_type>& new_global_node_num,
          size_t                               dis_nnod,
// modified:
          hack_array<geo_element_hack>&     gev)
{
    using namespace std;
    typedef geo_element::size_type size_type;
    //
    // vertices S[iloc] of new numbered element table still have ios numbering: fix it
    //
    // TODO: iterator loop on ge instead of gev
    size_type first_dis_igev = gev.ownership().first_index();
    for (size_type igev = 0, ngev = gev.size(); igev < ngev; igev++) {
      geo_element& S = gev [igev];
      for (size_type iloc = 0, nloc = S.n_node(); iloc < nloc; iloc++) {
        assert_macro (S[iloc] < dis_nnod, "node index "<<S[iloc]<<" out of range [0:"<<dis_nnod<<"[");
        assert_macro (new_global_node_num[S[iloc]] < dis_nnod, "new node index "<<new_global_node_num[S[iloc]] <<" out of range [0:"<<dis_nnod<<"[");
        S[iloc] = new_global_node_num [S[iloc]];
      }
    }
}

template <class T>
void
geo_rep<T,distributed>::build_external_entities ()
{
  distributor vertex_ownership = base::_geo_element[reference_element::p].ownership();
  size_type first_dis_iv = vertex_ownership.first_index();
  size_type  last_dis_iv = vertex_ownership.last_index();
  size_type       dis_nv = vertex_ownership.dis_size();
 
  distributor node_ownership = base::_node.ownership();
  size_type first_dis_inod = node_ownership.first_index();
  size_type  last_dis_inod = node_ownership.last_index();
  size_type       dis_nnod = node_ownership.dis_size();
 
  // 1a) list external vertex & nodes indexes from internal elements
  std::set<size_type> ext_vertex_set;
  std::set<size_type> ext_node_set;
  std::vector<size_type> dis_inod1;
  for (size_type dim = 1; dim <= base::_gs._map_dimension; dim++) {
    // loop on elements
    for (size_type ige = 0, nge = base::_gs.ownership_by_dimension[dim].size(); ige < nge; ige++) {
      const geo_element& K = get_geo_element(dim,ige);
        space_numbering::dis_idof (base::_piola_basis, base::_gs, K, dis_inod1);
        for (size_type loc_inod = 0, loc_nnod = dis_inod1.size(); loc_inod < loc_nnod; loc_inod++) {
          size_type dis_inod = dis_inod1 [loc_inod];
          assert_macro (dis_inod < dis_nnod, "node index "<< dis_inod <<" out of range [0:"<<dis_nnod<<"[");
          if (node_ownership.is_owned (dis_inod)) continue;  
          ext_node_set.insert (dis_inod);
          if (loc_inod >= K.size()) continue;
          // then reports vertex to node: same owner
          size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
          check_macro (dis_iv < dis_nv, "vertex index "<< dis_iv <<" out of range [0:"<<dis_nv<<"[");
          check_macro (!vertex_ownership.is_owned (dis_iv), "strange bug: not owned (nod) but is_owned(ver)");
          ext_vertex_set.insert (dis_iv);
        }
    }
  }
  // 1b) pull external vertices & nodes:
  base::_node.append_dis_indexes (ext_node_set);
  base::_geo_element [reference_element::p].append_dis_indexes (ext_vertex_set);
 
  // 2) K.edge(i) and K.face(i) have been completed by set set_element_side_index 
  //     but not propagated to external elements : update them by refreshing the external data
  for (size_type variant = reference_element::first_variant_by_dimension (1);
                 variant < reference_element:: last_variant_by_dimension (base::_gs._map_dimension); ++variant) {
    const std::map<size_type,geo_element_auto<>>& ext_gev = base::_geo_element [variant].get_dis_map_entries();
    index_set ext_dis_ie_set;
    for (auto i: ext_gev) { ext_dis_ie_set.insert (i.first); }
    base::_geo_element [variant].set_dis_indexes (ext_dis_ie_set);
  }
  // 3a) list external vertex indexes from external elements
  std::array<index_set,reference_element::max_variant> ext_dis_igev_set;
  for (size_type dim = base::_gs._map_dimension; dim >= 1; dim--) {
    for (size_type variant = reference_element::first_variant_by_dimension (dim);
                   variant < reference_element:: last_variant_by_dimension (dim); ++variant) {
      for (auto x : base::_geo_element[variant].get_dis_map_entries()) {
        const geo_element& K = x.second;
        for (size_type subgeo_dim = 0; subgeo_dim < K.dimension(); ++subgeo_dim) {
          for (size_type loc_is = 0, loc_ns = K.n_subgeo(subgeo_dim); loc_is < loc_ns; ++loc_is) {
            size_type dis_ige = (subgeo_dim == 0) ? base::dis_inod2dis_iv(K[loc_is]) : (subgeo_dim == 1) ? K.edge(loc_is) : K.face(loc_is);
            check_macro(dis_ige != std::numeric_limits<size_type>::max(), "invalid external subgeo dis_index");
            if (base::sizes().ownership_by_dimension [subgeo_dim].is_owned (dis_ige)) continue;
            // this subgeo is external: add it to ext_dis_igev_set[variant] :
            size_type dis_size = base::sizes().ownership_by_dimension [subgeo_dim].dis_size();
            check_macro(dis_ige < dis_size, "invalid external "<<subgeo_dim<<"d subgeo dis_index = "<<dis_ige<<" out of range [0:"<<dis_size<<"[");
            // convert dis_ige to dis_igev and get its variant
            size_type variant;
            size_type dis_igev = base::sizes().dis_ige2dis_igev_by_dimension (subgeo_dim, dis_ige, variant);
            ext_dis_igev_set [variant].insert (dis_igev);
          }
        }
      }
    }
  }
  // 3b) append external subgeo entities from all external elements
  for (size_type variant = reference_element::first_variant_by_dimension (1);
                 variant < reference_element:: last_variant_by_dimension (base::_gs._map_dimension); ++variant) {
    base::_geo_element [variant].append_dis_indexes (ext_dis_igev_set [variant]);
  }
  // 4a) list external nodes from all external entities, based on dis_idof(), for high order meshes 
  for (size_type dim = base::_gs._map_dimension; dim >= 1; dim--) {
    for (size_type variant = reference_element::first_variant_by_dimension (dim);
                   variant < reference_element:: last_variant_by_dimension (dim); ++variant) {
      for (auto x : base::_geo_element[variant].get_dis_map_entries()) {
        const geo_element& K = x.second;
        space_numbering::dis_idof (base::_piola_basis, base::_gs, K, dis_inod1);
        for (size_type loc_inod = 0, loc_nnod = dis_inod1.size(); loc_inod < loc_nnod; loc_inod++) {
          size_type dis_inod = dis_inod1 [loc_inod];
          assert_macro (dis_inod < dis_nnod, "node index "<< dis_inod <<" out of range [0:"<<dis_nnod<<"[");
          if (node_ownership.is_owned (dis_inod)) continue;  
          ext_node_set.insert (dis_inod);
        }
      }
    }
  }
  // 4b) pull external nodes:
  base::_node.set_dis_indexes (ext_node_set);
}

template <class T>
void
geo_rep<T,distributed>::set_element_side_index (size_type side_dim)
{
  distributor vertex_ownership = base::_gs.ownership_by_dimension [0];
  distributor   node_ownership = base::_node.ownership();
  size_type           nv   = vertex_ownership.size();
  size_type first_dis_iv   = vertex_ownership.first_index();
  size_type  last_dis_iv   = vertex_ownership. last_index();
  size_type     dis_nnod   = node_ownership.dis_size(); // for bound checks
  // ------------------------------------------------------------------------
  // 0) build external vertices index set
  // ------------------------------------------------------------------------
  index_set ext_iv_set; // size=O((N/nproc)^((d-1)/d))
  for (size_type dim = side_dim; dim <= base::_gs._map_dimension; dim++) {
    for (size_type var = reference_element::first_variant_by_dimension(dim);
                   var < reference_element:: last_variant_by_dimension(dim); var++) {
      for (size_type igev = 0, ngev = base::_geo_element[var].size(); igev < ngev; igev++) {
        const geo_element& K = base::_geo_element [var] [igev];
        for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
          size_type dis_inod = K[iloc];
          if (! node_ownership.is_owned(dis_inod)) { // vertex ownership follows node ownership
            size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
            ext_iv_set.insert (dis_iv);
          }
        }
      }
    }
  }
  // ------------------------------------------------------------------------
  // 1) ball(X) := { E; X is a vertex of E }
  // ------------------------------------------------------------------------
  index_set empty_set; // TODO: add a global allocator to empty_set
  disarray<index_set,distributed> ball [reference_element::max_variant];
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
    ball [variant].resize (vertex_ownership, empty_set);
    distributor gev_ownership  = base::_gs.ownership_by_variant [variant];
    size_type   first_dis_igev = gev_ownership.first_index();
    for (size_type igev = 0, ngev = base::_geo_element [variant].size(); igev < ngev; igev++) {
      const geo_element& S = base::_geo_element [variant] [igev];
      size_type dis_igev = first_dis_igev + igev;
      for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
        size_type dis_inod = S[iloc];
        assert_macro (dis_inod < dis_nnod, "S={"<<S<<"}: "<<iloc<<"-node index " << dis_inod << " out of range [0:"<<dis_nnod<<"[");
        size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
        if (vertex_ownership.is_owned(dis_iv)) {
          size_type iv = dis_iv - first_dis_iv;
          ball [variant][iv] += dis_igev;
        } else {
          index_set dis_igev_set;
          dis_igev_set += dis_igev;
          ball [variant].dis_entry (dis_iv) += dis_igev_set; // not so efficient: union with {dis_iegv}
        }
      }
    }
    ball [variant].dis_entry_assembly();
    // ------------------------------------------------------------------------
    // for all the dis_iv that are not handled by the current process
    // but are referenced by at least a side, get ext_ball[dis_iv]
    // ------------------------------------------------------------------------
    // Question: faudra-t'il inclure d'autres sommets dans ext_iv_set ?
    // -> ceux issus des triangles, tetra, ect, et externes a la partition ?
    // Question : est-ce que ca en ajouterait ? reponse : OUI (teste')
    // mais c'est pas utile pour l'instant : ball sert juste aux aretes
    ball [variant].set_dis_indexes (ext_iv_set);
  }
  // ------------------------------------------------------------------------
  // 2)  on parcourt chaque dis_E de ball(iv) et ext_ball(dis_iv) :
  //     si dis_E n'est pas local, on le met dans une liste ext_isid_set
  //     et on importe l'arete en tant que geo_element
  // ------------------------------------------------------------------------
  // scan ball(iv) for external sides
  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                 variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
 
    // 2.1) scan ball(iv) for external sides
    std::set<size_type> ext_igev_set; // size=O((N/nproc)^((d-1)/d))
    distributor gev_ownership  = base::_gs.ownership_by_variant [variant];
    size_type dis_ngev = gev_ownership.dis_size();
    for (size_type iv = 0, nv = vertex_ownership.size(); iv < nv; iv++) {
 
      const index_set& ball_iv = ball [variant] [iv];
      for (index_set::const_iterator iter = ball_iv.begin(), last = ball_iv.end(); iter != last; iter++) {
        size_type dis_igev = *iter;
        assert_macro (dis_igev < dis_ngev, "index "<<dis_igev<<" of element variant="<<variant<<" is out of range [0:"<<dis_ngev<<"[");
        if (! gev_ownership.is_owned (dis_igev)) {
          ext_igev_set.insert (dis_igev);
        }
      }
    }
    // 2.2) scan ball(dis_iv) for external sides
    typedef disarray<index_set>::scatter_map_type map_type;
    const map_type& ball_dis  = ball [variant].get_dis_map_entries();
    for (typename map_type::const_iterator iter_b = ball_dis.begin(), last_b = ball_dis.end(); iter_b != last_b; iter_b++) {
      size_type             dis_iv = (*iter_b).first;
      const index_set& ball_dis_iv = (*iter_b).second;
      for (index_set::const_iterator iter = ball_dis_iv.begin(), last = ball_dis_iv.end(); iter != last; iter++) {
        size_type dis_igev = *iter;
        assert_macro (dis_igev < dis_ngev, "index "<<dis_igev<<" of element variant="<<variant<<" is out of range [0:"<<dis_ngev<<"[");
        if (! gev_ownership.is_owned (dis_igev)) {
          ext_igev_set.insert (dis_igev);
        }
      }
    }
    // 2.3) import external sides
    base::_geo_element[variant].append_dis_indexes (ext_igev_set);
  }
  // ------------------------------------------------------------------------
  // 3) pour K dans partition(iproc)
  //      pour (dis_A,dis_B) arete de K
  //        set = dis_ball(dis_A) inter dis_ball(dis_B) = {dis_iedg}
  //        E = dis_edges(dis_iedg)
  //        => on numerote dis_iedg cette arete dans le geo_element K
  //        et on indique son orient en comparant a E, arete qui definit l'orient
  // ------------------------------------------------------------------------
  for (size_type dim = side_dim+1; dim <= base::_gs._map_dimension; dim++) {
    for (size_type var = reference_element::first_variant_by_dimension(dim);
                   var < reference_element:: last_variant_by_dimension(dim); var++) {
      for (size_type igev = 0, ngev = base::_geo_element [var].size(); igev < ngev; igev++) {
        geo_element& K = base::_geo_element [var] [igev];
        size_type S_iv  [4];
        size_type S_inod[4];
        for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
          size_type S_size = K.subgeo_size (side_dim, loc_isid);
          size_type S_variant = reference_element::variant (S_size, side_dim);
          size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
          S_inod [0] = K[loc_jv0];
          S_iv   [0] = base::dis_inod2dis_iv (S_inod[0]);
          const disarray<index_set,distributed>& ball_S_variant = ball [S_variant];
          index_set igev_set = ball_S_variant.dis_at (S_iv[0]);
          for (size_type sid_jloc = 1, sid_nloc = S_size; sid_jloc < sid_nloc; sid_jloc++) { 
            size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
            S_inod [sid_jloc] = K[loc_jv];
            S_iv   [sid_jloc] = base::dis_inod2dis_iv (S_inod[sid_jloc]);
            const index_set& ball_jv = ball_S_variant.dis_at (S_iv[sid_jloc]);
            igev_set.inplace_intersection (ball_jv);
          }
          check_macro (igev_set.size() >  0, "connectivity: "<<S_size<<"-side ("
                <<S_iv[0]<<","<<S_iv[1]<<","<<S_iv[2]<<","<<S_iv[3]<<") not found in the side set");
          check_macro (igev_set.size() == 1, "connectivity: the same side is multiply represented");
          size_type dis_igev = *(igev_set.begin());
          const geo_element& S = base::_geo_element[S_variant].dis_at(dis_igev);
          size_type dis_isid = S.dis_ie();
          if (side_dim == 1) {
            // side: edge
            geo_element::orientation_type orient = S.get_edge_orientation (S_inod[0], S_inod[1]);
            K.edge_indirect (loc_isid).set (orient, dis_isid);
          } else { // side_dim == 2 
            geo_element::orientation_type orient;
            geo_element::shift_type       shift;
            if (K.subgeo_size (side_dim, loc_isid) == 3) {
              // side: triangle
              S.get_orientation_and_shift (S_inod[0], S_inod[1], S_inod[2], orient, shift);
            } else {
              // side: quadrangle
              S.get_orientation_and_shift (S_inod[0], S_inod[1], S_inod[2], S_inod[3], orient, shift);
            }
            K.face_indirect (loc_isid).set (orient, dis_isid, shift);
          }
        }
      }
    }
  }
}
// --------------------------------------------------------------------------
// set permutation for nodes: ios_inod2dis_inod & inod2ios_dis_inod
// and redistribute ios_node[] into _node[]
// --------------------------------------------------------------------------
template <class T>
void
geo_rep<T,distributed>::set_ios_permutation (disarray<size_type,distributed>& idof2ios_dis_idof) const
{
  space_numbering::generic_set_ios_permutation (base::get_piola_basis(), map_dimension(), base::sizes(), _igev2ios_dis_igev, idof2ios_dis_idof);
}
template <class T>
void
geo_rep<T,distributed>::node_renumbering (const distributor& ios_node_ownership)
{
  std::array<size_type,reference_element::max_variant> loc_ndof_by_variant;
  reference_element::init_local_nnode_by_variant (base::order(), loc_ndof_by_variant);
  set_ios_permutation (_inod2ios_dis_inod);
}
// --------------------------------------------------------------------------
// get geo
// --------------------------------------------------------------------------
template <class T>
idiststream&
geo_rep<T,distributed>::get (idiststream& ips)
{
  using namespace std;
  check_macro (ips.good(), "bad input stream for geo.");
  communicator comm = base::_geo_element[reference_element::p].ownership().comm();
 
  size_type io_proc = idiststream::io_proc();
  size_type my_proc = comm.rank();
  // ------------------------------------------------------------------------
  // 1) read file
  // ------------------------------------------------------------------------
  //
  // 1.1) get header
  //
  check_macro (dis_scatch(ips, ips.comm(), "\nmesh"), "input stream does not contains a geo.");
  ips >> base::_version;
  check_macro (base::_version == 4, "geo version < 4 not supported (HINT: see geo -upgrade)");
 
  geo_header hdr;
  ips >> hdr;
  check_macro (! hdr.need_upgrade(), 
        "unsupported geo without connectivity in the distributed version: HINT: use geo_upgrade");
  base::_have_connectivity = true;
  base::_dimension         = hdr.dimension;
  base::_gs._map_dimension = hdr.map_dimension;
  base::_sys_coord         = hdr.sys_coord;
  base::_name              = "unnamed";
  base::_piola_basis.reset_family_index (hdr.order);
  size_type dis_nnod = hdr.dis_size_by_dimension [0];
  size_type dis_nedg = hdr.dis_size_by_dimension [1];
  size_type dis_nfac = hdr.dis_size_by_dimension [2];
  size_type dis_ne   = hdr.dis_size_by_dimension [base::_gs._map_dimension];
  //
  // 1.2) get node coordinates
  //
  size_type ios_size_by_variant [reference_element::max_variant];
  std::fill (ios_size_by_variant, ios_size_by_variant+reference_element::max_variant, 0);
  disarray<node_type> ios_node (dis_nnod);
  ios_node.get_values (ips, _point_get<T>(base::_dimension));
  check_macro (ips.good(), "bad input stream for geo.");
  //
  // 1.3) get elements
  //
  std::array<hack_array<geo_element_hack>, reference_element::max_variant> ios_geo_element;
  size_type ios_nv = 0;
  if (base::_gs._map_dimension == 0) {
    ios_nv = ios_node.size();
  } else {
    // set elt ios index and compute map_dim & ios_nv (ios_nv < ios_nnod when order > 1)
    size_type ios_ne = 0;
    size_type dis_ne = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
                   variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
      geo_element::parameter_type param (variant, 1);
      ios_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
      ios_geo_element [variant].get_values (ips);
      _ios_gs.first_by_variant   [variant] = distributor (distributor::decide, base::comm(), ios_ne);
      ios_ne  += ios_geo_element [variant].size();
      dis_ne  += ios_geo_element [variant].dis_size();
    }
    build_apparent_ios_ownership(
                ios_geo_element,
                base::_gs._map_dimension,
                _ios_gs.ownership_by_dimension [base::_gs._map_dimension],
                _ios_gs.ownership_by_variant);
 
    vector<size_type> local_is_vertex (dis_nnod, 0);
    for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
                   variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
      size_type first_ios_dis_v    = _ios_gs.first_by_variant [variant].dis_size();
      size_type first_ios_dis_igev = ios_geo_element [variant].ownership().first_index();
      size_type ios_igev = 0;
      for (iterator_by_variant iter = ios_geo_element [variant].begin(), last = ios_geo_element [variant].end();
                               iter != last; iter++, ios_igev++) {
        geo_element& K = *iter;
        size_type ios_dis_ie = first_ios_dis_v + first_ios_dis_igev + ios_igev;
        K.set_ios_dis_ie (ios_dis_ie); // OK
        ios_size_by_variant [K.variant()]++;
        if (base::order() > 1) {
          for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
            local_is_vertex [K[iloc]] = 1;
          }
        }
      }
    }
    if (base::order() == 1) {
      ios_nv = ios_node.size();
    } else {
      vector<size_type> global_is_vertex (dis_nnod, 0);
      mpi::all_reduce (
        comm,
        local_is_vertex.begin().operator->(),
        local_is_vertex.size(),
        global_is_vertex.begin().operator->(),
        mpi::maximum<size_type>());
      ios_nv = accumulate (global_is_vertex.begin() + ios_node.ownership().first_index(),
                           global_is_vertex.begin() + ios_node.ownership().last_index(),  0);
    }
  }
  distributor ios_vertex_ownership;
  if (base::order() == 1) {
    ios_vertex_ownership = ios_node.ownership();
  } else {
    ios_vertex_ownership = distributor (distributor::decide, base::comm(), ios_nv); // not sure...
  }
  size_type dis_nv = ios_vertex_ownership.dis_size();
  //
  // 1.4) create 0d vertex-elements
  //
  // set ios_dis_iv index as fisrt field of the idx_vertex pair:
  // # of node that are vertices:
  {
    geo_element::parameter_type param (reference_element::p, 1);
    ios_geo_element [reference_element::p].resize (ios_vertex_ownership, param);
    _ios_gs.ownership_by_variant   [reference_element::p] = ios_vertex_ownership;
    _ios_gs.ownership_by_dimension [0]                    = ios_vertex_ownership;
    size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
    for (size_type ios_iv = 0, ios_nv = ios_vertex_ownership.size(); ios_iv < ios_nv; ios_iv++) {
      geo_element& P = ios_geo_element [reference_element::p] [ios_iv];
      size_type ios_dis_iv = first_ios_dis_iv + ios_iv;
      P [0] = ios_dis_iv;
      P.set_ios_dis_ie (ios_dis_iv);
      ios_size_by_variant [P.variant()]++;
    }
  }
  //
  // 1.5) get faces & edges
  //
  if (base::_gs._map_dimension > 0) {
    for (size_type side_dim = base::_gs._map_dimension - 1; side_dim >= 1; side_dim--) {
      size_type ios_ngev = 0;
      size_type dis_ngev = 0;
      for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                     variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
        geo_element::parameter_type param (variant, 1);
        ios_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
        ios_geo_element [variant].get_values (ips); 
        _ios_gs.first_by_variant     [variant] = distributor (distributor::decide, base::comm(), ios_ngev);
        ios_ngev += ios_geo_element [variant].size();
        dis_ngev += ios_geo_element [variant].dis_size();
      }
      build_apparent_ios_ownership(
                ios_geo_element,
                side_dim,
                _ios_gs.ownership_by_dimension [side_dim],
                _ios_gs.ownership_by_variant);
 
      for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
                     variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
        size_type ios_igev = 0;
        size_type first_dis_v        = _ios_gs.first_by_variant [variant].dis_size();
        size_type first_ios_dis_igev = ios_geo_element [variant].ownership().first_index();
        for (iterator_by_variant iter = ios_geo_element [variant].begin(), last = ios_geo_element [variant].end();
                                 iter != last; iter++, ios_igev++) {
          geo_element& S = *iter;
          size_type ios_dis_igev = first_dis_v + first_ios_dis_igev + ios_igev;
          S.set_ios_dis_ie (ios_dis_igev); // OK
          ios_size_by_variant [S.variant()]++;
        }
      }
    }
  }
  // ------------------------------------------------------------------------
  // 2) mesh partition & element renumbering
  // ------------------------------------------------------------------------
  distributor true_ios_ownership_by_dimension [4];
  true_ios_ownership_by_dimension [base::_gs._map_dimension] 
   = build_true_ios_ge_ownership_by_dimension (ios_geo_element, base::_gs._map_dimension);
 
  disarray<size_type> partition = geo_mpi_partition (
    ios_geo_element,
    true_ios_ownership_by_dimension [base::_gs._map_dimension],
    base::_gs._map_dimension,
    dis_nv);
 
  // copy partition into partition_by_variant[]:
  std::array<disarray<size_t>, reference_element::max_variant> partition_by_variant;
  {
    typename disarray<size_type>::const_iterator iter = partition.begin();
    for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
                   variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
      partition_by_variant [variant].resize (ios_geo_element [variant].ownership());
      for (typename disarray<size_type>::iterator iter_by_var = partition_by_variant [variant].begin(),
                                               last_by_var = partition_by_variant [variant].end(); 
                                               iter_by_var != last_by_var; iter_by_var++, iter++) {
         *iter_by_var = *iter;
      }
    }
  }
  geo_element_renumbering_part2 (
        ios_geo_element,
        _ios_gs,
        dis_nv,
        base::_gs._map_dimension,
        base::_geo_element,
        base::_gs,
        _igev2ios_dis_igev,
        _ios_igev2dis_igev,
        _ios_ige2dis_ige,
        partition_by_variant);
 
 
  // copy partition into massive_partition_by_dimension[map_dim]:
  std::array<std::vector<size_t>, 4>  massive_partition_by_dimension;
  {
    std::vector<size_t>  tmp_massive_partition (dis_ne, 0);
    size_type true_first_dis_ige = true_ios_ownership_by_dimension [base::_gs._map_dimension].first_index();
    for (size_type ios_ie = 0, ios_ne = partition.size(); ios_ie < ios_ne; ios_ie++) {
      size_type ios_dis_ie = true_first_dis_ige + ios_ie;
      tmp_massive_partition [ios_dis_ie] = partition [ios_ie];
    } 
    massive_partition_by_dimension [base::_gs._map_dimension].resize (dis_ne);
    mpi::all_reduce (
        comm,
        tmp_massive_partition.begin().operator->(),
        tmp_massive_partition.size(),
        massive_partition_by_dimension[base::_gs._map_dimension].begin().operator->(),
        mpi::maximum<size_type>());
  }
  // propagate partition to subgeo : faces, edges & vertices
  for (size_type supergeo_dim = base::_gs.map_dimension(); supergeo_dim > 0; supergeo_dim--) {
    geo_element_renumbering_part1_new (
      ios_geo_element, 
      _ios_gs,
      supergeo_dim-1,
      massive_partition_by_dimension,
      partition_by_variant);
  }
  // ------------------------------------------------------------------------
  // 3) vertices renumbering 
  // ------------------------------------------------------------------------
  //
  // 3.1) redistribute the _geo_element[p] vertices
  //
  partition_by_variant [reference_element::p].resize (ios_vertex_ownership);
  size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
  for (size_type ios_iv = 0, ios_nv = ios_vertex_ownership.size(); ios_iv < ios_nv; ios_iv++) {
    size_type ios_dis_iv = first_ios_dis_iv + ios_iv;
    partition_by_variant [reference_element::p] [ios_iv] = massive_partition_by_dimension[0] [ios_dis_iv];
  }
  disarray<size_type> iv2ios_dis_iv;
  ios_geo_element [reference_element::p].repartition (
        partition_by_variant [reference_element::p],
        base::_geo_element [reference_element::p],
        iv2ios_dis_iv,
        _ios_ige2dis_ige[0]);
 
  distributor vertex_ownership = base::_geo_element[reference_element::p].ownership();
  base::_gs.ownership_by_dimension [0]                       = vertex_ownership;
  base::_gs.ownership_by_variant   [reference_element::p] = vertex_ownership;
  base::_gs.first_by_variant       [reference_element::p] = distributor (0, base::comm(), 0);
  //
  // 3.2) set the element[0] disarray 
  //
  size_type first_dis_iv = vertex_ownership.first_index();
  size_type  last_dis_iv = vertex_ownership.last_index();
  _igev2ios_dis_igev [reference_element::p].resize (vertex_ownership);
  for (size_type iv = 0, nv = base::_geo_element[reference_element::p].size(); iv < nv; iv++) {
    geo_element& P = base::_geo_element[reference_element::p] [iv];
    size_type dis_iv = first_dis_iv + iv;
    P[0]        = dis_iv;
    P.set_dis_ie (dis_iv);
    _igev2ios_dis_igev [reference_element::p] [iv] = P.ios_dis_ie();
  }
  // ------------------------------------------------------------------------
  // 4) edge & face renumbering 
  // ------------------------------------------------------------------------
  if (base::_gs._map_dimension > 0) {
    for (size_type side_dim = base::_gs._map_dimension-1; side_dim > 0; side_dim--) {
 
      geo_element_renumbering_part2 (
        ios_geo_element,
        _ios_gs,
        dis_nv,
        side_dim,
        base::_geo_element,
        base::_gs,
        _igev2ios_dis_igev,
        _ios_igev2dis_igev,
        _ios_ige2dis_ige,
        partition_by_variant);
    }
  }
  // ------------------------------------------------------------------------
  // 5) get domain, until end-of-file (requires ios_ige2dis_ige renumbering)
  // ------------------------------------------------------------------------
  do {
    domain_indirect_basic<distributed> dom;
    bool status = dom.get (ips, *this);
    if (!status) break;
    base::_domains.push_back (dom);
  } while (true);
  // ------------------------------------------------------------------------
  // 6) renumbering ios_nodes[]: set permutations inod2ios_dis_inod[]
  // ------------------------------------------------------------------------
  node_renumbering (ios_node.ownership());
  distributor node_ownership = _inod2ios_dis_inod.ownership();
  base::_gs.node_ownership = node_ownership;
  // ------------------------------------------------------------------------
  // 7) redistribute the nodes : from ios_node[] to node() disarrays
  // ------------------------------------------------------------------------
  _ios_inod2dis_inod.resize (ios_node.ownership(), std::numeric_limits<size_type>::max());
  _inod2ios_dis_inod.reverse_permutation (_ios_inod2dis_inod);
 
  // ------------------------------------------------------------------------
  // 8) propagate new node numbering
  // ------------------------------------------------------------------------
  // 8.1) build new node disarray:  base::_node
  base::_node.resize (node_ownership);
  ios_node.permutation_apply (_ios_inod2dis_inod, base::_node);
  // 
  // 8.2) global table of node renumbering
  //
  vector<size_type> new_local_node_num (dis_nnod, 0);
  size_type first_ios_inod = _ios_inod2dis_inod.ownership().first_index();
  size_type  last_ios_inod = _ios_inod2dis_inod.ownership(). last_index();
  for (size_type dis_ios_inod = first_ios_inod; dis_ios_inod < last_ios_inod; dis_ios_inod++) {
     size_type ios_inod = dis_ios_inod - first_ios_inod;
     new_local_node_num [dis_ios_inod] = _ios_inod2dis_inod[ios_inod];
  }
  vector<size_type> new_global_node_num (dis_nnod, 0);
  mpi::all_reduce (
        comm,
        new_local_node_num.begin().operator->(),
        new_local_node_num.size(),
        new_global_node_num.begin().operator->(),
        mpi::maximum<size_type>());
  //
  // 8.3) propagate the new node numbering into geo_elements:
  //      vertices K[iloc] of new numbered element table K still have ios numbering: fix it
  //
  for (size_type dim = base::_gs._map_dimension; dim > 0; dim--) {
    for (size_type variant = reference_element::first_variant_by_dimension(dim);
                   variant < reference_element:: last_variant_by_dimension(dim); variant++) {
      geo_element_renumbering_propagate (
        new_global_node_num,
        dis_nnod,
        base::_geo_element [variant]);
    }
  }
  // ------------------------------------------------------------------------
  // 9) set indexes on faces and edges of elements, for P2 approx
  // ------------------------------------------------------------------------
  for (size_type dim = base::_gs._map_dimension - 1; base::_gs._map_dimension > 0 && dim > 0; dim--) {
    set_element_side_index (dim);
  }
  // ------------------------------------------------------------------------
  // 10) set external entities, at partition boundaries
  // ------------------------------------------------------------------------
  build_external_entities ();
  // ------------------------------------------------------------------------
  // 11) bounding box: _xmin, _xmax
  // ------------------------------------------------------------------------ 
  base::compute_bbox();
  return ips;
}
// ----------------------------------------------------------------------------
// read from file
// ----------------------------------------------------------------------------
template <class T>
void
geo_rep<T,distributed>::load (
  std::string filename, 
  const communicator& comm)
{
  idiststream ips;
  ips.open (filename, "geo", comm);
  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
  get (ips);
  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
  std::string name = get_basename (root_name);
  base::_name = name;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_rep<Float,distributed>;
 
} // namespace rheolef
#endif // _RHEOLEF_HAVE_MPI