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

#ifndef _RHEO_GEO_H
#define _RHEO_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
/// 
/// =========================================================================

/*Class:geo
NAME: @code{geo} - the finite element mesh class
@cindex  mesh
@cindex RHEOPATH environment variable
@clindex geo
@fiindex @file{.geo} mesh
@fiindex @file{.gz} gzip
@toindex @code{gzip}
SYNOPSYS:       
@noindent
The @code{geo} class defines a container for a finite element
mesh. This describes the nodes coordinates and the connectivity.
@code{geo} can contains domains, usefull for boundary condition
setting.

EXAMPLE:
@noindent
 A sample usage of the geo class writes
 @example
   geo omega;
   cin >> omega;
   cout << mayavi << full << omega;
 @end example

DESCRIPTION:
@noindent
The empty constructor makes an empty geo.  A string argument, as
@example
   	geo g("square");
@end example
@noindent
will recursively look for a @file{square.geo[.gz]} file in the
directory mentionned by the RHEOPATH environment variable,
while @code{gzip} decompression is assumed.
If the file starts with @file{.} as @file{./square} or with a
@file{/}
as in @file{/home/oscar/square}, no search occurs.
Also, if the environment variable RHEOPATH is not set, the
default value is the current directory.

@noindent
Input and output on streams are available, and manipulators
works for text or graphic output (@pxref{geo command}).

@noindent
Finally, an STL-like interface provides efficient accesses to
nodes, elements and domains.

MESH ADAPTATION:
@noindent
The @code{geo_adapt} functions performs a mesh adaptation to
improve the @code{P1}-Lagrange interpolation of the @code{field}
given in argument (@pxref{field class}). 

AXISYMETRIC GEOMETRY:
@noindent
@cindex axisymetric geometry
The @code{coordinate_system} and @code{set_coordinate_system} members
supports both @code{cartesian}, @code{rz} and @code{zr} (axisymmetric)
coordinate systems. This information is used by the @code{form}
class (@pxref{form class}).

ACCESS TO CONNECTIVITY:
@noindent
The folowing code prints triangle vertex numbers
@example
  geo omega ("circle");
  for (geo::const_iterator i = g.begin(); i != i.end(); i++) @{
    const geo_element& K = *i;
    if (K.name() != 't') continue;
    for (geo::size_type j = 0; j < 3; j++)
      cout << K [j] << " ";
    cout << endl;
  @}
@end example
@xref{geo_element internal}.

ACCESS TO VERTICE COORDINATES:
The folowing code prints vertices coordinate
@example
  for (geo::const_iterator_vertex i = g.begin_vertex(); i != g.end_node(); i++) @{
    const point& xi = *i;
    for (geo::size_type j = 0; j < g.dimension(); j++)
      cout << xi [j] << " ";
    cout << endl;
  @}
@end example

ACCESS TO DOMAINS:
 The following code prints edges on domain:
@example
  for (geo::const_iterator_domain i = g.begin_domain(); i != i.end_domain(); i++) @{
    const domain& dom = *i;
    if (dom.dimension() != 2) continue;
    for (domain::const_iterator j = dom.begin(); j < dom.end(); j++) @{
      const geo_element& E = *j;
      cout << E [0] << " " << E[1] << endl;
    @}
    cout << endl;
  @}
@end example
@xref{domain class}.

ENVIRONS:
   RHEOPATH: search path for geo data file.
   Also @ref{form class}.
AUTHORS: Pierre.Saramito@imag.fr, Jocelyn.Etienne@imag.fr
DATE:   12 may 1997  update:  26 mar 2002
End:
*/


// ============================================================================
//  class geo
// ============================================================================

#include "rheolef/georep.h"
namespace rheolef { 

//<geo:
class geo : public smart_pointer<georep> {
public:
  
// typedefs:

  typedef georep::plot_options		plot_options;
  void write_gnuplot_postscript_options (std::ostream& plot, const plot_options& opt) const;

