common_tree.hpp

Go to the documentation of this file.

#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: public std::binary_function< T, T, T> {
        T operator()( const T a, const T b) const{
            return (a+b)/2.0;
        }
    }; //Compute_center
} //non-exported center computation.

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)));
    }
}

} //namspace common_tree

} // namespace moab

#endif //COMMON_TREE_HPP