///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is sequential in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

#include "rheolef/config.h"

#ifdef _RHEOLEF_HAVE_MPI
#include "rheolef/geo_domain.h"

namespace rheolef {

// ------------------------------------------------------------------------
// build edge and face connectivity from domain
// ------------------------------------------------------------------------
template <class T>
void
geo_rep<T,distributed>::domain_set_side_part1 (
    const domain_indirect_rep<distributed>&    indirect, 
    const geo_abstract_rep<T,distributed>&     bgd_omega,
    size_type                                  sid_dim,
    disarray<size_type>&                       bgd_isid2dom_dis_isid,
    disarray<size_type>&                       dom_isid2bgd_isid,
    disarray<size_type>&                       dom_isid2dom_ios_dis_isid,
    size_type                                  size_by_variant [reference_element::max_variant])
{
  if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
  // ------------------------------------------------------------------------
  // 1) side compact re-numbering
  // ------------------------------------------------------------------------
  // 1.1 loop on elements and mark used sides
  disarray<size_type> bgd_isid_is_on_domain     (bgd_omega.geo_element_ownership(sid_dim), 0); // logical, init to "false"
  disarray<size_type> bgd_ios_isid_is_on_domain (bgd_omega.geo_element_ios_ownership(sid_dim), 0);
  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
    size_type ige = indirect.oige (ioige).index();
    const geo_element& bgd_K = bgd_omega.get_geo_element (base::_gs._map_dimension, ige);
    for (size_type loc_isid = 0, loc_nsid = bgd_K.n_subgeo(sid_dim); loc_isid < loc_nsid; loc_isid++) {
      size_type bgd_dis_isid = 0; // TODO: = bgd_K.subgeo (sid_dim, loc_isid);
      switch (sid_dim) {
        case 0: {
	        size_type bgd_dis_inod = bgd_K[loc_isid];
	        bgd_dis_isid = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
	        break;
        }
        case 1: bgd_dis_isid = bgd_K.edge(loc_isid); break;
        case 2: bgd_dis_isid = bgd_K.face(loc_isid); break;
	default: error_macro ("domain: unexpected side dimension " << sid_dim);
      }
      size_type bgd_ios_dis_isid = bgd_omega.dis_ige2ios_dis_ige (sid_dim, bgd_dis_isid);
      bgd_isid_is_on_domain.dis_entry     (bgd_dis_isid)     += 1;
      bgd_ios_isid_is_on_domain.dis_entry (bgd_ios_dis_isid) += 1;
    }
  }
  bgd_isid_is_on_domain.dis_entry_assembly(); // logical "or" across processes
  bgd_ios_isid_is_on_domain.dis_entry_assembly();

