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

#ifndef _RHEO_GEOREP_H
#define _RHEO_GEOREP_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
/// 
/// =========================================================================

// ============================================================================
//  includes
// ============================================================================

#include <map>
#include "rheolef/tiny_matvec.h"
#include "rheolef/rheobase.h"
#include "rheolef/geo_element.h"
#include "rheolef/Vector.h"
#include "rheolef/point.h"
#include "rheolef/meshpoint.h"
#include "rheolef/domain.h"
#include "rheolef/reference_element.h"
#include "rheolef/rheostream.h"
#include "rheolef/adapt_option.h"
#include "rheolef/hazel.h"
#include "rheolef/hazel_1d.h"
#include "rheolef/fem_helper.h"
#include "rheolef/geo-visu.h"
#include "rheolef/interface.h"
#include "rheolef/geo_tree.h"
namespace rheolef { 

// ============================================================================
//  class georep
// ============================================================================


class georep : public occurence, public Vector<geo_element> {
public:
  
// typedefs:

  struct plot_options
   {
        bool bare;              // No comments shown
        bool empty;             // No scales shown
        bool silhouette;        // Show only the boundary domains, not the mesh
        bool meshonly;        	// Show only the mesh, not the boundary domains
        bool meshplot;        	// Show the mesh in field vector-plots
        bool xmirror;        	// Duplicate graph by symmetry about x=0 axis
        bool ymirror;        	// Duplicate graph by symmetry about y=0 axis
        point xmin;             // If xmin==xmax then use geo's
        point xmax;
	bool equal_scales;		
        Float border;           // In units of diagonal length
        bool postscript;        // Generate postscript output (gnuplot)
        bool noxplot;           // Generate postscript output and do not open
                                //  an X plot window (gnuplot)
        bool colour_postscript;
        Float linewidth;
	Float vscale;		// velocity scaling

	bool double_precision;	// Should the current mesh be saved in double precision
	
	bool execute;
	bool clean;
	bool verbose;

        plot_options(): bare(false), empty(false), silhouette(false), meshonly(false), meshplot(false),
		xmirror(false), ymirror(false),
                xmin(0,0,0), xmax(0,0,0), equal_scales(true), border(0),
                postscript(false), noxplot(false), colour_postscript(true),
                linewidth(0.1), vscale(0.2),
		execute(true), clean(true), verbose(false)
         {}
   };
  void write_gnuplot_postscript_options (std::ostream& plot, const plot_options& opt) const;

  typedef Vector<geo_element>       elemlist_type;
  typedef Vector<point>             nodelist_type;
  typedef Vector<point>::size_type  size_type;
  typedef Vector<domain>            domlist_type;
  
  typedef nodelist_type::const_iterator const_iterator_vertex;
  typedef nodelist_type::const_iterator const_iterator_node;
  typedef domlist_type::const_iterator  const_iterator_domain;
  typedef nodelist_type::iterator       iterator_vertex;
  typedef nodelist_type::iterator       iterator_node;
  typedef domlist_type::iterator        iterator_domain;

// allocators/deallocators:
  
  explicit georep (const std::string& filename, const std::string& coord_sys = "cartesian");
  georep(const georep& g, const domain& d);
  georep(const nodelist_type& p, const georep& g);
  georep(const georep&);
  georep();

  void label_interface(const domain& dom1, const domain& dom2, const std::string& name);
  void init_localizer (const domain& boundary, Float = -1, int = 0) const;
 
  // main adapt call
  friend void georep_adapt (const class field& criteria, georep& new_g,
	const adapt_option_type& = adapt_option_type(),
	bool reinterpolate_criteria = false);

  // with bamg gnerator:
  friend void georep_adapt_bamg (const class field& criteria, georep& new_g,
	const adapt_option_type& = adapt_option_type(),
	bool reinterpolate_criteria = false);

  friend void georep_adapt_bamg (const Vector<class field>& criteria, georep& new_g,
	const adapt_option_type& = adapt_option_type(),
	bool reinterpolate_criteria = false);

  friend void georep_metric_adapt_bamg (const class field& metric, georep& new_g,
	const adapt_option_type& = adapt_option_type());

  // with grummp gnerator:
  friend void georep_adapt_grummp (const class field& criteria, georep& new_g,
	const adapt_option_type& = adapt_option_type(),
	bool reinterpolate_criteria = false);

// modifiers:

  void build_subregion(const domain& start_from, const domain& dont_cross, 
  	std::string name, std::string name_of_complement);

// input/output:
  
  friend std::istream& operator >> (std::istream&, georep&);
  friend std::ostream& operator << (std::ostream&, const georep&);         
  void save () const;
  inline void use_double_precision_in_saving() { _save_double_precision=true; }
  std::ostream& dump (std::ostream& s = std::cerr) const;

// check consistency
  
