Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
element_tree.hpp
Go to the documentation of this file.
1 /**
2  * element_tree.hpp
3  * Ryan H. Lewis
4  * (C) 2012
5  *
6  * An element tree partitions a mesh composed of elements.
7  * We subdivide the bounding box of a mesh, and each element is
8  * either entirely on the left, entirely on the right, or crossing
9  * the diving line. We build a tree on the mesh with this property.
10  */
11 #include <vector>
12 #include <set>
13 #include <iostream>
14 #include <map>
15 #include <algorithm>
16 #include <bitset>
17 #include <numeric>
18 #include <cmath>
19 #include <tr1/unordered_map>
20 #include <limits>
21 
22 #include "common_tree.hpp"
23 
24 #ifndef ELEMENT_TREE_HPP
25 #define ELEMENT_TREE_HPP
26 namespace moab
27 {
28 // forward declarations
29 
30 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
31 class Element_tree;
32 
33 // non-exported functionality
34 namespace
35 {
36  namespace _element_tree
37  {
38  template < typename Iterator >
39  struct Iterator_comparator
40  {
41  typedef typename Iterator::value_type Value;
42  bool operator()( const Value& a, const Value& b )
43  {
44  return a->second.second.to_ulong() < b->second.second.to_ulong();
45  }
46  }; // Iterator_comparator
47 
48  template < typename Data >
49  struct Split_comparator
50  {
51  // we minimizes ||left| - |right|| + |middle|^2
52  double split_objective( const Data& a ) const
53  {
54  if( a.second.sizes[2] == 0 || a.second.sizes[0] == 0 )
55  {
56  return std::numeric_limits< std::size_t >::max();
57  }
58  const double total = a.second.sizes[0] + a.second.sizes[2];
59  const int max = 2 * ( a.second.sizes[2] > a.second.sizes[0] );
60 
61  return ( a.second.sizes[max] - a.second.sizes[2 * ( 1 - ( max == 2 ) )] ) / total;
62  }
63  bool operator()( const Data& a, const Data& b ) const
64  {
65  return split_objective( a ) < split_objective( b );
66  }
67  }; // Split_comparator
68 
69  template < typename Partition_data, typename Box >
70  void correct_bounding_box( const Partition_data& data, Box& box, const int child )
71  {
72  const int dim = data.dim;
73  switch( child )
74  {
75  case 0:
76  box.max[dim] = data.left_rightline;
77  break;
78  case 1:
79  box.max[dim] = data.right_line;
80  box.min[dim] = data.left_line;
81  break;
82  case 2:
83  box.min[dim] = data.right_leftline;
84  break;
85  }
86 #ifdef ELEMENT_TREE_DEBUG
87  print_vector( data.bounding_box.max );
88  print_vector( data.bounding_box.min );
89  print_vector( box.max );
90  print_vector( box.min );
91 #endif
92  }
93 
94  template < typename Box >
95  struct _Partition_data
96  {
97  typedef _Partition_data< Box > Self;
98  // default constructor
99  _Partition_data() : sizes( 3, 0 ), dim( 0 ) {}
100  _Partition_data( const Self& f )
101  {
102  *this = f;
103  }
104  _Partition_data( const Box& _box, int _dim )
105  : sizes( 3, 0 ), bounding_box( _box ), split( ( _box.max[_dim] + _box.min[_dim] ) / 2.0 ),
106  left_line( split ), right_line( split ), dim( _dim )
107  {
108  }
109  _Partition_data& operator=( const Self& f )
110  {
111  sizes = f.sizes;
112  bounding_box = f.bounding_box;
113  split = f.split;
114  left_line = f.left_line;
115  right_line = f.right_line;
116  right_leftline = f.right_leftline;
117  left_rightline = f.left_rightline;
118  dim = f.dim;
119  return *this;
120  }
121  std::vector< std::size_t > sizes;
123  double split;
124  double left_line;
125  double right_line;
128  int dim;
129  std::size_t left() const
130  {
131  return sizes[0];
132  }
133  std::size_t middle() const
134  {
135  return sizes[1];
136  }
137  std::size_t right() const
138  {
139  return sizes[2];
140  }
141  }; // Partition_data
142 
143  template < typename _Entity_handles, typename _Entities >
144  class _Node
145  {
146  // public types:
147  public:
148  typedef _Entity_handles Entity_handles;
149  typedef _Entities Entities;
150 
151  // private types:
152  private:
153  typedef _Node< _Entity_handles, _Entities > Self;
154 
155  // Constructors
156  public:
157  // Default constructor
158  _Node()
159  : children( 3, -1 ), left_( children[0] ), middle_( children[1] ), right_( children[2] ), dim( -1 ),
160  split( 0 ), left_line( 0 ), right_line( 0 ), entities( 0 )
161  {
162  }
163 
164  // Copy constructor
165  _Node( const Self& from )
166  : children( from.children ), left_( children[0] ), middle_( children[1] ), right_( children[2] ),
167  dim( from.dim ), split( from.split ), left_line( from.left_line ), right_line( from.right_line ),
168  entities( from.entities )
169  {
170  }
171 
172  public:
173  template < typename Iterator >
174  void assign_entities( const Iterator& begin, const Iterator& end )
175  {
176  entities.reserve( std::distance( begin, end ) );
177  for( Iterator i = begin; i != end; ++i )
178  {
179  entities.push_back( std::make_pair( ( *i )->second.first, ( *i )->first ) );
180  }
181  }
182 
183  // Functionality
184  public:
185  bool leaf() const
186  {
187  return children[0] == -1 && children[1] == -1 && children[2] == -1;
188  }
189  Self& operator=( const Self& from )
190  {
191  children = from.children;
192  dim = from.dim;
193  left_ = from.left_;
194  middle_ = from.middle_;
195  right_ = from.right_;
196  split = from.split;
197  left_line = from.left_line;
198  right_line = from.right_line;
199  entities = from.entities;
200  return *this;
201  }
202  template < typename Box >
203  Self& operator=( const _Partition_data< Box >& from )
204  {
205  dim = from.dim;
206  split = from.split;
207  left_line = from.left_line;
208  right_line = from.right_line;
209  return *this;
210  }
211  // private data members:
212  private:
213  // indices of children
214  std::vector< int > children;
216  int dim; // split dimension
217  double split; // split position
218  double left_line;
219  double right_line;
220  Entities entities;
221 
222  // Element_tree can touch my privates.
223  template < Entity_handles, typename B >
224  friend class moab::Element_tree;
225  }; // class Node
226 
227  } // namespace _element_tree
228 } // namespace
229 
230 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
232 {
233 
234  // public types
235  public:
236  typedef _Entity_handles Entity_handles;
237  typedef _Box Box;
238  typedef _Moab Moab;
239  typedef _Parametrizer Parametrizer;
240  typedef typename Entity_handles::value_type Entity_handle;
241 
242  // private types
243  private:
245  typedef std::pair< Box, Entity_handle > Leaf_element;
246  typedef _element_tree::_Node< Entity_handles, std::vector< Leaf_element > > Node;
247 // int is because we only need to store
248 #define MAX_ITERATIONS 2
250  typedef std::vector< Node > Nodes;
251  // TODO: we really want an unordered map here, make sure this is kosher..
252  typedef std::tr1::unordered_map< Entity_handle, Element_data > Element_map;
253  typedef typename std::vector< typename Element_map::iterator > Element_list;
254  typedef _element_tree::_Partition_data< Box > Partition_data;
255  // public methods
256  public:
257  // Constructor
258  Element_tree( Entity_handles& _entities, Moab& _moab, Box& _bounding_box, Parametrizer& _entity_contains )
259  : entity_handles_( _entities ), tree_(), moab( _moab ), bounding_box( _bounding_box ),
260  entity_contains( _entity_contains )
261  {
262  tree_.reserve( _entities.size() );
263  Element_map element_map( _entities.size() );
264  Partition_data _data;
265  common_tree::construct_element_map( entity_handles_, element_map, _data.bounding_box, moab );
266  bounding_box = _data.bounding_box;
267  _bounding_box = bounding_box;
268  Element_list element_ordering( element_map.size() );
269  std::size_t index = 0;
270  for( typename Element_map::iterator i = element_map.begin(); i != element_map.end(); ++i, ++index )
271  {
272  element_ordering[index] = i;
273  }
274  // We only build nonempty trees
275  if( element_ordering.size() )
276  {
277  // initially all bits are set
278  std::bitset< 3 > directions( 7 );
279  tree_.push_back( Node() );
280  int depth = 0;
281  build_tree( element_ordering.begin(), element_ordering.end(), 0, directions, _data, depth );
282  std::cout << "depth: " << depth << std::endl;
283  }
284  }
285 
286  // Copy constructor
289  {
290  }
291 
292  // private functionality
293  private:
294  template < typename Iterator, typename Split_data >
295  void compute_split( Iterator& begin, Iterator& end, Split_data& split_data, bool iteration = false )
296  {
297  typedef typename Iterator::value_type::value_type Map_value_type;
298  typedef typename Map_value_type::second_type::second_type Bitset;
299  // we will update the left/right line
300  double& left_line = split_data.left_line;
301  double& right_line = split_data.right_line;
302  double& split = split_data.split;
303  const int& dim = split_data.dim;
304 #ifdef ELEMENT_TREE_DEBUG
305  std::cout << std::endl;
306  std::cout << "-------------------" << std::endl;
307  std::cout << "dim: " << dim << " split: " << split << std::endl;
308  std::cout << "bounding_box min: ";
309  print_vector( split_data.bounding_box.min );
310  std::cout << "bounding_box max: ";
311  print_vector( split_data.bounding_box.max );
312 #endif
313  // for each elt determine if left/middle/right
314  for( Iterator i = begin; i != end; ++i )
315  {
316  const Box& box = ( *i )->second.first;
317  Bitset& bits = ( *i )->second.second;
318  // will be 0 if on left, will be 1 if in the middle
319  // and 2 if on the right;
320  const bool on_left = ( box.max[dim] < split );
321  const bool on_right = ( box.min[dim] > split );
322  const bool in_middle = !on_left && !on_right;
323  // set the corresponding bits in the bit vector
324  // looks like: [x_1 = 00 | x_2 = 00 | .. | z_1 = 00 | z_2 = 00]
325  // two bits, left = 00, middle = 01, right = 10
326  const int index = 4 * dim + 2 * iteration;
327  if( on_left )
328  {
329  split_data.sizes[0]++;
330  }
331  else if( in_middle )
332  {
333  split_data.sizes[1]++;
334  bits.set( index, 1 );
335  left_line = std::min( left_line, box.min[dim] );
336  right_line = std::max( right_line, box.max[dim] );
337  }
338  else if( on_right )
339  {
340  bits.set( index + 1, 1 );
341  split_data.sizes[2]++;
342  }
343  }
344 #ifdef ELEMENT_TREE_DEBUG
345  std::size_t _count = std::accumulate( split_data.sizes.begin(), split_data.sizes.end(), 0 );
346  std::size_t total = std::distance( begin, end );
347  if( total != _count )
348  {
349  std::cout << total << "vs. " << _count << std::endl;
350  }
351  std::cout << " left_line: " << left_line;
352  std::cout << " right_line: " << right_line << std::endl;
353  std::cout << "co/mputed partition size: ";
354  print_vector( split_data.sizes );
355  std::cout << "-------------------" << std::endl;
356 #endif
357  }
358 
359  template < typename Split_data >
360  bool update_split_line( Split_data& data ) const
361  {
362  const int max = 2 * ( data.sizes[2] > data.sizes[0] );
363  const int min = 2 * ( 1 - ( max == 2 ) );
364  bool one_side_empty = data.sizes[max] == 0 || data.sizes[min] == 0;
365  double balance_ratio = data.sizes[max] - data.sizes[min];
366  // if ( !one_side_empty && balance_ratio < .05*total){ return false; }
367  if( !one_side_empty )
368  {
369  // if we have some imbalance on left/right
370  // try to fix the situation
371  balance_ratio /= data.sizes[max];
372  data.split += ( max - 1 ) * balance_ratio * ( data.split / 2.0 );
373  }
374  else
375  {
376  // if the (left) side is empty move the split line just past the
377  // extent of the (left) line of the middle box.
378  // if middle encompasses everything then wiggle
379  // the split line a bit and hope for the best..
380  const double left_distance = std::abs( data.left_line - data.split );
381  const double right_distance = std::abs( data.right_line - data.split );
382  if( ( data.sizes[0] == 0 ) && ( data.sizes[2] != 0 ) )
383  {
384  data.split += right_distance;
385  }
386  else if( data.sizes[2] == 0 && data.sizes[0] != 0 )
387  {
388  data.split -= left_distance;
389  }
390  else
391  {
392  data.split *= 1.05;
393  }
394  }
395  data.left_line = data.right_line = data.split;
396  data.sizes.assign( data.sizes.size(), 0 );
397  return true;
398  }
399 
400  template < typename Iterator, typename Split_data, typename Directions >
401  void determine_split( Iterator& begin, Iterator& end, Split_data& data, const Directions& directions )
402  {
403  typedef typename Iterator::value_type Pair;
404  typedef typename Pair::value_type Map_value_type;
405  typedef typename Map_value_type::second_type::second_type Bitset;
406  typedef typename Map_value_type::second_type::first_type Box;
407  typedef typename std::map< std::size_t, Split_data > Splits;
408  typedef typename Splits::value_type Split;
409  typedef _element_tree::Split_comparator< Split > Comparator;
410  Splits splits;
411  for( std::size_t dir = 0; dir < directions.size(); ++dir )
412  {
413  if( directions.test( dir ) )
414  {
415  Split_data split_data( data.bounding_box, dir );
416  compute_split( begin, end, split_data );
417  splits.insert( std::make_pair( 2 * dir, split_data ) );
418  if( update_split_line( split_data ) )
419  {
420  compute_split( begin, end, split_data, true );
421  splits.insert( std::make_pair( 2 * dir + 1, split_data ) );
422  }
423  }
424  }
425  Split best = *std::min_element( splits.begin(), splits.end(), Comparator() );
426 #ifdef ELEMENT_TREE_DEBUG
427  std::cout << "best: " << Comparator().split_objective( best ) << " ";
428  print_vector( best.second.sizes );
429 #endif
430  const int dir = best.first / 2;
431  const int iter = best.first % 2;
432  double& left_rightline = best.second.left_rightline = best.second.bounding_box.min[dir];
433  double& right_leftline = best.second.right_leftline = best.second.bounding_box.max[dir];
434  Bitset mask( 0 );
435  mask.flip( 4 * dir + 2 * iter ).flip( 4 * dir + 2 * iter + 1 );
436  for( Iterator i = begin; i != end; ++i )
437  {
438  Bitset& bits = ( *i )->second.second;
439  const Box& box = ( *i )->second.first;
440  // replace 12 bits with just two.
441  bits &= mask;
442  bits >>= 4 * dir + 2 * iter;
443  // if box is labeled left/right but properly contained
444  // in the middle, move the element into the middle.
445  // we can shrink the size of left/right
446  switch( bits.to_ulong() )
447  {
448  case 0:
449  if( box.max[dir] > best.second.left_line )
450  {
451  left_rightline = std::max( left_rightline, box.max[dir] );
452  }
453  break;
454  case 2:
455  if( box.min[dir] < best.second.right_line )
456  {
457  right_leftline = std::min( right_leftline, box.max[dir] );
458  }
459  break;
460  }
461  }
462  data = best.second;
463  }
464 
465 // define here for now.
466 #define ELEMENTS_PER_LEAF 30
467 #define MAX_DEPTH 30
468 #define EPSILON 1e-1
469  template < typename Iterator, typename Node_index, typename Directions, typename Partition_data >
470  void build_tree( Iterator begin,
471  Iterator end,
472  const Node_index node,
473  const Directions& directions,
474  Partition_data& _data,
475  int& depth,
476  const bool is_middle = false )
477  {
478  std::size_t number_elements = std::distance( begin, end );
479  if( depth < MAX_DEPTH && number_elements > ELEMENTS_PER_LEAF && ( !is_middle || directions.any() ) )
480  {
481  determine_split( begin, end, _data, directions );
482  // count_sort( begin, end, _data);
483  std::sort( begin, end, _element_tree::Iterator_comparator< Iterator >() );
484  // update the tree
485  tree_[node] = _data;
486  Iterator middle_begin( begin + _data.left() );
487  Iterator middle_end( middle_begin + _data.middle() );
488  std::vector< int > depths( 3, depth );
489  // left subtree
490  if( _data.left() > 0 )
491  {
492  Partition_data data( _data );
493  tree_.push_back( Node() );
494  tree_[node].children[0] = tree_.size() - 1;
495  correct_bounding_box( _data, data.bounding_box, 0 );
496  Directions new_directions( directions );
497  const bool axis_is_very_small =
498  ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
499  new_directions.set( _data.dim, axis_is_very_small );
500  build_tree( begin, middle_begin, tree_[node].children[0], new_directions, data, ++depths[0],
501  is_middle );
502  }
503  // middle subtree
504  if( _data.middle() > 0 )
505  {
506  Partition_data data( _data );
507  tree_.push_back( Node() );
508  tree_[node].children[1] = tree_.size() - 1;
509  correct_bounding_box( _data, data.bounding_box, 1 );
510  // force the middle subtree to split
511  // in a different direction from this one
512  Directions new_directions( directions );
513  new_directions.flip( tree_[node].dim );
514  bool axis_is_very_small =
515  ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
516  new_directions.set( _data.dim, axis_is_very_small );
517  build_tree( middle_begin, middle_end, tree_[node].children[1], new_directions, data, ++depths[1],
518  true );
519  }
520  // right subtree
521  if( _data.right() > 0 )
522  {
523  Partition_data data( _data );
524  tree_.push_back( Node() );
525  tree_[node].children[2] = tree_.size() - 1;
526  correct_bounding_box( _data, data.bounding_box, 2 );
527  Directions new_directions( directions );
528  const bool axis_is_very_small =
529  ( data.bounding_box.max[_data.dim] - data.bounding_box.min[_data.dim] < EPSILON );
530  new_directions.set( _data.dim, axis_is_very_small );
531 
532  build_tree( middle_end, end, tree_[node].children[2], directions, data, ++depths[2], is_middle );
533  }
534  depth = *std::max_element( depths.begin(), depths.end() );
535  }
536  if( tree_[node].leaf() )
537  {
538  common_tree::assign_entities( tree_[node].entities, begin, end );
539  }
540  }
541 
542  template < typename Vector, typename Node_index, typename Result >
543  Result& _find_point( const Vector& point, const Node_index& index, Result& result ) const
544  {
545  typedef typename Node::Entities::const_iterator Entity_iterator;
546  typedef typename std::pair< bool, Vector > Return_type;
547  const Node& node = tree_[index];
548  if( node.leaf() )
549  {
550  // check each node
551  for( Entity_iterator i = node.entities.begin(); i != node.entities.end(); ++i )
552  {
553  if( common_tree::box_contains_point( i->first, point ) )
554  {
555  Return_type r = entity_contains( moab, i->second, point );
556  if( r.first )
557  {
558  result = std::make_pair( i->second, r.second );
559  }
560  return result;
561  }
562  }
563  return Result( 0, point );
564  }
565  if( point[node.dim] < node.left_line )
566  {
567  return _find_point( point, node.left_, result );
568  }
569  else if( point[node.dim] > node.right_line )
570  {
571  return _find_point( point, node.right_, result );
572  }
573  else
574  {
575  Entity_handle middle = _find_point( point, node.middle_, result );
576  if( middle != 0 )
577  {
578  return result;
579  }
580  if( point[node.dim] < node.split )
581  {
582  return _find_point( point, node.left_, result );
583  }
584  return _find_point( point, node.right_, result );
585  }
586  }
587 
588  // public functionality
589  public:
590  template < typename Vector, typename Result >
591  Result& find( const Vector& point, Result& result ) const
592  {
593  typedef typename Vector::const_iterator Point_iterator;
594  typedef typename Box::Pair Pair;
595  typedef typename Pair::first_type Box_iterator;
596  return _find_point( point, 0, result );
597  }
598 
599  // public accessor methods
600  public:
601  // private data members
602  private:
608 
609 }; // class Element_tree
610 
611 } // namespace moab
612 
613 #endif // ELEMENT_TREE_HPP