  // 1.2 counting & distribution for dom_isid 
  size_type dom_nsid = 0;
  for (size_type bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
    if (bgd_isid_is_on_domain[bgd_isid] != 0) dom_nsid++ ;
  }
  // 1.3 counting & distribution for dom_ios_isid 
  size_type dom_ios_nsid = 0;
  for (size_type bgd_ios_isid = 0, bgd_ios_nsid = bgd_omega.geo_element_ios_ownership(sid_dim).size(); bgd_ios_isid < bgd_ios_nsid; bgd_ios_isid++) {
    if (bgd_ios_isid_is_on_domain[bgd_ios_isid] != 0) dom_ios_nsid++ ;
  }
  // 1.4 numbering dom_isid & permutation: bgd_isid <--> dom_isid
  for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
                 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
    size_by_variant [variant] = 0;
  }
  distributor dom_isid_ownership (distributor::decide, comm, dom_nsid);
  size_type first_dom_dis_isid = dom_isid_ownership.first_index();
  bgd_isid2dom_dis_isid.resize (bgd_omega.geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
  dom_isid2bgd_isid.resize     (dom_isid_ownership,                       std::numeric_limits<size_type>::max());
  for (size_type dom_isid = 0, bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
    if (bgd_isid_is_on_domain[bgd_isid] == 0) continue;
    size_type dom_dis_isid = first_dom_dis_isid + dom_isid; // distributed case !
    bgd_isid2dom_dis_isid [bgd_isid] = dom_dis_isid;
    dom_isid2bgd_isid     [dom_isid] = bgd_isid;
    const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
    size_by_variant [bgd_S.variant()]++;
    dom_isid++;
  }
  // 1.5 numbering dom_ios_isid & permutation: bgd_ios_isid <--> dom_ios_isid
  distributor dom_ios_isid_ownership (distributor::decide, comm, dom_ios_nsid);
  size_type first_dom_ios_dis_isid = dom_ios_isid_ownership.first_index();
  disarray<size_type> bgd_ios_isid2dom_ios_dis_isid (bgd_omega.geo_element_ios_ownership(sid_dim), std::numeric_limits<size_type>::max());
  disarray<size_type> dom_ios_isid2bgd_ios_isid (dom_ios_isid_ownership, std::numeric_limits<size_type>::max());
  for (size_type dom_ios_isid = 0, bgd_ios_isid = 0, bgd_ios_nsid = bgd_omega.geo_element_ios_ownership(sid_dim).size(); bgd_ios_isid < bgd_ios_nsid; bgd_ios_isid++) {
    if (bgd_ios_isid_is_on_domain[bgd_ios_isid] == 0) continue;
    size_type dom_ios_dis_isid = first_dom_ios_dis_isid + dom_ios_isid;
    bgd_ios_isid2dom_ios_dis_isid [bgd_ios_isid] = dom_ios_dis_isid;
    dom_ios_isid2bgd_ios_isid [dom_ios_isid] = bgd_ios_isid;
    dom_ios_isid++;
  }
  // 1.6 permutation: bgd_isid --> dom_ios_isid
  disarray<size_type> bgd_isid2dom_ios_dis_isid (bgd_omega.geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
  for (size_type dom_ios_isid = 0; dom_ios_isid < dom_ios_nsid; dom_ios_isid++) {
    size_type bgd_ios_isid = dom_ios_isid2bgd_ios_isid [dom_ios_isid];
    size_type dom_ios_dis_isid = dom_ios_isid + first_dom_ios_dis_isid;
    size_type bgd_dis_isid = bgd_omega.ios_ige2dis_ige (sid_dim, bgd_ios_isid);
    bgd_isid2dom_ios_dis_isid.dis_entry (bgd_dis_isid) = dom_ios_dis_isid;
  }
  bgd_isid2dom_ios_dis_isid.dis_entry_assembly();

  // 1.7 permutations: dom_ios_isid <--> dom_isid
  dom_isid2dom_ios_dis_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
  // set reverse permutations for i/o:
  _ios_ige2dis_ige [sid_dim].resize (dom_ios_isid_ownership, std::numeric_limits<size_type>::max());
  for (size_type dom_isid = 0; dom_isid < dom_nsid; dom_isid++) {
    size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
    size_type dom_dis_isid = dom_isid + first_dom_dis_isid;
    size_type dom_ios_dis_isid = bgd_isid2dom_ios_dis_isid [bgd_isid];
    dom_isid2dom_ios_dis_isid [dom_isid] = dom_ios_dis_isid;
    _ios_ige2dis_ige [sid_dim].dis_entry (dom_ios_dis_isid) = dom_dis_isid;
  }
  _ios_ige2dis_ige [sid_dim].dis_entry_assembly();
  // ------------------------------------------------------------------------
  // 2. compute sizes and resize _geo_element[variant]
  // ------------------------------------------------------------------------
  size_type dom_nge = 0;
  size_type dom_dis_nge = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
                 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {

     distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);	
     geo_element::parameter_type param (variant, 1);
     base::_geo_element[variant].resize (dom_igev_ownership, param);
     base::_gs.ownership_by_variant[variant] = dom_igev_ownership;
     dom_nge     += dom_igev_ownership.size();
     dom_dis_nge += dom_igev_ownership.dis_size();
  }
  base::_gs.ownership_by_dimension[sid_dim] = distributor (dom_dis_nge, comm, dom_nge);

}
template <class T>
void
geo_rep<T,distributed>::domain_set_side_part2 (
    const domain_indirect_rep<distributed>&    indirect, 
    const geo_abstract_rep<T,distributed>&     bgd_omega,
    disarray<size_type>&                       bgd_iv2dom_dis_iv,
    size_type                                  sid_dim,
    disarray<size_type>&                       bgd_isid2dom_dis_isid,
    disarray<size_type>&                       dom_isid2bgd_isid,
    disarray<size_type>&                       dom_isid2dom_ios_dis_isid,
    size_type                                  size_by_variant [reference_element::max_variant])
{
  if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
  // ------------------------------------------------------------------------
  // 2) set _geo_element[variant] and S.set_ios_dis_ie
  //    also reperate external vertices
  // ------------------------------------------------------------------------
  std::set<size_type> ext_bgd_dis_iv_set;
  distributor dom_isid_ownership = dom_isid2bgd_isid.ownership();
  size_type first_dom_dis_isid = dom_isid_ownership.first_index();
  for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
    size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
    size_type dom_dis_isid = first_dom_dis_isid + dom_isid;
    size_type dom_ios_dis_isid = dom_isid2dom_ios_dis_isid [dom_isid];
    const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
    geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
    dom_S = bgd_S;
    dom_S.set_dis_ie     (dom_dis_isid);
    dom_S.set_ios_dis_ie (dom_ios_dis_isid);
    if (dom_S.dimension() == 0) continue;
    // set S face & edge : index, orientation & rotation will be set by propagate_numbering later
    for (size_type iloc = 0, nloc = dom_S.size(); iloc < nloc; iloc++) {
      size_type bgd_dis_inod = bgd_S[iloc];
      size_type bgd_dis_iv   = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
      if (! bgd_omega.geo_element_ownership(0).is_owned(bgd_dis_iv)) {
	ext_bgd_dis_iv_set.insert (bgd_dis_iv);
      }
    }
  }
  // ------------------------------------------------------------------------
  // 3) propagate new vertex numbering in all `dom_S' new sides
  // ------------------------------------------------------------------------
  if (sid_dim > 0) {
    bgd_iv2dom_dis_iv.append_dis_indexes (ext_bgd_dis_iv_set);
    for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
      geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
      for (size_type iloc = 0, nloc = dom_S.size(); iloc < nloc; iloc++) {
        size_type bgd_dis_inod = dom_S[iloc];
        size_type bgd_dis_iv   = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
        size_type dom_dis_iv   = bgd_iv2dom_dis_iv.dis_at (bgd_dis_iv);
        size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
        dom_S[iloc] = dom_dis_inod;
      }
    }
  }
}
/** build_from_domain : geo constructor

   Implementation notes:

   There are four vertices numbering:

    a) bgd_iv : background vertex index
             i.e. index of vertex as numbered in the backgroud geo omega

    b) bgd_ios_iv : background vertex index for i/o (nproc independent)

    c) dom_iv : index of a vertex as numbered in the current geo_domain
             this is a reduced set of vertices: vertices that not appears
             in the list of elements of the domain are skipped.
	     This numbering is constrained to follow the background
	     vertex distribution.

    d) dom_ios_iv : index of a vertex on the domain for i/o (nproc independent)
	     This numbering is constrained to follow the background
	     ios vertex distribution.

    The numbering bgd_iv and bgd_ios_iv are already defined.
    The numbering dom_iv and dom_ios_iv are defined here.
    The correspondances dom_iv <--> dom_ios_iv are also defined for i/o.

  Algorithm: for vertices renumbering

    1) scan the domain and mark vertex in the bgd numbering:
    bool bgd_iv_is_on_domain     [bgd_iv_ownership] = false
    bool bgd_ios_iv_is_on_domain [bgd_ios_iv_ownership] = false
    for K in omega.domain
      for iloc=0..K.size
        bgd_dis_iv = K[iloc]
        bgd_ios_dis_ie = omega.dis_ie2ios_dis_ie (0,bgd_iv)
        bgd_iv_is_on_domain     [[bgd_dis_iv]]     = true
        bgd_ios_iv_is_on_domain [[bgd_ios_dis_iv]] = true
    bgd_iv_is_on_domain.assembly
    bgd_ios_iv_is_on_domain.assembly
 
    2) counting & distribution for dom_iv 
    dom_nv = 0
    for bdg_iv = ...
      if bgd_iv_is_on_domain [bgd_iv] then
	dom_nv++
    dom_iv_ownership = distributor (-,comm,dom_nv)

    3) numbering dom_iv 
     & permutation: bgd_iv --> dom_iv
    disarray bgd_iv2dom_iv [bgd_iv_ownership]
    dom_iv = 0
    for bdg_iv = ...
      if bgd_iv_is_on_domain [bgd_iv] then
        bgd_iv2dom_iv [bgd_iv] = dom_iv
        dom_iv++

    4) counting & distribution for dom_ios_iv 
    dom_ios_nv = 0
    for bdg_ios_iv = ...
      if bgd_ios_iv_is_on_domain [bgd_ios_iv] then
	dom_ios_nv++
    dom_ios_iv_ownership = distributor (-,comm,dom_ios_nv)

    5) numbering dom_ios_iv 
     & permutation: bgd_ios_iv --> dom_ios_iv
    disarray bgd_ios_iv2dom_ios_iv [bgd_ios_iv_ownership]
    dom_ios_iv = 0
    for bdg_ios_iv = ...
      if bgd_ios_iv_is_on_domain [bgd_ios_iv] then
        bgd_ios_iv2dom_ios_iv [bgd_ios_iv] = dom_ios_iv
        dom_ios_iv++

    6) permutation: bgd_iv --> dom_ios_iv
    disarray bgd_iv2dom_ios_dis_iv   [bgd_iv_ownership]
    for bdg_ios_iv = ...
      if bgd_ios_iv_is_on_domain [bgd_ios_iv] then
        dom_ios_iv = bgd_ios_iv2dom_ios_iv [bgd_ios_iv]
	bgd_dis_iv = omega.ios_iv2dis_iv (bgd_ios_iv)
        bgd_iv2dom_ios_dis_iv [[bgd_dis_iv]] = dom_ios_dis_iv
    bgd_iv2dom_ios_dis_iv.assembly

    7) permutations: dom_ios_iv <--> dom_iv
    disarray dom_ios_iv2dom_dis_iv   [dom_iv_ownership]
    disarray dom_iv2dom_ios_dis_iv   [dom_ios_iv_ownership]
    for bdg_iv = ...
      if bgd_iv_is_on_domain [bgd_iv] then
        dom_iv = bgd_iv2dom_iv [bgd_iv]
        dom_ios_dis_iv = bgd_iv2dom_ios_dis_iv [bgd_iv]
        dom_iv2dom_ios_dis_iv [dom_iv] = dom_ios_dis_iv
        dom_dis_iv = first_dom_dis_iv + dom_iv
        dom_ios_iv2dom_dis_iv [[dom_ios_dis_iv]] = dom_dis_iv
    dom_ios_iv2dom_dis_iv.assembly

    8) loop on elements and copy from domain with new vertex numbering

  */ 