  void check()  const;
  void check (const geo_element&) const;

// accessors:
  
  std::string name() const;
  std::string basename() const;
  std::string familyname() const;
  size_type version() const;
  size_type dimension() const;
  size_type map_dimension() const;
  std::string coordinate_system () const; // "cartesian", "rz", "zr"
  fem_helper::coordinate_type coordinate_system_type() const { return _coord_sys; }
  size_type number() const;
  size_type serial_number() const;
  bool localizer_initialized () const;
  const domain& boundary () const;

  size_type size() const;
  size_type n_node() const;
  size_type n_vertex() const;
  size_type n_edge() const ;
  size_type n_face() const ;
  size_type n_triangle() const ;
  size_type n_quadrangle() const ;
  size_type n_volume() const ;
  size_type n_tetraedra() const ;
  size_type n_prism() const ;
  size_type n_hexaedra() const ;
  size_type n_subgeo(size_type d) const;
  size_type n_element(reference_element::enum_type t) const;

  const geo_element& element (size_type K_idx) const;
  const point&       vertex  (size_type i) const;
  Float measure (const geo_element& K) const;
  Float hmin() const;
  Float hmax() const;
  const point& xmin() const;
  const point& xmax() const;

  meshpoint hatter (const point& x, size_type K) const;
  point dehatter (const meshpoint& S, bool is_a_vector=false) const;

  const_iterator_node begin_vertex() const;
  const_iterator_node end_vertex() const;
  const_iterator_node begin_node() const;
  const_iterator_node end_node() const;

  // localizer:
  bool localize (const point& x, geo_element::size_type& element) const;
  void localize_nearest (const point& x, point& y, geo_element::size_type& element) const;
  bool trace (const point& x0, const point& v, point& x, Float& t, size_type& element) const;

  // access to domains:
  const domain& get_domain(size_type i) const;
  size_type n_domain() const;
  bool has_domain (const std::string& domname) const;
  const domain& operator[] (const std::string& domname) const;
  const_iterator_domain begin_domain() const;
  const_iterator_domain end_domain() const;
  
  void
  set_interface_orientation(const domain& d, const interface& bc) const;
  void 
  sort_interface(const domain&, const interface&) const;
  void
  interface_process (const domain& d, const interface& bc,
	geometric_event& event) const;
  class field
  get_normal (const class space& Nh, const domain& d,
  	bool axisymmetric_curvature, const interface& bc) const;
  class field
  tangent_spline (const space& Nh, const domain& d, const interface& bc) const;
  class field
  plane_curvature (const space& Nh, const domain& d, 
	const interface& bc) const;
  class field
  plane_curvature_spline (const space& Nh, const domain& d, 
	const interface& bc) const;
  class field
  plane_curvature_quadratic (const space& Nh, const domain& d, 
	const interface& bc) const;
  class field
  axi_curvature (const space& Nh, const domain& d, 
	const interface& bc) const;

  point normal (const geo_element& K, georep::size_type side) const;
  point normal (const geo_element& S) const;

  // construction of jump-interface domain data:
  void 
  jump_interface(const domain& interface, const domain& subgeo, 
	std::map<size_type,tiny_vector<size_type> >& special_elements, 
        std::map<size_type,size_type>& node_global_to_interface) const;

// comparator:

  bool operator == (const georep&) const;
  bool operator != (const georep&) const;
  
// modifiers

  template <class ElementContainer, class VertexContainer>
  void build (const ElementContainer&, const VertexContainer&);
 
  void set_name (const std::string&);
  void set_dimension (size_type d);
  void set_coordinate_system (const std::string&); // "cartesian", "rz", "zr"
  void set_serial_number (size_type);
  void set_number (size_type);
  void upgrade();

  // modify domains
  void erase_domain (const std::string& name);
  void insert_domain (const domain&);


  iterator_node begin_vertex();
  iterator_node end_vertex();
  iterator_node begin_node();
  iterator_node end_node();

  iterator_domain begin_domain();
  iterator_domain end_domain();


// counters by geo-dimension and element-type
protected:

  // n_point, n_edge, n_face, n_volume
  const size_type* get_geo_counter() const;

  // 	.. n_triangle ... n_tetraedra ..
  const size_type* get_element_counter() const;

  void reset_counters();
  friend void build_connectivity (georep&);

  void may_have_version_2() const; // fatal if not !

