geo.h

Go to the documentation of this file.

#ifndef _RHEOLEF_GEO_H
#define _RHEOLEF_GEO_H
//
// 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 distributed 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
//
// =========================================================================
// author: Pierre.Saramito@imag.fr
// date: 20 april 2013
 
namespace rheolef {
} // namespace rheolef
 
/*
  Implementation note: geo = mesh data structure
 
  geo is an alias to geo_basic<T,M>
        T: the current Float representation (e.g. double)
        M: the memory model (sequential or distributed)
           the default memory model has been chosen by the configure script
 
  geo_basic<T,M> is a smart_pointer_clone<geo_abstract_rep<T,M> >
        i.e. a pointer with shallow copy semantic on a pure virtual base class geo_abstract_rep<T,M>
        used by polymorphic hierarchy of classes.
 
  there are three concrete variants for this base class:
    geo_rep<T,M>                 : for usual meshes, as omega
    geo_domain_indirect_rep<T,M> : for usual boundary domains, as gamma=omega["boundary"]
    geo_domain_rep<T,M>          : for compacted boundary domains, as gamma=omega["boundary"]
                                   used when defining a space on a boundary domain, as
                                        space W (omega["boundary"], "P1");
                                   => elements and vertex are renumbered in a compact form
                                      for easy P1 dofs numbering.
 
   Most code is shared by these three classes and by the sequential/distributed variants.
   This leads to the following hierarchy of classes:
 
        geo_abstract_base_rep<T>             : virtual methods for M=seq
          geo_abstract_rep<T,M>              : M=seq,dis have separated impl; add methods for M=dis
            geo_domain_indirect_base_rep<T,M>: data + methods for domains
              geo_domain_indirect_rep<T,M>   : M=seq,dis have separated impl; add methods for M=dis
            geo_base_rep<T,M>                : data + methods for geometry
              geo_rep<T,M>                   : M=seq,dis have separated impl: add methods for M=dis
                geo_domain_rep<T,M>          : data + methods for compacted geometry on domain
*/
 
#include "rheolef/iorheo.h"
#include "rheolef/diststream.h"
#include "rheolef/linalg.h"
#include "rheolef/basis.h"
#include "rheolef/hack_array.h"
#include "rheolef/geo_size.h"
#include "rheolef/domain_indirect.h"
#include "rheolef/geo_header.h"
#include "rheolef/space_constant.h"
#include "rheolef/geo_locate.h"
#include "rheolef/geo_trace_ray_boundary.h"
#include "rheolef/geo_nearest.h"
#include "rheolef/rounder.h"
#include <unordered_map>
 
namespace rheolef {
 
namespace details {
 
struct zero_dimension {}; // used by space IR
 
} // namespace details
// =========================================================================
// point io helpers
// =========================================================================
template <class T>
struct _point_get {
  typedef typename point_basic<T>::size_type size_type;
  _point_get (size_type d1) : d(d1) {}
  size_type d;
  std::istream& operator() (std::istream& is, point_basic<T>& x) { return x.get (is, d); }
}; 

template <class T>
struct _point_put {
  typedef typename point_basic<T>::size_type size_type;
  _point_put (size_type d1) : d(d1) {}
  size_type d;
  std::ostream& operator() (std::ostream& os, const point_basic<T>& x) { return x.put (os, d); }
};

template <class T>
struct _round_point_put {
  typedef typename point_basic<T>::size_type size_type;
  _round_point_put (size_type d1, const T& eps1) : d(d1), round(eps1) {}
  size_type       d;
  rounder_type<T> round;
  std::ostream& operator() (std::ostream& os, const point_basic<T>& x) {
    for (size_t i = 0; i < d; ++i) {
      os << round(x[i]);
      if (i+1 != d) os << " ";
    }
    return os;
  }
};
 
// forward declaration:
template <class T, class M> class geo_basic;
template <class T, class M> class geo_domain_indirect_rep;
 
// =========================================================================
// =========================================================================
template<class T, class Ref, class Ptr, class IteratorByVariant>
struct geo_iterator {
  typedef geo_iterator<T,Ref,Ptr,IteratorByVariant>                  _self;
  typedef geo_iterator<T,T&,T*,typename hack_array<T>::iterator>  _nonconst_iterator;
 
// typedefs
 
  // see std::deque<T>::iterator : TODO: complete to a full random iterator
  typedef std::random_access_iterator_tag iterator_category;
  typedef T                               value_type;
  typedef Ptr                             pointer;
  typedef Ref                             reference;
  typedef typename T::size_type           size_type;
  typedef ptrdiff_t                       difference_type;
 
// allocators:
 
  template<class Geo>
  geo_iterator (size_type dim, size_type variant, IteratorByVariant iter, Geo& omega);
  geo_iterator (const _nonconst_iterator& y); // conversion from iter to const_iter
 
// accessors & modifiers:
 
  reference operator* () const { return *_iter_by_var; }
  pointer   operator->() const { return _iter_by_var.operator->(); }
 
  _self& operator++ () { 
    ++_iter_by_var;
    _reset_to_next_or_last();
    return *this;
  }
  void _reset_to_next_or_last () {
    while (_variant < _variant_max && _iter_by_var == _last_by_var [_variant]) {
      ++_variant; // then variant <= variant_max
      if (_variant < _variant_max) {
        _iter_by_var = _first_by_var [_variant];
      }
    }
  }
  _self operator++ (int) { _self tmp = *this; operator++(); return tmp; }
 
  bool operator== (const _self& y) const { return _variant == y._variant && _iter_by_var == y._iter_by_var; }
  bool operator!= (const _self& y) const { return ! operator== (y); }
 
// data:
 
  size_type                       _variant;
  size_type                       _variant_max;
  IteratorByVariant               _iter_by_var;
  IteratorByVariant               _first_by_var [reference_element::max_variant];
  IteratorByVariant                _last_by_var [reference_element::max_variant];
};
// =========================================================================
// =========================================================================
template <class T>
class geo_abstract_base_rep {
public:
 
// typedefs
 
    typedef enum { 
        geo                 = 0, 
        geo_domain          = 1, 
        geo_domain_indirect = 2, 
        max_variant         = 3 
    } geo_rep_variant_type; 
 
    typedef geo_element_hack::size_type                    size_type;
    typedef point_basic<T>                                 node_type;
    typedef reference_element::variant_type                variant_type;
    typedef geo_element&                                   reference;
    typedef const geo_element&                             const_reference;
    typedef space_constant::coordinate_type                coordinate_type;
    typedef std::map<size_type,geo_element_auto<>>         geo_element_map_type;
 
    typedef typename hack_array<geo_element_hack>::iterator
            iterator_by_variant;
    typedef geo_iterator<
              geo_element,
              geo_element&,
              geo_element*,
              iterator_by_variant>
            iterator;
 
    typedef typename hack_array<geo_element_hack>::const_iterator
            const_iterator_by_variant;
    typedef geo_iterator<
              geo_element,
              const geo_element&,
              const geo_element*,
              const_iterator_by_variant>
            const_iterator;
 
// allocators:
 
    geo_abstract_base_rep () {}
    virtual ~geo_abstract_base_rep () {}
 
// abstract accessors:
 
    virtual size_type variant() const = 0;
    virtual std::string name() const = 0;
    virtual std::string familyname() const = 0;
    virtual size_type dimension() const = 0;
    virtual size_type serial_number() const = 0;
    virtual size_type map_dimension() const = 0;
    virtual bool is_broken() const = 0;
    virtual coordinate_type coordinate_system() const = 0;
    virtual const basis_basic<T>& get_piola_basis() const = 0;
    virtual const node_type& xmin() const = 0;
    virtual const node_type& xmax() const = 0;
    virtual const T&         hmin() const = 0;
    virtual const T&         hmax() const = 0;
    virtual const geo_size&  sizes() const = 0;
    virtual const geo_size&  ios_sizes() const = 0;
 
