 |
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 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 /**\file OrientedBox.hpp
17 *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
18 *\date 2006-07-18
19 */
20
21 #include "moab/Interface.hpp"
22 #include "Internals.hpp"
23 #include "moab/OrientedBoxTreeTool.hpp"
24 #include "moab/Range.hpp"
25 #include "moab/CN.hpp"
26 #include "moab/GeomUtil.hpp"
27 #include "MBTagConventions.hpp"
28 #include <iostream>
29 #include <iomanip>
30 #include <algorithm>
31 #include <limits>
32 #include <cassert>
33 #include <cmath>
34
35 //#define MB_OBB_USE_VECTOR_QUERIES
36 //#define MB_OBB_USE_TYPE_QUERIES
37
38 namespace moab
39 {
40
41 #if defined( MB_OBB_USE_VECTOR_QUERIES ) && defined( MB_OBB_USE_TYPE_QUERIES )
42 #undef MB_OBB_USE_TYPE_QUERIES
43 #endif
44
45 const char DEFAULT_TAG_NAME[] = "OBB";
46
47 OrientedBoxTreeTool::Op::~Op() {}
48
49 OrientedBoxTreeTool::OrientedBoxTreeTool( Interface* i, const char* tag_name, bool destroy_created_trees )
50 : instance( i ), cleanUpTrees( destroy_created_trees )
51 {
52 if( !tag_name ) tag_name = DEFAULT_TAG_NAME;
53 ErrorCode rval = OrientedBox::tag_handle( tagHandle, instance, tag_name );
54 if( MB_SUCCESS != rval ) tagHandle = 0;
55 }
56
57 OrientedBoxTreeTool::~OrientedBoxTreeTool()
58 {
59 if( !cleanUpTrees ) return;
60
61 while( !createdTrees.empty() )
62 {
63 EntityHandle tree = createdTrees.back();
64 // make sure this is a tree (rather than some other, stale handle)
65 const void* data_ptr = 0;
66 ErrorCode rval = instance->tag_get_by_ptr( tagHandle, &tree, 1, &data_ptr );
67 if( MB_SUCCESS == rval ) rval = delete_tree( tree );
68 if( MB_SUCCESS != rval ) createdTrees.pop_back();
69 }
70 }
71
72 OrientedBoxTreeTool::Settings::Settings()
73 : max_leaf_entities( 8 ), max_depth( 0 ), worst_split_ratio( 0.95 ), best_split_ratio( 0.4 ),
74 set_options( MESHSET_SET )
75 {
76 }
77
78 bool OrientedBoxTreeTool::Settings::valid() const
79 {
80 return max_leaf_entities > 0 && max_depth >= 0 && worst_split_ratio <= 1.0 && best_split_ratio >= 0.0 &&
81 worst_split_ratio >= best_split_ratio;
82 }
83
84 /********************** Simple Tree Access Methods ****************************/
85
86 ErrorCode OrientedBoxTreeTool::box( EntityHandle set, OrientedBox& obb )
87 {
88 return instance->tag_get_data( tagHandle, &set, 1, &obb );
89 }
90
91 ErrorCode OrientedBoxTreeTool::box( EntityHandle set,
92 double center[3],
93 double axis1[3],
94 double axis2[3],
95 double axis3[3] )
96 {
97 OrientedBox obb;
98 ErrorCode rval = this->box( set, obb );
99 obb.center.get( center );
100 obb.scaled_axis( 0 ).get( axis1 );
101 obb.scaled_axis( 1 ).get( axis2 );
102 obb.scaled_axis( 2 ).get( axis3 );
103 return rval;
104 }
105
106 /********************** Tree Construction Code ****************************/
107
108 struct OrientedBoxTreeTool::SetData
109 {
110 EntityHandle handle;
111 OrientedBox::CovarienceData box_data;
112 // Range vertices;
113 };
114
115 ErrorCode OrientedBoxTreeTool::build( const Range& entities, EntityHandle& set_handle_out, const Settings* settings )
116 {
117 if( !entities.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
118 if( settings && !settings->valid() ) return MB_FAILURE;
119
120 return build_tree( entities, set_handle_out, 0, settings ? *settings : Settings() );
121 }
122
123 ErrorCode OrientedBoxTreeTool::join_trees( const Range& sets, EntityHandle& set_handle_out, const Settings* settings )
124 {
125 if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
126 if( settings && !settings->valid() ) return MB_FAILURE;
127
128 // Build initial set data list.
129 std::list< SetData > data;
130 for( Range::const_iterator i = sets.begin(); i != sets.end(); ++i )
131 {
132 Range elements;
133 ErrorCode rval = instance->get_entities_by_dimension( *i, 2, elements, true );
134 if( MB_SUCCESS != rval ) return rval;
135 if( elements.empty() ) continue;
136
137 data.push_back( SetData() );
138 SetData& set_data = data.back();
139 set_data.handle = *i;
140 rval = OrientedBox::covariance_data_from_tris( set_data.box_data, instance, elements );
141 if( MB_SUCCESS != rval ) return rval;
142 }
143
144 ErrorCode result = build_sets( data, set_handle_out, 0, settings ? *settings : Settings() );
145 if( MB_SUCCESS != result ) return result;
146
147 for( Range::reverse_iterator i = sets.rbegin(); i != sets.rend(); ++i )
148 {
149 createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), *i ), createdTrees.end() );
150 }
151 createdTrees.push_back( set_handle_out );
152 return MB_SUCCESS;
153 }
154
155 /**\brief Split triangles by which side of a plane they are on
156 *
157 * Given a plane specified as a bisecting plane normal to one
158 * of the axes of a box, split triangles based on which side
159 * of the plane they are on.
160 *\param instance MOAB instance
161 *\param box The oriented box containing all the entities
162 *\param axis The axis for which the split plane is orthogonal
163 *\param entities The entities intersecting plane
164 *\param left_list Output, entities to the left of the plane
165 *\param right_list Output, entities to the right of the plane
166 */
167 static ErrorCode split_box( Interface* instance,
168 const OrientedBox& box,
169 int axis,
170 const Range& entities,
171 Range& left_list,
172 Range& right_list )
173 {
174 ErrorCode rval;
175 left_list.clear();
176 right_list.clear();
177
178 std::vector< CartVect > coords;
179 for( Range::reverse_iterator i = entities.rbegin(); i != entities.rend(); ++i )
180 {
181 const EntityHandle* conn = NULL;
182 int conn_len = 0;
183 rval = instance->get_connectivity( *i, conn, conn_len );
184 if( MB_SUCCESS != rval ) return rval;
185
186 coords.resize( conn_len );
187 rval = instance->get_coords( conn, conn_len, coords[0].array() );
188 if( MB_SUCCESS != rval ) return rval;
189
190 CartVect centroid( 0.0 );
191 for( int j = 0; j < conn_len; ++j )
192 centroid += coords[j];
193 centroid /= conn_len;
194
195 if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
196 left_list.insert( *i );
197 else
198 right_list.insert( *i );
199 }
200
201 return MB_SUCCESS;
202 }
203
204 ErrorCode OrientedBoxTreeTool::build_tree( const Range& entities,
205 EntityHandle& set,
206 int depth,
207 const Settings& settings )
208 {
209 OrientedBox tmp_box;
210 ErrorCode rval;
211
212 if( entities.empty() )
213 {
214 Matrix3 axis;
215 tmp_box = OrientedBox( axis, CartVect( 0. ) );
216 }
217 else
218 {
219 rval = OrientedBox::compute_from_2d_cells( tmp_box, instance, entities );
220 if( MB_SUCCESS != rval ) return rval;
221 }
222
223 // create an entity set for the tree node
224 rval = instance->create_meshset( settings.set_options, set );
225 if( MB_SUCCESS != rval ) return rval;
226
227 rval = instance->tag_set_data( tagHandle, &set, 1, &tmp_box );
228 if( MB_SUCCESS != rval )
229 {
230 delete_tree( set );
231 return rval;
232 }
233
234 // check if should create children
235 bool leaf = true;
236 ++depth;
237 if( ( !settings.max_depth || depth < settings.max_depth ) &&
238 entities.size() > (unsigned)settings.max_leaf_entities )
239 {
240 // try splitting with planes normal to each axis of the box
241 // until we find an acceptable split
242 double best_ratio = settings.worst_split_ratio; // worst case ratio
243 Range best_left_list, best_right_list;
244 // Axes are sorted from shortest to longest, so search backwards
245 for( int axis = 2; best_ratio > settings.best_split_ratio && axis >= 0; --axis )
246 {
247 Range left_list, right_list;
248
249 rval = split_box( instance, tmp_box, axis, entities, left_list, right_list );
250 if( MB_SUCCESS != rval )
251 {
252 delete_tree( set );
253 return rval;
254 }
255
256 double ratio = fabs( (double)right_list.size() - left_list.size() ) / entities.size();
257
258 if( ratio < best_ratio )
259 {
260 best_ratio = ratio;
261 best_left_list.swap( left_list );
262 best_right_list.swap( right_list );
263 }
264 }
265
266 // create children
267 if( !best_left_list.empty() )
268 {
269 EntityHandle child = 0;
270
271 rval = build_tree( best_left_list, child, depth, settings );
272 if( MB_SUCCESS != rval )
273 {
274 delete_tree( set );
275 return rval;
276 }
277 rval = instance->add_child_meshset( set, child );
278 if( MB_SUCCESS != rval )
279 {
280 delete_tree( set );
281 delete_tree( child );
282 return rval;
283 }
284
285 rval = build_tree( best_right_list, child, depth, settings );
286 if( MB_SUCCESS != rval )
287 {
288 delete_tree( set );
289 return rval;
290 }
291 rval = instance->add_child_meshset( set, child );
292 if( MB_SUCCESS != rval )
293 {
294 delete_tree( set );
295 delete_tree( child );
296 return rval;
297 }
298
299 leaf = false;
300 }
301 }
302
303 if( leaf )
304 {
305 rval = instance->add_entities( set, entities );
306 if( MB_SUCCESS != rval )
307 {
308 delete_tree( set );
309 return rval;
310 }
311 }
312
313 createdTrees.push_back( set );
314 return MB_SUCCESS;
315 }
316
317 static ErrorCode split_sets( Interface*,
318 const OrientedBox& box,
319 int axis,
320 const std::list< OrientedBoxTreeTool::SetData >& sets,
321 std::list< OrientedBoxTreeTool::SetData >& left,
322 std::list< OrientedBoxTreeTool::SetData >& right )
323 {
324 left.clear();
325 right.clear();
326
327 std::list< OrientedBoxTreeTool::SetData >::const_iterator i;
328 for( i = sets.begin(); i != sets.end(); ++i )
329 {
330 CartVect centroid( i->box_data.center / i->box_data.area );
331 if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
332 left.push_back( *i );
333 else
334 right.push_back( *i );
335 }
336
337 return MB_SUCCESS;
338 }
339
340 ErrorCode OrientedBoxTreeTool::build_sets( std::list< SetData >& sets,
341 EntityHandle& node_set,
342 int depth,
343 const Settings& settings )
344 {
345 ErrorCode rval;
346 int count = sets.size();
347 if( 0 == count ) return MB_FAILURE;
348
349 // calculate box
350 OrientedBox obox;
351
352 // make vector go out of scope when done, so memory is released
353 {
354 Range elems;
355 std::vector< OrientedBox::CovarienceData > data( sets.size() );
356 data.clear();
357 for( std::list< SetData >::iterator i = sets.begin(); i != sets.end(); ++i )
358 {
359 data.push_back( i->box_data );
360 rval = instance->get_entities_by_dimension( i->handle, 2, elems, true );
361 if( MB_SUCCESS != rval ) return rval;
362 }
363
364 Range points;
365 rval = instance->get_adjacencies( elems, 0, false, points, Interface::UNION );
366 if( MB_SUCCESS != rval ) return rval;
367
368 rval = OrientedBox::compute_from_covariance_data( obox, instance, &data[0], data.size(), points );
369 if( MB_SUCCESS != rval ) return rval;
370 }
371
372 // If only one set in list...
373 if( count == 1 )
374 {
375 node_set = sets.front().handle;
376 return instance->tag_set_data( tagHandle, &node_set, 1, &obox );
377 }
378
379 // create an entity set for the tree node
380 rval = instance->create_meshset( settings.set_options, node_set );
381 if( MB_SUCCESS != rval ) return rval;
382
383 rval = instance->tag_set_data( tagHandle, &node_set, 1, &obox );
384 if( MB_SUCCESS != rval )
385 {
386 delete_tree( node_set );
387 return rval;
388 }
389
390 double best_ratio = 2.0;
391 std::list< SetData > best_left_list, best_right_list;
392 for( int axis = 0; axis < 2; ++axis )
393 {
394 std::list< SetData > left_list, right_list;
395 rval = split_sets( instance, obox, axis, sets, left_list, right_list );
396 if( MB_SUCCESS != rval )
397 {
398 delete_tree( node_set );
399 return rval;
400 }
401
402 double ratio = fabs( (double)right_list.size() - left_list.size() ) / sets.size();
403
404 if( ratio < best_ratio )
405 {
406 best_ratio = ratio;
407 best_left_list.swap( left_list );
408 best_right_list.swap( right_list );
409 }
410 }
411
412 // We must subdivide the list of sets, because we want to guarantee that
413 // there is a node in the tree corresponding to each of the sets. So if
414 // we couldn't find a usable split plane, just split them in an arbitrary
415 // manner.
416 if( best_left_list.empty() || best_right_list.empty() )
417 {
418 best_left_list.clear();
419 best_right_list.clear();
420 std::list< SetData >* lists[2] = { &best_left_list, &best_right_list };
421 int i = 0;
422 while( !sets.empty() )
423 {
424 lists[i]->push_back( sets.front() );
425 sets.pop_front();
426 i = 1 - i;
427 }
428 }
429 else
430 {
431 sets.clear(); // release memory before recursion
432 }
433
434 // Create child sets
435
436 EntityHandle child = 0;
437
438 rval = build_sets( best_left_list, child, depth + 1, settings );
439 if( MB_SUCCESS != rval )
440 {
441 delete_tree( node_set );
442 return rval;
443 }
444 rval = instance->add_child_meshset( node_set, child );
445 if( MB_SUCCESS != rval )
446 {
447 delete_tree( node_set );
448 delete_tree( child );
449 return rval;
450 }
451
452 rval = build_sets( best_right_list, child, depth + 1, settings );
453 if( MB_SUCCESS != rval )
454 {
455 delete_tree( node_set );
456 return rval;
457 }
458 rval = instance->add_child_meshset( node_set, child );
459 if( MB_SUCCESS != rval )
460 {
461 delete_tree( node_set );
462 delete_tree( child );
463 return rval;
464 }
465
466 return MB_SUCCESS;
467 }
468
469 ErrorCode OrientedBoxTreeTool::delete_tree( EntityHandle set )
470 {
471 std::vector< EntityHandle > children;
472 ErrorCode rval = instance->get_child_meshsets( set, children, 0 );
473 if( MB_SUCCESS != rval ) return rval;
474
475 createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), set ), createdTrees.end() );
476 children.insert( children.begin(), set );
477 return instance->delete_entities( &children[0], children.size() );
478 }
479
480 ErrorCode OrientedBoxTreeTool::remove_root( EntityHandle root )
481 {
482 std::vector< EntityHandle >::iterator i = find( createdTrees.begin(), createdTrees.end(), root );
483 if( i != createdTrees.end() )
484 {
485 createdTrees.erase( i );
486 return MB_SUCCESS;
487 }
488 else
489 return MB_ENTITY_NOT_FOUND;
490 }
491
492 /********************** Generic Tree Traversal ****************************/
493 struct Data
494 {
495 EntityHandle set;
496 int depth;
497 };
498 ErrorCode OrientedBoxTreeTool::preorder_traverse( EntityHandle set, Op& operation, TrvStats* accum )
499 {
500 ErrorCode rval;
501 std::vector< EntityHandle > children;
502 std::vector< Data > the_stack;
503 Data data = { set, 0 };
504 the_stack.push_back( data );
505 int max_depth = -1;
506
507 while( !the_stack.empty() )
508 {
509 data = the_stack.back();
510 the_stack.pop_back();
511
512 // increment traversal statistics
513 if( accum )
514 {
515 accum->increment( data.depth );
516 max_depth = std::max( max_depth, data.depth );
517 }
518
519 bool descend = true;
520 rval = operation.visit( data.set, data.depth, descend );
521 assert( MB_SUCCESS == rval );
522 if( MB_SUCCESS != rval ) return rval;
523
524 if( !descend ) continue;
525
526 children.clear();
527 rval = instance->get_child_meshsets( data.set, children );
528 assert( MB_SUCCESS == rval );
529 if( MB_SUCCESS != rval ) return rval;
530 if( children.empty() )
531 {
532 if( accum )
533 {
534 accum->increment_leaf( data.depth );
535 }
536 rval = operation.leaf( data.set );
537 assert( MB_SUCCESS == rval );
538 if( MB_SUCCESS != rval ) return rval;
539 }
540 else if( children.size() == 2 )
541 {
542 data.depth++;
543 data.set = children[0];
544 the_stack.push_back( data );
545 data.set = children[1];
546 the_stack.push_back( data );
547 }
548 else
549 return MB_MULTIPLE_ENTITIES_FOUND;
550 }
551
552 if( accum )
553 {
554 accum->end_traversal( max_depth );
555 }
556
557 return MB_SUCCESS;
558 }
559
560 /********************** General Sphere/Triangle Intersection ***************/
561
562 struct OBBTreeSITFrame
563 {
564 OBBTreeSITFrame( EntityHandle n, EntityHandle s, int dp ) : node( n ), surf( s ), depth( dp ) {}
565 EntityHandle node;
566 EntityHandle surf;
567 int depth;
568 };
569
570 ErrorCode OrientedBoxTreeTool::sphere_intersect_triangles( const double* center_v,
571 double radius,
572 EntityHandle tree_root,
573 std::vector< EntityHandle >& facets_out,
574 std::vector< EntityHandle >* sets_out,
575 TrvStats* accum )
576 {
577 const double radsqr = radius * radius;
578 OrientedBox b;
579 ErrorCode rval;
580 Range sets;
581 const CartVect center( center_v );
582 CartVect closest, coords[3];
583 const EntityHandle* conn;
584 int num_conn;
585 #ifndef MB_OBB_USE_VECTOR_QUERIES
586 Range tris;
587 Range::const_iterator t;
588 #else
589 std::vector< EntityHandle > tris;
590 std::vector< EntityHandle >::const_iterator t;
591 #endif
592
593 std::vector< OBBTreeSITFrame > stack;
594 std::vector< EntityHandle > children;
595 stack.reserve( 30 );
596 stack.push_back( OBBTreeSITFrame( tree_root, 0, 0 ) );
597
598 int max_depth = -1;
599
600 while( !stack.empty() )
601 {
602 EntityHandle surf = stack.back().surf;
603 EntityHandle node = stack.back().node;
604 int current_depth = stack.back().depth;
605 stack.pop_back();
606
607 // increment traversal statistics.
608 if( accum )
609 {
610 accum->increment( current_depth );
611 max_depth = std::max( max_depth, current_depth );
612 }
613
614 if( !surf && sets_out )
615 {
616 rval = get_moab_instance()->get_entities_by_type( node, MBENTITYSET, sets );
617 if( !sets.empty() ) surf = sets.front();
618 sets.clear();
619 }
620
621 // check if sphere intersects box
622 rval = box( node, b );
623 if( MB_SUCCESS != rval ) return rval;
624 b.closest_location_in_box( center, closest );
625 closest -= center;
626 if( ( closest % closest ) > radsqr ) continue;
627
628 // push child boxes on stack
629 children.clear();
630 rval = instance->get_child_meshsets( node, children );
631 if( MB_SUCCESS != rval ) return rval;
632 if( !children.empty() )
633 {
634 assert( children.size() == 2 );
635 stack.push_back( OBBTreeSITFrame( children[0], surf, current_depth + 1 ) );
636 stack.push_back( OBBTreeSITFrame( children[1], surf, current_depth + 1 ) );
637 continue;
638 }
639
640 if( accum )
641 {
642 accum->increment_leaf( current_depth );
643 }
644 // if leaf, intersect sphere with triangles
645 #ifndef MB_OBB_USE_VECTOR_QUERIES
646 #ifdef MB_OBB_USE_TYPE_QUERIES
647 rval = get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
648 #else
649 rval = get_moab_instance()->get_entities_by_handle( node, tris );
650 #endif
651 t = tris.begin();
652 #else
653 rval = get_moab_instance()->get_entities_by_handle( node, tris );
654 t = tris.lower_bound( MBTRI );
655 #endif
656 if( MB_SUCCESS != rval ) return rval;
657
658 for( t = tris.begin(); t != tris.end(); ++t )
659 {
660 #ifndef MB_OBB_USE_VECTOR_QUERIES
661 if( TYPE_FROM_HANDLE( *t ) != MBTRI ) break;
662 #elif !defined( MB_OBB_USE_TYPE_QUERIES )
663 if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
664 #endif
665 rval = get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
666 if( MB_SUCCESS != rval ) return rval;
667 if( num_conn != 3 ) continue;
668
669 rval = get_moab_instance()->get_coords( conn, num_conn, coords[0].array() );
670 if( MB_SUCCESS != rval ) return rval;
671
672 GeomUtil::closest_location_on_tri( center, coords, closest );
673 closest -= center;
674 if( ( closest % closest ) <= radsqr &&
675 std::find( facets_out.begin(), facets_out.end(), *t ) == facets_out.end() )
676 {
677 facets_out.push_back( *t );
678 if( sets_out ) sets_out->push_back( surf );
679 }
680 }
681 }
682
683 if( accum )
684 {
685 accum->end_traversal( max_depth );
686 }
687
688 return MB_SUCCESS;
689 }
690
691 /********************** General Ray/Tree and Ray/Triangle Intersection ***************/
692
693 class RayIntersector : public OrientedBoxTreeTool::Op
694 {
695 private:
696 OrientedBoxTreeTool* tool;
697 const CartVect b, m;
698 const double* len;
699 const double tol;
700 Range& boxes;
701
702 public:
703 RayIntersector( OrientedBoxTreeTool* tool_ptr,
704 const double* ray_point,
705 const double* unit_ray_dir,
706 const double* ray_length,
707 double tolerance,
708 Range& leaf_boxes )
709 : tool( tool_ptr ), b( ray_point ), m( unit_ray_dir ), len( ray_length ), tol( tolerance ), boxes( leaf_boxes )
710 {
711 }
712
713 virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
714 virtual ErrorCode leaf( EntityHandle node );
715 };
716
717 //#include <stdio.h>
718 // inline void dump_fragmentation( const Range& range ) {
719 // static FILE* file = fopen( "fragmentation", "w" );
720 // unsigned ranges = 0, entities = 0;
721 // for (Range::const_pair_iterator i = range.const_pair_begin(); i != range.const_pair_end(); ++i)
722 // {
723 // ++ranges;
724 // entities += i->second - i->first + 1;
725 // }
726 // fprintf( file, "%u %u\n", ranges, entities );
727 //}
728
729 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
730 std::vector< EntityHandle >& intersection_facets_out,
731 const Range& boxes,
732 double /*tolerance*/,
733 const double ray_point[3],
734 const double unit_ray_dir[3],
735 const double* ray_length,
736 unsigned int* raytri_test_count )
737 {
738 ErrorCode rval;
739 intersection_distances_out.clear();
740 #ifdef MB_OBB_USE_VECTOR_QUERIES
741 std::vector< EntityHandle > tris;
742 #endif
743
744 const CartVect point( ray_point );
745 const CartVect dir( unit_ray_dir );
746
747 for( Range::iterator b = boxes.begin(); b != boxes.end(); ++b )
748 {
749 #ifndef MB_OBB_USE_VECTOR_QUERIES
750 Range tris;
751 #ifdef MB_OBB_USE_TYPE_QUERIES
752 rval = instance->get_entities_by_type( *b, MBTRI, tris );
753 #else
754 rval = instance->get_entities_by_handle( *b, tris );
755 #endif
756 #else
757 tris.clear();
758 rval = instance->get_entities_by_handle( *b, tris );
759 #endif
760 if( MB_SUCCESS != rval ) return rval;
761 // dump_fragmentation( tris );
762
763 #ifndef MB_OBB_USE_VECTOR_QUERIES
764 for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
765 #else
766 for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
767 #endif
768 {
769 #ifndef MB_OBB_USE_TYPE_QUERIES
770 if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
771 #endif
772
773 const EntityHandle* conn = NULL;
774 int len = 0;
775 rval = instance->get_connectivity( *t, conn, len, true );
776 if( MB_SUCCESS != rval ) return rval;
777
778 CartVect coords[3];
779 rval = instance->get_coords( conn, 3, coords[0].array() );
780 if( MB_SUCCESS != rval ) return rval;
781
782 if( raytri_test_count ) *raytri_test_count += 1;
783
784 double td;
785 if( GeomUtil::plucker_ray_tri_intersect( coords, point, dir, td, ray_length ) )
786 {
787 intersection_distances_out.push_back( td );
788 intersection_facets_out.push_back( *t );
789 }
790 }
791 }
792
793 return MB_SUCCESS;
794 }
795
796 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
797 std::vector< EntityHandle >& intersection_facets_out,
798 EntityHandle root_set,
799 double tolerance,
800 const double ray_point[3],
801 const double unit_ray_dir[3],
802 const double* ray_length,
803 TrvStats* accum )
804 {
805 Range boxes;
806 ErrorCode rval;
807
808 rval = ray_intersect_boxes( boxes, root_set, tolerance, ray_point, unit_ray_dir, ray_length, accum );
809 if( MB_SUCCESS != rval ) return rval;
810
811 return ray_intersect_triangles( intersection_distances_out, intersection_facets_out, boxes, tolerance, ray_point,
812 unit_ray_dir, ray_length, accum ? &( accum->ray_tri_tests_count ) : NULL );
813 }
814
815 ErrorCode OrientedBoxTreeTool::ray_intersect_boxes( Range& boxes_out,
816 EntityHandle root_set,
817 double tolerance,
818 const double ray_point[3],
819 const double unit_ray_dir[3],
820 const double* ray_length,
821 TrvStats* accum )
822 {
823 RayIntersector op( this, ray_point, unit_ray_dir, ray_length, tolerance, boxes_out );
824 return preorder_traverse( root_set, op, accum );
825 }
826
827 ErrorCode RayIntersector::visit( EntityHandle node, int, bool& descend )
828 {
829 OrientedBox box;
830 ErrorCode rval = tool->box( node, box );
831 if( MB_SUCCESS != rval ) return rval;
832
833 descend = box.intersect_ray( b, m, tol, len );
834 return MB_SUCCESS;
835 }
836
837 ErrorCode RayIntersector::leaf( EntityHandle node )
838 {
839 boxes.insert( node );
840 return MB_SUCCESS;
841 }
842
843 /********************** Ray/Set Intersection ****************************/
844
845 ErrorCode OrientedBoxTreeTool::get_close_tris( CartVect int_pt,
846 double tol,
847 const EntityHandle* rootSet,
848 const EntityHandle* geomVol,
849 const Tag* senseTag,
850 std::vector< EntityHandle >& close_tris,
851 std::vector< int >& close_senses )
852 {
853 std::vector< EntityHandle > close_surfs;
854 ErrorCode rval = sphere_intersect_triangles( int_pt.array(), tol, *rootSet, close_tris, &close_surfs );
855 assert( MB_SUCCESS == rval );
856 if( MB_SUCCESS != rval ) return rval;
857
858 // for each surface, get the surf sense wrt parent volume
859 close_senses.resize( close_surfs.size() );
860 for( unsigned i = 0; i < close_surfs.size(); ++i )
861 {
862 EntityHandle vols[2];
863 rval = get_moab_instance()->tag_get_data( *senseTag, &( close_surfs[i] ), 1, vols );
864 assert( MB_SUCCESS == rval );
865 if( MB_SUCCESS != rval ) return rval;
866 if( vols[0] == vols[1] )
867 {
868 std::cerr << "error: surf has positive and negative sense wrt same volume" << std::endl;
869 return MB_FAILURE;
870 }
871 if( *geomVol == vols[0] )
872 {
873 close_senses[i] = 1;
874 }
875 else if( *geomVol == vols[1] )
876 {
877 close_senses[i] = -1;
878 }
879 else
880 {
881 return MB_FAILURE;
882 }
883 }
884
885 return MB_SUCCESS;
886 }
887
888 class RayIntersectSets : public OrientedBoxTreeTool::Op
889 {
890 private:
891 // Input
892 OrientedBoxTreeTool* tool;
893 const CartVect ray_origin;
894 const CartVect ray_direction;
895 OrientedBoxTreeTool::IntersectSearchWindow& search_win; /* length to search ahead/behind of ray origin */
896 const double tol; /* used for box.intersect_ray, radius of
897 neighborhood for adjacent triangles,
898 and old mode of add_intersection */
899 OrientedBoxTreeTool::IntRegCtxt& int_reg_callback;
900
901 // Optional Input - to screen RTIs by orientation and edge/node intersection
902 int* surfTriOrient; /* holds desired orientation of tri wrt surface */
903 int surfTriOrient_val;
904
905 // Other Variables
906 unsigned int* raytri_test_count;
907 EntityHandle lastSet;
908 int lastSetDepth;
909
910 public:
911 RayIntersectSets( OrientedBoxTreeTool* tool_ptr,
912 const double* ray_point,
913 const double* unit_ray_dir,
914 const double tolerance,
915 OrientedBoxTreeTool::IntersectSearchWindow& win,
916 unsigned int* ray_tri_test_count,
917 OrientedBoxTreeTool::IntRegCtxt& intRegCallback )
918 : tool( tool_ptr ), ray_origin( ray_point ), ray_direction( unit_ray_dir ), search_win( win ), tol( tolerance ),
919 int_reg_callback( intRegCallback ), surfTriOrient_val( 0 ), raytri_test_count( ray_tri_test_count ),
920 lastSet( 0 ), lastSetDepth( 0 )
921 {
922
923 // specified orientation should be 1 or -1, indicating ray and surface
924 // normal in the same or opposite directions, respectively.
925 if( int_reg_callback.getDesiredOrient() )
926 {
927 surfTriOrient = &surfTriOrient_val;
928 }
929 else
930 {
931 surfTriOrient = NULL;
932 }
933
934 // check the limits
935 if( search_win.first )
936 {
937 assert( 0 <= *( search_win.first ) );
938 }
939 if( search_win.second )
940 {
941 assert( 0 >= *( search_win.second ) );
942 }
943 };
944
945 virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
946 virtual ErrorCode leaf( EntityHandle node );
947 };
948
949 ErrorCode RayIntersectSets::visit( EntityHandle node, int depth, bool& descend )
950 {
951 OrientedBox box;
952 ErrorCode rval = tool->box( node, box );
953 assert( MB_SUCCESS == rval );
954 if( MB_SUCCESS != rval ) return rval;
955
956 descend = box.intersect_ray( ray_origin, ray_direction, tol, search_win.first, search_win.second );
957
958 if( lastSet && depth <= lastSetDepth ) lastSet = 0;
959
960 if( descend && !lastSet )
961 {
962 Range tmp_sets;
963 rval = tool->get_moab_instance()->get_entities_by_type( node, MBENTITYSET, tmp_sets );
964 assert( MB_SUCCESS == rval );
965 if( MB_SUCCESS != rval ) return rval;
966
967 if( !tmp_sets.empty() )
968 {
969 if( tmp_sets.size() > 1 ) return MB_FAILURE;
970 lastSet = *tmp_sets.begin();
971 lastSetDepth = depth;
972
973 rval = int_reg_callback.update_orient( lastSet, surfTriOrient );
974 if( MB_SUCCESS != rval ) return rval;
975 }
976 }
977
978 return MB_SUCCESS;
979 }
980
981 ErrorCode RayIntersectSets::leaf( EntityHandle node )
982 {
983 assert( lastSet );
984 if( !lastSet ) // if no surface has been visited yet, something's messed up.
985 return MB_FAILURE;
986
987 #ifndef MB_OBB_USE_VECTOR_QUERIES
988 Range tris;
989 #ifdef MB_OBB_USE_TYPE_QUERIES
990 ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
991 #else
992 ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
993 #endif
994 #else
995 std::vector< EntityHandle > tris;
996 ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
997 #endif
998 assert( MB_SUCCESS == rval );
999 if( MB_SUCCESS != rval ) return rval;
1000
1001 #ifndef MB_OBB_USE_VECTOR_QUERIES
1002 for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
1003 #else
1004 for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
1005 #endif
1006 {
1007 #ifndef MB_OBB_USE_TYPE_QUERIES
1008 if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
1009 #endif
1010
1011 const EntityHandle* conn;
1012 int num_conn;
1013 rval = tool->get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
1014 assert( MB_SUCCESS == rval );
1015 if( MB_SUCCESS != rval ) return rval;
1016
1017 CartVect coords[3];
1018 rval = tool->get_moab_instance()->get_coords( conn, 3, coords[0].array() );
1019 assert( MB_SUCCESS == rval );
1020 if( MB_SUCCESS != rval ) return rval;
1021
1022 if( raytri_test_count ) *raytri_test_count += 1;
1023
1024 double int_dist;
1025 GeomUtil::intersection_type int_type = GeomUtil::NONE;
1026
1027 if( GeomUtil::plucker_ray_tri_intersect( coords, ray_origin, ray_direction, int_dist, search_win.first,
1028 search_win.second, surfTriOrient, &int_type ) )
1029 {
1030 int_reg_callback.register_intersection( lastSet, *t, int_dist, search_win, int_type );
1031 }
1032 }
1033 return MB_SUCCESS;
1034 }
1035
1036 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
1037 std::vector< EntityHandle >& sets_out,
1038 std::vector< EntityHandle >& facets_out,
1039 EntityHandle root_set,
1040 const double tolerance,
1041 const double ray_point[3],
1042 const double unit_ray_dir[3],
1043 IntersectSearchWindow& search_win,
1044 IntRegCtxt& int_reg_callback,
1045 TrvStats* accum )
1046
1047 {
1048 RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1049 accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
1050 ErrorCode rval = preorder_traverse( root_set, op, accum );
1051
1052 distances_out = int_reg_callback.get_intersections();
1053 sets_out = int_reg_callback.get_sets();
1054 facets_out = int_reg_callback.get_facets();
1055
1056 return rval;
1057 }
1058
1059 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
1060 std::vector< EntityHandle >& sets_out,
1061 std::vector< EntityHandle >& facets_out,
1062 EntityHandle root_set,
1063 double tolerance,
1064 const double ray_point[3],
1065 const double unit_ray_dir[3],
1066 const double* ray_length,
1067 TrvStats* accum )
1068 {
1069 IntRegCtxt int_reg_ctxt;
1070
1071 OrientedBoxTreeTool::IntersectSearchWindow search_win( ray_length, (double*)0 );
1072
1073 RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1074 accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_ctxt );
1075
1076 ErrorCode rval = preorder_traverse( root_set, op, accum );
1077
1078 if( MB_SUCCESS != rval )
1079 {
1080 return rval;
1081 }
1082
1083 distances_out = int_reg_ctxt.get_intersections();
1084 sets_out = int_reg_ctxt.get_sets();
1085 facets_out = int_reg_ctxt.get_facets();
1086
1087 return MB_SUCCESS;
1088 }
1089
1090 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( EntityHandle root_set,
1091 const double tolerance,
1092 const double ray_point[3],
1093 const double unit_ray_dir[3],
1094 IntersectSearchWindow& search_win,
1095 IntRegCtxt& int_reg_callback,
1096 TrvStats* accum )
1097
1098 {
1099 RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1100 accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
1101 return preorder_traverse( root_set, op, accum );
1102 }
1103
1104 /********************** Closest Point code ***************/
1105
1106 struct OBBTreeCPFrame
1107 {
1108 OBBTreeCPFrame( double d, EntityHandle n, EntityHandle s, int dp )
1109 : dist_sqr( d ), node( n ), mset( s ), depth( dp )
1110 {
1111 }
1112 double dist_sqr;
1113 EntityHandle node;
1114 EntityHandle mset;
1115 int depth;
1116 };
1117
1118 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point,
1119 EntityHandle root,
1120 double* point_out,
1121 EntityHandle& facet_out,
1122 EntityHandle* set_out,
1123 TrvStats* accum )
1124 {
1125 ErrorCode rval;
1126 const CartVect loc( point );
1127 double smallest_dist_sqr = std::numeric_limits< double >::max();
1128
1129 EntityHandle current_set = 0;
1130 Range sets;
1131 std::vector< EntityHandle > children( 2 );
1132 std::vector< double > coords;
1133 std::vector< OBBTreeCPFrame > stack;
1134 int max_depth = -1;
1135
1136 stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
1137
1138 while( !stack.empty() )
1139 {
1140
1141 // pop from top of stack
1142 EntityHandle node = stack.back().node;
1143 double dist_sqr = stack.back().dist_sqr;
1144 current_set = stack.back().mset;
1145 int current_depth = stack.back().depth;
1146 stack.pop_back();
1147
1148 // If current best result is closer than the box, skip this tree node.
1149 if( dist_sqr > smallest_dist_sqr ) continue;
1150
1151 // increment traversal statistics.
1152 if( accum )
1153 {
1154 accum->increment( current_depth );
1155 max_depth = std::max( max_depth, current_depth );
1156 }
1157
1158 // Check if this node has a set associated with it
1159 if( set_out && !current_set )
1160 {
1161 sets.clear();
1162 rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
1163 if( MB_SUCCESS != rval ) return rval;
1164 if( !sets.empty() )
1165 {
1166 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
1167 current_set = sets.front();
1168 }
1169 }
1170
1171 // Get child boxes
1172 children.clear();
1173 rval = instance->get_child_meshsets( node, children );
1174 if( MB_SUCCESS != rval ) return rval;
1175
1176 // if not a leaf node
1177 if( !children.empty() )
1178 {
1179 if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1180
1181 // get boxes
1182 OrientedBox box1, box2;
1183 rval = box( children[0], box1 );
1184 if( MB_SUCCESS != rval ) return rval;
1185 rval = box( children[1], box2 );
1186 if( MB_SUCCESS != rval ) return rval;
1187
1188 // get distance from each box
1189 CartVect pt1, pt2;
1190 box1.closest_location_in_box( loc, pt1 );
1191 box2.closest_location_in_box( loc, pt2 );
1192 pt1 -= loc;
1193 pt2 -= loc;
1194 const double dsqr1 = pt1 % pt1;
1195 const double dsqr2 = pt2 % pt2;
1196
1197 // push children on tree such that closer one is on top
1198 if( dsqr1 < dsqr2 )
1199 {
1200 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1201 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1202 }
1203 else
1204 {
1205 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1206 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1207 }
1208 }
1209 else
1210 { // LEAF NODE
1211 if( accum )
1212 {
1213 accum->increment_leaf( current_depth );
1214 }
1215
1216 Range facets;
1217 rval = instance->get_entities_by_dimension( node, 2, facets );
1218 if( MB_SUCCESS != rval ) return rval;
1219
1220 const EntityHandle* conn = NULL;
1221 int len = 0;
1222 CartVect tmp, diff;
1223 for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
1224 {
1225 rval = instance->get_connectivity( *i, conn, len, true );
1226 if( MB_SUCCESS != rval ) return rval;
1227
1228 coords.resize( 3 * len );
1229 rval = instance->get_coords( conn, len, &coords[0] );
1230 if( MB_SUCCESS != rval ) return rval;
1231
1232 if( len == 3 )
1233 GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
1234 else
1235 GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
1236
1237 diff = tmp - loc;
1238 dist_sqr = diff % diff;
1239 if( dist_sqr < smallest_dist_sqr )
1240 {
1241 smallest_dist_sqr = dist_sqr;
1242 facet_out = *i;
1243 tmp.get( point_out );
1244 if( set_out ) *set_out = current_set;
1245 }
1246 }
1247 } // LEAF NODE
1248 }
1249
1250 if( accum )
1251 {
1252 accum->end_traversal( max_depth );
1253 }
1254
1255 return MB_SUCCESS;
1256 }
1257
1258 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point,
1259 EntityHandle root,
1260 double tolerance,
1261 std::vector< EntityHandle >& facets_out,
1262 std::vector< EntityHandle >* sets_out,
1263 TrvStats* accum )
1264 {
1265 ErrorCode rval;
1266 const CartVect loc( point );
1267 double smallest_dist_sqr = std::numeric_limits< double >::max();
1268 double smallest_dist = smallest_dist_sqr;
1269
1270 EntityHandle current_set = 0;
1271 Range sets;
1272 std::vector< EntityHandle > children( 2 );
1273 std::vector< double > coords;
1274 std::vector< OBBTreeCPFrame > stack;
1275 int max_depth = -1;
1276
1277 stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
1278
1279 while( !stack.empty() )
1280 {
1281
1282 // pop from top of stack
1283 EntityHandle node = stack.back().node;
1284 double dist_sqr = stack.back().dist_sqr;
1285 current_set = stack.back().mset;
1286 int current_depth = stack.back().depth;
1287 stack.pop_back();
1288
1289 // If current best result is closer than the box, skip this tree node.
1290 if( dist_sqr > smallest_dist_sqr + tolerance ) continue;
1291
1292 // increment traversal statistics.
1293 if( accum )
1294 {
1295 accum->increment( current_depth );
1296 max_depth = std::max( max_depth, current_depth );
1297 }
1298
1299 // Check if this node has a set associated with it
1300 if( sets_out && !current_set )
1301 {
1302 sets.clear();
1303 rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
1304 if( MB_SUCCESS != rval ) return rval;
1305 if( !sets.empty() )
1306 {
1307 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
1308 current_set = *sets.begin();
1309 }
1310 }
1311
1312 // Get child boxes
1313 children.clear();
1314 rval = instance->get_child_meshsets( node, children );
1315 if( MB_SUCCESS != rval ) return rval;
1316
1317 // if not a leaf node
1318 if( !children.empty() )
1319 {
1320 if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1321
1322 // get boxes
1323 OrientedBox box1, box2;
1324 rval = box( children[0], box1 );
1325 if( MB_SUCCESS != rval ) return rval;
1326 rval = box( children[1], box2 );
1327 if( MB_SUCCESS != rval ) return rval;
1328
1329 // get distance from each box
1330 CartVect pt1, pt2;
1331 box1.closest_location_in_box( loc, pt1 );
1332 box2.closest_location_in_box( loc, pt2 );
1333 pt1 -= loc;
1334 pt2 -= loc;
1335 const double dsqr1 = pt1 % pt1;
1336 const double dsqr2 = pt2 % pt2;
1337
1338 // push children on tree such that closer one is on top
1339 if( dsqr1 < dsqr2 )
1340 {
1341 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1342 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1343 }
1344 else
1345 {
1346 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1347 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1348 }
1349 }
1350 else
1351 { // LEAF NODE
1352 if( accum )
1353 {
1354 accum->increment_leaf( current_depth );
1355 }
1356
1357 Range facets;
1358 rval = instance->get_entities_by_dimension( node, 2, facets );
1359 if( MB_SUCCESS != rval ) return rval;
1360
1361 const EntityHandle* conn = NULL;
1362 int len = 0;
1363 CartVect tmp, diff;
1364 for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
1365 {
1366 rval = instance->get_connectivity( *i, conn, len, true );
1367 if( MB_SUCCESS != rval ) return rval;
1368
1369 coords.resize( 3 * len );
1370 rval = instance->get_coords( conn, len, &coords[0] );
1371 if( MB_SUCCESS != rval ) return rval;
1372
1373 if( len == 3 )
1374 GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
1375 else
1376 GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
1377
1378 diff = tmp - loc;
1379 dist_sqr = diff % diff;
1380 if( dist_sqr < smallest_dist_sqr )
1381 {
1382 if( 0.5 * dist_sqr < 0.5 * smallest_dist_sqr + tolerance * ( 0.5 * tolerance - smallest_dist ) )
1383 {
1384 facets_out.clear();
1385 if( sets_out ) sets_out->clear();
1386 }
1387 smallest_dist_sqr = dist_sqr;
1388 smallest_dist = sqrt( smallest_dist_sqr );
1389 // put closest result at start of lists
1390 facets_out.push_back( *i );
1391 std::swap( facets_out.front(), facets_out.back() );
1392 if( sets_out )
1393 {
1394 sets_out->push_back( current_set );
1395 std::swap( sets_out->front(), sets_out->back() );
1396 }
1397 }
1398 else if( dist_sqr <= smallest_dist_sqr + tolerance * ( tolerance + 2 * smallest_dist ) )
1399 {
1400 facets_out.push_back( *i );
1401 if( sets_out ) sets_out->push_back( current_set );
1402 }
1403 }
1404 } // LEAF NODE
1405 }
1406
1407 if( accum )
1408 {
1409 accum->end_traversal( max_depth );
1410 }
1411
1412 return MB_SUCCESS;
1413 }
1414
1415 /********************** Tree Printing Code ****************************/
1416
1417 class TreeLayoutPrinter : public OrientedBoxTreeTool::Op
1418 {
1419 public:
1420 TreeLayoutPrinter( std::ostream& stream, Interface* instance );
1421
1422 virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
1423 virtual ErrorCode leaf( EntityHandle node );
1424
1425 private:
1426 Interface* instance;
1427 std::ostream& outputStream;
1428 std::vector< bool > path;
1429 };
1430
1431 TreeLayoutPrinter::TreeLayoutPrinter( std::ostream& stream, Interface* interface )
1432 : instance( interface ), outputStream( stream )
1433 {
1434 }
1435
1436 ErrorCode TreeLayoutPrinter::visit( EntityHandle node, int depth, bool& descend )
1437 {
1438 descend = true;
1439
1440 if( (unsigned)depth > path.size() )
1441 {
1442 // assert(depth+1 == path.size); // preorder traversal
1443 path.push_back( true );
1444 }
1445 else
1446 {
1447 path.resize( depth );
1448 if( depth ) path.back() = false;
1449 }
1450
1451 for( unsigned i = 0; i + 1 < path.size(); ++i )
1452 {
1453 if( path[i] )
1454 outputStream << "| ";
1455 else
1456 outputStream << " ";
1457 }
1458 if( depth )
1459 {
1460 if( path.back() )
1461 outputStream << "+---";
1462 else
1463 outputStream << "\\---";
1464 }
1465 outputStream << instance->id_from_handle( node ) << std::endl;
1466 return MB_SUCCESS;
1467 }
1468
1469 ErrorCode TreeLayoutPrinter::leaf( EntityHandle )
1470 {
1471 return MB_SUCCESS;
1472 }
1473
1474 class TreeNodePrinter : public OrientedBoxTreeTool::Op
1475 {
1476 public:
1477 TreeNodePrinter( std::ostream& stream,
1478 bool list_contents,
1479 bool list_box,
1480 const char* id_tag_name,
1481 OrientedBoxTreeTool* tool_ptr );
1482
1483 virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
1484
1485 virtual ErrorCode leaf( EntityHandle )
1486 {
1487 return MB_SUCCESS;
1488 }
1489
1490 private:
1491 ErrorCode print_geometry( EntityHandle node );
1492 ErrorCode print_contents( EntityHandle node );
1493 ErrorCode print_counts( EntityHandle node );
1494
1495 bool printContents;
1496 bool printGeometry;
1497 bool haveTag;
1498 Tag tag, gidTag, geomTag;
1499 Interface* instance;
1500 OrientedBoxTreeTool* tool;
1501 std::ostream& outputStream;
1502 };
1503
1504 TreeNodePrinter::TreeNodePrinter( std::ostream& stream,
1505 bool list_contents,
1506 bool list_box,
1507 const char* id_tag_name,
1508 OrientedBoxTreeTool* tool_ptr )
1509 : printContents( list_contents ), printGeometry( list_box ), haveTag( false ), tag( 0 ), gidTag( 0 ), geomTag( 0 ),
1510 instance( tool_ptr->get_moab_instance() ), tool( tool_ptr ), outputStream( stream )
1511 {
1512 ErrorCode rval;
1513 if( id_tag_name )
1514 {
1515 rval = instance->tag_get_handle( id_tag_name, 1, MB_TYPE_INTEGER, tag );
1516 if( !rval )
1517 {
1518 std::cerr << "Could not get tag \"" << id_tag_name << "\"\n";
1519 stream << "Could not get tag \"" << id_tag_name << "\"\n";
1520 }
1521 else
1522 {
1523 haveTag = true;
1524 }
1525 }
1526
1527 gidTag = instance->globalId_tag();
1528
1529 rval = instance->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );
1530 if( MB_SUCCESS != rval ) geomTag = 0;
1531 }
1532
1533 ErrorCode TreeNodePrinter::visit( EntityHandle node, int, bool& descend )
1534 {
1535 descend = true;
1536 EntityHandle setid = instance->id_from_handle( node );
1537 outputStream << setid << ":" << std::endl;
1538
1539 Range surfs;
1540 ErrorCode r3 = MB_SUCCESS;
1541 if( geomTag )
1542 {
1543 const int two = 2;
1544 const void* tagdata[] = { &two };
1545 r3 = instance->get_entities_by_type_and_tag( node, MBENTITYSET, &geomTag, tagdata, 1, surfs );
1546
1547 if( MB_SUCCESS == r3 && surfs.size() == 1 )
1548 {
1549 EntityHandle surf = *surfs.begin();
1550 int id;
1551 if( gidTag && MB_SUCCESS == instance->tag_get_data( gidTag, &surf, 1, &id ) )
1552 outputStream << " Surface " << id << std::endl;
1553 else
1554 outputStream << " Surface w/ unknown ID (" << surf << ")" << std::endl;
1555 }
1556 }
1557
1558 ErrorCode r1 = printGeometry ? print_geometry( node ) : MB_SUCCESS;
1559 ErrorCode r2 = printContents ? print_contents( node ) : print_counts( node );
1560 outputStream << std::endl;
1561
1562 if( MB_SUCCESS != r1 )
1563 return r1;
1564 else if( MB_SUCCESS != r2 )
1565 return r2;
1566 else
1567 return r3;
1568 }
1569
1570 ErrorCode TreeNodePrinter::print_geometry( EntityHandle node )
1571 {
1572 OrientedBox box;
1573 ErrorCode rval = tool->box( node, box );
1574 if( MB_SUCCESS != rval ) return rval;
1575
1576 CartVect length = box.dimensions();
1577
1578 outputStream << box.center << " Radius: " << box.inner_radius() << " - " << box.outer_radius() << std::endl
1579 << '+' << box.axis( 0 ) << " : " << length[0] << std::endl
1580 << 'x' << box.axis( 1 ) << " : " << length[1] << std::endl
1581 << 'x' << box.axis( 2 ) << " : " << length[2] << std::endl;
1582 return MB_SUCCESS;
1583 }
1584
1585 ErrorCode TreeNodePrinter::print_counts( EntityHandle node )
1586 {
1587 for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
1588 {
1589 int count = 0;
1590 ErrorCode rval = instance->get_number_entities_by_type( node, type, count );
1591 if( MB_SUCCESS != rval ) return rval;
1592 if( count > 0 ) outputStream << " " << count << " " << CN::EntityTypeName( type ) << std::endl;
1593 }
1594 return MB_SUCCESS;
1595 }
1596
1597 ErrorCode TreeNodePrinter::print_contents( EntityHandle node )
1598 {
1599 // list contents
1600 for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
1601 {
1602 Range range;
1603 ErrorCode rval = instance->get_entities_by_type( node, type, range );
1604 if( MB_SUCCESS != rval ) return rval;
1605 if( range.empty() ) continue;
1606 outputStream << " " << CN::EntityTypeName( type ) << " ";
1607 std::vector< int > ids( range.size() );
1608 if( haveTag )
1609 {
1610 rval = instance->tag_get_data( tag, range, &ids[0] );
1611 std::sort( ids.begin(), ids.end() );
1612 }
1613 else
1614 {
1615 Range::iterator ri = range.begin();
1616 std::vector< int >::iterator vi = ids.begin();
1617 while( ri != range.end() )
1618 {
1619 *vi = instance->id_from_handle( *ri );
1620 ++ri;
1621 ++vi;
1622 }
1623 }
1624
1625 unsigned i = 0;
1626 for( ;; )
1627 {
1628 unsigned beg = i, end;
1629 do
1630 {
1631 end = i++;
1632 } while( i < ids.size() && ids[end] + 1 == ids[i] );
1633 if( end == beg )
1634 outputStream << ids[end];
1635 else if( end == beg + 1 )
1636 outputStream << ids[beg] << ", " << ids[end];
1637 else
1638 outputStream << ids[beg] << "-" << ids[end];
1639
1640 if( i == ids.size() )
1641 {
1642 outputStream << std::endl;
1643 break;
1644 }
1645 else
1646 outputStream << ", ";
1647 }
1648 }
1649
1650 return MB_SUCCESS;
1651 }
1652
1653 void OrientedBoxTreeTool::print( EntityHandle set, std::ostream& str, bool list, const char* tag )
1654 {
1655 TreeLayoutPrinter op1( str, instance );
1656 TreeNodePrinter op2( str, list, true, tag, this );
1657 ErrorCode r1 = preorder_traverse( set, op1 );
1658 str << std::endl;
1659 ErrorCode r2 = preorder_traverse( set, op2 );
1660 if( r1 != MB_SUCCESS || r2 != MB_SUCCESS )
1661 {
1662 std::cerr << "Errors encountered while printing tree\n";
1663 str << "Errors encountered while printing tree\n";
1664 }
1665 }
1666
1667 /********************* Traversal Metrics Code **************************/
1668
1669 void OrientedBoxTreeTool::TrvStats::reset()
1670 {
1671 nodes_visited_count.clear();
1672 leaves_visited_count.clear();
1673 traversals_ended_count.clear();
1674 ray_tri_tests_count = 0;
1675 }
1676
1677 void OrientedBoxTreeTool::TrvStats::increment( unsigned depth )
1678 {
1679
1680 while( nodes_visited_count.size() <= depth )
1681 {
1682 nodes_visited_count.push_back( 0 );
1683 leaves_visited_count.push_back( 0 );
1684 traversals_ended_count.push_back( 0 );
1685 }
1686 nodes_visited_count[depth] += 1;
1687 }
1688
1689 void OrientedBoxTreeTool::TrvStats::increment_leaf( unsigned depth )
1690 {
1691 // assume safe depth, because increment is called first
1692 leaves_visited_count[depth] += 1;
1693 }
1694
1695 void OrientedBoxTreeTool::TrvStats::end_traversal( unsigned depth )
1696 {
1697 // assume safe depth, because increment is always called on a given
1698 // tree level first
1699 traversals_ended_count[depth] += 1;
1700 }
1701
1702 void OrientedBoxTreeTool::TrvStats::print( std::ostream& str ) const
1703 {
1704
1705 const std::string h1 = "OBBTree Depth";
1706 const std::string h2 = " - NodesVisited";
1707 const std::string h3 = " - LeavesVisited";
1708 const std::string h4 = " - TraversalsEnded";
1709
1710 str << h1 << h2 << h3 << h4 << std::endl;
1711
1712 unsigned num_visited = 0, num_leaves = 0, num_traversals = 0;
1713 for( unsigned i = 0; i < traversals_ended_count.size(); ++i )
1714 {
1715
1716 num_visited += nodes_visited_count[i];
1717 num_leaves += leaves_visited_count[i];
1718 num_traversals += traversals_ended_count[i];
1719
1720 str << std::setw( h1.length() ) << i << std::setw( h2.length() ) << nodes_visited_count[i]
1721 << std::setw( h3.length() ) << leaves_visited_count[i] << std::setw( h4.length() )
1722 << traversals_ended_count[i] << std::endl;
1723 }
1724
1725 str << std::setw( h1.length() ) << "---- Totals:" << std::setw( h2.length() ) << num_visited
1726 << std::setw( h3.length() ) << num_leaves << std::setw( h4.length() ) << num_traversals << std::endl;
1727
1728 if( ray_tri_tests_count )
1729 {
1730 str << std::setw( h1.length() ) << "---- Total ray-tri tests: " << ray_tri_tests_count << std::endl;
1731 }
1732 }
1733
1734 /********************** Tree Statistics Code ****************************/
1735
1736 class StatData
1737 {
1738 public:
1739 struct Ratio
1740 {
1741 double min, max, sum, sqr;
1742 int hist[10];
1743 Ratio()
1744 : min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() ), sum( 0.0 ),
1745 sqr( 0.0 )
1746 {
1747 hist[0] = hist[1] = hist[2] = hist[3] = hist[4] = hist[5] = hist[6] = hist[7] = hist[8] = hist[9] = 0;
1748 }
1749 void accum( double v )
1750 {
1751 if( v < min ) min = v;
1752 if( v > max ) max = v;
1753 sum += v;
1754 sqr += v * v;
1755 int i = (int)( 10 * v );
1756 if( i < 0 )
1757 i = 0;
1758 else if( i > 9 )
1759 i = 9;
1760 ++hist[i];
1761 }
1762 };
1763
1764 template < typename T >
1765 struct Stat
1766 {
1767 T min, max;
1768 double sum, sqr;
1769 Stat() : sum( 0.0 ), sqr( 0.0 )
1770 {
1771 std::numeric_limits< T > lim;
1772 min = lim.max();
1773 if( lim.is_integer )
1774 max = lim.min();
1775 else
1776 max = -lim.max();
1777 }
1778 void accum( T v )
1779 {
1780 if( v < min ) min = v;
1781 if( v > max ) max = v;
1782 sum += v;
1783 sqr += (double)v * v;
1784 }
1785 };
1786
1787 StatData() : count( 0 ) {}
1788
1789 Ratio volume;
1790 Ratio entities;
1791 Ratio radius;
1792 Stat< unsigned > leaf_ent;
1793 Stat< double > vol;
1794 Stat< double > area;
1795 std::vector< unsigned > leaf_depth;
1796 unsigned count;
1797 };
1798
1799 static int measure( const CartVect& v, double& result )
1800 {
1801 const double tol = 1e-6;
1802 int dims = 0;
1803 result = 1;
1804 for( int i = 0; i < 3; ++i )
1805 if( v[i] > tol )
1806 {
1807 ++dims;
1808 result *= v[i];
1809 }
1810 return dims;
1811 }
1812
1813 ErrorCode OrientedBoxTreeTool::recursive_stats( OrientedBoxTreeTool* tool,
1814 Interface* inst,
1815 EntityHandle set,
1816 int depth,
1817 StatData& data,
1818 unsigned& count_out,
1819 CartVect& dimensions_out )
1820 {
1821 ErrorCode rval;
1822 OrientedBox tmp_box;
1823 std::vector< EntityHandle > children( 2 );
1824 unsigned counts[2];
1825 bool isleaf;
1826
1827 ++data.count;
1828
1829 rval = tool->box( set, tmp_box );
1830 if( MB_SUCCESS != rval ) return rval;
1831 children.clear();
1832 rval = inst->get_child_meshsets( set, children );
1833 if( MB_SUCCESS != rval ) return rval;
1834 isleaf = children.empty();
1835 if( !isleaf && children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1836
1837 dimensions_out = tmp_box.dimensions();
1838 data.radius.accum( tmp_box.inner_radius() / tmp_box.outer_radius() );
1839 data.vol.accum( tmp_box.volume() );
1840 data.area.accum( tmp_box.area() );
1841
1842 if( isleaf )
1843 {
1844 if( data.leaf_depth.size() <= (unsigned)depth ) data.leaf_depth.resize( depth + 1, 0 );
1845 ++data.leaf_depth[depth];
1846
1847 int count;
1848 rval = inst->get_number_entities_by_handle( set, count );
1849 if( MB_SUCCESS != rval ) return rval;
1850 count_out = count;
1851 data.leaf_ent.accum( count_out );
1852 }
1853 else
1854 {
1855 for( int i = 0; i < 2; ++i )
1856 {
1857 CartVect dims;
1858 rval = recursive_stats( tool, inst, children[i], depth + 1, data, counts[i], dims );
1859 if( MB_SUCCESS != rval ) return rval;
1860 double this_measure, chld_measure;
1861 int this_dim = measure( dimensions_out, this_measure );
1862 int chld_dim = measure( dims, chld_measure );
1863 double ratio;
1864 if( chld_dim < this_dim )
1865 ratio = 0;
1866 else
1867 ratio = chld_measure / this_measure;
1868
1869 data.volume.accum( ratio );
1870 }
1871 count_out = counts[0] + counts[1];
1872 data.entities.accum( (double)counts[0] / count_out );
1873 data.entities.accum( (double)counts[1] / count_out );
1874 }
1875 return MB_SUCCESS;
1876 }
1877
1878 static inline double std_dev( double sqr, double sum, double count )
1879 {
1880 sum /= count;
1881 sqr /= count;
1882 return sqrt( sqr - sum * sum );
1883 }
1884
1885 //#define WW <<std::setw(10)<<std::fixed<<
1886 #define WE << std::setw( 10 ) <<
1887 #define WW WE
1888 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set,
1889 unsigned& total_entities,
1890 double& rv,
1891 double& tot_node_volume,
1892 double& tot_to_root_volume,
1893 unsigned& tree_height,
1894 unsigned& node_count,
1895 unsigned& num_leaves )
1896 {
1897 StatData d;
1898 ErrorCode rval;
1899 unsigned i;
1900 CartVect total_dim;
1901
1902 rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
1903 if( MB_SUCCESS != rval ) return rval;
1904
1905 tree_height = d.leaf_depth.size();
1906 unsigned min_leaf_depth = tree_height;
1907 num_leaves = 0;
1908 unsigned max_leaf_per_depth = 0;
1909 // double sum_leaf_depth = 0, sqr_leaf_depth = 0;
1910 for( i = 0; i < d.leaf_depth.size(); ++i )
1911 {
1912 unsigned val = d.leaf_depth[i];
1913 num_leaves += val;
1914 // sum_leaf_depth += (double)val*i;
1915 // sqr_leaf_depth += (double)val*i*i;
1916 if( val && i < min_leaf_depth ) min_leaf_depth = i;
1917 if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
1918 }
1919 rv = total_dim[0] * total_dim[1] * total_dim[2];
1920 tot_node_volume = d.vol.sum;
1921 tot_to_root_volume = d.vol.sum / rv;
1922 node_count = d.count;
1923
1924 return MB_SUCCESS;
1925 }
1926
1927 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, std::ostream& s )
1928 {
1929 StatData d;
1930 ErrorCode rval;
1931 unsigned total_entities, i;
1932 CartVect total_dim;
1933
1934 rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
1935 if( MB_SUCCESS != rval ) return rval;
1936
1937 unsigned tree_height = d.leaf_depth.size();
1938 unsigned min_leaf_depth = tree_height, num_leaves = 0;
1939 unsigned max_leaf_per_depth = 0;
1940 double sum_leaf_depth = 0, sqr_leaf_depth = 0;
1941 for( i = 0; i < d.leaf_depth.size(); ++i )
1942 {
1943 unsigned val = d.leaf_depth[i];
1944 num_leaves += val;
1945 sum_leaf_depth += (double)val * i;
1946 sqr_leaf_depth += (double)val * i * i;
1947 if( val && i < min_leaf_depth ) min_leaf_depth = i;
1948 if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
1949 }
1950 unsigned num_non_leaf = d.count - num_leaves;
1951
1952 double rv = total_dim[0] * total_dim[1] * total_dim[2];
1953 s << "entities in tree: " << total_entities << std::endl
1954 << "root volume: " << rv << std::endl
1955 << "total node volume: " << d.vol.sum << std::endl
1956 << "total/root volume: " << d.vol.sum / rv << std::endl
1957 << "tree height: " << tree_height << std::endl
1958 << "node count: " << d.count << std::endl
1959 << "leaf count: " << num_leaves << std::endl
1960 << std::endl;
1961
1962 double avg_leaf_depth = sum_leaf_depth / num_leaves;
1963 double rms_leaf_depth = sqrt( sqr_leaf_depth / num_leaves );
1964 double std_leaf_depth = std_dev( sqr_leaf_depth, sum_leaf_depth, num_leaves );
1965
1966 double avg_leaf_ent = d.leaf_ent.sum / num_leaves;
1967 double rms_leaf_ent = sqrt( d.leaf_ent.sqr / num_leaves );
1968 double std_leaf_ent = std_dev( d.leaf_ent.sqr, d.leaf_ent.sum, num_leaves );
1969
1970 unsigned num_child = 2 * num_non_leaf;
1971
1972 double avg_vol_ratio = d.volume.sum / num_child;
1973 double rms_vol_ratio = sqrt( d.volume.sqr / num_child );
1974 double std_vol_ratio = std_dev( d.volume.sqr, d.volume.sum, num_child );
1975
1976 double avg_ent_ratio = d.entities.sum / num_child;
1977 double rms_ent_ratio = sqrt( d.entities.sqr / num_child );
1978 double std_ent_ratio = std_dev( d.entities.sqr, d.entities.sum, num_child );
1979
1980 double avg_rad_ratio = d.radius.sum / d.count;
1981 double rms_rad_ratio = sqrt( d.radius.sqr / d.count );
1982 double std_rad_ratio = std_dev( d.radius.sqr, d.radius.sum, d.count );
1983
1984 double avg_vol = d.vol.sum / d.count;
1985 double rms_vol = sqrt( d.vol.sqr / d.count );
1986 double std_vol = std_dev( d.vol.sqr, d.vol.sum, d.count );
1987
1988 double avg_area = d.area.sum / d.count;
1989 double rms_area = sqrt( d.area.sqr / d.count );
1990 double std_area = std_dev( d.area.sqr, d.area.sum, d.count );
1991
1992 int prec = s.precision();
1993 s << " " WW "Minimum" WW "Average" WW "RMS" WW "Maximum" WW "Std.Dev." << std::endl;
1994 s << std::setprecision( 1 )
1995 << "Leaf Depth " WW min_leaf_depth WW avg_leaf_depth WW rms_leaf_depth WW d.leaf_depth.size() -
1996 1 WW std_leaf_depth
1997 << std::endl;
1998 s << std::setprecision( 0 )
1999 << "Entities/Leaf " WW d.leaf_ent.min WW avg_leaf_ent WW rms_leaf_ent WW d.leaf_ent.max WW std_leaf_ent
2000 << std::endl;
2001 s << std::setprecision( 3 )
2002 << "Child Volume Ratio " WW d.volume.min WW avg_vol_ratio WW rms_vol_ratio WW d.volume.max WW std_vol_ratio
2003 << std::endl;
2004 s << std::setprecision( 3 )
2005 << "Child Entity Ratio " WW d.entities.min WW avg_ent_ratio WW rms_ent_ratio WW d.entities.max WW std_ent_ratio
2006 << std::endl;
2007 s << std::setprecision( 3 )
2008 << "Box Radius Ratio " WW d.radius.min WW avg_rad_ratio WW rms_rad_ratio WW d.radius.max WW std_rad_ratio
2009 << std::endl;
2010 s << std::setprecision( 0 ) << "Box volume " WE d.vol.min WE avg_vol WE rms_vol WE d.vol.max WE std_vol
2011 << std::endl;
2012 s << std::setprecision( 0 ) << "Largest side area " WE d.area.min WE avg_area WE rms_area WE d.area.max WE std_area
2013 << std::endl;
2014 s << std::setprecision( prec ) << std::endl;
2015
2016 s << "Leaf Depth Histogram (Root depth is 0)" << std::endl;
2017 double f = 60.0 / max_leaf_per_depth;
2018 for( i = min_leaf_depth; i < d.leaf_depth.size(); ++i )
2019 s << std::setw( 2 ) << i << " " << std::setw( 5 ) << d.leaf_depth[i] << " |" << std::setfill( '*' )
2020 << std::setw( (int)floor( f * d.leaf_depth[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2021 s << std::endl;
2022
2023 s << "Child/Parent Volume Ratio Histogram" << std::endl;
2024 f = 60.0 / *( std::max_element( d.volume.hist, d.volume.hist + 10 ) );
2025 for( i = 0; i < 10u; ++i )
2026 s << "0." << i << " " << std::setw( 5 ) << d.volume.hist[i] << " |" << std::setfill( '*' )
2027 << std::setw( (int)floor( f * d.volume.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2028 s << std::endl;
2029
2030 s << "Child/Parent Entity Count Ratio Histogram" << std::endl;
2031 f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 10 ) );
2032 for( i = 0; i < 10u; ++i )
2033 s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
2034 << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2035 s << std::endl;
2036
2037 s << "Inner/Outer Radius Ratio Histogram (~0.70 for cube)" << std::endl;
2038 // max radius ratio for a box is about 0.7071. Move any boxes
2039 // in the .7 bucket into .6 and print .0 to .6.
2040 d.radius.hist[6] += d.radius.hist[7];
2041 f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 7 ) );
2042 for( i = 0; i < 7u; ++i )
2043 s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
2044 << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2045 s << std::endl;
2046
2047 return MB_SUCCESS;
2048 }
2049
2050 } // namespace moab
1.9.1.