Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
element_tree.hpp
Go to the documentation of this file.
00001 /**
00002  * element_tree.hpp
00003  * Ryan H. Lewis
00004  * (C) 2012
00005  *
00006  * An element tree partitions a mesh composed of elements.
00007  * We subdivide the bounding box of a mesh, and each element is
00008  * either entirely on the left, entirely on the right, or crossing
00009  * the diving line. We build a tree on the mesh with this property.
00010  */
00011 #include <vector>
00012 #include <set>
00013 #include <iostream>
00014 #include <map>
00015 #include <algorithm>
00016 #include <bitset>
00017 #include <numeric>
00018 #include <cmath>
00019 #include <tr1/unordered_map>
00020 #include <limits>
00021 
00022 #include "common_tree.hpp"
00023 
00024 #ifndef ELEMENT_TREE_HPP
00025 #define ELEMENT_TREE_HPP
00026 namespace moab
00027 {
00028 // forward declarations
00029 
00030 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
00031 class Element_tree;
00032 
00033 // non-exported functionality
00034 namespace
00035 {
00036     namespace _element_tree
00037     {
00038         template < typename Iterator >
00039         struct Iterator_comparator
00040         {
00041             typedef typename Iterator::value_type Value;
00042             bool operator()( const Value& a, const Value& b )
00043             {
00044                 return a->second.second.to_ulong() < b->second.second.to_ulong();
00045             }
00046         };  // Iterator_comparator
00047 
00048         template < typename Data >
00049         struct Split_comparator
00050         {
00051             // we minimizes ||left| - |right|| + |middle|^2
00052             double split_objective( const Data& a ) const
00053             {
00054                 if( a.second.sizes[2] == 0 || a.second.sizes[0] == 0 )
00055                 {
00056                     return std::numeric_limits< std::size_t >::max();
00057                 }
00058                 const double total = a.second.sizes[0] + a.second.sizes[2];
00059                 const int max      = 2 * ( a.second.sizes[2] > a.second.sizes[0] );
00060 
00061                 return ( a.second.sizes[max] - a.second.sizes[2 * ( 1 - ( max == 2 ) )] ) / total;
00062             }
00063             bool operator()( const Data& a, const Data& b ) const
00064             {
00065                 return split_objective( a ) < split_objective( b );
00066             }
00067         };  // Split_comparator
00068 
00069         template < typename Partition_data, typename Box >
00070         void correct_bounding_box( const Partition_data& data, Box& box, const int child )
00071         {
00072             const int dim = data.dim;
00073             switch( child )
00074             {
00075                 case 0:
00076                     box.max[dim] = data.left_rightline;
00077                     break;
00078                 case 1:
00079                     box.max[dim] = data.right_line;
00080                     box.min[dim] = data.left_line;
00081                     break;
00082                 case 2:
00083                     box.min[dim] = data.right_leftline;
00084                     break;
00085             }
00086 #ifdef ELEMENT_TREE_DEBUG
00087             print_vector( data.bounding_box.max );
00088             print_vector( data.bounding_box.min );
00089             print_vector( box.max );
00090             print_vector( box.min );
00091 #endif
00092         }
00093 
00094         template < typename Box >
00095         struct _Partition_data
00096         {
00097             typedef _Partition_data< Box > Self;
00098             // default constructor
00099             _Partition_data() : sizes( 3, 0 ), dim( 0 ) {}
00100             _Partition_data( const Self& f )
00101             {
00102                 *this = f;
00103             }
00104             _Partition_data( const Box& _box, int _dim )
00105                 : sizes( 3, 0 ), bounding_box( _box ), split( ( _box.max[_dim] + _box.min[_dim] ) / 2.0 ),
00106                   left_line( split ), right_line( split ), dim( _dim )
00107             {
00108             }
00109             _Partition_data& operator=( const Self& f )
00110             {
00111                 sizes          = f.sizes;
00112                 bounding_box   = f.bounding_box;
00113                 split          = f.split;
00114                 left_line      = f.left_line;
00115                 right_line     = f.right_line;
00116                 right_leftline = f.right_leftline;
00117                 left_rightline = f.left_rightline;
00118                 dim            = f.dim;
00119                 return *this;
00120             }
00121             std::vector< std::size_t > sizes;
00122             Box bounding_box;
00123             double split;
00124             double left_line;
00125             double right_line;
00126             double right_leftline;
00127             double left_rightline;
00128             int dim;
00129             std::size_t left() const
00130             {
00131                 return sizes[0];
00132             }
00133             std::size_t middle() const
00134             {
00135                 return sizes[1];
00136             }
00137             std::size_t right() const
00138             {
00139                 return sizes[2];
00140             }
00141         };  // Partition_data
00142 
00143         template < typename _Entity_handles, typename _Entities >
00144         class _Node
00145         {
00146             // public types:
00147           public:
00148             typedef _Entity_handles Entity_handles;
00149             typedef _Entities Entities;
00150 
00151             // private types:
00152           private:
00153             typedef _Node< _Entity_handles, _Entities > Self;
00154 
00155             // Constructors
00156           public:
00157             // Default constructor
00158             _Node()
00159                 : children( 3, -1 ), left_( children[0] ), middle_( children[1] ), right_( children[2] ), dim( -1 ),
00160                   split( 0 ), left_line( 0 ), right_line( 0 ), entities( 0 )
00161             {
00162             }
00163 
00164             // Copy constructor
00165             _Node( const Self& from )
00166                 : children( from.children ), left_( children[0] ), middle_( children[1] ), right_( children[2] ),
00167                   dim( from.dim ), split( from.split ), left_line( from.left_line ), right_line( from.right_line ),
00168                   entities( from.entities )
00169             {
00170             }
00171 
00172           public:
00173             template < typename Iterator >
00174             void assign_entities( const Iterator& begin, const Iterator& end )
00175             {
00176                 entities.reserve( std::distance( begin, end ) );
00177                 for( Iterator i = begin; i != end; ++i )
00178                 {
00179                     entities.push_back( std::make_pair( ( *i )->second.first, ( *i )->first ) );
00180                 }
00181             }
00182 
00183             // Functionality
00184           public:
00185             bool leaf() const
00186             {
00187                 return children[0] == -1 && children[1] == -1 && children[2] == -1;
00188             }
00189             Self& operator=( const Self& from )
00190             {
00191                 children   = from.children;
00192                 dim        = from.dim;
00193                 left_      = from.left_;
00194                 middle_    = from.middle_;
00195                 right_     = from.right_;
00196                 split      = from.split;
00197                 left_line  = from.left_line;
00198                 right_line = from.right_line;
00199                 entities   = from.entities;
00200                 return *this;
00201             }
00202             template < typename Box >
00203             Self& operator=( const _Partition_data< Box >& from )
00204             {
00205                 dim        = from.dim;
00206                 split      = from.split;
00207                 left_line  = from.left_line;
00208                 right_line = from.right_line;
00209                 return *this;
00210             }
00211             // private data members:
00212           private:
00213             // indices of children
00214             std::vector< int > children;
00215             int &left_, middle_, right_;
00216             int dim;       // split dimension
00217             double split;  // split position
00218             double left_line;
00219             double right_line;
00220             Entities entities;
00221 
00222             // Element_tree can touch my privates.
00223             template < Entity_handles, typename B >
00224             friend class moab::Element_tree;
00225         };  // class Node
00226 
00227     }  // namespace _element_tree
00228 }  // namespace
00229 
00230 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
00231 class Element_tree
00232 {
00233 
00234     // public types
00235   public:
00236     typedef _Entity_handles Entity_handles;
00237     typedef _Box Box;
00238     typedef _Moab Moab;
00239     typedef _Parametrizer Parametrizer;
00240     typedef typename Entity_handles::value_type Entity_handle;
00241 
00242     // private types
00243   private:
00244     typedef Element_tree< _Entity_handles, _Box, _Moab, _Parametrizer > Self;
00245     typedef std::pair< Box, Entity_handle > Leaf_element;
00246     typedef _element_tree::_Node< Entity_handles, std::vector< Leaf_element > > Node;
00247 // int is because we only need to store
00248 #define MAX_ITERATIONS 2
00249     typedef common_tree::_Element_data< Box, std::bitset< NUM_DIM * MAX_ITERATIONS * 2 > > Element_data;
00250     typedef std::vector< Node > Nodes;
00251     // TODO: we really want an unordered map here, make sure this is kosher..
00252     typedef std::tr1::unordered_map< Entity_handle, Element_data > Element_map;
00253     typedef typename std::vector< typename Element_map::iterator > Element_list;
00254     typedef _element_tree::_Partition_data< Box > Partition_data;
00255     // public methods
00256   public:
00257     // Constructor
00258     Element_tree( Entity_handles& _entities, Moab& _moab, Box& _bounding_box, Parametrizer& _entity_contains )
00259         : entity_handles_( _entities ), tree_(), moab( _moab ), bounding_box( _bounding_box ),
00260           entity_contains( _entity_contains )
00261     {
00262         tree_.reserve( _entities.size() );
00263         Element_map element_map( _entities.size() );
00264         Partition_data _data;
00265         common_tree::construct_element_map( entity_handles_, element_map, _data.bounding_box, moab );
00266         bounding_box  = _data.bounding_box;
00267         _bounding_box = bounding_box;
00268         Element_list element_ordering( element_map.size() );
00269         std::size_t index = 0;
00270         for( typename Element_map::iterator i = element_map.begin(); i != element_map.end(); ++i, ++index )
00271         {
00272             element_ordering[index] = i;
00273         }
00274         // We only build nonempty trees
00275         if( element_ordering.size() )
00276         {
00277             // initially all bits are set
00278             std::bitset< 3 > directions( 7 );
00279             tree_.push_back( Node() );
00280             int depth = 0;
00281             build_tree( element_ordering.begin(), element_ordering.end(), 0, directions, _data, depth );
00282             std::cout << "depth: " << depth << std::endl;
00283         }
00284     }
00285 
00286     // Copy constructor
00287     Element_tree( Self& s )
00288         : entity_handles_( s.entity_handles_ ), tree_( s.tree_ ), moab( s.moab ), bounding_box( s.bounding_box )
00289     {
00290     }
00291 
00292     // private functionality
00293   private:
00294     template < typename Iterator, typename Split_data >
00295     void compute_split( Iterator& begin, Iterator& end, Split_data& split_data, bool iteration = false )
00296     {
00297         typedef typename Iterator::value_type::value_type Map_value_type;
00298         typedef typename Map_value_type::second_type::second_type Bitset;
00299         // we will update the left/right line
00300         double& left_line  = split_data.left_line;
00301         double& right_line = split_data.right_line;
00302         double& split      = split_data.split;
00303         const int& dim     = split_data.dim;
00304 #ifdef ELEMENT_TREE_DEBUG
00305         std::cout << std::endl;
00306         std::cout << "-------------------" << std::endl;
00307         std::cout << "dim: " << dim << " split: " << split << std::endl;
00308         std::cout << "bounding_box min: ";
00309         print_vector( split_data.bounding_box.min );
00310         std::cout << "bounding_box max: ";
00311         print_vector( split_data.bounding_box.max );
00312 #endif
00313         // for each elt determine if left/middle/right
00314         for( Iterator i = begin; i != end; ++i )
00315         {
00316             const Box& box = ( *i )->second.first;
00317             Bitset& bits   = ( *i )->second.second;
00318             // will be 0 if on left, will be 1 if in the middle
00319             // and 2 if on the right;
00320             const bool on_left   = ( box.max[dim] < split );
00321             const bool on_right  = ( box.min[dim] > split );
00322             const bool in_middle = !on_left && !on_right;
00323             // set the corresponding bits in the bit vector
00324             // looks like: [x_1 = 00 | x_2 = 00 | .. | z_1 = 00 | z_2 = 00]
00325             // two bits, left = 00, middle = 01, right = 10
00326             const int index = 4 * dim + 2 * iteration;
00327             if( on_left )
00328             {
00329                 split_data.sizes[0]++;
00330             }
00331             else if( in_middle )
00332             {
00333                 split_data.sizes[1]++;
00334                 bits.set( index, 1 );
00335                 left_line  = std::min( left_line, box.min[dim] );
00336                 right_line = std::max( right_line, box.max[dim] );
00337             }
00338             else if( on_right )
00339             {
00340                 bits.set( index + 1, 1 );
00341                 split_data.sizes[2]++;
00342             }
00343         }
00344 #ifdef ELEMENT_TREE_DEBUG
00345         std::size_t _count = std::accumulate( split_data.sizes.begin(), split_data.sizes.end(), 0 );
00346         std::size_t total  = std::distance( begin, end );
00347         if( total != _count )
00348         {
00349             std::cout << total << "vs. " << _count << std::endl;
00350         }
00351         std::cout << " left_line: " << left_line;
00352         std::cout << " right_line: " << right_line << std::endl;
00353         std::cout << "co/mputed partition size: ";
00354         print_vector( split_data.sizes );
00355         std::cout << "-------------------" << std::endl;
00356 #endif
00357     }
00358 
00359     template < typename Split_data >
00360     bool update_split_line( Split_data& data ) const
00361     {
00362         const int max        = 2 * ( data.sizes[2] > data.sizes[0] );
00363         const int min        = 2 * ( 1 - ( max == 2 ) );
00364         bool one_side_empty  = data.sizes[max] == 0 || data.sizes[min] == 0;
00365         double balance_ratio = data.sizes[max] - data.sizes[min];
00366         // if ( !one_side_empty && balance_ratio < .05*total){ return false; }
00367         if( !one_side_empty )
00368         {
00369             // if we have some imbalance on left/right
00370             // try to fix the situation
00371             balance_ratio /= data.sizes[max];
00372             data.split += ( max - 1 ) * balance_ratio * ( data.split / 2.0 );
00373         }
00374         else
00375         {
00376             // if the (left) side is empty move the split line just past the
00377             // extent of the (left) line of the middle box.
00378             // if middle encompasses everything then wiggle
00379             // the split line a bit and hope for the best..
00380             const double left_distance  = std::abs( data.left_line - data.split );
00381             const double right_distance = std::abs( data.right_line - data.split );
00382             if( ( data.sizes[0] == 0 ) && ( data.sizes[2] != 0 ) )
00383             {
00384                 data.split += right_distance;
00385             }
00386             else if( data.sizes[2] == 0 && data.sizes[0] != 0 )
00387             {
00388                 data.split -= left_distance;
00389             }
00390             else
00391             {
00392                 data.split *= 1.05;
00393             }
00394         }
00395         data.left_line = data.right_line = data.split;
00396         data.sizes.assign( data.sizes.size(), 0 );
00397         return true;
00398     }
00399 
00400     template < typename Iterator, typename Split_data, typename Directions >
00401     void determine_split( Iterator& begin, Iterator& end, Split_data& data, const Directions& directions )
00402     {
00403         typedef typename Iterator::value_type Pair;
00404         typedef typename Pair::value_type Map_value_type;
00405         typedef typename Map_value_type::second_type::second_type Bitset;
00406         typedef typename Map_value_type::second_type::first_type Box;
00407         typedef typename std::map< std::size_t, Split_data > Splits;
00408         typedef typename Splits::value_type Split;
00409         typedef _element_tree::Split_comparator< Split > Comparator;
00410         Splits splits;
00411         for( std::size_t dir = 0; dir < directions.size(); ++dir )
00412         {
00413             if( directions.test( dir ) )
00414             {
00415                 Split_data split_data( data.bounding_box, dir );
00416                 compute_split( begin, end, split_data );
00417                 splits.insert( std::make_pair( 2 * dir, split_data ) );
00418                 if( update_split_line( split_data ) )
00419                 {
00420                     compute_split( begin, end, split_data, true );
00421                     splits.insert( std::make_pair( 2 * dir + 1, split_data ) );
00422                 }
00423             }
00424         }
00425         Split best = *std::min_element( splits.begin(), splits.end(), Comparator() );
00426 #ifdef ELEMENT_TREE_DEBUG
00427         std::cout << "best: " << Comparator().split_objective( best ) << " ";
00428         print_vector( best.second.sizes );
00429 #endif
00430         const int dir          = best.first / 2;
00431         const int iter         = best.first % 2;
00432         double& left_rightline = best.second.left_rightline = best.second.bounding_box.min[dir];
00433         double& right_leftline = best.second.right_leftline = best.second.bounding_box.max[dir];
00434         Bitset mask( 0 );
00435         mask.flip( 4 * dir + 2 * iter ).flip( 4 * dir + 2 * iter + 1 );
00436         for( Iterator i = begin; i != end; ++i )
00437         {
00438             Bitset& bits   = ( *i )->second.second;
00439             const Box& box = ( *i )->second.first;
00440             // replace 12 bits with just two.
00441             bits &= mask;
00442             bits >>= 4 * dir + 2 * iter;
00443             // if box is labeled left/right but properly contained
00444             // in the middle, move the element into the middle.
00445             // we can shrink the size of left/right
00446             switch( bits.to_ulong() )
00447             {
00448                 case 0:
00449                     if( box.max[dir] > best.second.left_line )
00450                     {
00451                         left_rightline = std::max( left_rightline, box.max[dir] );
00452                     }
00453                     break;
00454                 case 2:
00455                     if( box.min[dir] < best.second.right_line )
00456                     {
00457                         right_leftline = std::min( right_leftline, box.max[dir] );
00458                     }
00459                     break;
00460             }
00461         }
00462         data = best.second;
00463     }
00464 
00465 // define here for now.
00466 #define ELEMENTS_PER_LEAF 30
00467 #define MAX_DEPTH         30
00468 #define EPSILON           1e-1
00469     template < typename Iterator, typename Node_index, typename Directions, typename Partition_data >
00470     void build_tree( Iterator begin,
00471                      Iterator end,
00472                      const Node_index node,
00473                      const Directions& directions,
00474                      Partition_data& _data,
00475                      int& depth,
00476                      const bool is_middle = false )
00477     {
00478         std::size_t number_elements = std::distance( begin, end );
00479         if( depth < MAX_DEPTH && number_elements > ELEMENTS_PER_LEAF && ( !is_middle || directions.any() ) )
00480         {
00481             determine_split( begin, end, _data, directions );
00482             // count_sort( begin, end, _data);
00483             std::sort( begin, end, _element_tree::Iterator_comparator< Iterator >() );
00484             // update the tree
00485             tree_[node] = _data;
00486             Iterator middle_begin( begin + _data.left() );
00487             Iterator middle_end( middle_begin + _data.middle() );
00488             std::vector< int > depths( 3, depth );
00489             // left subtree
00490             if( _data.left() > 0 )
00491             {
00492                 Partition_data data( _data );
00493                 tree_.push_back( Node() );
00494                 tree_[node].children[0] = tree_.size() - 1;
00495                 correct_bounding_box( _data, data.bounding_box, 0 );
00496                 Directions new_directions( directions );
00497                 const bool axis_is_very_small =
00498                     ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
00499                 new_directions.set( _data.dim, axis_is_very_small );
00500                 build_tree( begin, middle_begin, tree_[node].children[0], new_directions, data, ++depths[0],
00501                             is_middle );
00502             }
00503             // middle subtree
00504             if( _data.middle() > 0 )
00505             {
00506                 Partition_data data( _data );
00507                 tree_.push_back( Node() );
00508                 tree_[node].children[1] = tree_.size() - 1;
00509                 correct_bounding_box( _data, data.bounding_box, 1 );
00510                 // force the middle subtree to split
00511                 // in a different direction from this one
00512                 Directions new_directions( directions );
00513                 new_directions.flip( tree_[node].dim );
00514                 bool axis_is_very_small =
00515                     ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
00516                 new_directions.set( _data.dim, axis_is_very_small );
00517                 build_tree( middle_begin, middle_end, tree_[node].children[1], new_directions, data, ++depths[1],
00518                             true );
00519             }
00520             // right subtree
00521             if( _data.right() > 0 )
00522             {
00523                 Partition_data data( _data );
00524                 tree_.push_back( Node() );
00525                 tree_[node].children[2] = tree_.size() - 1;
00526                 correct_bounding_box( _data, data.bounding_box, 2 );
00527                 Directions new_directions( directions );
00528                 const bool axis_is_very_small =
00529                     ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
00530                 new_directions.set( _data.dim, axis_is_very_small );
00531 
00532                 build_tree( middle_end, end, tree_[node].children[2], directions, data, ++depths[2], is_middle );
00533             }
00534             depth = *std::max_element( depths.begin(), depths.end() );
00535         }
00536         if( tree_[node].leaf() )
00537         {
00538             common_tree::assign_entities( tree_[node].entities, begin, end );
00539         }
00540     }
00541 
00542     template < typename Vector, typename Node_index, typename Result >
00543     Result& _find_point( const Vector& point, const Node_index& index, Result& result ) const
00544     {
00545         typedef typename Node::Entities::const_iterator Entity_iterator;
00546         typedef typename std::pair< bool, Vector > Return_type;
00547         const Node& node = tree_[index];
00548         if( node.leaf() )
00549         {
00550             // check each node
00551             for( Entity_iterator i = node.entities.begin(); i != node.entities.end(); ++i )
00552             {
00553                 if( common_tree::box_contains_point( i->first, point ) )
00554                 {
00555                     Return_type r = entity_contains( moab, i->second, point );
00556                     if( r.first )
00557                     {
00558                         result = std::make_pair( i->second, r.second );
00559                     }
00560                     return result;
00561                 }
00562             }
00563             return Result( 0, point );
00564         }
00565         if( point[node.dim] < node.left_line )
00566         {
00567             return _find_point( point, node.left_, result );
00568         }
00569         else if( point[node.dim] > node.right_line )
00570         {
00571             return _find_point( point, node.right_, result );
00572         }
00573         else
00574         {
00575             Entity_handle middle = _find_point( point, node.middle_, result );
00576             if( middle != 0 )
00577             {
00578                 return result;
00579             }
00580             if( point[node.dim] < node.split )
00581             {
00582                 return _find_point( point, node.left_, result );
00583             }
00584             return _find_point( point, node.right_, result );
00585         }
00586     }
00587 
00588     // public functionality
00589   public:
00590     template < typename Vector, typename Result >
00591     Result& find( const Vector& point, Result& result ) const
00592     {
00593         typedef typename Vector::const_iterator Point_iterator;
00594         typedef typename Box::Pair Pair;
00595         typedef typename Pair::first_type Box_iterator;
00596         return _find_point( point, 0, result );
00597     }
00598 
00599     // public accessor methods
00600   public:
00601     // private data members
00602   private:
00603     const Entity_handles& entity_handles_;
00604     Nodes tree_;
00605     Moab& moab;
00606     Box bounding_box;
00607     Parametrizer entity_contains;
00608 
00609 };  // class Element_tree
00610 
00611 }  // namespace moab
00612 
00613 #endif  // ELEMENT_TREE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines