Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
obb_analysis.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 
5 #include "moab/Core.hpp"
6 #include "moab/Interface.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 
17 {
18 
19  public:
20  int count;
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,
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 
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 
235 
236  double tot_vol;
237 
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 
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 }