 |
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 * common_tree.hpp
3 * Ryan H. Lewis
4 * (C) 2012
5 * Functionality common to all trees.
6 */
7 #include <vector>
8 #include <set>
9 #include <iostream>
10 #include <map>
11 #include <algorithm>
12 #include <bitset>
13 #include <numeric>
14 #include <cmath>
15 #include <tr1/unordered_map>
16 #include <limits>
17
18 #ifndef COMMON_TREE_HPP
19 #define COMMON_TREE_HPP
20 #define NUM_DIM 3
21 #define TREE_DEBUG
22 namespace moab
23 {
24 namespace common_tree
25 {
26
27 template < typename T, typename Stream >
28 void print_vector( const T& v, Stream& out )
29 {
30 typedef typename T::const_iterator Iterator;
31 out << "[ ";
32 for( Iterator i = v.begin(); i != v.end(); ++i )
33 {
34 out << *i;
35 if( i + 1 != v.end() )
36 {
37 out << ", ";
38 }
39 }
40 out << " ]" << std::endl;
41 }
42
43 #ifdef TREE_DEBUG
44 template < typename T >
45 void print_vector( const T& begin, const T& end )
46 {
47 std::cout << "[ ";
48 for( T i = begin; i != end; ++i )
49 {
50 std::cout << ( *i )->second.second.to_ulong();
51 if( i + 1 != end )
52 {
53 std::cout << ", ";
54 }
55 }
56 std::cout << " ]" << std::endl;
57 }
58 #endif
59
60 template < typename _Box, typename _Point >
61 bool box_contains_point( const _Box& box, const _Point& p, const double tol )
62 {
63 for( std::size_t i = 0; i < box.min.size(); ++i )
64 {
65 if( p[i] < ( box.min[i] - tol ) || p[i] > ( box.max[i] ) + tol )
66 {
67 return false;
68 }
69 }
70 return true;
71 }
72
73 template < typename _Box >
74 bool box_contains_box( const _Box& a, const _Box& b, const double tol )
75 {
76 for( std::size_t i = 0; i < a.min.size(); ++i )
77 {
78 if( b.min[i] < ( a.min[i] - tol ) )
79 {
80 return false;
81 }
82 if( b.max[i] > ( a.max[i] + tol ) )
83 {
84 return false;
85 }
86 }
87 return true;
88 }
89
90 namespace
91 {
92 template < typename T >
93 struct Compute_center
94 {
95 // deprecation of binary_function
96 typedef T first_argument_type;
97 typedef T second_argument_type;
98 typedef T result_type;
99
100 T operator()( const T a, const T b ) const
101 {
102 return ( a + b ) / 2.0;
103 }
104 }; // Compute_center
105 } // namespace
106
107 template < typename Vector >
108 inline void compute_box_center( Vector& max, Vector& min, Vector& center )
109 {
110 typedef typename Vector::value_type Unit;
111 center = min;
112 std::transform( max.begin(), max.end(), center.begin(), center.begin(), Compute_center< Unit >() );
113 }
114
115 template < typename Box >
116 inline typename Box::value_type compute_box_center( const Box& box, const int dim )
117 {
118 return ( box.max[dim] + box.min[dim] ) / 2.0;
119 }
120
121 template < typename T = float >
122 class Box
123 {
124 public:
125 typedef T value_type;
126 typedef std::vector< T > Vector;
127 Box() : max( 3, 0.0 ), min( 3, 0.0 ) {}
128 Box( const Box& from ) : max( from.max ), min( from.min ) {}
129 template < typename Iterator >
130 Box( const Iterator begin, const Iterator end ) : max( begin, end ), min( begin, end )
131 {
132 }
133 Box& operator=( const Box& from )
134 {
135 max = from.max;
136 min = from.min;
137 return *this;
138 }
139 Vector max;
140 Vector min;
141 }; // Box
142
143 template < typename T >
144 std::ostream& operator<<( std::ostream& out, const Box< T >& box )
145 {
146 out << "Max: ";
147 print_vector( box.max, out );
148 out << "Min: ";
149 print_vector( box.min, out );
150 return out;
151 }
152
153 // essentially a pair, but with an added constructor.
154 template < typename T1, typename T2 >
155 struct _Element_data
156 {
157 typedef T1 first_type;
158 typedef T2 second_type;
159 T1 first;
160 T2 second;
161 _Element_data() : first( T1() ), second( T2() ) {}
162 _Element_data( const T1& x ) : first( x ), second( T2() ) {}
163 _Element_data( const T1& x, T2& y ) : first( x ), second( y ) {}
164 template < typename U, typename V >
165 _Element_data( const _Element_data< U, V >& p ) : first( p.first ), second( p.second )
166 {
167 }
168 }; // Element_data
169
170 template < typename Entities, typename Iterator >
171 void assign_entities( Entities& entities, const Iterator& begin, const Iterator& end )
172 {
173 entities.reserve( std::distance( begin, end ) );
174 for( Iterator i = begin; i != end; ++i )
175 {
176 entities.push_back( std::make_pair( ( *i )->second.first, ( *i )->first ) );
177 }
178 }
179
180 template < typename Coordinate, typename Coordinate_iterator >
181 void update_bounding_max( Coordinate& max, Coordinate_iterator j )
182 {
183 typedef typename Coordinate::iterator Iterator;
184 for( Iterator i = max.begin(); i != max.end(); ++i, ++j )
185 {
186 *i = std::max( *i, *j );
187 }
188 }
189
190 template < typename Coordinate, typename Coordinate_iterator >
191 void update_bounding_min( Coordinate& min, Coordinate_iterator j )
192 {
193 typedef typename Coordinate::iterator Iterator;
194 for( Iterator i = min.begin(); i != min.end(); ++i, ++j )
195 {
196 *i = std::min( *i, *j );
197 }
198 }
199
200 template < typename Box >
201 void update_bounding_box( Box& a, const Box& b )
202 {
203 update_bounding_max( a.max, b.max.begin() );
204 update_bounding_min( a.min, b.min.begin() );
205 #ifdef COMMON_TREE_DEBUG
206 if( !box_contains_box( a, b ) )
207 {
208 std::cout << a << b << std::endl;
209 }
210 #endif
211 }
212
213 template < typename Entity_map, typename Ordering >
214 void construct_ordering( Entity_map& entity_map, Ordering& entity_ordering )
215 {
216 entity_ordering.reserve( entity_map.size() );
217 typedef typename Entity_map::iterator Map_iterator;
218 for( Map_iterator i = entity_map.begin(); i != entity_map.end(); ++i )
219 {
220 entity_ordering.push_back( i );
221 }
222 }
223
224 // Input: A bunch of entity handles
225 // Output: A map from handle -> Data
226 // Requirements: Data contains at least a bounding box.
227 // And a non-default constructor which takes only a Box&
228 template < typename Entity_handles, typename Element_map, typename Bounding_box, typename Moab >
229 void construct_element_map( const Entity_handles& elements,
230 Element_map& map,
231 Bounding_box& bounding_box,
232 Moab& moab )
233 {
234 typedef typename Element_map::mapped_type Box_data;
235 typedef typename Entity_handles::value_type Entity_handle;
236 typedef typename Entity_handles::iterator Entity_handles_iterator;
237 typedef typename Box_data::first_type::value_type Unit;
238 typedef typename std::vector< Unit > Coordinates;
239 typedef typename Coordinates::iterator Coordinate_iterator;
240
241 for( Entity_handles_iterator i = elements.begin(); i != elements.end(); ++i )
242 {
243 // TODO: not generic enough. Why dim != 3
244 const int DIM = 3;
245 int num_vertices = 0;
246 // Commence un-necessary deep copying.
247 const Entity_handle* vertex_handle;
248 moab.get_connectivity( *i, vertex_handle, num_vertices );
249 Coordinates coordinate( DIM * num_vertices, 0.0 );
250 moab.get_coords( vertex_handle, num_vertices, &coordinate[0] );
251 Bounding_box box( coordinate.begin(), coordinate.begin() + 3 );
252 if( i == elements.begin() )
253 {
254 bounding_box = box;
255 }
256 for( Coordinate_iterator j = coordinate.begin() + DIM; j != coordinate.end(); j += DIM )
257 {
258 update_bounding_max( box.max, j );
259 update_bounding_min( box.min, j );
260 }
261 update_bounding_box( bounding_box, box );
262 map.insert( std::make_pair( *i, Box_data( box ) ) );
263 }
264 }
265
266 } // namespace common_tree
267
268 } // namespace moab
269
270 #endif // COMMON_TREE_HPP
1.9.1.