 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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;
122 Box bounding_box;
123 double split;
124 double left_line;
125 double right_line;
126 double right_leftline;
127 double left_rightline;
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;
215 int &left_, middle_, right_;
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 >
231 class Element_tree
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:
244 typedef Element_tree< _Entity_handles, _Box, _Moab, _Parametrizer > Self;
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
249 typedef common_tree::_Element_data< Box, std::bitset< NUM_DIM * MAX_ITERATIONS * 2 > > Element_data;
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
287 Element_tree( Self& s )
288 : entity_handles_( s.entity_handles_ ), tree_( s.tree_ ), moab( s.moab ), bounding_box( s.bounding_box )
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:
603 const Entity_handles& entity_handles_;
604 Nodes tree_;
605 Moab& moab;
606 Box bounding_box;
607 Parametrizer entity_contains;
608
609 }; // class Element_tree
610
611 } // namespace moab
612
613 #endif // ELEMENT_TREE_HPP
1.9.1.