    virtual const distributor&  geo_element_ownership (size_type dim) const = 0;
    virtual const_reference get_geo_element           (size_type dim, size_type ige) const = 0;
    virtual const geo_element& bgd2dom_geo_element (const geo_element& bgd_K) const { return bgd_K; }
    virtual const geo_element& dom2bgd_geo_element (const geo_element& dom_K) const { return dom_K; }
    virtual size_type neighbour (size_type ie, size_type loc_isid) const = 0;
    virtual void neighbour_guard() const = 0;
 
    virtual const_iterator_by_variant begin_by_variant (variant_type variant) const = 0;
    virtual const_iterator_by_variant   end_by_variant (variant_type variant) const = 0;
    virtual const geo_element_map_type& get_external_geo_element_map (size_type variant) const = 0;
 
    const_iterator begin (size_type dim) const;
    const_iterator end   (size_type dim) const;
 
    virtual size_type n_node() const = 0;
    virtual const node_type& node   (size_type inod) const = 0;
    virtual const node_type& dis_node (size_type dis_inod) const = 0;
    virtual void dis_inod (const geo_element& K, std::vector<size_type>& dis_inod) const = 0;
    virtual size_type dis_inod2dis_iv (size_type dis_inod) const = 0;
 
    virtual size_type n_domain_indirect () const = 0;
    virtual bool have_domain_indirect (const std::string& name) const = 0;
    virtual void reset_order (size_type order) = 0;
 
