Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
common_tree.hpp
Go to the documentation of this file.
00001 /**
00002  * common_tree.hpp
00003  * Ryan H. Lewis
00004  * (C) 2012
00005  * Functionality common to all trees.
00006  */
00007 #include <vector>
00008 #include <set>
00009 #include <iostream>
00010 #include <map>
00011 #include <algorithm>
00012 #include <bitset>
00013 #include <numeric>
00014 #include <cmath>
00015 #include <tr1/unordered_map>
00016 #include <limits>
00017 
00018 #ifndef COMMON_TREE_HPP
00019 #define COMMON_TREE_HPP
00020 #define NUM_DIM 3
00021 #define TREE_DEBUG
00022 namespace moab
00023 {
00024 namespace common_tree
00025 {
00026 
00027     template < typename T, typename Stream >
00028     void print_vector( const T& v, Stream& out )
00029     {
00030         typedef typename T::const_iterator Iterator;
00031         out << "[ ";
00032         for( Iterator i = v.begin(); i != v.end(); ++i )
00033         {
00034             out << *i;
00035             if( i + 1 != v.end() )
00036             {
00037                 out << ", ";
00038             }
00039         }
00040         out << " ]" << std::endl;
00041     }
00042 
00043 #ifdef TREE_DEBUG
00044     template < typename T >
00045     void print_vector( const T& begin, const T& end )
00046     {
00047         std::cout << "[ ";
00048         for( T i = begin; i != end; ++i )
00049         {
00050             std::cout << ( *i )->second.second.to_ulong();
00051             if( i + 1 != end )
00052             {
00053                 std::cout << ", ";
00054             }
00055         }
00056         std::cout << " ]" << std::endl;
00057     }
00058 #endif
00059 
00060     template < typename _Box, typename _Point >
00061     bool box_contains_point( const _Box& box, const _Point& p, const double tol )
00062     {
00063         for( std::size_t i = 0; i < box.min.size(); ++i )
00064         {
00065             if( p[i] < ( box.min[i] - tol ) || p[i] > ( box.max[i] ) + tol )
00066             {
00067                 return false;
00068             }
00069         }
00070         return true;
00071     }
00072 
00073     template < typename _Box >
00074     bool box_contains_box( const _Box& a, const _Box& b, const double tol )
00075     {
00076         for( std::size_t i = 0; i < a.min.size(); ++i )
00077         {
00078             if( b.min[i] < ( a.min[i] - tol ) )
00079             {
00080                 return false;
00081             }
00082             if( b.max[i] > ( a.max[i] + tol ) )
00083             {
00084                 return false;
00085             }
00086         }
00087         return true;
00088     }
00089 
00090     namespace
00091     {
00092         template < typename T >
00093         struct Compute_center
00094         {
00095             // deprecation of binary_function
00096             typedef T first_argument_type;
00097             typedef T second_argument_type;
00098             typedef T result_type;
00099 
00100             T operator()( const T a, const T b ) const
00101             {
00102                 return ( a + b ) / 2.0;
00103             }
00104         };  // Compute_center
00105     }       // namespace
00106 
00107     template < typename Vector >
00108     inline void compute_box_center( Vector& max, Vector& min, Vector& center )
00109     {
00110         typedef typename Vector::value_type Unit;
00111         center = min;
00112         std::transform( max.begin(), max.end(), center.begin(), center.begin(), Compute_center< Unit >() );
00113     }
00114 
00115     template < typename Box >
00116     inline typename Box::value_type compute_box_center( const Box& box, const int dim )
00117     {
00118         return ( box.max[dim] + box.min[dim] ) / 2.0;
00119     }
00120 
00121     template < typename T = float >
00122     class Box
00123     {
00124       public:
00125         typedef T value_type;
00126         typedef std::vector< T > Vector;
00127         Box() : max( 3, 0.0 ), min( 3, 0.0 ) {}
00128         Box( const Box& from ) : max( from.max ), min( from.min ) {}
00129         template < typename Iterator >
00130         Box( const Iterator begin, const Iterator end ) : max( begin, end ), min( begin, end )
00131         {
00132         }
00133         Box& operator=( const Box& from )
00134         {
00135             max = from.max;
00136             min = from.min;
00137             return *this;
00138         }
00139         Vector max;
00140         Vector min;
00141     };  // Box
00142 
00143     template < typename T >
00144     std::ostream& operator<<( std::ostream& out, const Box< T >& box )
00145     {
00146         out << "Max: ";
00147         print_vector( box.max, out );
00148         out << "Min: ";
00149         print_vector( box.min, out );
00150         return out;
00151     }
00152 
00153     // essentially a pair, but with an added constructor.
00154     template < typename T1, typename T2 >
00155     struct _Element_data
00156     {
00157         typedef T1 first_type;
00158         typedef T2 second_type;
00159         T1 first;
00160         T2 second;
00161         _Element_data() : first( T1() ), second( T2() ) {}
00162         _Element_data( const T1& x ) : first( x ), second( T2() ) {}
00163         _Element_data( const T1& x, T2& y ) : first( x ), second( y ) {}
00164         template < typename U, typename V >
00165         _Element_data( const _Element_data< U, V >& p ) : first( p.first ), second( p.second )
00166         {
00167         }
00168     };  // Element_data
00169 
00170     template < typename Entities, typename Iterator >
00171     void assign_entities( Entities& entities, const Iterator& begin, const Iterator& end )
00172     {
00173         entities.reserve( std::distance( begin, end ) );
00174         for( Iterator i = begin; i != end; ++i )
00175         {
00176             entities.push_back( std::make_pair( ( *i )->second.first, ( *i )->first ) );
00177         }
00178     }
00179 
00180     template < typename Coordinate, typename Coordinate_iterator >
00181     void update_bounding_max( Coordinate& max, Coordinate_iterator j )
00182     {
00183         typedef typename Coordinate::iterator Iterator;
00184         for( Iterator i = max.begin(); i != max.end(); ++i, ++j )
00185         {
00186             *i = std::max( *i, *j );
00187         }
00188     }
00189 
00190     template < typename Coordinate, typename Coordinate_iterator >
00191     void update_bounding_min( Coordinate& min, Coordinate_iterator j )
00192     {
00193         typedef typename Coordinate::iterator Iterator;
00194         for( Iterator i = min.begin(); i != min.end(); ++i, ++j )
00195         {
00196             *i = std::min( *i, *j );
00197         }
00198     }
00199 
00200     template < typename Box >
00201     void update_bounding_box( Box& a, const Box& b )
00202     {
00203         update_bounding_max( a.max, b.max.begin() );
00204         update_bounding_min( a.min, b.min.begin() );
00205 #ifdef COMMON_TREE_DEBUG
00206         if( !box_contains_box( a, b ) )
00207         {
00208             std::cout << a << b << std::endl;
00209         }
00210 #endif
00211     }
00212 
00213     template < typename Entity_map, typename Ordering >
00214     void construct_ordering( Entity_map& entity_map, Ordering& entity_ordering )
00215     {
00216         entity_ordering.reserve( entity_map.size() );
00217         typedef typename Entity_map::iterator Map_iterator;
00218         for( Map_iterator i = entity_map.begin(); i != entity_map.end(); ++i )
00219         {
00220             entity_ordering.push_back( i );
00221         }
00222     }
00223 
00224     // Input: A bunch of entity handles
00225     // Output: A map from handle -> Data
00226     // Requirements: Data contains at least a bounding box.
00227     // And a non-default constructor which takes only a Box&
00228     template < typename Entity_handles, typename Element_map, typename Bounding_box, typename Moab >
00229     void construct_element_map( const Entity_handles& elements,
00230                                 Element_map& map,
00231                                 Bounding_box& bounding_box,
00232                                 Moab& moab )
00233     {
00234         typedef typename Element_map::mapped_type Box_data;
00235         typedef typename Entity_handles::value_type Entity_handle;
00236         typedef typename Entity_handles::iterator Entity_handles_iterator;
00237         typedef typename Box_data::first_type::value_type Unit;
00238         typedef typename std::vector< Unit > Coordinates;
00239         typedef typename Coordinates::iterator Coordinate_iterator;
00240 
00241         for( Entity_handles_iterator i = elements.begin(); i != elements.end(); ++i )
00242         {
00243             // TODO: not generic enough. Why dim != 3
00244             const int DIM    = 3;
00245             int num_vertices = 0;
00246             // Commence un-necessary deep copying.
00247             const Entity_handle* vertex_handle;
00248             moab.get_connectivity( *i, vertex_handle, num_vertices );
00249             Coordinates coordinate( DIM * num_vertices, 0.0 );
00250             moab.get_coords( vertex_handle, num_vertices, &coordinate[0] );
00251             Bounding_box box( coordinate.begin(), coordinate.begin() + 3 );
00252             if( i == elements.begin() )
00253             {
00254                 bounding_box = box;
00255             }
00256             for( Coordinate_iterator j = coordinate.begin() + DIM; j != coordinate.end(); j += DIM )
00257             {
00258                 update_bounding_max( box.max, j );
00259                 update_bounding_min( box.min, j );
00260             }
00261             update_bounding_box( bounding_box, box );
00262             map.insert( std::make_pair( *i, Box_data( box ) ) );
00263         }
00264     }
00265 
00266 }  // namespace common_tree
00267 
00268 }  // namespace moab
00269 
00270 #endif  // COMMON_TREE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines