librheolef-dev 5.93-2 / usr / include / rheolef / geo_element_v2.h

#ifndef _RHEOLEF_GEO_ELEMENT_V2_H
#define _RHEOLEF_GEO_ELEMENT_V2_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
///
/// =========================================================================
//
// the geo_element hierarchy of classes
//
// author: Pierre.Saramito@imag.fr
//
// date: 19 december 2010
//
#include "rheolef/distributed.h"
#include "rheolef/reference_element.h"
#include "rheolef/polymorphic_traits.h"

#include <boost/serialization/serialization.hpp>
#include <boost/array.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/mpl/size.hpp>

namespace rheolef {

/*Class:geo_element
NAME:  @code{geo_element} - element of a mesh
@cindex  geometrical element
@clindex geo_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 internal}).
  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}).

IMPLEMENTATION NOTE:
  @code{geo_element} is a pure virtual class, basis of a hierarchy
  of classes such as @code{geo_element_t} for triangles, etc.
SEE ALSO: "geo"(3)
DATE:   20 dec 2010
End:
*/
//<verbatim:
class geo_element {
public:
    typedef std::vector<int>::size_type  size_type;
    typedef reference_element::enum_type enum_type; 
    typedef short int                    orientation_type; // for sign (+1,-1)

    static char tab_name [reference_element::max_size];

    virtual ~geo_element () {}
    geo_element& operator= (const geo_element& K) {
        _dis_ie     = K._dis_ie;
        _ios_dis_ie = K._ios_dis_ie;
        return *this;
    }
    template<class Reference>
    static void intrusive_set (Reference lhs, const geo_element& K);
    virtual enum_type variant() const = 0;
    enum_type type() const { return variant(); } // backward compatibility
    virtual char      name() const = 0;
    virtual size_type size() const = 0;
    virtual size_type dimension() const = 0;
    operator reference_element () const { return reference_element(variant()); }
    virtual size_type  operator[] (size_type i) const = 0;
    virtual size_type& operator[] (size_type i)       = 0;
    size_type dis_ie() const { return _dis_ie; }
    size_type ios_dis_ie() const { return _ios_dis_ie; } 
    void set_dis_ie    (size_type dis_ie1)     { _dis_ie = dis_ie1; } 
    void set_ios_dis_ie(size_type ios_dis_ie1) { _ios_dis_ie = ios_dis_ie1; } 
    size_type index() const { return dis_ie(); } // for backward compatibility
    virtual size_type edge (size_type i_edge) const { error_macro ("not yet"); return 0; } 

    bool get_orientation_and_shift (const geo_element& S, 
	  orientation_type& orient, size_type& shift) const;

    template <class Pool> struct get;
    template<class Archive>
    void serialize (Archive& ar, const unsigned int version) {
	ar & _dis_ie;
	ar & _ios_dis_ie;
    }
// data:
protected:
    size_type _dis_ie;
    size_type _ios_dis_ie;
};
//>verbatim:
// =====================================================================
class geo_element_p : public geo_element {
    static const enum_type _variant =  reference_element::p;
    static const char      _name = 'p';
    static const size_type _size =  1;
    static const size_type _dim  =  0;
public:
    typedef boost::array<size_type,_size> array_type;

    geo_element_p () : geo_element() {
	geo_element::_dis_ie = std::numeric_limits<size_type>::max();
    }
    geo_element_p (size_type a) : geo_element() {
	geo_element::_dis_ie = a;
    } 
    ~geo_element_p () {}
    geo_element_p& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_p& K1 = dynamic_cast<const geo_element_p&>(K);
	geo_element::_dis_ie = K.geo_element::dis_ie();
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type) const { return geo_element::_dis_ie; }
    size_type& operator[] (size_type)       { return geo_element::_dis_ie; }
    template<class Archive>
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
    }
};
// =====================================================================
class geo_element_e : public geo_element {
    static const enum_type _variant =  reference_element::e;
    static const char      _name = 'e';
    static const size_type _size =  2;
    static const size_type _dim  =  1;
public:
    typedef boost::array<size_type,_size> array_type;

    geo_element_e () : geo_element() {
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_e (size_type a, size_type b) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
    } 
    geo_element_e (const array_type& K) {  _vertex = K; }
    ~geo_element_e () {}
    geo_element_e& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_e& K1 = dynamic_cast<const geo_element_e&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive>
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// =====================================================================
class geo_element_t : public geo_element {
    static const enum_type _variant =  reference_element::t;
    static const char      _name = 't';
    static const size_type _size =  3;
    static const size_type _dim  =  2;
public:
    typedef boost::array<size_type,_size> array_type;

