geo_element_v4.h

Go to the documentation of this file.

#ifndef _RHEOLEF_GEO_ELEMENT_V4_H
#define _RHEOLEF_GEO_ELEMENT_V4_H
//      => la memoire est contigue globalement, l'acces plus rapide avec
//    => petit gaspillage de memoire
#include "rheolef/reference_element.h"
#include "rheolef/geo_element_indirect.h"
#include "rheolef/heap_allocator.h"
#include "rheolef/array.h"
#include <boost/serialization/serialization.hpp>
#include <boost/serialization/base_object.hpp>

namespace rheolef {

template <class A> class geo_element_auto;

class geo_element {
public:


  enum {
        _variant_offset    = 0, // i.e. type, as triangle(t) or tetra(T), etc
        _order_offset      = 1, // i.e. k, when Pk curved element
        _dis_ie_offset     = 2, // internal numbering, depend upon partitionand nproc
        _ios_dis_ie_offset = 3, // i/o numbering, independent of parition and nproc
        _master_offset     = 4, // (d-1)-side has one or two master d-element that contains it
        _last_offset       = 6  // here starts node indexes, face indexes, etc
  };
  // => waste a lot of place
  typedef reference_element::size_type    size_type;
  typedef reference_element::variant_type variant_type;
  typedef size_type*                         iterator;
  typedef const size_type*                   const_iterator;
  typedef size_type                          raw_type;

  typedef geo_element                                  generic_type;
  typedef geo_element_auto<heap_allocator<size_type> > automatic_type;

  typedef geo_element_indirect::orientation_type orientation_type; // for sign (+1,-1)
  typedef geo_element_indirect::shift_type       shift_type;       // for 0..3 face shift
  struct  parameter_type {
    variant_type variant;
    size_type    order;
    parameter_type (variant_type v = reference_element::max_variant, size_type o = 0)
      : variant(v), order(o) {} 
  };


  geo_element& operator= (const geo_element& K)
  {
    reset (K.variant(), K.order()); // resize auto, nothing for hack
    std::copy (K._data_begin(), K._data_begin() + _data_size(), _data_begin());
    reset (K.variant(), K.order()); // reset order=1 for hack, resize nothing for auto
    return *this;
  }
  virtual ~geo_element() {}
  virtual void reset (variant_type variant, size_type order) = 0;


  operator reference_element () const { return reference_element(variant()); }


  variant_type variant()    const { return variant_type( *(_data_begin() +  _variant_offset)); }
  size_type order()         const { return *(_data_begin() + _order_offset); }
  size_type dis_ie()        const { return *(_data_begin() + _dis_ie_offset); }
  size_type ios_dis_ie()    const { return *(_data_begin() + _ios_dis_ie_offset); }
  size_type master (bool i) const { return *(_data_begin() + _master_offset + i); }

  size_type dimension()  const { return reference_element::dimension (variant()); }
  size_type size()       const { return reference_element::n_vertex (variant()); }
  char      name()       const { return reference_element::name     (variant()); }
  size_type n_node()     const { return reference_element::n_node (variant(), order()); }

  void set_dis_ie         (size_type dis_ie) { *(_data_begin() + _dis_ie_offset)     = dis_ie; }
  void set_ios_dis_ie (size_type ios_dis_ie) { *(_data_begin() + _ios_dis_ie_offset) = ios_dis_ie; }
  void set_master (bool i, size_type dis_ie) const {
    const_iterator p = _data_begin() + _master_offset + i; // mutable member fct
    *(const_cast<iterator>(p)) = dis_ie;
  }

  iterator       begin()       { return _data_begin() +  _node_offset (variant(), order()); }
  const_iterator begin() const { return _data_begin() +  _node_offset (variant(), order()); }
  iterator       end()         { return begin() + size(); }
  const_iterator end()   const { return begin() + size(); }
  size_type& operator[] (size_type loc_inod)       { return *(begin() + loc_inod); }
  size_type  operator[] (size_type loc_inod) const { return *(begin() + loc_inod); }
  size_type& node       (size_type loc_inod)       { return operator[] (loc_inod); }
  size_type  node       (size_type loc_inod) const { return operator[] (loc_inod); }

  iterator       begin(size_type node_subgeo_dim)       { return begin() + first_inod (node_subgeo_dim); }
  const_iterator begin(size_type node_subgeo_dim) const { return begin() + first_inod (node_subgeo_dim); }
  iterator       end  (size_type node_subgeo_dim)       { return begin() +  last_inod (node_subgeo_dim); }
  const_iterator end  (size_type node_subgeo_dim) const { return begin() +  last_inod (node_subgeo_dim); }

  const geo_element_indirect&  edge_indirect (size_type i) const { 
    const_iterator p = _data_begin() +  _edge_offset (variant(), order()) + i;
    return *(reinterpret_cast<const geo_element_indirect*>(p));
  }
  const geo_element_indirect&  face_indirect (size_type i) const { 
    const_iterator p = _data_begin() +  _face_offset (variant(), order()) + i;
    return *(reinterpret_cast<const geo_element_indirect*>(p));
  }
  geo_element_indirect&  edge_indirect (size_type i)       {
    iterator p = _data_begin() +  _edge_offset (variant(), order()) + i;
    return *(reinterpret_cast<geo_element_indirect*>(p));
  }
  geo_element_indirect&  face_indirect (size_type i)       {
    iterator p = _data_begin() +  _face_offset (variant(), order()) + i;
    return *(reinterpret_cast<geo_element_indirect*>(p));
  }
  size_type edge (size_type i) const { return (dimension() <= 1) ? dis_ie() : edge_indirect(i).index(); }
  size_type face (size_type i) const { return (dimension() <= 2) ? dis_ie() : face_indirect(i).index(); }

  size_type n_subgeo (size_type subgeo_dim) const {
     return reference_element::n_subgeo (variant(), subgeo_dim); }
  size_type subgeo_n_node (size_type subgeo_dim, size_type loc_isid) const {
     return reference_element::subgeo_n_node (variant(), order(), subgeo_dim, loc_isid); }
  size_type subgeo_local_node (size_type subgeo_dim, size_type loc_isid, size_type loc_jsidnod) const {
     return reference_element::subgeo_local_node (variant(), order(), subgeo_dim, loc_isid, loc_jsidnod); }
  size_type subgeo_size (size_type subgeo_dim, size_type loc_isid) const {
     return reference_element::subgeo_n_node (variant(), 1, subgeo_dim, loc_isid); }
  size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const { 
     return reference_element::subgeo_local_node (variant(), 1, subgeo_dim, i_subgeo, i_subgeo_vertex); }
  size_type first_inod (size_type subgeo_dim) const {
     return reference_element::first_inod (variant(), order(), subgeo_dim); }
  size_type  last_inod (size_type subgeo_dim) const {
     return reference_element::last_inod (variant(), order(), subgeo_dim); }

  size_type n_edge () const { return n_subgeo (1); }
  size_type n_face () const { return n_subgeo (2); }


    orientation_type get_side_informations (
                const geo_element& S,
                size_type& loc_isid,
                size_type& shift) const;

    orientation_type get_side_orientation (const geo_element& S) const;

    bool get_orientation_and_shift (const geo_element& S, 
          orientation_type& orient, shift_type& shift) const;
    orientation_type get_edge_orientation (size_type dis_iv0, size_type dis_iv1) const;
    void get_orientation_and_shift (
      size_type dis_iv0, size_type dis_iv1, size_type dis_iv2,
      orientation_type&  orient,
      shift_type&        shift) const;
    void get_orientation_and_shift (
      size_type dis_iv0, size_type dis_iv1, size_type dis_iv2, size_type dis_iv3,
      orientation_type&  orient,
      shift_type&        shift) const;


    template <class T, class M>
    bool contains (const array<point_basic<T>,M>& node, const point_basic<T>& x) const;

// i/o;

  friend std::istream& operator>> (std::istream& is, geo_element& K);
  friend std::ostream& operator<< (std::ostream& os, const geo_element& K);


  static size_type fix_edge_indirect (
      const geo_element& K,
      size_type          loc_iedg,
      size_type          loc_iedg_j, 
      size_type          order);

  static size_type fix_edge_indirect (
      orientation_type  orient,
      size_type         order,
      size_type         loc_iedg_j);