    virtual size_type seq_locate (
        const point_basic<T>& x,
        size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
    virtual size_type dis_locate (
        const point_basic<T>& x,
        size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
 
    virtual size_type seq_trace_move (
        const point_basic<T>&     x,
        const point_basic<T>&     v,
              point_basic<T>&     y) const = 0;
    virtual size_type dis_trace_move (
        const point_basic<T>&     x,
        const point_basic<T>&     v,
              point_basic<T>&     y) const = 0;
 
    virtual size_type seq_nearest (
        const point_basic<T>&    x,
              point_basic<T>&    x_nearest) const = 0;
    virtual size_type dis_nearest (
        const point_basic<T>&    x,
              point_basic<T>&    x_nearest) const = 0;
 
// virtual i/o:
 
    virtual odiststream& put         (odiststream& ops) const = 0;
    virtual bool check(bool verbose) const = 0;
 
// deduced comparator:
 
    bool operator== (const geo_abstract_base_rep<T>& omega2) const {
        return name() == omega2.name(); }
};
template<class T, class Ref, class Ptr, class IteratorByVariant>
template<class Geo>
geo_iterator<T,Ref,Ptr,IteratorByVariant>::geo_iterator (
  size_type                       dim,
  size_type                       variant,
  IteratorByVariant               iter,
  Geo&                            omega)
 : _variant     (variant),
   _variant_max (reference_element::last_variant_by_dimension(dim)),
   _iter_by_var (iter),
   _first_by_var(),
    _last_by_var()
{
   for (size_type variant = reference_element::first_variant_by_dimension(dim);
                  variant < reference_element::last_variant_by_dimension(dim); variant++) {
     _first_by_var [variant] = omega.begin_by_variant (variant);
      _last_by_var [variant] = omega.  end_by_variant (variant);
   }
   _reset_to_next_or_last();
}
template<class T, class Ref, class Ptr, class IteratorByVariant>
geo_iterator<T,Ref,Ptr,IteratorByVariant>::geo_iterator (const _nonconst_iterator& y)
 : _variant     (y._variant),
   _variant_max (y._variant_max),
   _iter_by_var (y._iter_by_var),
   _first_by_var(),
    _last_by_var()
{
   std::copy (y._first_by_var, y._first_by_var + reference_element::max_variant, _first_by_var);
   std::copy (y. _last_by_var, y. _last_by_var + reference_element::max_variant,  _last_by_var);
}
// =========================================================================
// =========================================================================
template <class T, class M>
class geo_abstract_rep {};
 
template <class T>
class geo_abstract_rep<T,sequential> : public geo_abstract_base_rep<T> {
public:
 
// typedefs
 
    typedef geo_abstract_base_rep<T>                  base;
    typedef typename base::size_type                  size_type;
    typedef typename base::node_type                  node_type;
    typedef typename base::variant_type               variant_type;
    typedef typename base::iterator                   iterator;
    typedef typename base::const_iterator             const_iterator;
    typedef typename base::iterator_by_variant        iterator_by_variant;
    typedef typename base::const_iterator_by_variant  const_iterator_by_variant;
    typedef typename base::reference                  reference;
    typedef typename base::const_reference            const_reference;
 
// allocators:
 
    geo_abstract_rep () {}
    virtual geo_abstract_rep<T,sequential>* clone() const = 0;
    virtual ~geo_abstract_rep () {}
 
// abstract accessors:
 
    virtual const domain_indirect_basic<sequential>& get_domain_indirect (size_type i) const = 0;
    virtual const domain_indirect_basic<sequential>& get_domain_indirect (const std::string& name) const = 0;
    virtual void  insert_domain_indirect (const domain_indirect_basic<sequential>& dom) const = 0;
 
    virtual const disarray<node_type,sequential>& get_nodes() const = 0;
    virtual void set_nodes (const disarray<node_type,sequential>&) = 0;
 
    virtual void locate (
                const disarray<point_basic<T>, sequential>&    x,
                      disarray<size_type, sequential>&         dis_ie,
                bool  do_check = false) const = 0;
    virtual void trace_ray_boundary (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y,
                bool  do_check = false) const = 0;
    virtual void trace_move (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y) const = 0;
    virtual void nearest (
                const disarray<point_basic<T>,sequential>&     x,
                      disarray<point_basic<T>,sequential>&     x_nearest,
                      disarray<size_type, sequential>&         dis_ie) const = 0;
};
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
class geo_abstract_rep<T,distributed> : public geo_abstract_base_rep<T> {
public:
 
// typedefs
 
    typedef geo_abstract_base_rep<T>       base;
    typedef typename base::size_type       size_type;
    typedef typename base::node_type       node_type;
    typedef typename base::const_reference const_reference;
    typedef std::map <size_type, node_type, std::less<size_type>,
            heap_allocator<std::pair<size_type,node_type> > >    node_map_type;
 
// allocators:
 
    geo_abstract_rep () {}
    virtual geo_abstract_rep<T,distributed>* clone() const = 0;
    virtual ~geo_abstract_rep () {}
 
// abstract accessors:
 
    virtual distributor geo_element_ios_ownership (size_type dim) const = 0;
    virtual const_reference   dis_get_geo_element (size_type dim, size_type dis_ige) const = 0;
    virtual size_type             ige2ios_dis_ige (size_type dim, size_type     ige) const = 0;
    virtual size_type         dis_ige2ios_dis_ige (size_type dim, size_type dis_ige) const = 0;
    virtual size_type             ios_ige2dis_ige (size_type dim, size_type ios_ige) const = 0;
 
    virtual const domain_indirect_basic<distributed>& get_domain_indirect (size_type i) const = 0;
    virtual const domain_indirect_basic<distributed>& get_domain_indirect (const std::string& name) const = 0;
    virtual void insert_domain_indirect (const domain_indirect_basic<distributed>& dom) const = 0;
 
    virtual const disarray<node_type,distributed>& get_nodes() const = 0;
    virtual void set_nodes (const disarray<node_type,distributed>&) = 0;
 
    virtual void locate (
                const disarray<point_basic<T>,distributed>&     x,
                      disarray<size_type, distributed>&         dis_ie,
                bool  do_check = true) const = 0;
    virtual void trace_ray_boundary (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y,
                bool  do_check = false) const = 0;
    virtual void trace_move (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y) const = 0;
    virtual void nearest (
                const disarray<point_basic<T>,distributed>&     x,
                      disarray<point_basic<T>,distributed>&     x_nearest,
                      disarray<size_type, distributed>&         dis_ie) const = 0;
 
// utility:
 
    virtual void set_ios_permutation (disarray<size_type,distributed>& idof2ios_dis_idof) const = 0;
 
    // used by space_constritution for ios numbering
    virtual const std::array<disarray<size_type,distributed>,reference_element::max_variant>&
    get_igev2ios_dis_igev() const = 0;
 
};
#endif // _RHEOLEF_HAVE_MPI
// =========================================================================
// =========================================================================
// NOTE: since geo_rep<seq> contains sequential arrays for vertices and elts,
//  the geo_rep<mpi> cannot derive from geo_rep<seq>. The solution is to
//  derive both geo_rep<seq> and geo_rep<mpi> classes from a generic base class
//  named geo_base_rep that takes the memory model (seq or mpi) as template argument.
template <class T, class M>
class geo_base_rep : public geo_abstract_rep<T,M> {
public:
// typedefs:
 
    typedef geo_abstract_rep<T,M>                     base;
    typedef typename base::size_type                  size_type;
    typedef typename base::node_type                  node_type;
    typedef typename base::variant_type               variant_type;
    typedef typename base::iterator                   iterator;
    typedef typename base::const_iterator             const_iterator;
    typedef typename base::iterator_by_variant        iterator_by_variant;
    typedef typename base::const_iterator_by_variant  const_iterator_by_variant;
    typedef typename base::reference                  reference;
    typedef typename base::const_reference            const_reference;
    typedef typename base::coordinate_type            coordinate_type;
 
// allocators:
 
    geo_base_rep ();
    geo_base_rep (const geo_base_rep<T,M>&);
 
    void build_from_list (
      const geo_basic<T,M>&                               lambda,
      const disarray<point_basic<T>,M>&                   node_list,
      const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,M>,
                         reference_element::max_variant>& elt_list);
 
    ~geo_base_rep();
 
// abstract accessors defined:
  
    size_type          variant() const { return geo_abstract_base_rep<T>::geo; }
    std::string     familyname() const { return _name; }
    std::string           name() const;
    size_type    serial_number() const { return _serial_number; }
    size_type        dimension() const { return _dimension; }
    size_type    map_dimension() const { return _gs.map_dimension(); }
    bool             is_broken() const { return false; }
    coordinate_type coordinate_system() const { return _sys_coord; }
    void set_coordinate_system (coordinate_type sys_coord) { _sys_coord = sys_coord; }
    void set_name (std::string name) { _name = name; }
    void set_dimension (size_type dim) { _dimension = dim; }
    void set_serial_number(size_type i) { _serial_number = i; }
    const basis_basic<T>& get_piola_basis() const { return _piola_basis; }
    const node_type&      xmin() const { return _xmin; }
    const node_type&      xmax() const { return _xmax; }
    const T&              hmin() const { return _hmin; }
    const T&              hmax() const { return _hmax; }
    const geo_size&      sizes() const { return _gs; }
    const geo_size&  ios_sizes() const { return _gs; }
 
    const distributor& geo_element_ownership (size_type dim) const { return _gs.ownership_by_dimension[dim]; }
    const_reference     get_geo_element (size_type dim, size_type ige) const;
    const_reference dis_get_geo_element (size_type dim, size_type dis_ige) const;
 
    const_iterator_by_variant begin_by_variant (variant_type variant) const;
          iterator_by_variant begin_by_variant (variant_type variant);
    const_iterator_by_variant   end_by_variant (variant_type variant) const;
          iterator_by_variant   end_by_variant (variant_type variant);
 
    const node_type&     node (size_type     inod) const { return _node [inod]; }
    const node_type& dis_node (size_type dis_inod) const { return _node.dis_at (dis_inod); }
    const node_type&     node (const geo_element& K, size_type loc_inod) const;
    void             dis_inod (const geo_element& K, std::vector<size_type>& dis_inod) const;
    const disarray<node_type,M>& get_nodes() const { return _node; }
 
    void set_nodes (const disarray<node_type,M>& x) { _node = x; _node.reset_dis_indexes(); compute_bbox(); }
    size_type n_domain_indirect () const { return _domains.size(); }
    bool have_domain_indirect (const std::string& name) const;
    const domain_indirect_basic<M>& get_domain_indirect (size_type i) const { return _domains[i]; }
    const domain_indirect_basic<M>& get_domain_indirect (const std::string& name) const;
    void  insert_domain_indirect (const domain_indirect_basic<M>& dom) const;
 
    size_type seq_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const;
    size_type dis_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const;
    size_type seq_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const;
    size_type dis_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const;
    size_type seq_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const;
    size_type dis_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const;
 
    size_type neighbour (size_type ie, size_type loc_isid) const;
    void neighbour_guard() const;
 
// additional accessors & modifier:
 
    reference get_geo_element (size_type dim, size_type ige);
    iterator begin (size_type dim);
    iterator end   (size_type dim);
 
// deduced accessors:
 
    size_type     size(size_type dim) const;
    size_type dis_size(size_type dim) const;
    const distributor&        ownership() const { return geo_element_ownership (map_dimension()); }
    const distributor& vertex_ownership() const { return geo_element_ownership (0); }
    const communicator&            comm() const { return ownership().comm(); }
    size_type            order() const { return get_piola_basis().degree(); }
 
    size_type        n_node() const { return _node.    size(); }
    size_type    dis_n_node() const { return _node.dis_size(); }
    size_type      n_vertex() const { return     size (0); }
    size_type          size() const { return     size (map_dimension()); }
    size_type  dis_n_vertex() const { return dis_size (0); }
    size_type      dis_size() const { return dis_size (map_dimension()); }
    size_type    dis_n_edge() const { return dis_size (1); }
    size_type    dis_n_face() const { return dis_size (2); }
    size_type dis_inod2dis_iv (size_type dis_inod) const;
    size_type dis_iv2dis_inod (size_type dis_iv)   const;
 
    const_reference operator[] (size_type ie) const { return get_geo_element (map_dimension(), ie); }
    reference       operator[] (size_type ie)       { return get_geo_element (map_dimension(), ie); }
 
    const_iterator begin (size_type dim) const { return base::begin (dim); }
    const_iterator end   (size_type dim) const { return base::end   (dim); }
 
    const_iterator begin()               const { return base::begin (map_dimension()); }
    const_iterator end()                 const { return base::end   (map_dimension()); }
 
    const_iterator begin_edge()          const { return base::begin (1); }
    const_iterator end_edge()            const { return base::end   (1); }
    const_iterator begin_face()          const { return base::begin (2); }
    const_iterator end_face()            const { return base::end   (2); }
 
protected:
    void compute_bbox();
    void init_neighbour() const;
    template<class U> friend void add_ball_externals (const geo_base_rep<U,M>&, const disarray<index_set,M>&);
 
    friend class geo_rep<T,M>;
 
// data:
// 0) header:
    std::string                              _name;
    size_type                                _version;
    size_type                                _serial_number;
// 1) connectivity:
    std::array<hack_array<geo_element_hack,M>, reference_element::max_variant> _geo_element;
    geo_size                                 _gs; // counters by geo_element dimension: 0,1,2,3
    mutable std::vector<domain_indirect_basic<M> > _domains;
    bool                                     _have_connectivity; // e.g.list of  edges in a 2d triangular mesh
    mutable bool                             _have_neighbour; // inter-element connectivity
// 2) coordinates:
    disarray<node_type, M>                   _node;
    size_type                                _dimension;
    coordinate_type                          _sys_coord;
    node_type                                _xmin; // bounding box
    node_type                                _xmax;
    T                                        _hmin;
    T                                        _hmax;
    basis_basic<T>                           _piola_basis;
    geo_locate<T,M>                          _locator;
    geo_trace_ray_boundary<T,M>              _tracer_ray_boundary;
    geo_nearest<T,M>                         _nearestor;
 
// static data member: loaded geo table
public:
    typedef std::unordered_map<std::string,void*> loaded_map_t;
    static loaded_map_t& loaded_map();
protected:
    static loaded_map_t _loaded_map;
};
template <class T, class M>
inline
void
geo_base_rep<T,M>::neighbour_guard() const
{
  if (_have_neighbour) return;
  _have_neighbour = true;
  init_neighbour();
}
template <class T, class M>
typename geo_base_rep<T,M>::iterator_by_variant
geo_base_rep<T,M>::begin_by_variant (variant_type variant)
{
  return _geo_element [variant].begin(); 
}
template <class T, class M>
typename geo_base_rep<T,M>::const_iterator_by_variant
geo_base_rep<T,M>::begin_by_variant (variant_type variant) const
{
  return _geo_element [variant].begin(); 
}
template <class T, class M>
typename geo_base_rep<T,M>::iterator_by_variant
geo_base_rep<T,M>::end_by_variant (variant_type variant)
{
  return _geo_element [variant].end(); 
}
template <class T, class M>
typename geo_base_rep<T,M>::const_iterator_by_variant
geo_base_rep<T,M>::end_by_variant (variant_type variant) const
{
  return _geo_element [variant].end(); 
}

template <class T, class M>
inline
typename geo_base_rep<T,M>::iterator
geo_base_rep<T,M>::begin (size_type dim)
{
  variant_type variant = reference_element::first_variant_by_dimension(dim);
  iterator_by_variant iter = begin_by_variant (variant);
  iterator res = iterator (dim, variant, iter, *this);
  return res;
}
template <class T, class M>
inline
typename geo_base_rep<T,M>::iterator
geo_base_rep<T,M>::end (size_type dim) 
{
  variant_type variant = reference_element::last_variant_by_dimension(dim);
  iterator_by_variant iter = end_by_variant (variant - 1);
  iterator res = iterator (dim, variant, iter, *this);
  return res;
}
template <class T>
inline
typename geo_abstract_base_rep<T>::const_iterator
geo_abstract_base_rep<T>::begin (size_type dim) const
{
  variant_type variant = reference_element::first_variant_by_dimension(dim);
  const_iterator_by_variant iter = begin_by_variant (variant);
  return const_iterator (dim, variant, iter, *this);
}
template <class T>
inline
typename geo_abstract_base_rep<T>::const_iterator
geo_abstract_base_rep<T>::end (size_type dim) const
{
  variant_type variant = reference_element::last_variant_by_dimension(dim);
  const_iterator_by_variant iter = end_by_variant (variant - 1);
  return const_iterator (dim, variant, iter, *this);
}
// =========================================================================
// =========================================================================
template <class T, class M> class geo_rep {};
template <class T, class M>
void geo_build_by_subdividing (
            geo_rep  <T,M>&            new_omega,
      const geo_basic<T,M>&            old_omega,
      typename geo_rep<T,M>::size_type k);
 
template <class T>
class geo_rep<T,sequential> : public geo_base_rep<T,sequential> {
public:
// typedefs:
 
    typedef geo_base_rep<T,sequential>               base;
    typedef typename base::size_type                 size_type;
    typedef typename base::node_type                 node_type;
    typedef typename base::variant_type              variant_type;
    typedef typename base::reference                 reference;
    typedef typename base::const_reference           const_reference;
    typedef typename base::iterator                  iterator;
    typedef typename base::const_iterator            const_iterator;
    typedef typename base::iterator_by_variant       iterator_by_variant;
    typedef typename base::const_iterator_by_variant const_iterator_by_variant;
    typedef typename base::coordinate_type           coordinate_type;
    typedef typename base::geo_element_map_type      geo_element_map_type;
 
// allocators:
 
    geo_rep();
    geo_rep (const geo_rep<T,sequential>&);
    geo_abstract_rep<T,sequential>* clone() const;
 
    // build from_list (for level set)
    geo_rep (
      const geo_basic<T,sequential>&                      lambda,
      const disarray<point_basic<T>,sequential>&          node_list,
      const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,sequential>,
                         reference_element::max_variant>& elt_list);
 
    void build_from_domain (
        const domain_indirect_rep<sequential>&              indirect,
        const geo_abstract_rep<T,sequential>&               omega,
              std::map<size_type,size_type>&                bgd_ie2dom_ie,
              std::map<size_type,size_type>&                dis_bgd_ie2dis_dom_ie);
 
// abstract accessors redefined:
 
    void locate (
                const disarray<point_basic<T>, sequential>&    x,
                      disarray<size_type, sequential>&         dis_ie,
                bool do_check = false) const;
    void trace_ray_boundary (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y,
                bool  do_check = false) const;
    void trace_move (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y) const;
    void nearest (
                const disarray<point_basic<T>,sequential>&     x,
                      disarray<point_basic<T>,sequential>&     x_nearest,
                      disarray<size_type, sequential>&         dis_ie) const;
 
// herited accessors:
 
    size_type map_dimension () const { return base::map_dimension(); }
    bool           is_broken() const { return base::is_broken(); }
    const distributor& geo_element_ownership(size_type dim) const { return base::geo_element_ownership (dim); }
    const_reference get_geo_element (size_type dim, size_type ige) const { return base::get_geo_element (dim, ige); }
    reference       get_geo_element (size_type dim, size_type ige)       { return base::get_geo_element (dim, ige); }
    iterator       begin (size_type dim)       { return base::begin(dim); }
    iterator       end   (size_type dim)       { return base::end  (dim); }
    const_iterator begin (size_type dim) const { return base::begin(dim); }
    const_iterator end   (size_type dim) const { return base::end  (dim); }
    const geo_element_map_type& get_external_geo_element_map (size_type variant) const;
 
    const domain_indirect_basic<sequential>& get_domain_indirect (size_type idom) const { return base::get_domain_indirect (idom); }
    const domain_indirect_basic<sequential>& get_domain_indirect (const std::string& name) const { return base::get_domain_indirect (name); }
 
// deduced accessors:
 
    const distributor& vertex_ownership()        const { return geo_element_ownership(0); }
    const_reference operator[] (size_type ie) const { return get_geo_element (map_dimension(), ie); }
    reference       operator[] (size_type ie)       { return get_geo_element (map_dimension(), ie); }
 
// i/o:
 
    idiststream& get         (idiststream&);
    odiststream& put_geo     (odiststream&) const;
    odiststream& put         (odiststream& ops) const { return put_geo(ops); }
 
    void dump (std::string name) const;
    void load (std::string name, const communicator& = communicator());
    bool check(bool verbose) const;
 
// modifier:
 
    void reset_order (size_type order);
    void build_by_subdividing (const geo_basic<T,sequential>& omega, size_type k);
    void build_from_data (
        const geo_header&                                               hdr,
        const disarray<node_type, sequential>&                          node,
              std::array<disarray<geo_element_auto<>,sequential>, reference_element::max_variant>&
                                                                        tmp_geo_element,
        bool                                                            do_upgrade);
 
// internal:
protected:
    idiststream& get_standard (idiststream&, const geo_header&);
    idiststream& get_upgrade  (idiststream&, const geo_header&);
    void build_connectivity (
        std::array<disarray<geo_element_auto<>,sequential>, reference_element::max_variant>& tmp_geo_element);
 
    void build_connectivity_sides (
                size_type side_dim, 
                std::array<disarray<geo_element_auto<>,sequential>, reference_element::max_variant>& tmp_geo_element);
    void set_element_side_index (size_type side_dim);
    void domain_set_side_part1 (
                const domain_indirect_rep<sequential>&    indirect,
                const geo_abstract_rep<T,sequential>&     bgd_omega,
                size_type                                 sid_dim,
                disarray<size_type,sequential>&           bgd_isid2dom_dis_isid,
                disarray<size_type,sequential>&           dom_isid2bgd_isid,
                disarray<size_type,sequential>&           dom_isid2dom_ios_dis_isid,
                size_type                                 size_by_variant [reference_element::max_variant]);
    void domain_set_side_part2 (
                const domain_indirect_rep<sequential>&    indirect,
                const geo_abstract_rep<T,sequential>&     bgd_omega,
                disarray<size_type,sequential>&           bgd_iv2dom_dis_iv,
                size_type                                 sid_dim,
                disarray<size_type,sequential>&           bgd_isid2dom_dis_isid,
                disarray<size_type,sequential>&           dom_isid2bgd_isid,
                disarray<size_type,sequential>&           dom_isid2dom_ios_dis_isid,
                size_type                                 size_by_variant [reference_element::max_variant]);
    void build_external_entities () {} // for distributed compat
 
// friends:
 
  friend void geo_build_by_subdividing<> (
            geo_rep<T,sequential>&              new_omega,
      const geo_basic<T,sequential>&            old_omega,
      typename geo_rep<T,sequential>::size_type k);
};
#ifdef _RHEOLEF_HAVE_MPI
// =========================================================================
// =========================================================================
template <class T>
class geo_rep<T,distributed> : public geo_base_rep<T,distributed> {
public:
// typedefs:
 
    typedef geo_base_rep<T,distributed>              base;
    typedef typename base::size_type                 size_type;
    typedef typename base::node_type                 node_type;
    typedef typename base::variant_type              variant_type;
    typedef typename base::node_map_type             node_map_type;
    typedef typename base::reference                 reference;
    typedef typename base::const_reference           const_reference;
    typedef typename base::iterator                  iterator;
    typedef typename base::const_iterator            const_iterator;
    typedef typename base::iterator_by_variant       iterator_by_variant;
    typedef typename base::const_iterator_by_variant const_iterator_by_variant;
    typedef typename base::coordinate_type           coordinate_type;
    typedef typename base::geo_element_map_type      geo_element_map_type;
 
// allocators:
 
    geo_rep ();
    geo_rep (const geo_rep<T,distributed>&);
    geo_abstract_rep<T,distributed>* clone() const;
 
    void build_from_domain (
        const domain_indirect_rep<distributed>&             indirect,
        const geo_abstract_rep<T,distributed>&              omega,
              std::map<size_type,size_type>&                bgd_ie2dom_ie,
              std::map<size_type,size_type>&                dis_bgd_ie2dis_dom_ie);
 
    // build from_list (for level set)
    geo_rep (
      const geo_basic<T,distributed>&                     lambda,
      const disarray<point_basic<T>,distributed>&         node_list,
      const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,distributed>,
                         reference_element::max_variant>& elt_list);
 
// abstract accessors defined:
 
    distributor geo_element_ios_ownership (size_type dim) const;
 
    size_type             ige2ios_dis_ige (size_type dim, size_type ige) const;
    size_type         dis_ige2ios_dis_ige (size_type dim, size_type dis_ige) const;
    size_type             ios_ige2dis_ige (size_type dim, size_type ios_ige) const;
    const geo_size&  ios_sizes() const { return _ios_gs; }
    void locate (
                const disarray<point_basic<T>,distributed>&     x,
                      disarray<size_type, distributed>&         dis_ie,
                      bool do_check = false) const;
    void trace_ray_boundary (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y,
                bool  do_check = false) const;
    void trace_move (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y) const;
    void nearest (
                const disarray<point_basic<T>,distributed>&     x,
                      disarray<point_basic<T>,distributed>&     x_nearest,
                      disarray<size_type, distributed>&         dis_ie) const;
 
// herited accessors:
 
    size_type map_dimension () const { return base::map_dimension(); }
    bool           is_broken() const { return base::is_broken(); }
    size_type size (size_type dim) const { return base::size(dim); }
    const distributor& geo_element_ownership(size_type dim) const { return base::geo_element_ownership (dim); }
    const_reference get_geo_element (size_type dim, size_type ige) const { return base::get_geo_element (dim, ige); }
    reference       get_geo_element (size_type dim, size_type ige)       { return base::get_geo_element (dim, ige); }
 
    const_iterator begin (size_type dim) const { return base::begin(dim); }
    const_iterator end   (size_type dim) const { return base::end  (dim); }
    iterator       begin (size_type dim)       { return base::begin(dim); }
    iterator       end   (size_type dim)       { return base::end  (dim); }
    const geo_element_map_type& get_external_geo_element_map (size_type variant) const 
        { return base::_geo_element[variant].get_dis_map_entries(); }
 
    const domain_indirect_basic<distributed>& get_domain_indirect (size_type idom) const { return base::get_domain_indirect (idom); }
    const domain_indirect_basic<distributed>& get_domain_indirect (const std::string& name) const { return base::get_domain_indirect (name); }
 
// deduced accessors:
 
    size_type          size ()                   const { return size (map_dimension()); }
    const distributor& vertex_ownership()        const { return geo_element_ownership(0); }
    const_reference operator[] (size_type ie) const { return get_geo_element (map_dimension(), ie); }
#ifdef TODO
          reference operator[] (size_type ie)       { return get_geo_element (map_dimension(), ie); }
#endif // TODO
 
// modifier:
 
    void reset_order (size_type order);
    void build_by_subdividing (const geo_basic<T,distributed>& omega, size_type k);
 
// i/o:
 
    idiststream& get (idiststream&);
    odiststream& put (odiststream&) const;
    void dump (std::string name) const;
    void load (std::string name, const communicator& comm);
    bool check(bool verbose) const;
 
// utilities:
 
    void set_ios_permutation (disarray<size_type,distributed>& idof2ios_dis_idof) const;
 
protected:
// internal:
    void build_external_entities ();
    void set_element_side_index (size_type side_dim);
    void 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]);
    void 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]);
 
    void node_renumbering (const distributor& ios_node_ownership);
 