  typedef georep::iterator              iterator;
  typedef georep::const_iterator        const_iterator;
  typedef georep::elemlist_type         elemlist_type;
  typedef georep::nodelist_type         nodelist_type;
  typedef georep::size_type             size_type;
  typedef georep::domlist_type          domlist_type;
  
  typedef georep::const_iterator_node   const_iterator_node;
  typedef georep::const_iterator_vertex const_iterator_vertex;
  typedef georep::const_iterator_domain const_iterator_domain;
  typedef georep::iterator_node         iterator_node;
  typedef georep::iterator_vertex       iterator_vertex;
  typedef georep::iterator_domain       iterator_domain;

// allocators/deallocators:
  
  explicit geo (const std::string& filename, const std::string& coord_sys = "cartesian");
  geo(const geo& g, const domain& d);
  geo(const nodelist_type& p, const geo& g);
  geo();
  
  friend geo geo_adapt (const class field& criteria, const Float& hcoef, 
	bool reinterpolate_criteria = false);

  friend geo geo_adapt (const class field& criteria, 
	const adapt_option_type& = adapt_option_type(),
	bool reinterpolate_criteria = false);
  
  friend geo geo_metric_adapt (const field& mh, 
  	const adapt_option_type& = adapt_option_type());

// input/output:
  
  friend std::istream& operator >> (std::istream&, geo&);
  friend std::ostream& operator << (std::ostream&, const geo&);         
  void save () const;
  void use_double_precision_in_saving();
  std::ostream& dump (std::ostream& s = std::cerr) const;
  std::ostream& put_mtv_domains (std::ostream&, size_type=0) const;

// accessors:
  
  const point&       vertex  (size_type i) const;
  const geo_element& element (size_type K_idx) const;
  Float              measure (const geo_element& K) const;

  std::string name() const;
  //! Family name plus number
  std::string basename() const;
  //! For moving boundary problems
  std::string familyname() const;
  //! Number of moving boundary mesh
  size_type number() const;
  //! Refinment iteration for the current mesh number
  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;

  size_type serial_number() const;

  const domain& boundary() const;
  void build_subregion(const domain& start_from, const domain& dont_cross, 
  	std::string name, std::string name_of_complement="");
  //! Builds a new domain on the side of domain `dont_cross' on which `start_from'
  //! lies. These must have no intersection.

  size_type size() const;
  size_type n_element() const;	// same as size()
  size_type n_vertex() const;
  size_type n_vertice() const;
  size_type n_node() const;	// same as n_vertex()
  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;

  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) const;
  point dehatter (const point& x_hat, size_type e) const;

  const_iterator begin() const;
  const_iterator end() const;
  const_iterator_node begin_node() const;
  const_iterator_node end_node() const;
  const_iterator_vertex begin_vertex() const;
  const_iterator_vertex end_vertex() 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;

  point normal (const geo_element& S) const;
  point normal (const geo_element& K, georep::size_type side_idx) const;
  void 
  sort_interface(const domain&, const interface&) const;
  class field
  normal (const class space& Nh, const domain& d, 
	const std::string& region="") const;
  class field
  normal (const class space& Nh, const domain& d, const interface& bc) const;
  class field
  tangent (const class space& Nh, const domain& d, 
	const std::string& region="") const;
  class field
  tangent (const class space& Nh, const domain& d, const interface& bc) const;
   //! Gives normal to d in discontinuous space Nh (P0 or P1d) and can
   //! initialize orientation of domain through `bc' data structure.
   //! Currently limited to straight-sided elements.
  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;
   //! Gives curvature of d in the (x,y) or (r,z) plane.
   //! Currently limited to straight-sided elements.
  class field
  plane_curvature_spline (const space& Nh, const domain& d, 
	const interface& bc) const;
   //! Gives curvature of d in the (x,y) or (r,z) plane based on a spline interpolation.
  class field
  plane_curvature_quadratic (const space& Nh, const domain& d, 
	const interface& bc) const;
   //! Gives curvature of d in the (x,y) or (r,z) plane based on a local quadratic interpolation.
  
