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

#ifndef _RHEOLEF_GEO_TREE_H
#define _RHEOLEF_GEO_TREE_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
///
/// =========================================================================

#include "rheolef/point.h"
#include "rheolef/geo_element.h"
#include "rheolef/heap_object.h"
namespace rheolef { 

#undef _RHEOLEF_ENABLE_GMP

#if defined(RHEOLEF_ENABLE_GMP) && defined(_RHEOLEF_HAVE_GMPXX_H)
#include <gmpxx.h>
#define _RHEOLEF_USE_GMP 1
#endif // _RHEOLEF_HAVE_GMPXX_H

class geo_tree {
public:
// typedefs and related constants:

    typedef size_t                  size_type;

#if defined(_RHEOLEF_USE_GMP)
    typedef mpz_class               big_integer_type;
    static const size_type deep_max = std::numeric_limits<Float>::digits;
#elif defined(HAVE_UNSIGNED_LONG_LONG_INT)
    typedef unsigned long long int  big_integer_type; 
    static const size_type deep_max = 8*sizeof(big_integer_type) - 2;
#else 
    typedef unsigned long int       big_integer_type; 
    static const size_type deep_max = 8*sizeof(big_integer_type) - 2;
#endif

    typedef basic_point<big_integer_type> ipoint;
    static const big_integer_type         coord_max;
   
// allocator:

    geo_tree ();
    geo_tree (size_type dim, const point& xmin, const point& xmax,
	 size_type n_vertex, std::vector<point>::const_iterator vertices,
	 size_type n_elt,    std::vector<geo_element>::const_iterator elements);
    void initialize (size_type dim, const point& xmin, const point& xmax,
	 size_type n_vertex, std::vector<point>::const_iterator vertices,
	 size_type n_elt,    std::vector<geo_element>::const_iterator elements);

// accessors:

    size_type search_element_containing (std::vector<geo_element>::const_iterator elements, 
        const point& x) const;
    size_type find_close_element (std::vector<geo_element>::const_iterator elements,
		const point& x, point& y) const;
    size_type find_close_vertex   (const point& x) const;

// inputs/outputs:

    void put_gnuplot (std::string basename = "output") const;

// internal sub-class:
    
    class tree {
    public:
      // allocator:

	tree();
	~tree();

      // accessors/modifiers:

        bool is_leaf() const;
        bool is_node() const;
	size_type  n_leaf() const;
	size_type  leaf(size_type i) const;
	size_type& leaf(size_type i);
	tree*      node(size_type i) const;
	tree*&     node(size_type i);
        void set_leaf();
        void incr_node_level();
        void add_leaf(size_type vertex_idx);
	void insert_data (size_type data_idx);
	const std::set<size_type>& get_element_list() const;

      // data:
      protected:
	union node_or_leaf {
	  // array [d_pow_2] of tree* or _vertex_index
          tree*     *_node;  
	  size_type *_leaf;
        };
	int                   _n;
	node_or_leaf          _u;
	std::set<size_type>   _l;
	void _init(bool is_node, size_type two_pow_dim);
	friend class geo_tree;
    };
    void set_node(tree* t) const;
    tree* new_node(bool is_node = true);
    tree* new_leaf(bool is_leaf = true);
    size_type  n_leaf_max() const;

// internal utilities:

    size_type find_close_vertex_i (const ipoint& ix) const;

    void bbox_element_i (const geo_element& K, ipoint& xmin, ipoint& xmax) const;
    bool belongs_to_bounding_box (const point& x) const;
    point project_into_bounding_box (const point& x) const;
    ipoint bbox_component_i (const ipoint& x, size_type sect, const big_integer_type& h) const;
    size_type sector_i (const ipoint& x, const big_integer_type& h) const;
    void insert (size_type x_idx);
    void insert (const geo_element& K);

// i/o for debuging purpose:

    void put_marker (std::ostream& out, size_type x_idx) const;
    void put_box_i  (std::ostream& out, const ipoint& x, const big_integer_type& h) const;
    void put_cross_i (std::ostream& out, const ipoint& x, const big_integer_type& h) const;
    void dump_reccursive (std::ostream& out, tree* p, size_type level) const;
    void dump (std::ostream& out = std::cerr) const;
    point  to_fcoord (const ipoint& x) const;
    ipoint to_icoord (const point& x)  const;
    ipoint round_point_i (const ipoint& x) const;
    template <class T>
    static Float bigfloat_to_float (T x);

// internal predicates on floats:

    size_type search_element_containing (std::vector<geo_element>::const_iterator elements, 
        const std::set<size_type>& l, const point& x) const;

    static bool belongs_to_segment_1d (
        const Float& x,
	const Float& a, const Float& b);

    static int sgn_mixt_2d (const point& a,
        const point& b, const point& c);

    static int sgn_mixt_3d (const point& a, const point& b,
        const point& c, const point& d);

    static bool belongs_to_triangle_2d (
	 const point& x, 
	 const point& a, const point& b, const point& c);

    static bool belongs_to_tetra_3d (
	 const point& x, 
	 const point& a, const point& b, const point& c, const point& d);

    bool belongs_to_element (const point& a, const geo_element& K) const;

// internal predicates on integers:

    bool box_intersect_box_i (
	const ipoint& x1, const big_integer_type& h1,
	const ipoint& x2, const big_integer_type& h2) const;

    bool box_intersect_box_i (
    	const ipoint& xmin, const ipoint& xmax,
    	const ipoint& x, const big_integer_type& d) const;

    bool segment_intersect_box_i (const ipoint & a, const ipoint & b,
	 const ipoint& x, const big_integer_type& h) const;
    bool segment_subset_box_i (const ipoint & a, const ipoint & b,
	 const ipoint& x, const big_integer_type& h) const;
    bool triangle_subset_box_i (const ipoint& a, const ipoint& b,
	 const ipoint& c, const ipoint& x, const big_integer_type& h) const;

    bool belongs_to_box_i (const ipoint& a,
         const ipoint& x, const big_integer_type& h) const;

    bool element_subset_box_i (const geo_element& K,
	 const ipoint& x, const big_integer_type& h) const;
    bool bbox_element_intersect_box_i (const geo_element& K,
	 const ipoint& x, const big_integer_type& h) const;
    bool element_intersect_box_i (const geo_element& K,
	 const ipoint& x, const big_integer_type& h) const;

    void put_element_list (std::vector<geo_element>::const_iterator elements,
        const std::set<size_type>& l) const;
// data:
protected:
    size_type           _dim;	    // 1, 2 or 3
    size_type           _two_pow_dim; // = 2^dim
    std::vector<ipoint> _ivertex;
    point               _xmin;
    point               _xmax;
    point               _xmin_exact;
    point               _xmax_exact;
    basic_point<Float>  _scale;
    tree*               _root;
    heap_object<tree>   _heap;

    // DEBUG purpose: internal copy of geo data
    std::vector<geo_element>  _element;
    std::vector<point>   _fvertex;
    bool                 _initialized;
public:
    // should not be used
    geo_tree (const geo_tree&);
private:
    geo_tree& operator= (const geo_tree&);
};
// -------------------------------------------------------------
// inlined
// -------------------------------------------------------------

inline
geo_tree::tree*
geo_tree::new_leaf(bool is_leaf)
{
    return new_node (!is_leaf);
}
inline
bool
geo_tree::tree::is_leaf() const
{
    return (_n >= 0);
}
inline
bool
geo_tree::tree::is_node() const
{
    return ! is_leaf();
}
inline
geo_tree::size_type
geo_tree::tree::leaf(size_type i) const
{
    return _u._leaf[i];
}
inline
geo_tree::size_type&
geo_tree::tree::leaf(size_type i)
{
    return _u._leaf[i];
}
inline
geo_tree::tree*
geo_tree::tree::node(size_type i) const
{
    return _u._node[i];
}
inline
geo_tree::tree*&
geo_tree::tree::node(size_type i)
{
    return _u._node[i];
}
inline
geo_tree::size_type
geo_tree::tree::n_leaf() const
{
    return _n;
}
inline
geo_tree::size_type
geo_tree::n_leaf_max() const
{
    return _two_pow_dim;
}
inline
void
geo_tree::set_node(tree* t) const
{
    t->_n = - _two_pow_dim;
    for (size_type i = 0; i < _two_pow_dim; i++)  
        t->_u._node[i] = 0;
}
inline
void
geo_tree::tree::incr_node_level()
{
    --_n;
}
inline
void
geo_tree::tree::add_leaf(size_type vertex_idx)
{
    _u._leaf [_n] = vertex_idx;
    _n++;
}
inline
void
geo_tree::tree::insert_data (size_type data_idx)
{
    _l.insert (data_idx);
}
inline
const std::set<geo_tree::size_type>&
geo_tree::tree::get_element_list() const
{
    return _l;
}
#ifdef _RHEOLEF_USE_GMP
template <class T>
inline
Float
geo_tree::bigfloat_to_float (T x) { return mpf_class(x).get_d(); }
#else
template <class T>
inline
Float 
geo_tree::bigfloat_to_float (T x) { return Float(x); }
#endif
}// namespace rheolef
#endif // _RHEOLEF_GEO_TREE_H