 |
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 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4
5 #include "moab/Core.hpp"
6 #include "moab/Interface.hpp"
7 #include "moab/OrientedBoxTreeTool.hpp"
8 #include "moab/CartVect.hpp"
9 #include "moab/OrientedBox.hpp"
10 #include "moab/GeomTopoTool.hpp"
11
12 #include "cgm2moab.hpp"
13
14 using namespace moab;
15
16 class TriCounter : public OrientedBoxTreeTool::Op
17 {
18
19 public:
20 int count;
21 Interface* mbi;
22 OrientedBoxTreeTool* tool;
23 const CartVect& pt;
24
25 TriCounter( Interface* mbi_p, OrientedBoxTreeTool* tool_p, const CartVect& p )
26 : OrientedBoxTreeTool::Op(), count( 0 ), mbi( mbi_p ), tool( tool_p ), pt( p )
27 {
28 }
29
30 virtual ErrorCode visit( EntityHandle node, int, bool& descend )
31 {
32 OrientedBox box;
33 ErrorCode rval = tool->box( node, box );
34
35 descend = box.contained( pt, 1e-6 );
36
37 return rval;
38 }
39
40 virtual ErrorCode leaf( EntityHandle node )
41 {
42
43 int numtris;
44 ErrorCode rval = tool->get_moab_instance()->get_number_entities_by_type( node, MBTRI, numtris );
45 count += numtris;
46 return rval;
47 }
48 };
49
50 ErrorCode obbvis_create( GeomTopoTool& gtt, std::vector< int >& volumes, int grid, std::string& filename )
51 {
52 OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
53
54 CartVect min, max;
55 EntityHandle vol = gtt.entity_by_id( 3, volumes.front() );
56 double middle[3];
57 double axis1[3], axis2[3], axis3[3];
58 double minPt[3], maxPt[3];
59 ErrorCode rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
60 // compute min and max verticies
61 for( int i = 0; i < 3; i++ )
62 {
63 double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
64 minPt[i] = middle[i] - sum;
65 maxPt[i] = middle[i] + sum;
66 }
67 CHECKERR( gtt, rval );
68
69 /* Compute an axis-aligned bounding box of all the requested volumes */
70 for( std::vector< int >::iterator i = volumes.begin() + 1; i != volumes.end(); ++i )
71 {
72 CartVect i_min, i_max;
73 vol = gtt.entity_by_id( 3, *i );
74 rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
75 // compute min and max verticies
76 for( int j = 0; j < 3; j++ )
77 {
78 double sum = fabs( axis1[j] ) + fabs( axis2[j] ) + fabs( axis3[j] );
79 minPt[j] = middle[j] - sum;
80 maxPt[j] = middle[j] + sum;
81 }
82
83 for( int j = 0; j < 3; ++j )
84 {
85 min[j] = std::min( min[j], i_min[j] );
86 max[j] = std::max( max[j], i_max[j] );
87 }
88 }
89
90 // These vectors could be repurposed to describe an OBB without changing the loops below
91 CartVect center( middle );
92 CartVect v1( axis1 );
93 CartVect v2( axis2 );
94 CartVect v3( axis3 );
95
96 /* Compute the vertices of the visualization grid. Calculation points are at the center
97 of each cell in this grid, so make grid+1 vertices in each direction. */
98 int numpoints = pow( (double)( grid + 1 ), 3 );
99 double* pgrid = new double[numpoints * 3];
100 int idx = 0;
101
102 for( int i = 0; i < numpoints * 3; ++i )
103 pgrid[i] = 0.0;
104
105 for( int i = 0; i <= grid; ++i )
106 {
107 CartVect x = -v1 + ( ( v1 * 2.0 ) * ( i / (double)grid ) );
108
109 for( int j = 0; j <= grid; ++j )
110 {
111 CartVect y = -v2 + ( ( v2 * 2.0 ) * ( j / (double)grid ) );
112
113 for( int k = 0; k <= grid; ++k )
114 {
115 CartVect z = -v3 + ( ( v3 * 2.0 ) * ( k / (double)grid ) );
116
117 CartVect p = center + x + y + z;
118 for( int d = 0; d < 3; ++d )
119 {
120 pgrid[idx++] = p[d];
121 }
122 }
123 }
124 }
125
126 /* Create a new MOAB to use for output, and build the vertex grid */
127 Core mb2;
128 Range r;
129 rval = mb2.create_vertices( pgrid, numpoints, r );
130 CHECKERR( mb2, rval );
131
132 Tag lttag;
133 rval = mb2.tag_get_handle( "LEAFTRIS", sizeof( int ), MB_TYPE_INTEGER, lttag,
134 MB_TAG_EXCL | MB_TAG_CREAT | MB_TAG_BYTES | MB_TAG_DENSE, 0 );
135 CHECKERR( mb2, rval );
136
137 int row = grid + 1;
138 int side = row * row;
139 EntityHandle connect[8];
140 EntityHandle hex;
141
142 // offset from grid corner to grid center
143 CartVect grid_hex_center_offset = ( v1 + v2 + v3 ) * 2 * ( 1.0 / grid );
144
145 // collect all the surfaces from the requested volumes to iterate over --
146 // this prevents checking a shared surface more than once.
147 Range surfs;
148 for( std::vector< int >::iterator it = volumes.begin(); it != volumes.end(); ++it )
149 {
150
151 vol = gtt.entity_by_id( 3, *it );
152 Range it_surfs;
153 rval = gtt.get_moab_instance()->get_child_meshsets( vol, it_surfs );
154 CHECKERR( gtt, rval );
155 surfs.merge( it_surfs );
156 }
157 std::cout << "visualizing " << surfs.size() << " surfaces." << std::endl;
158
159 /* Build hexes for any point with 1 or more leaftris */
160 for( int i = 0; i < grid; ++i )
161 {
162 for( int j = 0; j < grid; ++j )
163 {
164 for( int k = 0; k < grid; ++k )
165 {
166
167 idx = ( side * k ) + ( row * j ) + i;
168 assert( idx + 1 + row + side > numpoints - 1 );
169
170 CartVect loc = CartVect( ( pgrid + ( idx * 3 ) ) ) + grid_hex_center_offset;
171 TriCounter tc( gtt.get_moab_instance(), &obbtool, loc );
172
173 for( Range::iterator it = surfs.begin(); it != surfs.end(); ++it )
174 {
175
176 EntityHandle surf_tree;
177 rval = gtt.get_root( *it, surf_tree );
178 CHECKERR( gtt, rval );
179
180 rval = obbtool.preorder_traverse( surf_tree, tc );
181 CHECKERR( gtt, rval );
182 }
183
184 if( tc.count == 0 ) continue;
185
186 connect[0] = r[idx];
187 connect[1] = r[idx + 1];
188 connect[2] = r[idx + 1 + row];
189 connect[3] = r[idx + row];
190 connect[4] = r[idx + side];
191 connect[5] = r[idx + 1 + side];
192 connect[6] = r[idx + 1 + row + side];
193 connect[7] = r[idx + row + side];
194
195 rval = mb2.create_element( MBHEX, connect, 8, hex );
196 CHECKERR( mb2, rval );
197
198 rval = mb2.tag_set_data( lttag, &hex, 1, &( tc.count ) );
199 CHECKERR( mb2, rval );
200 }
201 }
202 }
203
204 if( verbose )
205 {
206 std::cout << "Writing " << filename << std::endl;
207 }
208 rval = mb2.write_file( filename.c_str() );
209 CHECKERR( mb2, rval );
210
211 return rval;
212 }
213
214 // stats code borrowed from OrientedBoxTreeTool.cpp
215 static inline double std_dev( double sqr, double sum, double count )
216 {
217 sum /= count;
218 sqr /= count;
219 return sqrt( sqr - sum * sum );
220 }
221
222 class TriStats : public OrientedBoxTreeTool::Op
223 {
224
225 public:
226 unsigned min, max, sum, leaves;
227 double sqr;
228
229 const static unsigned ten_buckets_max = 5;
230 unsigned ten_buckets[ten_buckets_max];
231 double ten_buckets_vol[ten_buckets_max];
232
233 Interface* mbi;
234 OrientedBoxTreeTool* tool;
235
236 double tot_vol;
237
238 TriStats( Interface* mbi_p, OrientedBoxTreeTool* tool_p, EntityHandle root )
239 : OrientedBoxTreeTool::Op(), sum( 0 ), leaves( 0 ), sqr( 0 ), mbi( mbi_p ), tool( tool_p )
240 {
241 min = std::numeric_limits< unsigned >::max();
242 max = std::numeric_limits< unsigned >::min();
243
244 for( unsigned i = 0; i < ten_buckets_max; ++i )
245 {
246 ten_buckets[i] = 0;
247 ten_buckets_vol[i] = 0.;
248 }
249
250 OrientedBox box;
251 ErrorCode rval = tool->box( root, box );
252 CHECKERR( mbi, rval );
253 tot_vol = box.volume();
254 }
255
256 virtual ErrorCode visit( EntityHandle, int, bool& descend )
257 {
258
259 descend = true;
260 return MB_SUCCESS;
261 }
262
263 virtual ErrorCode leaf( EntityHandle node )
264 {
265
266 Range tris;
267 ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
268 unsigned count = tris.size();
269
270 sum += count;
271 sqr += ( count * count );
272 if( min > count ) min = count;
273 if( max < count ) max = count;
274
275 for( unsigned i = 0; i < ten_buckets_max; ++i )
276 {
277 if( count > std::pow( (double)10, (int)( i + 1 ) ) )
278 {
279 ten_buckets[i] += 1;
280 OrientedBox box;
281 rval = tool->box( node, box );
282 CHECKERR( mbi, rval );
283 ten_buckets_vol[i] += box.volume();
284 }
285 }
286
287 leaves++;
288
289 return rval;
290 }
291
292 std::string commafy( int num )
293 {
294 std::stringstream str;
295 str << num;
296 std::string s = str.str();
297
298 int n = s.size();
299 for( int i = n - 3; i >= 1; i -= 3 )
300 {
301 s.insert( i, 1, ',' );
302 n++;
303 }
304
305 return s;
306 }
307
308 void write_results( std::ostream& out )
309 {
310
311 out << commafy( sum ) << " triangles in " << commafy( leaves ) << " leaves." << std::endl;
312
313 double avg = sum / (double)leaves;
314 double stddev = std_dev( sqr, sum, leaves );
315
316 out << "Tris per leaf: Min " << min << ", Max " << max << ", avg " << avg << ", stddev " << stddev << std::endl;
317
318 for( unsigned i = 0; i < ten_buckets_max; ++i )
319 {
320 if( ten_buckets[i] )
321 {
322 out << "Leaves exceeding " << std::pow( (double)10, (int)( i + 1 ) )
323 << " triangles: " << ten_buckets[i];
324
325 double frac_total_vol = ten_buckets_vol[i] / tot_vol;
326 double avg_ftv = frac_total_vol / ten_buckets[i];
327
328 out << " (avg " << avg_ftv * 100.0 << "% of OBB volume)" << std::endl;
329 }
330 }
331 }
332 };
333
334 /* no longer supported
335 static std::string make_property_string( DagMC& dag, EntityHandle eh, std::vector<std::string>
336 &properties )
337 {
338 ErrorCode ret;
339 std::string propstring;
340 for( std::vector<std::string>::iterator p = properties.begin();
341 p != properties.end(); ++p )
342 {
343 if( dag.has_prop( eh, *p ) ){
344 std::vector< std::string> vals;
345 ret = dag.prop_values( eh, *p, vals );
346 CHECKERR(dag,ret);
347 propstring += *p;
348 if( vals.size() == 1 ){
349 propstring += "=";
350 propstring += vals[0];
351 }
352 else if( vals.size() > 1 ){
353 // this property has multiple values, list within brackets
354 propstring += "=[";
355 for( std::vector<std::string>::iterator i = vals.begin();
356 i != vals.end(); ++i )
357 {
358 propstring += *i;
359 propstring += ",";
360 }
361 // replace the last trailing comma with a close braket
362 propstring[ propstring.length()-1 ] = ']';
363 }
364 propstring += ", ";
365 }
366 }
367 if( propstring.length() ){
368 propstring.resize( propstring.length() - 2 ); // drop trailing comma
369 }
370 return propstring;
371 }
372 */
373
374 ErrorCode obbstat_write( GeomTopoTool& gtt,
375 std::vector< int >& volumes,
376 std::vector< std::string >& /*properties*/,
377 std::ostream& out )
378 {
379
380 ErrorCode ret = MB_SUCCESS;
381 OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
382
383 // can assume that volume numbers are valid.
384 for( std::vector< int >::iterator i = volumes.begin(); i != volumes.end(); ++i )
385 {
386 EntityHandle vol_root;
387 EntityHandle vol = gtt.entity_by_id( 3, *i );
388 CHECKERR( gtt, ret );
389
390 if( vol == 0 )
391 {
392 std::cerr << "ERROR: volume " << *i << " has no entity." << std::endl;
393 continue;
394 }
395
396 ret = gtt.get_root( vol, vol_root );
397 CHECKERR( gtt, ret );
398
399 out << "\nVolume " << *i << " " << std::flush;
400
401 if( gtt.is_implicit_complement( vol ) ) out << "(implicit complement) ";
402 out << std::endl;
403
404 // get all surfaces in volume
405 Range surfs;
406 ret = gtt.get_moab_instance()->get_child_meshsets( vol, surfs );
407 CHECKERR( gtt, ret );
408
409 out << " with " << surfs.size() << " surfaces" << std::endl;
410
411 TriStats ts( gtt.get_moab_instance(), &obbtool, vol_root );
412 ret = obbtool.preorder_traverse( vol_root, ts );
413 CHECKERR( gtt, ret );
414 ts.write_results( out );
415
416 if( verbose )
417 {
418 out << "Surface list: " << std::flush;
419 for( Range::iterator j = surfs.begin(); j != surfs.end(); ++j )
420 {
421 out << gtt.global_id( *j );
422 if( j + 1 != surfs.end() ) out << ",";
423 }
424 out << std::endl;
425 ret = obbtool.stats( vol_root, out );
426 CHECKERR( gtt, ret );
427 }
428
429 out << "\n ------------ " << std::endl;
430 }
431
432 return ret;
433 }
1.9.1.