    geo_element_t () : geo_element() {
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_t (size_type a, size_type b, size_type c) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
	_vertex[2] = c;
    } 
    geo_element_t (const array_type& K) {  _vertex = K; }
    ~geo_element_t () {}
    geo_element_t& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_t& K1 = dynamic_cast<const geo_element_t&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive>
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// =====================================================================
class geo_element_q : public geo_element {
    static const enum_type _variant =  reference_element::q;
    static const char      _name = 'q';
    static const size_type _size =  4;
    static const size_type _dim  =  2;
public:
    typedef boost::array<size_type,_size> array_type;
    geo_element_q () : geo_element() { 
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_q (size_type a, size_type b, size_type c, size_type d) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
	_vertex[2] = c;
	_vertex[3] = d;
    } 
    geo_element_q (const array_type& K) {  _vertex = K; }
    ~geo_element_q () {}
    geo_element_q& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_q& K1 = dynamic_cast<const geo_element_q&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive> 
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// =====================================================================
class geo_element_T : public geo_element {
    static const enum_type _variant =  reference_element::T;
    static const char      _name = 'T';
    static const size_type _size =  4;
    static const size_type _dim  =  3;
public:
    typedef boost::array<size_type,_size> array_type;
    geo_element_T () : geo_element() { 
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_T (size_type a, size_type b, size_type c, size_type d) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
	_vertex[2] = c;
	_vertex[3] = d;
    } 
    geo_element_T (const array_type& K) {  _vertex = K; }
    ~geo_element_T () {}
    geo_element_T& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_T& K1 = dynamic_cast<const geo_element_T&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive> 
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// =====================================================================
class geo_element_P : public geo_element {
    static const enum_type _variant =  reference_element::P;
    static const char      _name = 'P';
    static const size_type _size =  6;
    static const size_type _dim  =  3;
public:
    typedef boost::array<size_type,_size> array_type;
    geo_element_P () : geo_element() { 
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_P (size_type a, size_type b, size_type c, size_type d,
                   size_type e, size_type f) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
	_vertex[2] = c;
	_vertex[3] = d;
	_vertex[4] = e;
	_vertex[5] = f;
    } 
    geo_element_P (const array_type& K) {  _vertex = K; }
    ~geo_element_P () {}
    geo_element_P& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_P& K1 = dynamic_cast<const geo_element_P&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive> 
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// =====================================================================
class geo_element_H : public geo_element {
    static const enum_type _variant =  reference_element::H;
    static const char      _name = 'H';
    static const size_type _size =  8;
    static const size_type _dim  =  3;
public:
    typedef boost::array<size_type,_size> array_type;
    geo_element_H () : geo_element() { 
	std::fill (_vertex.begin(), _vertex.end(), std::numeric_limits<size_type>::max());
    }
    geo_element_H (size_type a, size_type b, size_type c, size_type d,
                   size_type e, size_type f, size_type g, size_type h) : geo_element() {
	_vertex[0] = a;
	_vertex[1] = b;
	_vertex[2] = c;
	_vertex[3] = d;
	_vertex[4] = e;
	_vertex[5] = f;
	_vertex[6] = g;
	_vertex[7] = h;
    } 
    geo_element_H (const array_type& K) {  _vertex = K; }
    ~geo_element_H () {}
    geo_element_H& operator= (const geo_element& K) {
        geo_element::operator= (K);
	check_macro (K.variant() == variant(), "incompatible assignment from `"<<K.name()<<"' to `"<<name()<<"' geo_element");	
	const geo_element_H& K1 = dynamic_cast<const geo_element_H&>(K);
        std::copy (K1._vertex.begin(), K1._vertex.end(), _vertex.begin());
        return *this;
    }
    enum_type variant() const   { return _variant; };
    char      name() const      { return _name; };
    size_type size() const      { return _size; }
    size_type dimension() const { return _dim; };
    size_type  operator[] (size_type i) const { return _vertex[i]; }
    size_type& operator[] (size_type i)       { return _vertex[i]; }
    template<class Archive> 
    void serialize (Archive& ar, const unsigned int version) {
        ar & boost::serialization::base_object<geo_element>(*this);
	ar & _vertex;
    }
// data:
protected:
    array_type _vertex;
};
// -------------------------------------------------------------------
// input associated to pool allocation
// -------------------------------------------------------------------
template <class Pool>
struct geo_element::get {
  geo_element* operator() (std::istream& is, Pool& pool) {
    extern char skip_blancs_and_tabs (std::istream&);
    char c;
    is >> std::ws >> c;
    if (!isdigit(c)) {
      geo_element* ptr = 0;
      switch (c) {
        case 'p': ptr = pool.template new_entry<geo_element_p>(); break;
        case 'e': ptr = pool.template new_entry<geo_element_e>(); break;
        case 't': ptr = pool.template new_entry<geo_element_t>(); break;
        case 'q': ptr = pool.template new_entry<geo_element_q>(); break;
        case 'T': ptr = pool.template new_entry<geo_element_T>(); break;
        case 'P': ptr = pool.template new_entry<geo_element_P>(); break;
        case 'H': ptr = pool.template new_entry<geo_element_H>(); break;
        default : error_macro ("unexpected element variant `"<<c<<"'");
      }
      geo_element& K = *ptr;
      for (size_t i = 0, n = K.size(); i < n; i++) {
	is >> K[i];
      }
      return ptr;
    }
    // here, starts with a digit: triangle (e.g. "21 22 23"), edge (e.g. "12 13") or point (e.g. "11")
    // end of element input: with end-of-line or non-digit input char
    // Note: with 4 vertices, there is an ambiguity between quadrangle and thetraedra
    size_type tmp [3];
    size_type nvertex = 0;
    while (is.good() && isdigit(c) && nvertex < 3) {
        is.unget(); // or is.putback(c); ?
        is >> tmp [nvertex];
        nvertex++;
        c = skip_blancs_and_tabs (is);
    }
    geo_element* ptr = 0;
    switch (nvertex) {
      case 1: ptr = pool.template new_entry<geo_element_p>(); break;
     case 2: ptr = pool.template new_entry<geo_element_e>(); break;
     case 3: ptr = pool.template new_entry<geo_element_t>(); break;
     default: error_macro ("unexpected element with "<<nvertex<<" vertices");
    }
    geo_element& K = *ptr;
    for (size_t i = 0, n = K.size(); i < n; i++) {
      K[i] = tmp[i];
    }
    return ptr;
  }
};
// -------------------------------------------------------------------
// output
// -------------------------------------------------------------------
inline
std::ostream&
operator<< (std::ostream& os, const geo_element& K) {
    if (K.size() > 2) { os << K.name() << "\t"; }
    for (geo_element::size_type iloc = 0; iloc < K.size(); iloc++) {
      os << K[iloc];
      if (iloc < K.size() - 1) os << " ";
    }
    return os;
}
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) {
    if (K.size() > 2) { os << K.name() << "\t"; }
    for (geo_element::size_type iloc = 0; iloc < K.size(); iloc++) {
      os << perm [K[iloc]];
      if (iloc < K.size() - 1) os << " ";
    }
    return os;
  }
  const std::vector<size_type>& perm;
};
template <class Array>
void
dump (const geo_element& K1, Array& K2) {
  	check_macro (K1.size() == K2.size(), "invalid element sizes: K1.size="<<K1.size()<<" and K2.size="<<K2.size());
        for (geo_element::size_type iloc = 0; iloc < K1.size(); iloc++) {
       	    K2[iloc] = K1[iloc];
	}
} 
// -------------------------------------------------------------------
// helper for polymorphic_array and polymorphic_map
// -------------------------------------------------------------------
template<class Reference>
inline
void
geo_element::intrusive_set (Reference lhs, const geo_element& K) 
{
  switch (K.name()) {
    case 'p': lhs = dynamic_cast<const geo_element_p&>(K); break;
    case 'e': lhs = dynamic_cast<const geo_element_e&>(K); break;
    case 't': lhs = dynamic_cast<const geo_element_t&>(K); break;
    case 'q': lhs = dynamic_cast<const geo_element_q&>(K); break;
    case 'T': lhs = dynamic_cast<const geo_element_T&>(K); break;
    case 'P': lhs = dynamic_cast<const geo_element_P&>(K); break;
    case 'H': lhs = dynamic_cast<const geo_element_H&>(K); break;
    default: error_macro ("unexpected geo_element type `" << K.name() << "'");
  }
}
// -------------------------------------------------------------------
// helper for long typedef list:
// -------------------------------------------------------------------

namespace mpl = boost::mpl;

// default polymorphic array derived type allows the 3d full family:
template <>
struct polymorphic_traits<geo_element> {
    typedef mpl::vector7<
	geo_element_p,
	geo_element_e,
	geo_element_t,
	geo_element_q,
	geo_element_T,
	geo_element_P,
	geo_element_H>
    derived_type;
    static const size_t n_derived = mpl::size<derived_type>::value;
};

} // namespace rheolef

#ifdef _RHEOLEF_HAVE_MPI
// =====================================================================
// Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
// When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous 
// copy operations.
// To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
namespace boost {
 namespace mpi {
  template <> struct is_mpi_datatype<rheolef::geo_element_p> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_e> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_t> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_q> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_T> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_P> : mpl::true_ { };
  template <> struct is_mpi_datatype<rheolef::geo_element_H> : mpl::true_ { };
 } // namespace mpi
} // namespace boost
#endif // _RHEOLEF_HAVE_MPI

#endif // _RHEOLEF_GEO_ELEMENT_V2_H