common_tree.hpp

Go to the documentation of this file.

/**
 * common_tree.hpp
 * Ryan H. Lewis
 * (C) 2012
 * Functionality common to all trees.
 */
#include <vector>
#include <set>
#include <iostream>
#include <map>
#include <algorithm>
#include <bitset>
#include <numeric>
#include <cmath>
#include <tr1/unordered_map>
#include <limits>
 
#ifndef COMMON_TREE_HPP
#define COMMON_TREE_HPP
#define NUM_DIM 3
#define TREE_DEBUG
namespace moab
{
namespace common_tree
{
 
 template < typename T, typename Stream >
 void print_vector( const T& v, Stream& out )
 {
 typedef typename T::const_iterator Iterator;
 out << "[ ";
 for( Iterator i = v.begin(); i != v.end(); ++i )
 {
 out << *i;
 if( i + 1 != v.end() )
 {
 out << ", ";
 }
 }
 out << " ]" << std::endl;
 }
 
#ifdef TREE_DEBUG
 template < typename T >
 void print_vector( const T& begin, const T& end )
 {
 std::cout << "[ ";
 for( T i = begin; i != end; ++i )
 {
 std::cout << ( *i )->second.second.to_ulong();
 if( i + 1 != end )
 {
 std::cout << ", ";
 }
 }
 std::cout << " ]" << std::endl;
 }
#endif
 
 template < typename _Box, typename _Point >
 bool box_contains_point( const _Box& box, const _Point& p, const double tol )
 {
 for( std::size_t i = 0; i < box.min.size(); ++i )
 {
 if( p[i] < ( box.min[i] - tol ) || p[i] > ( box.max[i] ) + tol )
 {
 return false;
 }
 }
 return true;
 }
 
 template < typename _Box >
 bool box_contains_box( const _Box& a, const _Box& b, const double tol )
 {
 for( std::size_t i = 0; i < a.min.size(); ++i )
 {
 if( b.min[i] < ( a.min[i] - tol ) )
 {
 return false;
 }
 if( b.max[i] > ( a.max[i] + tol ) )
 {
 return false;
 }
 }
 return true;
 }
 
 namespace
 {
 template < typename T >
 struct Compute_center
 {
 // deprecation of binary_function
 typedef T first_argument_type;
 typedef T second_argument_type;
 typedef T result_type;
 
 T operator()( const T a, const T b ) const
 {
 return ( a + b ) / 2.0;
 }
 }; // Compute_center
 } // namespace
 
 template < typename Vector >
 inline void compute_box_center( Vector& max, Vector& min, Vector& center )
 {
 typedef typename Vector::value_type Unit;
 center = min;
 std::transform( max.begin(), max.end(), center.begin(), center.begin(), Compute_center< Unit >() );
 }
 
 template < typename Box >
 inline typename Box::value_type compute_box_center( const Box& box, const int dim )
 {
 return ( box.max[dim] + box.min[dim] ) / 2.0;
 }
 
 template < typename T = float >
 class Box
 {
 public:
 typedef T value_type;
 typedef std::vector< T > Vector;
 Box() : max( 3, 0.0 ), min( 3, 0.0 ) {}
 Box( const Box& from ) : max( from.max ), min( from.min ) {}
 template < typename Iterator >
 Box( const Iterator begin, const Iterator end ) : max( begin, end ), min( begin, end )
 {
 }
 Box& operator=( const Box& from )
 {
 max = from.max;
 min = from.min;
 return *this;
 }
 Vector max;
 Vector min;
 }; // Box
 
 template < typename T >
 std::ostream& operator<<( std::ostream& out, const Box< T >& box )
 {
 out << "Max: ";
 print_vector( box.max, out );
 out << "Min: ";
 print_vector( box.min, out );
 return out;
 }
 