  class field
  axisymmetric_curvature (const class space& Nh, const domain& d) const;
  class field
  axisymmetric_curvature (const class space& Nh, const domain& d, 
  	const interface& bc) const;
   //! Specific to "rz" and "zr" coordinate systems:
   //! Gives curvature of d in a plane orthogonal to the (r,z) plane.
   //! Can also initialize orientation of domain through `bc' data structure.
   //! Currently limited to straight-sided elements.
  void
  interface_process (const domain& d, const interface& bc,
	geometric_event& event) const;
   //! Detects event along domain d sorted according to bc's.

  // 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 geo&) const;
  bool operator != (const geo&) const;
  
// modifiers

  //! build from a list of vertex and elements:
  template <class ElementContainer, class VertexContainer>
  void build (const ElementContainer&, const VertexContainer&);

  void set_name (const std::string&);
  void set_coordinate_system (const std::string&); // "cartesian", "rz", "zr"
  void upgrade();
  void label_interface(const domain& dom1, const domain& dom2, const std::string& name);

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

// accessors to internals:
  bool localizer_initialized () const;
  void init_localizer (const domain& boundary, Float tolerance=-1, int list_size=0) const;
  const size_type* get_geo_counter() const;
  const size_type* get_element_counter() const;
  
  // TO BE REMOVED
  int gnuplot2d (const std::string& basename, 
		 plot_options& opt) const;

// check consistency
  
  void check()  const;

protected:

  friend class spacerep;
  friend class field;

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

// modifiers to internals:

  void set_dimension (size_type d);
  iterator begin();
  iterator end();
  iterator_node begin_node();
  iterator_node end_node();
  iterator_vertex begin_vertex();
  iterator_vertex end_vertex();
};
//>geo:

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

inline
geo::geo()
  : smart_pointer<georep> (new_macro(georep)) {}

inline
geo::geo (const std::string& filename, const std::string& coord_sys)
  : smart_pointer<georep> (new_macro(georep(filename,coord_sys))) {}

inline
geo::geo(const geo& g, const domain& d)
  : smart_pointer<georep> (new_macro(georep(g.data(),d))) {}

inline
geo::geo(const nodelist_type& p, const geo& g)
  : smart_pointer<georep> (new_macro(georep(p,g.data()))) {}