// data:
    disarray<size_type>          _inod2ios_dis_inod;  // permutation for node
    disarray<size_type>          _ios_inod2dis_inod;  // reverse permutation for node
    std::array<disarray<size_type>, 4> _ios_ige2dis_ige; // reverse permutation for geo_element[dim]
    geo_size                     _ios_gs;
    std::array<disarray<size_type,distributed>,reference_element::max_variant> _igev2ios_dis_igev;
    std::array<disarray<size_type,distributed>,reference_element::max_variant> _ios_igev2dis_igev;
 
public:
    // used by space_constritution for ios numbering
    const std::array<disarray<size_type,distributed>,reference_element::max_variant>&
    get_igev2ios_dis_igev() const { return _igev2ios_dis_igev; }
 
// friends:
 
  friend void geo_build_by_subdividing<> (
            geo_rep<T,distributed>&              new_omega,
      const geo_basic<T,distributed>&            old_omega,
      typename geo_rep<T,distributed>::size_type k);
};
#endif // _RHEOLEF_HAVE_MPI
 
// =========================================================================
// =========================================================================
template <class T, class M = rheo_default_memory_model>
class geo_basic {
public:
    typedef M memory_type;
};
// =========================================================================
// =========================================================================
// handler for complex geo names as "square[boundary]"
template<class T, class M> geo_basic<T,M> geo_load (const std::string& name);
 