template <class T>
void
geo_rep<T,distributed>::build_from_domain (
        const domain_indirect_rep<distributed>&              indirect,
        const geo_abstract_rep<T,distributed>&               bgd_omega,
              std::map<size_type,size_type>&                 bgd_ie2dom_ie,
              std::map<size_type,size_type>&                 bgd_dis_ie2dom_dis_ie)
{
  base::_name = bgd_omega.name() + "[" + indirect.name() + "]";
  base::_version       = 4;
  base::_sys_coord = bgd_omega.coordinate_system();
  base::_dimension   = bgd_omega.dimension();
  base::_piola_basis = bgd_omega.get_piola_basis();
  base::_gs._map_dimension = indirect.map_dimension();
  size_type map_dim = base::_gs._map_dimension;
  size_type size_by_variant [reference_element::max_variant];
  std::fill (size_by_variant, size_by_variant+reference_element::max_variant, 0);
  // ------------------------------------------------------------------------
  // 1) _geo_element[0]: compact renumbering
  // ------------------------------------------------------------------------
  std::array<disarray<size_type,distributed>,4>  bgd_ige2dom_dis_ige;
  std::array<disarray<size_type,distributed>,4>  dom_ige2bgd_ige;
  disarray<size_type> dom_isid2dom_ios_dis_isid [4];
  domain_set_side_part1 (indirect, bgd_omega, 0,
	bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
        dom_isid2dom_ios_dis_isid [0], size_by_variant);
  domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], 0,
	bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
        dom_isid2dom_ios_dis_isid [0], size_by_variant);
  // ------------------------------------------------------------------------
  // 2) detect extern vertices & count elements by variants
  // ------------------------------------------------------------------------
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
    size_by_variant [variant] = 0;
  }
  std::set<size_type> ext_bgd_dis_iv_set;
  distributor ioige_ownership = indirect.ownership();
  size_type first_dis_ioige = ioige_ownership.first_index();
  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
    size_type ige = indirect.oige (ioige).index();
    bgd_ie2dom_ie [ige] = ioige;
    const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
    size_by_variant [bgd_K.variant()]++;
    for (size_type iloc = 0; iloc < bgd_K.size(); iloc++) {
      size_type bgd_dis_inod = bgd_K[iloc];
      size_type bgd_dis_iv   = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
      if (! bgd_omega.geo_element_ownership(0).is_owned(bgd_dis_iv)) {
	ext_bgd_dis_iv_set.insert (bgd_dis_iv);
      }
    }
  }
  bgd_ige2dom_dis_ige[0].append_dis_indexes (ext_bgd_dis_iv_set);
  // ------------------------------------------------------------------------
  // 3) _geo_element[map_dim]: compact vertex numbering
  // ------------------------------------------------------------------------
  // resize _geo_element[]
  size_type dis_nge = 0;
  size_type nge     = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {

     distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);	
     geo_element::parameter_type param (variant, 1);
     base::_geo_element[variant].resize (dom_igev_ownership, param);
     base::_gs.ownership_by_variant [variant] = dom_igev_ownership;
     dis_nge += dom_igev_ownership.dis_size();
     nge     += dom_igev_ownership.size();
  }
  base::_gs.ownership_by_dimension [map_dim] = distributor (dis_nge, base::comm(), nge);
  // ------------------------------------------------------------------------
  // 4) count all geo_element[] and set   base::_gs.ownership_by_variant[]
  //   => then can determine the node count (order > 1)
  // ------------------------------------------------------------------------
  for (size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
    domain_set_side_part1 (indirect, bgd_omega, sid_dim,
	bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
        dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
  }
  // ------------------------------------------------------------------------
  // 5) count nodes & set   base::_gs.node_ownership
  //  => then dis_iv2dis_inod works (used at step 6)
  // ------------------------------------------------------------------------
  std::array<size_type,reference_element::max_variant> loc_nnod_by_variant ;
  reference_element::init_local_nnode_by_variant (base::order(), loc_nnod_by_variant);
  size_type nnod = 0;
  for (size_type variant = 0;
                 variant < reference_element::last_variant_by_dimension(base::map_dimension());
                 variant++) {
    nnod += base::_gs.ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
  }
  distributor dom_node_ownership (distributor::decide, base::comm(), nnod);
  base::_gs.node_ownership = dom_node_ownership;
  // ------------------------------------------------------------------------
  // 6) set _geo_element[] values
  // ------------------------------------------------------------------------
  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
    size_type dis_ioige = first_dis_ioige + ioige;
    size_type ige = indirect.oige (ioige).index();
          geo_element& dom_K =           get_geo_element (map_dim, ioige);
    const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim,   ige);
    dom_K = bgd_K;
    dom_K.set_dis_ie (dis_ioige);
    size_type ini_dis_ioige = indirect.ioige2ini_dis_ioige (ioige);
    dom_K.set_ios_dis_ie (ini_dis_ioige);
    // set K face & edge : index, orientation & rotation will be set by propagate_numbering later
    for (size_type iloc = 0, nloc = dom_K.size(); iloc < nloc; iloc++) {
      size_type bgd_dis_inod = bgd_K[iloc];
      size_type bgd_dis_iv   = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
      size_type dom_dis_iv   = bgd_ige2dom_dis_ige[0].dis_at (bgd_dis_iv);
      size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
      dom_K[iloc] = dom_dis_inod;
    }
  }
  // reset also dom_ige2bgd_ige[map_dim] : used by _node[] for order > 1 geometries
  if (base::_gs._map_dimension > 0) {
    dom_ige2bgd_ige    [base::_gs._map_dimension].resize(indirect.ownership());
    bgd_ige2dom_dis_ige[base::_gs._map_dimension].resize(bgd_omega.geo_element_ownership(base::_gs._map_dimension), std::numeric_limits<size_type>::max());
    size_type first_dis_ioige = indirect.ownership().first_index();
    for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
      size_type bgd_ige = indirect.oige (ioige).index();
      size_type dis_ioige = first_dis_ioige + ioige;
      dom_ige2bgd_ige    [base::_gs._map_dimension] [ioige]   = bgd_ige; 
      bgd_ige2dom_dis_ige[base::_gs._map_dimension] [bgd_ige] = dis_ioige;
    }
  }
  // ------------------------------------------------------------------------
  // 7) _geo_element[1&2]: compact renumbering
  //    _ios_ige2dis_ige[1&2]: idem
  // ------------------------------------------------------------------------
  for (size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
    domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], sid_dim,
	bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
        dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
  }
  // ------------------------------------------------------------------------
  // 8) set reverse permutations for i/o
  // ------------------------------------------------------------------------
  _ios_ige2dis_ige [map_dim] = indirect._ini_ioige2dis_ioige;
  // ------------------------------------------------------------------------
  // 9) set geo_size
  // ------------------------------------------------------------------------
  /* There is still a problem with ios_ownership_by_variant[] :
      - first, its formaly "ini_ownership_by_variant" for domain_indirect
        since i/o are associated to "ini" numbering for domain_indirect.
      - physically, geo_elements of the background geo are never distributed
        with the "ini" ownership so, the associated number of variants are unavailable
     Notice that this "ini" domain_indirect numbering becomes the "ios" for the geo_domain.
     Also, the new numbering is denoted by "dom" for geo_domain.
     Thus, there are at least two solutions for geo_domain :
     1) on non-mixed domains, the variant disarray is not necessary, and assembly communications also
        can be avoided. Just identify the used variant related to map_dimension and store its
        number from dom_ige_domain.size in _ios_ownership_by_variant[].
     2) on mixed domains & for sid_dim=2,3 :
        via the new "dom_ige" numbering, access physically to elements and store variants in:
		disarray<size_t> dom_variant (dom_ige_ownership)
	then, permut it using ios_ige2dom_dis_ige and get
		disarray<size_t> ios_variant (ios_ige_ownership)
	finaly, count variants and store it in
  		_ios_ownership_by_variant[]
   */
  size_type  dis_size_by_variant [reference_element::max_variant];
  std::fill (dis_size_by_variant, dis_size_by_variant+reference_element::max_variant, 0);
  mpi::all_reduce (base::comm(), size_by_variant, reference_element::max_variant, dis_size_by_variant, std::plus<size_type>());

  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);
  ios_size_by_variant [reference_element::p] = _ios_ige2dis_ige[0].ownership().size();
  ios_size_by_variant [reference_element::e] = _ios_ige2dis_ige[1].ownership().size();

  std::vector<size_type> loc_ios_size_by_variant_by_proc [reference_element::max_variant];
  std::vector<size_type>     ios_size_by_variant_by_proc [reference_element::max_variant];
  for (size_type dim = 2; dim <= map_dim; dim++) {
    bool is_mixed = ((dim == 2) &&
                        (dis_size_by_variant[reference_element::t] != 0 &&
                         dis_size_by_variant[reference_element::q] != 0))
                      ||
                     ((dim == 3) &&
                       ((dis_size_by_variant[reference_element::T] != 0 &&
	                 dis_size_by_variant[reference_element::P] != 0) ||
	                (dis_size_by_variant[reference_element::P] != 0 &&
		         dis_size_by_variant[reference_element::H] != 0) ||
	                (dis_size_by_variant[reference_element::H] != 0 &&
	                 dis_size_by_variant[reference_element::T] != 0)));
    if (!is_mixed) {
      switch (dim) {
        case 2:
          if (        dis_size_by_variant[reference_element::t] != 0) {
                      ios_size_by_variant[reference_element::t] = _ios_ige2dis_ige[2].ownership().size();
          } else {
                      ios_size_by_variant[reference_element::q] = _ios_ige2dis_ige[2].ownership().size();
          }
	  break;
        case 3:
	default:
          if (        dis_size_by_variant[reference_element::T] != 0) {
                      ios_size_by_variant[reference_element::T] = _ios_ige2dis_ige[3].ownership().size();
          } else 
          if (        dis_size_by_variant[reference_element::P] != 0) {
                      ios_size_by_variant[reference_element::P] = _ios_ige2dis_ige[3].ownership().size();
          } else {
                      ios_size_by_variant[reference_element::H] = _ios_ige2dis_ige[3].ownership().size();
          }
	  break;
      }
      continue;
    }
    // here we have a mixed domain:
    size_type   nproc = base::comm().size();
    size_type my_proc = base::comm().rank();
    switch (dim) {
      case 2:
            ios_size_by_variant_by_proc [reference_element::t].resize(nproc, 0);
            ios_size_by_variant_by_proc [reference_element::q].resize(nproc, 0);
        loc_ios_size_by_variant_by_proc [reference_element::t].resize(nproc, 0);
        loc_ios_size_by_variant_by_proc [reference_element::q].resize(nproc, 0);
	break;
      case 3:
      default:
            ios_size_by_variant_by_proc [reference_element::T].resize(nproc, 0);
            ios_size_by_variant_by_proc [reference_element::P].resize(nproc, 0);
            ios_size_by_variant_by_proc [reference_element::H].resize(nproc, 0);
        loc_ios_size_by_variant_by_proc [reference_element::T].resize(nproc, 0);
        loc_ios_size_by_variant_by_proc [reference_element::P].resize(nproc, 0);
        loc_ios_size_by_variant_by_proc [reference_element::H].resize(nproc, 0);
	break;
    } 
    const distributor& ios_ige_ownership_dim = _ios_ige2dis_ige[dim].ownership(); // has been set by domain_set_sides
    for (size_type ie = 0, ne = base::sizes().ownership_by_dimension[dim].size(); ie < ne; ie++) {
        const geo_element& K = get_geo_element (dim, ie);
	size_type ios_dis_ie = K.ios_dis_ie();
	size_type iproc = ios_ige_ownership_dim.find_owner(ios_dis_ie);
        loc_ios_size_by_variant_by_proc [K.variant()][iproc]++;
    }
    switch (dim) {
      case 2:
	mpi::all_reduce (base::comm(),
	  loc_ios_size_by_variant_by_proc [reference_element::t].begin().operator->(), nproc,
	      ios_size_by_variant_by_proc [reference_element::t].begin().operator->(), std::plus<size_type>());
	mpi::all_reduce (base::comm(),
	  loc_ios_size_by_variant_by_proc [reference_element::q].begin().operator->(), nproc,
	      ios_size_by_variant_by_proc [reference_element::q].begin().operator->(), std::plus<size_type>());
	ios_size_by_variant[reference_element::t] = ios_size_by_variant_by_proc [reference_element::t][my_proc];
	ios_size_by_variant[reference_element::q] = ios_size_by_variant_by_proc [reference_element::q][my_proc];
	break;
      case 3:
      default:
	// TODO: TPH
        error_macro ("3D domain \""<<base::_name<<"\" with mixed element variants: not yet");
	break;
    }
    // end of domain with mixed elements
  }
  // then, ios_size_by_variant[] is completed and we set ios_sizes:
  for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
    size_type first_ios_v = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(dim);
                   variant < reference_element:: last_variant_by_dimension(dim); variant++) {

      _ios_gs.ownership_by_variant [variant] = distributor (distributor::decide, base::comm(), ios_size_by_variant [variant]);
      _ios_gs.first_by_variant     [variant] = distributor (distributor::decide, base::comm(), first_ios_v);
      first_ios_v += ios_size_by_variant [variant];
    }
    size_type ios_nge = first_ios_v;
    _ios_gs.ownership_by_dimension [dim] = distributor (distributor::decide, base::comm(), ios_nge);
  }
  // ------------------------------------------------------------------------
  // 10) resize _igev2ios_dis_igev : used by Pk_numbering
  // ------------------------------------------------------------------------
  for (size_type variant = 0;
                 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension);
                 variant++) {
    _igev2ios_dis_igev [variant].resize (
			base::_gs.ownership_by_variant [variant], 
			std::numeric_limits<size_type>::max());
  }
  // for dim=0,1 : no variants and igev2ios_dis_igev[] = ige2ios_dis_ige[]
  for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
    size_type ige = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(dim);
                   variant < reference_element:: last_variant_by_dimension(dim);
                   variant++) {
      for (size_type igev = 0, ngev = _igev2ios_dis_igev [variant].size(); igev < ngev; igev++, ige++) {
        const geo_element& K = base::_geo_element [variant][igev];
        size_type ios_dis_ige       = K.ios_dis_ie();
        size_type iproc             = _ios_gs.ownership_by_dimension [dim].find_owner(ios_dis_ige);
        size_type first_ios_dis_ige = _ios_gs.ownership_by_dimension [dim].first_index (iproc);
	assert_macro (ios_dis_ige >= first_ios_dis_ige, "invalid index");
        size_type ios_ige = ios_dis_ige - first_ios_dis_ige;
        size_type first_v  = _ios_gs.first_by_variant [variant].size (iproc);
	assert_macro (ios_ige >= first_v, "invalid index");
        size_type ios_igev = ios_ige - first_v;
        size_type first_ios_dis_igev = _ios_gs.ownership_by_variant [variant].first_index (iproc);
        size_type ios_dis_igev = first_ios_dis_igev + ios_igev;
        _igev2ios_dis_igev [variant][igev] = ios_dis_igev;
      }
#ifdef TODO
      _ios_igev2dis_igev [variant].resize (
	_ios_gs.ownership_by_variant [variant], // this ownership is not yet set
	std::numeric_limits<size_type>::max());
      _igev2ios_dis_igev [variant].reverse_permutation (_ios_igev2dis_igev [variant]); // not used ?
#endif // TODO
    }
  }
  // ------------------------------------------------------------------------
  // 11) raw copy _node
  // TODO: use shallow copy & indirection :
  //	node(inod) { return bgd_omega.node(inod2bgd_inod(inod))}
  // ------------------------------------------------------------------------
  // 11.a) resize _node[]
  set_ios_permutation (_inod2ios_dis_inod);
  distributor ios_dom_node_ownership (dom_node_ownership.dis_size(), base::comm(), distributor::decide);
  _ios_inod2dis_inod.resize (ios_dom_node_ownership);
  _inod2ios_dis_inod.reverse_permutation (_ios_inod2dis_inod);

  // 11.b) copy _node[] from bgd_omega.node[]
  base::_node.resize (dom_node_ownership);
  disarray<size_type> dom_inod2bgd_inod (dom_node_ownership, std::numeric_limits<size_type>::max());
  size_type dom_inod = 0;
  size_type first_bgd_inod_v = 0;
  for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
    size_type dom_ige = 0;
    for (size_type variant = reference_element::first_variant_by_dimension(dim);
                   variant < reference_element:: last_variant_by_dimension(dim);
                   variant++) {
      if (loc_nnod_by_variant [variant] == 0) continue;
      size_type first_bgd_v = bgd_omega.sizes().first_by_variant [variant].size();
      for (size_type dom_igev = 0, dom_ngev = base::_geo_element [variant].size(); dom_igev < dom_ngev; dom_igev++, dom_ige++) {
        const geo_element& K = base::_geo_element [variant][dom_igev];
        size_type bgd_ige = dom_ige2bgd_ige [dim][dom_ige];
        assert_macro (bgd_ige >= first_bgd_v, "invalid index");
        size_type bgd_igev = bgd_ige - first_bgd_v;
  	for (size_type loc_inod = 0, loc_nnod = loc_nnod_by_variant [variant]; loc_inod < loc_nnod; loc_inod++, dom_inod++) {
          size_type bgd_inod = first_bgd_inod_v + bgd_igev * loc_nnod_by_variant [variant] + loc_inod;
          dom_inod2bgd_inod [dom_inod] = bgd_inod;
	  base::_node [dom_inod] = bgd_omega.node (bgd_inod);
        }
      }
      first_bgd_inod_v += bgd_omega.sizes().ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
    }
  }
  // ------------------------------------------------------------------------
  // 12) propagate new edge & face numbering to all `dom_S' elements
  // ------------------------------------------------------------------------
  for (size_type dim = 1; dim < base::_gs._map_dimension; dim++) {
    set_element_side_index (dim);
  }
  // ------------------------------------------------------------------------
  // 13) reset vertex P[0] value (useful for band zero vertex domain) 
  // note: should be before "build_external_entities" as some vertex will be exported
  // ------------------------------------------------------------------------
  size_type first_dom_dis_iv = base::_geo_element[reference_element::p].ownership().first_index();
  for (size_type dom_iv = 0, dom_nv = base::_geo_element[reference_element::p].size(); dom_iv < dom_nv; dom_iv++) {
    geo_element& P = base::_geo_element[reference_element::p][dom_iv];
    size_type dom_dis_iv = first_dom_dis_iv + dom_iv;
    size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
    P[0] = dom_dis_inod;
  }
  // ------------------------------------------------------------------------
  // 17) bgd_dis_ie2dom_dis_ie, used by space("square[sides]",Pk) for HDG lambda multiplier
  // ------------------------------------------------------------------------
  // step 0: get direct access to bgd_omega _geo_element table
  const geo_base_rep<T,distributed>* ptr_bgd_omega = dynamic_cast<const geo_base_rep<T,distributed>*>(&bgd_omega);
  check_macro (ptr_bgd_omega != 0, "invalid bgd_omega");
  const geo_base_rep<T,distributed>& bgd_omega1 = *ptr_bgd_omega;
  
  // step 1: list exported elements from bgd_omega
  size_type map_d = base::_gs._map_dimension;
  // size_type first_bgd_igev_by_variant = 0;
  index_set ext_bgd_dis_ie_set;
  for (size_type variant = reference_element::first_variant_by_dimension (map_d);
                 variant < reference_element:: last_variant_by_dimension (map_d); ++variant) {
    const std::map<size_type,geo_element_auto<>>& ext_bgd_gev = bgd_omega1._geo_element [variant].get_dis_map_entries();
    for (auto x: ext_bgd_gev) {
      size_type bgd_dis_igev          = x.first;
      const geo_element_auto<>& bgd_K = x.second;
      size_type bgd_dis_ie            = bgd_K.dis_ie(); // = bgd_dis_igev + first_bgd_igev_by_variant;
      ext_bgd_dis_ie_set.insert (bgd_dis_ie);
    }
    // first_bgd_igev_by_variant += bgd_omega1._geo_element [variant].dis_size();
  }
  bgd_ige2dom_dis_ige[map_d].append_dis_indexes (ext_bgd_dis_ie_set);

  // step 2: copy bgd_ige2dom_dis_ige[map_d].externals since this array is a temporary
  // also: build ext_dom_dis_igev_set [variant], for each variant for map_d
  index_set ext_dom_dis_ie_set;
  const std::map <size_type, size_type>& bgd_dis_ie2dom_dis_ie_tmp = bgd_ige2dom_dis_ige[map_d].get_dis_map_entries();
  std::array<index_set,reference_element::max_variant> ext_dom_dis_igev_set;
  for (auto x : bgd_dis_ie2dom_dis_ie_tmp) {
    size_type bgd_dis_ie = x.first;
    size_type dom_dis_ie = x.second;
    if (dom_dis_ie == std::numeric_limits<size_type>::max()) {
      // bgd_dis_ie is not part of the present domain 
      continue;
    }
    bgd_dis_ie2dom_dis_ie [bgd_dis_ie] = dom_dis_ie;
    ext_dom_dis_ie_set.insert (dom_dis_ie);
    // convert dom_dis_ie to dom_dis_igev and get its variant:
    size_type variant;
    size_type dom_dis_igev = base::sizes().dis_ige2dis_igev_by_dimension (map_d, dom_dis_ie, variant);
    ext_dom_dis_igev_set [variant].insert (dom_dis_igev);
  }
  // step 3: communicate domain external elements 
  for (size_type variant = reference_element::first_variant_by_dimension (map_d);
                 variant < reference_element:: last_variant_by_dimension (map_d); ++variant) {
    base::_geo_element [variant].append_dis_indexes (ext_dom_dis_igev_set[variant]);
  }
  // ------------------------------------------------------------------------
  // 14) external entities
  // ------------------------------------------------------------------------
  build_external_entities();
  // ------------------------------------------------------------------------
  // 15) domains : intersection with the current domain
  // ------------------------------------------------------------------------
  const geo_base_rep<T,distributed> *ptr = dynamic_cast<const geo_base_rep<T,distributed>*> (&bgd_omega);
  check_macro (ptr != 0, "cannot build domains on \""<<base::_name<<"\"");
  const geo_base_rep<T,distributed>& bgd_omega2 = *ptr;
  for (size_type idom = 0, ndom = bgd_omega2.n_domain_indirect(); idom < ndom; ++idom) {
    const domain_indirect_basic<distributed>& bgd_dom = bgd_omega2.get_domain_indirect(idom);
    if (bgd_dom.name() == indirect.name()) continue; // myself
    size_type dom_map_dim = bgd_dom.map_dimension();
    if (dom_map_dim > base::_gs._map_dimension) continue; // dom has upper dimension
    std::vector<size_type> ie_list;
    size_type first_dom_dis_ige = base::_gs.ownership_by_dimension[dom_map_dim].first_index();
    for (domain_indirect_basic<distributed>::const_iterator_ioige iter = bgd_dom.ioige_begin(), last = bgd_dom.ioige_end(); iter != last; ++iter) {
      const geo_element_indirect& ioige = *iter;
      size_type bgd_ige = ioige.index();
      size_type dom_dis_ige = bgd_ige2dom_dis_ige [dom_map_dim][bgd_ige];
      if (dom_dis_ige == std::numeric_limits<size_type>::max()) continue; // side do not belongs to dom
      check_macro (dom_dis_ige >= first_dom_dis_ige, "invalid index");
      size_type dom_ige = dom_dis_ige - first_dom_dis_ige;
      ie_list.push_back(dom_ige);
    }
    size_type ie_list_dis_size = mpi::all_reduce (base::comm(), ie_list.size(), std::plus<size_type>());
    if (ie_list_dis_size == 0) {
      continue; // empty intersection
    }
    domain_indirect_basic<distributed> dom (*this, bgd_dom.name(), bgd_dom.map_dimension(), base::comm(), ie_list);
    base::_domains.push_back (dom);
  }
  // ------------------------------------------------------------------------
  // 16) some post treatments (after build_external_entities and so one!)
  // ------------------------------------------------------------------------
  base::compute_bbox();
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_rep<Float,distributed>;

} // namespace rheolef
#endif // _RHEOLEF_HAVE_MPI