inline  
geo
geo_adapt (const class field& c,
	const adapt_option_type& opt,
        bool reinterpolate_criteria)
{ 
  geo new_g;
  georep_adapt (c, new_g.data(), opt, reinterpolate_criteria);
  return new_g;
}
inline
geo
geo_adapt (const class field& c,
	const Float& hcoef, 
	bool reinterpolate_criteria)
{
  adapt_option_type opt;
  opt.hcoef = hcoef;
  return geo_adapt (c, opt, reinterpolate_criteria);
}
inline
geo
geo_metric_adapt (const class field& mh, 
  	const adapt_option_type& opt)
{
	geo new_g;
	georep_metric_adapt_bamg (mh, new_g.data(), opt);
	return new_g;
}
inline
std::istream&
operator >> (std::istream& is, geo& g)
{
    return is >> g.data();
}
inline
std::ostream&
operator << (std::ostream& os, const geo& g)
{
    return os << g.data();
}
inline
void
geo::save() const
{
    data().save();
}
inline
void 
geo::use_double_precision_in_saving()
{
    data().use_double_precision_in_saving();
}
inline
std::ostream&
geo::dump(std::ostream& os) const
{
    return data().dump();
}
inline
std::ostream&
geo::put_mtv_domains (std::ostream& os, size_type max_dim) const
{
    return data().put_mtv_domains(os, (max_dim==0)?dimension():max_dim);
}
inline
void
geo::check() const
{
    data().check();
}
inline
std::string
geo::basename() const 
{
    return data().basename();  
}
inline
std::string
geo::familyname() const 
{
    return data().familyname();  
}
inline
std::string
geo::name() const 
{
    return data().name();  
}
inline
geo::size_type
geo::version() const 
{
    return data().version();  
}
inline
geo::size_type
geo::dimension() const
{
    return data().dimension();
}
inline
geo::size_type
geo::map_dimension() const
{
    return data().map_dimension();
}
inline
std::string
geo::coordinate_system() const
{
    return data().coordinate_system();
}
inline
fem_helper::coordinate_type
geo::coordinate_system_type() const
{
    return data().coordinate_system_type();
}
inline
geo::size_type
geo::number() const 
{
    return data().number();  
}
inline
geo::size_type
geo::serial_number() const 
{
    return data().serial_number();  
}
inline
const domain&
geo::boundary() const
{
    return data().boundary();
}
inline
void  
geo::build_subregion(const domain& start_from, const domain& dont_cross, 
  	std::string name, std::string name_of_complement)
{
    return data().build_subregion(start_from, dont_cross, name, name_of_complement);
}
inline
geo::size_type
geo::size() const
{
    return data().size();
}
inline
geo::size_type
geo::n_element() const
{
    return size();
}
inline
geo::size_type
geo::n_node() const
{
    return data().n_node();
}
inline
geo::size_type
geo::n_vertex() const
{
    return n_node();
}
//! TO BE REMOVED: Kept for backward compatibility
inline
geo::size_type
geo::n_vertice() const
{
    return n_node();
}
inline 
geo::size_type
geo::n_edge() const
{
  return data().n_edge();
}
inline
geo::size_type
geo::n_face() const
{
  return data().n_face();
}
inline
geo::size_type
geo::n_triangle() const
{
  return data().n_triangle();
}
inline
geo::size_type
geo::n_quadrangle() const
{
  return data().n_quadrangle();
}
inline
geo::size_type
geo::n_volume() const
{
  return data().n_volume();
}
inline
geo::size_type
geo::n_tetraedra() const
{
  return data().n_tetraedra();
}
inline
geo::size_type
geo::n_prism() const
{
  return data().n_prism();
}
inline
geo::size_type
geo::n_hexaedra() const
{
  return data().n_hexaedra();
}
inline
geo::size_type
geo::n_subgeo(georep::size_type d) const
{
  return data().n_subgeo(d) ;
}
inline
geo::size_type 
geo::n_element(reference_element::enum_type t) const
{
  return data().n_element(t);
}
inline
Float
geo::hmin() const
{
  return data().hmin();
}
inline
Float
geo::hmax() const
{
  return data().hmax();
}
inline
const point&
geo::xmin() const
{
  return data().xmin();
}
inline
const point&
geo::xmax() const
{
  return data().xmax();
}
inline
geo::const_iterator
geo::begin() const
{
  return data().begin();
}
inline
geo::const_iterator
geo::end() const
{
  return data().end();
}
inline
geo::const_iterator_node
geo::begin_node() const
{
  return data().begin_node();
}
inline
geo::const_iterator_node
geo::end_node() const
{
  return data().end_node();
}
inline
geo::const_iterator_vertex
geo::begin_vertex() const
{
  return data().begin_vertex();
}
inline
geo::const_iterator_vertex
geo::end_vertex() const
{
  return data().end_vertex();
}
inline
geo::iterator
geo::begin()
{
  return data().begin();
}
inline
geo::iterator
geo::end()
{
  return data().end();
}
inline
geo::iterator_vertex
geo::begin_vertex()
{
  return data().begin_vertex();
}
inline
geo::iterator_vertex
geo::end_vertex()
{
  return data().end_vertex();
}
inline
geo::iterator_node
geo::begin_node()
{
  return data().begin_node();
}
inline
geo::iterator_node
geo::end_node()
{
  return data().end_node();
}
inline
const domain& 
geo::get_domain(size_type i) const
{
  return data().get_domain(i);
}
inline
geo::size_type 
geo::n_domain() const
{
  return data().n_domain();
}
inline
bool
geo::has_domain (const std::string& domname) const
{
  return data().has_domain(domname);
}
inline
const domain& 
geo::operator[] (const std::string& domname) const
{
 return data().operator[](domname);
}
inline
geo::const_iterator_domain
geo::begin_domain() const
{
  return data().begin_domain();
}
inline
geo::const_iterator_domain
geo::end_domain() const
{
  return data().end_domain();
}
inline
bool
geo::operator == (const geo& x) const
{
    return data() == x.data();
}
inline
bool
geo::operator != (const geo& x) const
{
    return data() != x.data();
}
template <class ElementContainer, class VertexContainer>
inline
void
geo::build (const ElementContainer& e, const VertexContainer& v) 
{
    data().build(e,v);  
}
inline
void
geo::set_name(const std::string& s)
{
    data().set_name(s);  
}
inline
void
geo::set_dimension(size_type d)
{
    data().set_dimension(d);  
}
inline
void
geo::set_coordinate_system(const std::string& cs)
{
    data().set_coordinate_system(cs);  
}
inline
void
geo::upgrade()
{
    data().upgrade();  
}
inline
void
geo::erase_domain(const std::string& name)
{
    data().erase_domain(name);  
}
inline
void
geo::insert_domain (const domain& d)
{
    data().insert_domain(d);  
}
inline
void 
geo::label_interface(const domain& dom1, const domain& dom2, const std::string& name)
{
    data().label_interface(dom1, dom2, name);
}
inline
const geo::size_type*
geo::get_geo_counter() const
{
    return data().get_geo_counter();  
}
inline
const geo::size_type* 
geo::get_element_counter() const
{
    return data().get_element_counter();  
}
inline
void
geo::may_have_version_2() const
{
    data().may_have_version_2();
}
inline
meshpoint 
geo::hatter (const point& x, size_type K) const
{
    return data().hatter(x,K);
}
inline
point 
geo::dehatter (const meshpoint& S) const
{
    return data().dehatter(S);
}
inline
point
geo::dehatter (const point& x_hat, size_type e) const
{
    return data().dehatter(meshpoint(e,x_hat));
}
inline
bool 
geo::localize (const point& x, geo_element::size_type& element) const
{
    return data().localize (x,element);
}
inline
void
geo::localize_nearest    (const point& x, point& y, geo_element::size_type& element) const
{
    return data().localize_nearest (x,y,element);
}
inline
bool 
geo::trace (const point& x0, const point& v, point& x, Float& t, size_type& element) const
{
    return data().trace (x0, v, x, t, element); 
}
inline
void
geo::init_localizer (const domain& boundary, Float tolerance, int list_size) const
{
    data().init_localizer(boundary, tolerance, list_size);
}
inline
bool
geo::localizer_initialized () const
{
    return data().localizer_initialized();
}
inline
void 
geo::write_gnuplot_postscript_options (std::ostream& plot, const plot_options& opt) const
 {
    data().write_gnuplot_postscript_options(plot,opt);
 }
inline
int geo::gnuplot2d (const std::string& basename, 
     plot_options& opt) const
 {
    return data().gnuplot2d (basename, opt);
 }
inline
void
geo::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
{
    data().jump_interface(interface, subgeo, special_elements, node_global_to_interface);
} 
inline
void 
geo::sort_interface(const domain& d, const interface& bc) const
{ 
    data().sort_interface(d, bc);
}
inline
const point&
geo::vertex (size_type i) const
{
  return *(begin_vertex() + i);
}
inline
const geo_element&
geo::element (size_type K_idx) const
{
  return *(begin() + K_idx);
}
inline
Float
geo::measure (const geo_element& K) const
{
  return data().measure(K);
}
inline
point
geo::normal (const geo_element& S) const
{
  return data().normal (S);
}
inline
point
geo::normal (const geo_element& K, georep::size_type side) const
{
  return data().normal (K, side);
}
}// namespace rheolef
#endif // _RHEO_GEO_H