// guards for omega.boundary(), omega.internal_sides() and omega.sides()
template<class T, class M> void boundary_guard       (const geo_basic<T,M>&);
template<class T, class M> void internal_sides_guard (const geo_basic<T,M>&);
template<class T, class M> void sides_guard          (const geo_basic<T,M>&);
 
// [verbatim_geo_basic]
template <class T>
class geo_basic<T,sequential> : public smart_pointer_clone<geo_abstract_rep<T,sequential> > {
public:
 
// typedefs:
 
    typedef sequential                              memory_type;
    typedef geo_abstract_rep<T,sequential>          rep;
    typedef geo_rep<T,sequential>                   rep_geo_rep;
    typedef smart_pointer_clone<rep>                base;
    typedef typename rep::size_type                 size_type;
    typedef typename rep::node_type                 node_type;
    typedef typename rep::variant_type              variant_type;
    typedef typename rep::reference                 reference;
    typedef typename rep::const_reference           const_reference;
    typedef typename rep::iterator                  iterator;
    typedef typename rep::const_iterator            const_iterator;
    typedef typename rep::iterator_by_variant       iterator_by_variant;
    typedef typename rep::const_iterator_by_variant const_iterator_by_variant;
    typedef typename rep::coordinate_type           coordinate_type;
    typedef typename rep::geo_element_map_type      geo_element_map_type;
 
// allocators:
 