  class field get_normal_P0 (const space& Nh, const domain& d, 
			     bool axisymmetric_curvature, const interface& bc) const;
  class field get_normal_P1P2 (const space& Nh, const domain& d, 
			       const interface& bc) const;
// low level i/o primitives:
public:
  std::istream& get_rheo (std::istream&);
  std::ostream& put(std::ostream&) const;
  std::ostream& put_vtk (std::ostream&, bool lattice = false) const;
  std::ostream& put_vtk_P2_iso_P1 (std::ostream&) const;
  std::ostream& put_vtk_P1d (std::ostream&) const;
  std::istream& get_vtk_polydata (std::istream&);
  std::ostream& put_vtk_polydata (std::ostream&) const;
  std::ostream& put_bamg (std::ostream&) const;
  std::istream& get_bamg (std::istream&, bool check_optional_domain_name = true);
  std::ostream& put_gmshcad (std::ostream&, std::string bdry_name = "") const;
  std::ostream& put_gmsh(std::ostream&, bool only_boundary = false) const;
  std::istream& get_gmsh(std::istream&);
  std::istream& get_grummp(std::istream&);
  std::ostream& put_qmg (std::ostream&) const;
  std::istream& get_qmg (std::istream&);
  std::ostream& put_stl (std::ostream&, const domain&) const;

  std::istream& get_cemagref (std::istream&);
  std::istream& get_cemagref (std::istream&, class field& z, bool load_z = true);
  std::ostream& put_cemagref (std::ostream&) const;
  std::ostream& put_cemagref (std::ostream&, const class field& z, bool store_z = true) const;

  std::ostream& put_mmg3d (std::ostream&, bool only_boundary = false) const;
  std::istream& get_mmg3d (std::istream&, bool check_optional_domain_name = true);

  std::ostream& put_mtv_domains (std::ostream&, size_type max_dim) const;

  std::istream& get_tetgen (std::istream&, bool check_optional_domain_name = true);
  std::ostream& put_tetgen (std::ostream&, bool only_boundary = false) const;

  int gnuplot1d (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int gnuplot2d (const std::string& basename, 
		 plot_options& opt) const;
  int gnuplot2d (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int gnuplot3d (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int plotmtv   (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int atom      (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int x2d       (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true) const;
  int vtk_tcl   (const std::string& basename, 
		 bool execute = true, bool clean = true, 
		 bool verbose = true, bool fill = false,
		 bool shrink  = false,bool tube = false,
		 bool ball    = false,bool full = false,
		 bool cut     = false,bool split= false,
		 const point& origin = point(std::numeric_limits<Float>::max(),0,0),
		 const point& normal = point()) const;
  int visu_mayavi (const geo_visu_option_type& opt) const;
protected:

  friend class geo;

// data:
protected:

  std::string 	      _name;
  size_type	      _number;
  size_type	      _serial_number;
  size_type           _dim;        	// 0, 1, 2 or 3 for elements
  fem_helper::coordinate_type _coord_sys;  // cartesian, axisymmetric
  nodelist_type       _x;		// coordinates
  domlist_type 	      _domlist;    	// domains (boundaries, etc..)
  point 	      _xmin;		// bounding box
  point 	      _xmax;
  bool		      _save_double_precision;
 
  size_type           _count_geo [4];
  size_type           _count_element [reference_element::max_size];
  size_type           _version;
  mutable geo_tree    _tree;
  mutable hazel	      _h;
  mutable hazel_1d    _h1d;
  mutable bool        _localizer_initialized;
  mutable domain      _boundary;
  mutable bool        _boundary_initialized;
};
void georep_calculate_metric_bamg (const Vector<class field>& criteria,
	class field& metric, const adapt_option_type& opt);
  
void georep_calculate_metric_bamg (const class field& u0,
	class field& metric, const adapt_option_type& opt);

// ================== [ inlined ] ==================================================

inline
georep::georep()

 : occurence(),
   Vector<geo_element>(),
   _name(),
   _number(std::numeric_limits<size_type>::max()),
   _serial_number(std::numeric_limits<size_type>::max()),
   _dim(0), 
   _coord_sys(fem_helper::cartesian),
   _x(), 
   _domlist(), 
   _xmin(), 
   _xmax(),
   _save_double_precision(false),
   _version(0),
   _tree(),
   _h(),
   _h1d(),
   _localizer_initialized(false),
   _boundary(),
   _boundary_initialized(false)
{
   reset_counters();
}
inline
std::string
georep::basename() const 
{
  if (_number == std::numeric_limits<size_type>::max())
        return _name;
  else
     return _name + "+" + itos(_number);
}
inline
std::string
georep::familyname() const 
{
  return _name ;
}
inline
std::string
georep::name() const 
{
  if ( (_number == std::numeric_limits<size_type>::max())
	 &&(_serial_number == std::numeric_limits<size_type>::max()) )
		return _name;
  else if (_serial_number == std::numeric_limits<size_type>::max())
        return _name + "+" + itos(_number);
  else if (_number == std::numeric_limits<size_type>::max())
        return _name + "-" + itos(_serial_number);
  else
        return _name + "+" + itos(_number) + "-" + itos(_serial_number);
}
inline
void
georep::set_name(const std::string& s)
{
  _name = s;
  _serial_number = std::numeric_limits<size_type>::max();
  _number = std::numeric_limits<size_type>::max();
}
inline
void
georep::set_dimension (size_type d)
{
  _dim = d;
}
inline
void
georep::set_number(size_type i)
{
  _number = i;
}
inline
void
georep::set_serial_number(size_type i)
{
  _serial_number = i;
}
inline
bool
georep::operator != (const georep& x) const
{
    return !(*this == x);
}
inline
georep::size_type
georep::dimension() const
{
  return _dim;
}
inline
georep::size_type
georep::version() const
{
  return _version;
}
inline
georep::size_type
georep::number() const
{
  return _number;
}
inline
georep::size_type
georep::serial_number() const
{
  return _serial_number;
}
inline
georep::size_type
georep::size() const
{
  return Vector<geo_element>::size();
}
inline
georep::size_type
georep::n_node() const
{
  return _x.size();
}
inline
georep::size_type
georep::n_vertex() const
{
  return _x.size();
}
inline 
georep::size_type
georep::n_edge() const
{
  return _count_geo[1];
}
inline
georep::size_type
georep::n_face() const
{
  return _count_geo[2] ;
}
inline
georep::size_type
georep::n_triangle() const
{
  return _count_element[reference_element::t] ;
}
inline
georep::size_type
georep::n_quadrangle() const
{
  return _count_element[reference_element::q] ;
}
inline
georep::size_type
georep::n_volume() const
{
  return _count_geo[3] ;
}
inline
georep::size_type
georep::n_tetraedra() const
{
  return _count_element[reference_element::T] ;
}
inline
georep::size_type
georep::n_prism() const
{
  return _count_element[reference_element::P] ;
}
inline
georep::size_type
georep::n_hexaedra() const
{
  return _count_element[reference_element::H] ;
}
inline
georep::size_type
georep::n_subgeo(georep::size_type d) const
{
  return _count_geo[d] ;
}
inline
georep::size_type
georep::n_element(reference_element::enum_type t) const
{
  return _count_element[t];
}
inline
georep::const_iterator_node
georep::begin_node() const
{
  return _x.begin();
}
inline
georep::const_iterator_node
georep::end_node() const
{
  return _x.end();
}
inline
georep::const_iterator_vertex
georep::begin_vertex() const
{
  return _x.begin();
}
inline
georep::const_iterator_vertex
georep::end_vertex() const
{
  return _x.end();
}
inline
georep::iterator_vertex 
georep::begin_vertex()
{
  return _x.begin();
}
inline
georep::iterator_vertex 
georep::end_vertex()
{
  return _x.end();
}
inline
georep::iterator_node 
georep::begin_node()
{
  return _x.begin();
}
inline
georep::iterator_node 
georep::end_node()
{
  return _x.end();
}
inline
const point&
georep::xmin() const
{
  return _xmin;
}
inline
const point&
georep::xmax() const
{
  return _xmax;
}
inline
const domain& 
georep::get_domain(size_type i) const
{
    return _domlist.at(i);
}
inline
georep::const_iterator_domain
georep::begin_domain() const
{
  return _domlist.begin();
}
inline
georep::const_iterator_domain
georep::end_domain() const
{
  return _domlist.end();
}
inline
georep::iterator_domain 
georep::begin_domain()
{
  return _domlist.begin();
}
inline
georep::iterator_domain 
georep::end_domain()
{
  return _domlist.end();
}
inline
georep::size_type 
georep::n_domain() const
{
  return _domlist.size();
}
inline
const georep::size_type *
georep::get_geo_counter() const
{
  return _count_geo;
}
inline
const georep::size_type *
georep::get_element_counter() const
{
  return _count_element;
}
inline
bool
georep::localizer_initialized () const
{
    return _localizer_initialized;
}
std::istream&
get_nodes (
        std::istream& s,
        georep::iterator_node iter_p,     
        georep::iterator_node last_p,   
        point& xmin, 
        point& xmax,
        georep::size_type dim);

std::istream&
get_nodes_bamg (
        std::istream& s,
        georep::iterator_node iter_p,
        georep::iterator_node last_p,
        point& xmin,
        point& xmax,
        std::vector<domain>& domains,
        std::vector<size_t>& vertice_domain_bamg_number);

inline
void
georep::write_gnuplot_postscript_options (std::ostream& plot, const plot_options& opt) const
{
	plot << "set terminal postscript portrait " << (opt.colour_postscript?("colour "):("mono "))
	     << "solid " << "linewidth " << opt.linewidth << std::endl;

}
inline
const geo_element&
georep::element (size_type K_idx) const 
{
  return *(begin()+K_idx);
}
inline
const point&
georep::vertex (size_type i) const 
{
  return *(begin_vertex()+i);
}
}// namespace rheolef
#endif // _RHEO_GEOREP_H