librheolef-dev 6.7-6 / usr / include / rheolef / geo_element_v4.h

#ifndef _RHEOLEF_GEO_ELEMENT_V4_H
#define _RHEOLEF_GEO_ELEMENT_V4_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
///
/// =========================================================================
//
// geo_element class with arbitrarily order
// 
// solution with a massive scratch vector
// for all geo_elements of the same variant & order
// An accessors as operator[] returns a true reference
// to the base class
//
//  - avantages :
//    * tableau massif de geo_element_e qui ont tous la meme taille
//      => la memoire est contigue globalement, l'acces plus rapide avec
//         les memoires caches actuelles
//    * la memoire n'est plus geree par les objets eux-memes mais en amont
//
//  - inconvenient : duplication de vptr et de order dans la table
//    => petit gaspillage de memoire
//
#include "rheolef/reference_element.h"
#include "rheolef/geo_element_indirect.h"
#include "rheolef/heap_allocator.h"
#include "rheolef/disarray.h"
#include "rheolef/reference_element_face_transformation.h"

#include <boost/serialization/serialization.hpp>
#include <boost/serialization/base_object.hpp>

namespace rheolef {

// --------------------------------------------------------------------
// geo_element: abstract class
// --------------------------------------------------------------------
template <class A> class geo_element_auto;

/*Class:geo_element
NAME:  @code{geo_element} - element of a mesh
@cindex  geometrical element
@clindex reference element
@clindex geo_element
@clindex reference_element
@clindex geo
DESCRIPTION:
  Defines geometrical elements and sides
  as a set of vertice and edge indexes.
  This element is obtained after a Piola transformation
  from a reference element (@pxref{reference_element iclass}).
  Indexes are related to arrays of edges and vertices.
  These arrays are included in the description of the mesh.
  Thus, this class is related of a given mesh instance
  (@pxref{geo class}).
EXAMPLE:
  This is the test of geo_element:
  @example
    geo_element_auto<> K;
    K.set_name('t') ;
    cout << "n_vertices: " << K.size()      << endl
         << "n_edges   : " << K.n_edges()   << endl
         << "dimension : " << K.dimension() << endl << endl;
    for(geo_element::size_type i = 0; i < K.size(); i++) 
        K[i] = i*10 ;
    for(geo_element::size_type i = 0; i < K.n_edges(); i++)
        K.set_edge(i, i*10+5) ;
    cout << "vertices: local -> global" << endl;
    for (geo_element::size_type vloc = 0; vloc < K.size(); vloc++)
        cout << vloc << "-> " << K[vloc] << endl;
    cout << endl 
         << "edges: local -> global" << endl;
    for (geo_element::size_type eloc = 0; eloc < K.n_edges(); eloc++) @{
        geo_element::size_type vloc1 = subgeo_local_vertex(1, eloc, 0);
        geo_element::size_type vloc2 = subgeo_local_vertex(1, eloc, 1);
        cout << eloc << "-> " << K.edge(eloc) << endl
             << "local_vertex_from_edge(" << eloc 
             << ") -> (" << vloc1 << ", " << vloc2 << ")" << endl;
    @}
  @end example
SEE ALSO: "geo"(3)
AUTHOR: Pierre.Saramito@imag.fr
DATE:   10 oct 2001, initial: 10 jan 1998
METHODS: @geo_element
End:
*/
class geo_element {
public:

// typedefs:

  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
  };
  // Implementation note: _master_offset  reserve 2 size_type but is used only for sides,
  // i.e. tri or quad in 3d mesh, edge in 2d mesh, or point in 1d
  // => waste a lot of place
  // it would be better with a polymorphic class
  // and the geo class would define an array of smart_pointers on this class
  // so, we could define edge with or without the 2 size_type for master elements
  // then, hack_array will become obsolete (good thing)
  // and reference_element could be also polymorphic, avoiding large swich (variant)
  // in all internal loops. This change of implementation will be considered
  // in the future.
  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) {} 
  };

// affectation:

  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;

// implicit conversion:

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

// accessors & modifiers:

  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 accessors:

    // seach S in all sides of K
    orientation_type get_side_informations (
  		const geo_element& S,
  		size_type& loc_isid,
  		size_type& shift) const;

    void get_side_informations (
  		const geo_element& S,
		side_information_type& sid) const;

    orientation_type get_side_orientation (const geo_element& S) const;

    // compare two sides: S and *this
    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;

// geometric predicate

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

// i/o;

  void put (std::ostream& is) const;
  void get (std::istream& os);

// static: fix orientation & shift helpers, for 2d edges & 3d faces:

  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);

// internals:
//protected:

  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
};
// specialization, from disarray.h (because of a g++ bug when T=float128)
template<>
struct _disarray_put_element_type<geo_element> {
  std::ostream& operator() (std::ostream& os, const geo_element& K) const { K.put(os); return os; }
};
template<>
struct _disarray_get_element_type<geo_element> {
  std::istream& operator() (std::istream& is, geo_element& K) const { K.get(is); return is; }
};
inline
std::istream&
operator>> (std::istream& is, geo_element& K)
{
  K.get (is);
  return is;
}
inline
std::ostream&
operator<< (std::ostream& os, const geo_element& K)
{
  K.put (os);
  return os;
}
// --------------------------------------------------------------------
// geo_element_auto: generic dynamically allocated class
// --------------------------------------------------------------------
template<class A = std::allocator<std::vector<int>::size_type> >
class geo_element_auto : public geo_element {
public:

// typedefs:

  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;

// allocators:

  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); }

// internals:

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

//protected:

        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;

// data:

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

// constants & typedefs:

  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;

// allocators:

  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());
  }

// accesors & modifiers

  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);
  }
// internals:
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;
};
// -------------------------------------------------------------------
// permuted io helper
// --------------------------------------------------------------------
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