    geo_basic ();
    geo_basic (std::string name, const communicator& comm = communicator());
    void load (std::string name, const communicator& comm = communicator());
    geo_basic (const domain_indirect_basic<sequential>& dom, const geo_basic<T,sequential>& omega);
 
    // build from_list (for level set)
    geo_basic (
      const geo_basic<T,sequential>&                      lambda,
      const disarray<point_basic<T>,sequential>&          node_list,
      const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,sequential>,
                         reference_element::max_variant>& elt_list)
    : base (new_macro(rep_geo_rep(lambda,node_list,elt_list))) {}
 
// accessors:
 
    std::string                    name() const { return base::data().name(); }
    std::string              familyname() const { return base::data().familyname(); }
    size_type                 dimension() const { return base::data().dimension(); }
    size_type             map_dimension() const { return base::data().map_dimension(); }
    bool                      is_broken() const { return base::data().is_broken(); }
    size_type             serial_number() const { return base::data().serial_number(); }
    size_type                   variant() const { return base::data().variant(); }
    coordinate_type   coordinate_system() const { return base::data().coordinate_system(); }
    std::string  coordinate_system_name() const { return space_constant::coordinate_system_name(coordinate_system()); }
    const basis_basic<T>& get_piola_basis() const { return base::data().get_piola_basis(); }
    size_type                     order() const { return base::data().get_piola_basis().degree(); }
    const node_type&               xmin() const { return base::data().xmin(); }
    const node_type&               xmax() const { return base::data().xmax(); }
    const T&                       hmin() const { return base::data().hmin(); }
    const T&                       hmax() const { return base::data().hmax(); }
    const distributor& geo_element_ownership(size_type dim) const { return base::data().geo_element_ownership(dim); }
    const geo_size&      sizes()             const { return base::data().sizes(); }
    const geo_size&  ios_sizes()             const { return base::data().ios_sizes(); }
    const_reference get_geo_element (size_type dim, size_type ige) const { return base::data().get_geo_element (dim, ige); }
    const_reference dis_get_geo_element (size_type dim, size_type dis_ige) const
                { return get_geo_element (dim, dis_ige); }
    const geo_element& bgd2dom_geo_element (const geo_element& bgd_K) const { return base::data().bgd2dom_geo_element (bgd_K); }
    const geo_element& dom2bgd_geo_element (const geo_element& dom_K) const { return base::data().dom2bgd_geo_element (dom_K); }
    size_type neighbour (size_type ie, size_type loc_isid) const {
                          return base::data().neighbour (ie, loc_isid); }
    void neighbour_guard() const { base::data().neighbour_guard(); }
    size_type        n_node()   const { return base::data().n_node(); }
    const node_type&     node(size_type     inod) const { return base::data().node(inod); }
    const node_type& dis_node(size_type dis_inod) const { return base::data().dis_node(dis_inod); }
    void dis_inod (const geo_element& K, std::vector<size_type>& dis_inod) const {
                return base::data().dis_inod(K,dis_inod); }
    const disarray<node_type,sequential>& get_nodes() const { return base::data().get_nodes(); }
    size_type dis_inod2dis_iv (size_type dis_inod) const { return base::data().dis_inod2dis_iv(dis_inod); }
 
    size_type n_domain_indirect () const { return base::data().n_domain_indirect (); }
    bool have_domain_indirect (const std::string& name) const { return base::data().have_domain_indirect (name); }
    const domain_indirect_basic<sequential>& get_domain_indirect (size_type i) const {
          return base::data().get_domain_indirect (i); }
    const domain_indirect_basic<sequential>& get_domain_indirect (const std::string& name) const {
          return base::data().get_domain_indirect (name); }
    void  insert_domain_indirect (const domain_indirect_basic<sequential>& dom) const {
          base::data().insert_domain_indirect (dom); }
 
    size_type n_domain () const { return base::data().n_domain_indirect (); }
    geo_basic<T,sequential> get_domain (size_type i) const;
    geo_basic<T,sequential> operator[] (const std::string& name) const;
    geo_basic<T,sequential> boundary() const;
    geo_basic<T,sequential> internal_sides() const;
    geo_basic<T,sequential> sides() const;
 
// modifiers:
 
    void set_name (std::string name);
    void set_dimension (size_type dim);
    void set_serial_number (size_type i);
    void reset_order (size_type order);
    void set_coordinate_system (coordinate_type sys_coord);
    void set_coordinate_system (std::string sys_coord_name) { set_coordinate_system (space_constant::coordinate_system(sys_coord_name)); }
    void set_nodes (const disarray<node_type,sequential>& x);
 
// extended accessors:
 
    const communicator& comm()        const { return geo_element_ownership (0).comm(); }
    size_type     size(size_type dim) const { return base::data().geo_element_ownership(dim).size(); }
    size_type dis_size(size_type dim) const { return base::data().geo_element_ownership(dim).dis_size(); }
    size_type     size()              const { return size     (map_dimension()); }
    size_type dis_size()              const { return dis_size (map_dimension()); }
    size_type     n_vertex()          const { return size     (0); }
    size_type dis_n_vertex()          const { return dis_size (0); }
    const_reference operator[] (size_type ie) const { return get_geo_element (map_dimension(), ie); }
    const_iterator begin (size_type dim) const { return base::data().begin(dim); }
    const_iterator end   (size_type dim) const { return base::data().end  (dim); }
    const_iterator begin ()              const { return begin(map_dimension()); }
    const_iterator end   ()              const { return end  (map_dimension()); }
 
// comparator:
 
    bool operator== (const geo_basic<T,sequential>& omega2) const { return base::data().operator== (omega2.data()); }
 
// i/o:
 
    void save (std::string filename = "") const;
// [verbatim_geo_basic]
 
// internals:
    geo_basic (details::zero_dimension, const communicator& comm = communicator());
    bool check (bool verbose = true) const { return base::data().check(verbose); }
    idiststream& get (idiststream& ips);
    odiststream& put (odiststream& ops) const;
 
    void build_by_subdividing (const geo_basic<T,sequential>& omega, size_type k);
    void build_from_data (
        const geo_header&                                               hdr,
        const disarray<node_type, sequential>&                          node,
              std::array<disarray<geo_element_auto<>,sequential>, reference_element::max_variant>&
                                                                        tmp_geo_element,
        bool                                                            do_upgrade);
 
    // locators:
    size_type seq_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const
                { return base::data().seq_locate (x, dis_ie_guest); }
    size_type dis_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const
                { return base::data().dis_locate (x, dis_ie_guest); }
    void locate (       
                const disarray<point_basic<T>, sequential>& x,
                disarray<size_type, sequential>& dis_ie) const
                { return base::data().locate (x, dis_ie); }
    size_type seq_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
                                        { return base::data().seq_trace_move (x,v,y); }
    size_type dis_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
                                        { return base::data().dis_trace_move (x,v,y); }
    void trace_ray_boundary (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y) const
                                        { return base::data().trace_ray_boundary (x,v,dis_ie,y); }
    void trace_move (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y) const
                                        { return base::data().trace_move (x,v,dis_ie,y); }
    size_type seq_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const
                                        { return base::data().seq_nearest (x, x_nearest); }
    size_type dis_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const
                                        { return base::data().dis_nearest (x, x_nearest); }
    void nearest (
                const disarray<point_basic<T>,sequential>&     x,
                      disarray<point_basic<T>,sequential>&     x_nearest,
                      disarray<size_type, sequential>&         dis_ie) const
                                        { base::data().nearest (x, x_nearest, dis_ie); }
 
    const_iterator_by_variant begin_by_variant (variant_type variant) const 
        { return base::data().begin_by_variant (variant); }
    const_iterator_by_variant   end_by_variant (variant_type variant) const
        { return base::data().  end_by_variant (variant); }
    const geo_element_map_type& get_external_geo_element_map (size_type variant) const 
        { return base::data().get_external_geo_element_map(variant); }
 
    const geo_basic<T,sequential>& get_background_geo() const; // code in geo_domain.h
          geo_basic<T,sequential>  get_background_domain() const;
 
    // for compatibility with distributed interface: 
    size_type ige2ios_dis_ige (size_type dim, size_type ige) const { return ige; }
    size_type dis_ige2ios_dis_ige (size_type dim, size_type dis_ige) const { return dis_ige; }
    size_type ios_ige2dis_ige (size_type dim, size_type ios_ige) const { return ios_ige; }
    void set_ios_permutation (disarray<size_type,sequential>& idof2ios_dis_idof) const {}
 
// [verbatim_geo_basic_cont]
};
template <class T, class M>
idiststream& operator>> (idiststream& ips, geo_basic<T,M>& omega);
 
template <class T, class M>
odiststream& operator<< (odiststream& ops, const geo_basic<T,M>& omega);
// [verbatim_geo_basic_cont]
 
template <class T>
inline
geo_basic<T,sequential>::geo_basic()
  : base (new_macro((geo_rep<T,sequential>)))
{
}
template <class T>
inline
geo_basic<T,sequential>::geo_basic (std::string name, const communicator& comm)
  : base (0)
{
    base::operator= (geo_load<T,sequential>(name));
}
template <class T>
inline
void
geo_basic<T,sequential>::load (std::string name, const communicator& comm)
{
    base::operator= (geo_load<T,sequential>(name));
}
template <class T>
inline
geo_basic<T,sequential>
geo_basic<T,sequential>::boundary() const
{
   boundary_guard (*this);
   return operator[] ("boundary");
}
template <class T>
inline
geo_basic<T,sequential>
geo_basic<T,sequential>::internal_sides() const
{
   internal_sides_guard (*this);
   return operator[] ("internal_sides");
}
template <class T>
inline
geo_basic<T,sequential>
geo_basic<T,sequential>::sides() const
{
   sides_guard (*this);
   return operator[] ("sides");
}
#ifdef _RHEOLEF_HAVE_MPI
// =========================================================================
// =========================================================================
template <class T>
class geo_basic<T,distributed> : public smart_pointer_clone<geo_abstract_rep<T,distributed> > {
public:
 
// typedefs:
 
    typedef distributed                             memory_type;
    typedef geo_abstract_rep<T,distributed>         rep;
    typedef geo_rep<T,distributed>                  rep_geo_rep;
    typedef smart_pointer_clone<rep>                base;
    typedef typename rep::size_type                 size_type;
    typedef typename rep::node_type                 node_type;
    typedef typename rep::variant_type              variant_type;
    typedef typename rep::node_map_type             node_map_type;
    typedef typename rep::reference                 reference;
    typedef typename rep::const_reference           const_reference;
    typedef typename rep::iterator                  iterator;
    typedef typename rep::const_iterator            const_iterator;
    typedef typename rep::iterator_by_variant       iterator_by_variant;
    typedef typename rep::const_iterator_by_variant const_iterator_by_variant;
    typedef typename rep::coordinate_type           coordinate_type;
    typedef typename rep::geo_element_map_type      geo_element_map_type;
 
// allocators:
 
    geo_basic ();
    geo_basic (std::string name, const communicator& comm = communicator());
    void load (std::string name, const communicator& comm = communicator());
    geo_basic (const domain_indirect_basic<distributed>& dom, const geo_basic<T,distributed>& omega);
    geo_basic (details::zero_dimension, const communicator& comm = communicator());
 
    // build from_list (for level set)
    geo_basic (
      const geo_basic<T,distributed>&                     lambda,
      const disarray<point_basic<T>,distributed>&                   node_list,
      const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,distributed>,
                         reference_element::max_variant>& elt_list)
    : base (new_macro(rep_geo_rep(lambda,node_list,elt_list))) {}
 
// accessors:
 
    std::string                    name() const { return base::data().name(); }
    std::string              familyname() const { return base::data().familyname(); }
    size_type                 dimension() const { return base::data().dimension(); }
    size_type             map_dimension() const { return base::data().map_dimension(); }
    bool                      is_broken() const { return base::data().is_broken(); }
    size_type             serial_number() const { return base::data().serial_number(); }
    size_type                   variant() const { return base::data().variant(); }
    coordinate_type   coordinate_system() const { return base::data().coordinate_system(); }
    std::string  coordinate_system_name() const { return space_constant::coordinate_system_name(coordinate_system()); }
    const basis_basic<T>& get_piola_basis() const { return base::data().get_piola_basis(); }
    size_type                     order() const { return base::data().get_piola_basis().degree(); }
    const node_type&               xmin() const { return base::data().xmin(); }
    const node_type&               xmax() const { return base::data().xmax(); }
    const T&                       hmin() const { return base::data().hmin(); }
    const T&                       hmax() const { return base::data().hmax(); }
    const distributor& geo_element_ownership(size_type dim) const
                        { return base::data().geo_element_ownership (dim); }
    const geo_size&      sizes()             const { return base::data().sizes(); }
    const geo_size&  ios_sizes()             const { return base::data().ios_sizes(); }
    const_reference get_geo_element (size_type dim, size_type ige) const
                { return base::data().get_geo_element (dim, ige); }
    const_reference dis_get_geo_element (size_type dim, size_type dis_ige) const
                { return base::data().dis_get_geo_element (dim, dis_ige); }
    const geo_element& bgd2dom_geo_element (const geo_element& bgd_K) const
                { return base::data().bgd2dom_geo_element (bgd_K); }
    const geo_element& dom2bgd_geo_element (const geo_element& dom_K) const
                { return base::data().dom2bgd_geo_element (dom_K); }
    size_type neighbour (size_type ie, size_type loc_isid) const {
                          return base::data().neighbour (ie, loc_isid); }
    void neighbour_guard() const { base::data().neighbour_guard(); }
    distributor geo_element_ios_ownership (size_type dim) const {
        return base::data().geo_element_ios_ownership (dim); }
    size_type ige2ios_dis_ige (size_type dim, size_type ige) const {
        return base::data().ige2ios_dis_ige (dim,ige); }
    size_type dis_ige2ios_dis_ige (size_type dim, size_type dis_ige) const {
        return base::data().dis_ige2ios_dis_ige (dim,dis_ige); }
    size_type ios_ige2dis_ige (size_type dim, size_type ios_ige) const {
        return base::data().ios_ige2dis_ige (dim, ios_ige); }
    size_type        n_node() const { return base::data().n_node(); }
    const node_type&     node(size_type     inod) const { return base::data().node(inod); }
    const node_type& dis_node(size_type dis_inod) const { return base::data().dis_node(dis_inod); }
    void dis_inod (const geo_element& K, std::vector<size_type>& dis_inod) const {
                return base::data().dis_inod(K,dis_inod); }
    const disarray<node_type,distributed>& get_nodes() const { return base::data().get_nodes(); }
 
    size_type n_domain_indirect () const { return base::data().n_domain_indirect (); }
    bool have_domain_indirect (const std::string& name) const { return base::data().have_domain_indirect (name); }
    const domain_indirect_basic<distributed>& get_domain_indirect (size_type i) const {
          return base::data().get_domain_indirect (i); }
    const domain_indirect_basic<distributed>& get_domain_indirect (const std::string& name) const {
          return base::data().get_domain_indirect (name); }
    void  insert_domain_indirect (const domain_indirect_basic<distributed>& dom) const {
          base::data().insert_domain_indirect (dom); }
 
