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