  static void loc_tri_inod2lattice (
      size_type               loc_tri_inod,
      size_type               order,
      point_basic<size_type>& ij_lattice);

  static void loc_qua_inod2lattice (
      size_type               loc_qua_inod,
      size_type               order,
      point_basic<size_type>& ij_lattice);

  static size_type fix_triangle_indirect (
      const geo_element& K,
      size_type          loc_itri,
      size_type          loc_itri_j, 
      size_type          order);

  static size_type fix_triangle_indirect (
      orientation_type  orient,
      shift_type        shift,
      size_type         order,
      size_type         loc_itri_j);

  static size_type fix_quadrangle_indirect (
      const geo_element& K,
      size_type          loc_iqua,
      size_type          loc_iqua_j, 
      size_type          order);

  static size_type fix_quadrangle_indirect (
      orientation_type  orient,
      shift_type        shift,
      size_type         order,
      size_type         loc_iqua_j);


  static size_type _edge_offset (variant_type variant, size_type order) { return _last_offset; }
  static size_type _face_offset (variant_type variant, size_type order) { return _edge_offset(variant,order) + reference_element::n_sub_edge(variant); }
  static size_type _node_offset (variant_type variant, size_type order) { return _face_offset(variant,order) + reference_element::n_sub_face(variant); }
  static size_type _data_size   (variant_type variant, size_type order) { return _node_offset(variant,order) + reference_element::n_node(variant,order); }
  static size_type _data_size   (const parameter_type& p) { return _data_size (p.variant,p.order); }

  size_type _data_size() const { return _data_size (variant(),order()); }

  virtual       iterator _data_begin() = 0;
  virtual       iterator _data_end()   = 0;
  virtual const_iterator _data_begin() const = 0;
  virtual const_iterator _data_end()   const = 0;

  template<class Archive>
  void serialize (Archive& ar, const unsigned int version) {
  }
#ifdef TO_CLEAN
  template<class A = std::allocator<std::vector<int>::size_type> > class geo_element_auto;
#endif // TO_CLEAN
};
template<class A = std::allocator<std::vector<int>::size_type> >
class geo_element_auto : public geo_element {
public:


  typedef A                               allocator_type;
  typedef reference_element::size_type    size_type;
  typedef reference_element::variant_type variant_type;
  typedef geo_element::iterator           iterator;
  typedef geo_element::const_iterator     const_iterator;
  typedef geo_element::parameter_type     parameter_type;
  typedef geo_element                     generic_type;
  typedef geo_element::automatic_type     automatic_type;


  explicit geo_element_auto (const A& alloc = A())
    : _data (_last_offset, std::numeric_limits<size_type>::max(), alloc)
    {
        _data [_variant_offset] = reference_element::max_variant;
        _data [_order_offset]   = 0;
    }
  explicit geo_element_auto (variant_type variant, size_type order = 1, const A& alloc = A())
    : _data (_data_size(variant,order), std::numeric_limits<size_type>::max(), alloc)
    {
        _data [_variant_offset] = variant;
        _data [_order_offset]   = order;
    }
  explicit geo_element_auto (parameter_type p, const A& alloc = A())
    : _data (_data_size(p), std::numeric_limits<size_type>::max(), alloc)
    {
        _data [_variant_offset] = p.variant;
        _data [_order_offset]   = p.order;
    }
  geo_element_auto (const geo_element& K)
    : _data (K._data_size(), size_type(0), A()) // cree un nouvel allocateur
    { std::copy (K._data_begin(), K._data_end(), _data.begin()); }
  geo_element_auto (const geo_element_auto<A>& K)
    : _data (K._data.size(), size_type(0), K._data.get_allocator()) // re-utilise l'allocateur precedent
    { std::copy (K._data.begin(), K._data.end(), _data.begin()); }
  template <class A2>
  geo_element_auto (const geo_element_auto<A2>& K)
    : _data (K._data.size(), size_type(0), A())         // cree un nouvel allocateur
    { std::copy (K._data.begin(), K._data.end(), _data.begin()); }
  const geo_element_auto<A>& operator= (const geo_element& K)
  {
    _data.resize(K._data_size());
    std::copy (K._data_begin(), K._data_end(), _data.begin());
    return *this;
  }
  void reset (variant_type variant, size_type order) {
        _data.resize (_data_size(variant,order), std::numeric_limits<size_type>::max());
        _data [_variant_offset] = variant;
        _data [_order_offset]   = order;
  }
  void reset (const parameter_type& param) { reset (param.variant, param.order); }