 // essentially a pair, but with an added constructor.
 template < typename T1, typename T2 >
 struct _Element_data
 {
 typedef T1 first_type;
 typedef T2 second_type;
 T1 first;
 T2 second;
 _Element_data() : first( T1() ), second( T2() ) {}
 _Element_data( const T1& x ) : first( x ), second( T2() ) {}
 _Element_data( const T1& x, T2& y ) : first( x ), second( y ) {}
 template < typename U, typename V >
 _Element_data( const _Element_data< U, V >& p ) : first( p.first ), second( p.second )
 {
 }
 }; // Element_data
 
 template < typename Entities, typename Iterator >
 void assign_entities( Entities& entities, const Iterator& begin, const Iterator& end )
 {
 entities.reserve( std::distance( begin, end ) );
 for( Iterator i = begin; i != end; ++i )
 {
 entities.push_back( std::make_pair( ( *i )->second.first, ( *i )->first ) );
 }
 }
 
 template < typename Coordinate, typename Coordinate_iterator >
 void update_bounding_max( Coordinate& max, Coordinate_iterator j )
 {
 typedef typename Coordinate::iterator Iterator;
 for( Iterator i = max.begin(); i != max.end(); ++i, ++j )
 {
 *i = std::max( *i, *j );
 }
 }
 
 template < typename Coordinate, typename Coordinate_iterator >
 void update_bounding_min( Coordinate& min, Coordinate_iterator j )
 {
 typedef typename Coordinate::iterator Iterator;
 for( Iterator i = min.begin(); i != min.end(); ++i, ++j )
 {
 *i = std::min( *i, *j );
 }
 }
 
 template < typename Box >
 void update_bounding_box( Box& a, const Box& b )
 {
 update_bounding_max( a.max, b.max.begin() );
 update_bounding_min( a.min, b.min.begin() );
#ifdef COMMON_TREE_DEBUG
 if( !box_contains_box( a, b ) )
 {
 std::cout << a << b << std::endl;
 }
#endif
 }
 
 template < typename Entity_map, typename Ordering >
 void construct_ordering( Entity_map& entity_map, Ordering& entity_ordering )
 {
 entity_ordering.reserve( entity_map.size() );
 typedef typename Entity_map::iterator Map_iterator;
 for( Map_iterator i = entity_map.begin(); i != entity_map.end(); ++i )
 {
 entity_ordering.push_back( i );
 }
 }
 
 // Input: A bunch of entity handles
 // Output: A map from handle -> Data
 // Requirements: Data contains at least a bounding box.
 // And a non-default constructor which takes only a Box&
 template < typename Entity_handles, typename Element_map, typename Bounding_box, typename Moab >
 void construct_element_map( const Entity_handles& elements,
 Element_map& map,
 Bounding_box& bounding_box,
 Moab& moab )
 {
 typedef typename Element_map::mapped_type Box_data;
 typedef typename Entity_handles::value_type Entity_handle;
 typedef typename Entity_handles::iterator Entity_handles_iterator;
 typedef typename Box_data::first_type::value_type Unit;
 typedef typename std::vector< Unit > Coordinates;
 typedef typename Coordinates::iterator Coordinate_iterator;
 
 for( Entity_handles_iterator i = elements.begin(); i != elements.end(); ++i )
 {
 // TODO: not generic enough. Why dim != 3
 const int DIM = 3;
 int num_vertices = 0;
 // Commence un-necessary deep copying.
 const Entity_handle* vertex_handle;
 moab.get_connectivity( *i, vertex_handle, num_vertices );
 Coordinates coordinate( DIM * num_vertices, 0.0 );
 moab.get_coords( vertex_handle, num_vertices, &coordinate[0] );
 Bounding_box box( coordinate.begin(), coordinate.begin() + 3 );
 if( i == elements.begin() )
 {
 bounding_box = box;
 }
 for( Coordinate_iterator j = coordinate.begin() + DIM; j != coordinate.end(); j += DIM )
 {
 update_bounding_max( box.max, j );
 update_bounding_min( box.min, j );
 }
 update_bounding_box( bounding_box, box );
 map.insert( std::make_pair( *i, Box_data( box ) ) );
 }
 }
 
} // namespace common_tree
 
} // namespace moab
 
#endif // COMMON_TREE_HPP