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