    size_type n_domain () const { return base::data().n_domain_indirect (); }
    geo_basic<T,distributed> get_domain (size_type i) const;
    geo_basic<T,distributed> operator[] (const std::string& name) const;
    geo_basic<T,distributed> boundary() const;
    geo_basic<T,distributed> internal_sides() const;
    geo_basic<T,distributed> sides() const;
 
    size_type seq_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const
                { return base::data().seq_locate (x, dis_ie_guest); }
    size_type dis_locate (
                const point_basic<T>& x,
                size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const
                { return base::data().dis_locate (x, dis_ie_guest); }
    void locate (const disarray<point_basic<T>, distributed>& x, disarray<size_type, distributed>& dis_ie) const
                { return base::data().locate (x, dis_ie); }
    size_type seq_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
                                        { return base::data().seq_trace_move (x,v,y); }
    size_type dis_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
                                        { return base::data().dis_trace_move (x,v,y); }
    void trace_ray_boundary (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y) const
                                        { return base::data().trace_ray_boundary (x,v,dis_ie,y); }
    void trace_move (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie,
                      disarray<point_basic<T>,distributed>&     y) const
                                        { return base::data().trace_move (x,v,dis_ie,y); }
    size_type seq_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const
                                        { return base::data().seq_nearest (x, x_nearest); }
    size_type dis_nearest (
                const point_basic<T>&    x,
                      point_basic<T>&    x_nearest) const
                                        { return base::data().dis_nearest (x, x_nearest); }
    void nearest (
                const disarray<point_basic<T>,distributed>&     x,
                      disarray<point_basic<T>,distributed>&     x_nearest,
                      disarray<size_type, distributed>&         dis_ie) const
                                        { base::data().nearest (x, x_nearest, dis_ie); }
// modifiers:
 
    void set_nodes (const disarray<node_type,distributed>& x);
    void reset_order (size_type order);
    size_type dis_inod2dis_iv (size_type dis_inod) const { return base::data().dis_inod2dis_iv(dis_inod); }
    void set_coordinate_system (coordinate_type sys_coord);
    void set_coordinate_system (std::string sys_coord_name) { set_coordinate_system (space_constant::coordinate_system(sys_coord_name)); }
    void set_dimension (size_type dim);
    void set_serial_number (size_type i);
    void set_name (std::string name);
    void build_by_subdividing (const geo_basic<T,distributed>& omega, size_type k);
 
// extended accessors:
 
    size_type     size(size_type dim) const { return base::data().geo_element_ownership(dim).size(); }
    size_type dis_size(size_type dim) const { return base::data().geo_element_ownership(dim).dis_size(); }
    const communicator& comm()        const { return geo_element_ownership (0).comm(); }
    size_type     size()              const { return size     (map_dimension()); }
    size_type dis_size()              const { return dis_size (map_dimension()); }
    size_type     n_vertex()          const { return size     (0); }
    size_type dis_n_vertex()          const { return dis_size (0); }
    const_reference operator[] (size_type ie) const     
                { return get_geo_element (map_dimension(), ie); }
 
    const_iterator begin (size_type dim) const { return base::data().begin(dim); }
    const_iterator end   (size_type dim) const { return base::data().end  (dim); }
    const_iterator begin ()              const { return begin(map_dimension()); }
    const_iterator end   ()              const { return end  (map_dimension()); }
 
    const_iterator_by_variant begin_by_variant (variant_type variant) const 
        { return base::data().begin_by_variant (variant); }
    const_iterator_by_variant   end_by_variant (variant_type variant) const
        { return base::data().  end_by_variant (variant); }
    const geo_element_map_type& get_external_geo_element_map (size_type variant) const 
        { return base::data().get_external_geo_element_map(variant); }
 
    const geo_basic<T,distributed>& get_background_geo() const; // code in geo_domain.h
          geo_basic<T,distributed>  get_background_domain() const;
 
// comparator:
 
    bool operator== (const geo_basic<T,distributed>& omega2) const { return base::data().operator== (omega2.data()); }
 
// i/o:
 
    odiststream& put (odiststream& ops) const { return base::data().put (ops); }
    idiststream& get (idiststream& ips);
    void save (std::string filename = "") const;
    bool check (bool verbose = true) const { return base::data().check(verbose); }
 
// utilities:
 
    void set_ios_permutation (disarray<size_type,distributed>& idof2ios_dis_idof) const
     { base::data().set_ios_permutation (idof2ios_dis_idof); }
 
    // used by space_constritution for ios numbering
    const std::array<disarray<size_type,distributed>,reference_element::max_variant>&
    get_igev2ios_dis_igev() const { return base::data().get_igev2ios_dis_igev(); }
};
#endif // _RHEOLEF_HAVE_MPI
 
// [verbatim_geo]
typedef geo_basic<Float,rheo_default_memory_model> geo;
// [verbatim_geo]
 
// ==============================================================================
// inlined: geo<T,distributed>
// ==============================================================================
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
inline
geo_basic<T,distributed>::geo_basic()
  : base (new_macro((geo_rep<T,distributed>)))
{
}
template <class T>
inline
geo_basic<T,distributed>::geo_basic (std::string name, const communicator& comm)
  : base (0)
{
    base::operator= (geo_load<T,distributed>(name));
}
template <class T>
inline
void
geo_basic<T,distributed>::load (std::string name, const communicator& comm)
{
    base::operator= (geo_load<T,distributed>(name));
}
template <class T>
inline
idiststream&
geo_basic<T,distributed>::get (idiststream& ips)
{
    // allocate a new geo_rep object (TODO: do a dynamic_cast ?)
    geo_rep<T,distributed>* ptr = new_macro((geo_rep<T,distributed>));
    ptr->get (ips);
    base::operator= (ptr);
    return ips;
}
template <class T>
inline
geo_basic<T,distributed>
geo_basic<T,distributed>::boundary() const
{
   boundary_guard (*this);
   return operator[] ("boundary");
}
template <class T>
inline
geo_basic<T,distributed>
geo_basic<T,distributed>::internal_sides() const
{
   internal_sides_guard (*this);
   return operator[] ("internal_sides");
}
template <class T>
inline
geo_basic<T,distributed>
geo_basic<T,distributed>::sides() const
{
   sides_guard (*this);
   return operator[] ("sides");
}
#endif // _RHEOLEF_HAVE_MPI
 
// ==============================================================================
// inlined: geo<T,M>
// ==============================================================================
template <class T, class M>
inline 
idiststream&
operator>> (idiststream& ips, geo_basic<T,M>& omega)
{
    return omega.get (ips);
}
template <class T, class M>
inline 
odiststream&
operator<< (odiststream& ops, const geo_basic<T,M>& omega)
{
    return omega.put (ops);
}
 
} // namespace rheolef
#endif // _RHEOLEF_GEO_H
 
#ifdef _RHEOLEF_GEO_DEMO_TST_CC
// demo for the geo class documentation
#include "rheolef/environment.h"
using namespace rheolef;
using namespace std;
int main(int argc, char**argv) {
  environment rheolef (argc, argv);
  geo omega (argv[1]);
  // [verbatim_geo_demo_tst]
  cout << omega.size() << " " << omega.n_node() << endl;
  for (size_t i = 0, n = omega.size(); i < n; ++i) {
    const geo_element& K = omega[i];
    cout << K.name();
    for (size_t j = 0, m = K.size(); j < m; ++j)
      cout << " " << K[j];
    cout << endl;
  }
  for (size_t jv = 0, nv = omega.n_node(); jv < nv; ++jv)
    cout << omega.node(jv) << endl;
  // [verbatim_geo_demo_tst]
}
#endif // _RHEOLEF_GEO_DEMO_TST_CC