    template<class Archive>
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
        ar & _data;
    }


        iterator _data_begin()       { return _data.begin().operator->(); }
  const_iterator _data_begin() const { return _data.begin().operator->(); }
        iterator _data_end()         { return _data.end().operator->(); }
  const_iterator _data_end()   const { return _data.end().operator->(); }
 
  template <class A2> friend class geo_element_auto;


  std::vector<size_type,A> _data;
};
class geo_element_hack : public geo_element {
public:


  enum {
        _vtable_size = 1        /* = sizeof(geo_element_X_hack)/sizeof(size_type) */
  };

  typedef geo_element::size_type          size_type;
  typedef geo_element::variant_type       variant_type;
  typedef geo_element::iterator           iterator;
  typedef geo_element::const_iterator     const_iterator;
  typedef geo_element::parameter_type     parameter_type;
  typedef size_type                          raw_type;
  typedef geo_element                     generic_type;
  typedef geo_element::automatic_type     automatic_type;


  geo_element_hack () : geo_element() {}
  template <class A>
  geo_element_hack (const geo_element_auto<A>& K) : geo_element()
  {
    check_macro (K.variant() == variant(), "incompatible conversion");
#ifdef TO_CLEAN
    check_macro (K.order() == order(),     "incompatible conversion");
#endif // TO_CLEAN
    std::copy (K._data.begin(), K._data.begin() + _data_size(), _data_begin());
  }


  void reset (variant_type variant1, size_type order1) {
    check_macro (variant1 == variant(), "cannot change variant from "<<variant()<<" to "<<variant1<<" in a raw element");
#ifdef TO_CLEAN
    check_macro (order1 == 1,           "cannot change order "<<order1<< " > 1 in a raw element");
#endif // TO_CLEAN
    _set_data (_order_offset, 1);
  }
protected:

  static size_type _size_of     (const parameter_type& p) { return _vtable_size + _data_size(p); }

  iterator       _data_begin()       { return reinterpret_cast<iterator>      (this) + _vtable_size; }
  const_iterator _data_begin() const { return reinterpret_cast<const_iterator>(this) + _vtable_size; }
  iterator       _data_end()         { return _data_begin() + _data_size(); }
  const_iterator _data_end()   const { return _data_begin() + _data_size(); }

  size_type  _get_data    (size_type i) const              { return *(_data_begin() + i); }
  size_type& _get_data_ref(size_type i)                    { return *(_data_begin() + i); }
  void       _set_data    (size_type i, size_type value)   { _get_data_ref(i) = value; }

  void _reset (variant_type variant, size_type order) {
        check_macro (order == 1, "cannot set order "<<order<< " > 1 in a raw element");
        _set_data (_variant_offset, variant);
        _set_data (_order_offset,   order);
        for (size_type i = _order_offset+1, n = _data_size (variant,order); i < n; i++) {
           _set_data (i, std::numeric_limits<size_type>::max());
        }
  }
  void _set_parameter (const parameter_type& p) { _reset (p.variant, p.order); }

  template <class T, class A> friend class hack_array_seq_rep;
};
struct geo_element_permuted_put {
  typedef geo_element::size_type size_type;
  geo_element_permuted_put (const std::vector<size_type>& perm1) : perm(perm1) {}
  std::ostream& operator() (std::ostream& os, const geo_element& K) {
    static const bool do_verbose = true;
    if (do_verbose || K.size() > 2 || K.order() > 1) { os << K.name() << "\t"; }
    if (do_verbose || K.order() > 1) { os << "p" << K.order() << " "; }
    for (geo_element::size_type iloc = 0; iloc < K.n_node(); iloc++) {
      os << perm [K[iloc]];
      if (iloc < K.n_node() - 1) os << " ";
    }
    return os;
  }
  const std::vector<size_type>& perm;
};

}// namespace rheolef
#endif // _RHEOLEF_GEO_ELEMENT_V4_H