1 #define IS_BUILDING_MB
2 #include "moab/Core.hpp"
3 #include "moab/CartVect.hpp"
4 #include "moab/OrientedBoxTreeTool.hpp"
5 #include "moab/OrientedBox.hpp"
6 #include "Internals.hpp"
7 #include "moab/Range.hpp"
8
9 #include <iostream>
10 #include <iomanip>
11 #include <cstdio>
12 #include <limits>
13 #include <stdlib.h>
14 #include <time.h>
15 #if !defined(_MSC_VER) && !defined(__MINGW32__)
16 # include <unistd.h>
17 # include <sys/stat.h>
18 # include <fcntl.h>
19 #endif
20
21 using namespace moab;
22
23 std::string clock_to_string( clock_t t );
24 std::string mem_to_string( unsigned long mem );
25
26 const int MAX_TAG_VALUE = 20;
27 const char* const TAG_NAME = "OBB_ID";
28 const char* const TREE_TAG = "OBB_ROOT";
29 const char* root_tag = TREE_TAG;
30
31 ErrorCode get_root( Interface* moab, EntityHandle& root );
32 EntityHandle build_tree( Interface* interface,
33 OrientedBoxTreeTool::Settings settings );
34 void delete_existing_tree( Interface* interface );
35 void print_stats( Interface* interface );
36 void tag_triangles( Interface* interface );
37 void tag_vertices( Interface* interface );
38 void write_tree_blocks( Interface* interface, const char* file );
39
usage(bool err=true)40 static void usage( bool err = true )
41 {
42 std::ostream& s = err ? std::cerr : std::cout;
43 s << "obb_tree_tool [-s|-S] [-d <int>] [-n <int>] <input file> <output file>" << std::endl
44 << "obb_tree_tool [-h]" << std::endl;
45 if (!err) {
46 OrientedBoxTreeTool::Settings st;
47 s << "Tool to build adaptive kd-Tree from triangles" << std::endl;
48 s << " -s Use range-based sets for tree nodes" << std::endl
49 << " -S Use vector-based sets for tree nodes" << std::endl
50 << " -d <int> Specify maximum depth for tree. Default: " << st.max_depth << std::endl
51 << " -n <int> Specify maximum entities per leaf. Default: " << st.max_leaf_entities << std::endl
52 << " -m <real> Specify worst split ratio. Default: " << st.worst_split_ratio << std::endl
53 << " -M <real> Specify best split ratio. Default: " << st.best_split_ratio << std::endl
54 << " -t Tag triangles will tree cell number." << std::endl
55 << " -T Write tree boxes to file." << std::endl
56 << " -N Specify mesh tag containing tree root. Default: \"" << TREE_TAG << '"' << std::endl
57 << std::endl;
58 }
59 exit( err );
60 }
61
62 #if defined(_MSC_VER) || defined(__MINGW32__)
memory_use(unsigned long long & vsize,unsigned long long & rss)63 static void memory_use( unsigned long long& vsize, unsigned long long& rss )
64 { vsize = rss = 0; }
65 #else
memory_use(unsigned long long & vsize,unsigned long long & rss)66 static void memory_use( unsigned long long& vsize, unsigned long long& rss )
67 {
68 char buffer[512];
69 unsigned long lvsize;
70 long lrss;
71 int filp = open( "/proc/self/stat", O_RDONLY );
72 ssize_t r = read( filp, buffer, sizeof(buffer)-1 );
73 close( filp );
74 if (r < 0) r = 0;
75 lvsize = lrss = 0;
76 buffer[r] = '\0';
77 sscanf( buffer, "%*d %*s %*c " // pid command state
78 "%*d %*d " // ppid pgrp
79 "%*d %*d %*d " // session tty_nr tpgid
80 "%*u " // flags
81 "%*u %*u %*u %*u " // minflt cminflt majflt cmajflt
82 "%*u %*u %*d %*d " // utime stime cutime cstime
83 "%*d %*d %*d " // priority nice (unused)
84 "%*d %*u " // itrealval starttime
85 "%lu %ld", &lvsize, &lrss );
86 rss = lrss*getpagesize();
87 vsize = lvsize;
88 }
89 #endif
90
parseint(int & i,int argc,char * argv[])91 static int parseint( int& i, int argc, char* argv[] )
92 {
93 char* end;
94 ++i;
95 if (i == argc) {
96 std::cerr << "Expected value following '" << argv[i-1] << "'" << std::endl;
97 usage();
98 }
99
100 int result = strtol( argv[i], &end, 0 );
101 if (result < 0 || *end) {
102 std::cerr << "Expected positive integer following '" << argv[i-1] << "'" << std::endl;
103 usage();
104 }
105
106 return result;
107 }
108
parsedouble(int & i,int argc,char * argv[])109 static double parsedouble( int& i, int argc, char* argv[] )
110 {
111 char* end;
112 ++i;
113 if (i == argc) {
114 std::cerr << "Expected value following '" << argv[i-1] << "'" << std::endl;
115 usage();
116 }
117
118 double result = strtod( argv[i], &end );
119 if (result < 0 || *end) {
120 std::cerr << "Expected positive real number following '" << argv[i-1] << "'" << std::endl;
121 usage();
122 }
123
124 return result;
125 }
126
127
main(int argc,char * argv[])128 int main( int argc, char* argv[] )
129 {
130 const char* input_file = 0;
131 const char* output_file = 0;
132 const char* tree_file = 0;
133 OrientedBoxTreeTool::Settings settings;
134 bool tag_tris = false;
135 clock_t load_time, build_time, stat_time, tag_time, write_time, block_time;
136
137 for (int i = 1; i < argc; ++i) {
138 if (argv[i][0] != '-') {
139 if (!input_file)
140 input_file = argv[i];
141 else if (!output_file)
142 output_file = argv[i];
143 else
144 usage();
145 continue;
146 }
147
148 if (!argv[i][1] || argv[i][2])
149 usage();
150
151 switch (argv[i][1]) {
152 case 's': settings.set_options = MESHSET_SET; break;
153 case 'S': settings.set_options = MESHSET_ORDERED; break;
154 case 'd': settings.max_depth = parseint( i, argc, argv ); break;
155 case 'n': settings.max_leaf_entities = parseint( i, argc, argv ); break;
156 case 'm': settings.worst_split_ratio = parsedouble( i, argc, argv ); break;
157 case 'M': settings.best_split_ratio = parsedouble( i, argc, argv ); break;
158 case 't': tag_tris = true; break;
159 case 'T': if (++i == argc) usage(); tree_file = argv[i]; break;
160 case 'N': if (++i == argc) usage(); root_tag = argv[i]; break;
161 case 'h': usage(false);
162 default: usage();
163 }
164 }
165
166 if (!output_file)
167 usage();
168
169 ErrorCode rval;
170 Core moab_core;
171 Interface* interface = &moab_core;
172
173 load_time = clock();
174 rval = interface->load_mesh( input_file );
175 if (MB_SUCCESS != rval) {
176 std::cerr << "Error reading file: " << input_file << std::endl;
177 exit(2);
178 }
179 load_time = clock() - load_time;
180
181
182 delete_existing_tree( interface );
183
184 std::cout << "Building tree..." << std::endl;
185 build_time = clock();
186 build_tree( interface, settings );
187 build_time = clock() - build_time;
188
189 std::cout << "Calculating stats..." << std::endl;
190 print_stats( interface );
191 stat_time = clock() - build_time;
192
193 if (tag_tris) {
194 std::cout << "Tagging tree..." << std::endl;
195 tag_triangles( interface );
196 tag_vertices( interface );
197 }
198 tag_time = clock() - stat_time;
199
200
201 std::cout << "Writing file... "; std::cout.flush();
202 rval = interface->write_mesh( output_file );
203 if (MB_SUCCESS != rval) {
204 std::cerr << "Error writing file: " << output_file << std::endl;
205 exit(3);
206 }
207 write_time = clock() - tag_time;
208 std::cout << "Wrote " << output_file << std::endl;
209
210 if (tree_file) {
211 std::cout << "Writing tree block rep..."; std::cout.flush();
212 write_tree_blocks( interface, tree_file );
213 std::cout << "Wrote " << tree_file << std::endl;
214 }
215 block_time = clock() - write_time;
216
217 std::cout << "Times: "
218 << " Load"
219 << " Build"
220 << " Stats"
221 << " Write";
222 if (tag_tris)
223 std::cout << "Tag Sets";
224 if (tree_file)
225 std::cout << "Block ";
226 std::cout << std::endl;
227
228 std::cout << " "
229 << std::setw(8) << clock_to_string(load_time)
230 << std::setw(8) << clock_to_string(build_time)
231 << std::setw(8) << clock_to_string(stat_time)
232 << std::setw(8) << clock_to_string(write_time);
233 if (tag_tris)
234 std::cout << std::setw(8) << clock_to_string(tag_time);
235 if (tree_file)
236 std::cout << std::setw(8) << clock_to_string(block_time);
237 std::cout << std::endl;
238
239 return 0;
240 }
241
get_root(Interface * moab,EntityHandle & root)242 ErrorCode get_root( Interface* moab, EntityHandle& root )
243 {
244 Tag tag;
245 ErrorCode rval;
246
247 rval = moab->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, tag );
248 if (MB_SUCCESS != rval)
249 return rval;
250
251 const EntityHandle mesh = 0;
252 return moab->tag_get_data( tag, &mesh, 1, &root );
253 }
254
255
delete_existing_tree(Interface * interface)256 void delete_existing_tree( Interface* interface )
257 {
258 EntityHandle root;
259 ErrorCode rval = get_root(interface, root);
260 if (MB_SUCCESS == rval) {
261 OrientedBoxTreeTool tool(interface);
262 rval = tool.delete_tree( root );
263 if (MB_SUCCESS != rval) {
264 std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
265 exit( 5 );
266 }
267 }
268 }
269
build_tree(Interface * interface,OrientedBoxTreeTool::Settings settings)270 EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings )
271 {
272 ErrorCode rval;
273 EntityHandle root = 0;
274 Range triangles;
275
276 rval = interface->get_entities_by_type( 0, MBTRI, triangles );
277 if (MB_SUCCESS != rval || triangles.empty()) {
278 std::cerr << "No triangles from which to build tree." << std::endl;
279 exit(4);
280 }
281
282 OrientedBoxTreeTool tool( interface );
283 rval = tool.build( triangles, root, &settings );
284 if (MB_SUCCESS != rval || !root) {
285 std::cerr << "Tree construction failed." << std::endl;
286 exit(4);
287 }
288
289 // store tree root
290 Tag roottag;
291 rval = interface->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, roottag, MB_TAG_CREAT|MB_TAG_SPARSE );
292 if (MB_SUCCESS != rval) {
293 std::cout << "Failed to create root tag: \"" << root_tag << '"' << std::endl;
294 exit(2);
295 }
296 const EntityHandle mesh = 0;
297 rval = interface->tag_set_data( roottag, &mesh, 1, &root );
298 if (MB_SUCCESS != rval) {
299 std::cout << "Failed to set root tag: \"" << root_tag << '"' << std::endl;
300 exit(2);
301 }
302
303 return root;
304 }
305
clock_to_string(clock_t t)306 std::string clock_to_string( clock_t t )
307 {
308 char unit[5] = "s";
309 char buffer[256];
310 double dt = t / (double)CLOCKS_PER_SEC;
311 //if (dt > 300) {
312 // dt /= 60;
313 // strcpy( unit, "min" );
314 //}
315 //if (dt > 300) {
316 // dt /= 60;
317 // strcpy( unit, "hr" );
318 //}
319 //if (dt > 100) {
320 // dt /= 24;
321 // strcpy( unit, "days" );
322 //}
323 sprintf(buffer,"%0.2f%s",dt,unit);
324 return buffer;
325 }
326
mem_to_string(unsigned long mem)327 std::string mem_to_string( unsigned long mem )
328 {
329 char unit[3] = "B";
330 if (mem > 9*1024) {
331 mem = (mem + 512) / 1024;
332 strcpy( unit, "kB" );
333 }
334 if (mem > 9*1024) {
335 mem = (mem + 512) / 1024;
336 strcpy( unit, "MB" );
337 }
338 if (mem > 9*1024) {
339 mem = (mem + 512) / 1024;
340 strcpy( unit, "GB" );
341 }
342 char buffer[256];
343 sprintf(buffer, "%lu %s", mem, unit );
344 return buffer;
345 }
346
347 template <typename T>
348 struct SimpleStat
349 {
350 T min, max, sum, sqr;
351 size_t count;
352 SimpleStat();
353 void add( T value );
avgSimpleStat354 double avg() const { return (double)sum / count; }
rmsSimpleStat355 double rms() const { return sqrt( (double)sqr / count ); }
devSimpleStat356 double dev() const { return sqrt( (count * (double)sqr - (double)sum * (double)sum) / ((double)count * (count - 1) ) ); }
357 };
358
SimpleStat()359 template <typename T> SimpleStat<T>::SimpleStat()
360 : min( std::numeric_limits<T>::max() ),
361 max( std::numeric_limits<T>::min() ),
362 sum( 0 ), sqr( 0 ), count( 0 )
363 {}
364
print_stats(Interface * interface)365 void print_stats( Interface* interface )
366 {
367 EntityHandle root;
368 Range range;
369 ErrorCode rval;
370 rval = get_root( interface, root );
371 if (MB_SUCCESS != rval) {
372 std::cerr << "Internal error: Failed to retrieve root." << std::endl;
373 exit(5);
374 }
375 OrientedBoxTreeTool tool(interface);
376
377 Range tree_sets, triangles, verts;
378 //interface->get_child_meshsets( root, tree_sets, 0 );
379 interface->get_entities_by_type( 0, MBENTITYSET, tree_sets );
380 tree_sets.erase( tree_sets.begin(), Range::lower_bound( tree_sets.begin(), tree_sets.end(), root ) );
381 interface->get_entities_by_type( 0, MBTRI, triangles );
382 interface->get_entities_by_type( 0, MBVERTEX, verts );
383 triangles.merge( verts );
384 tree_sets.insert( root );
385 unsigned long long set_used, set_amortized, set_store_used, set_store_amortized,
386 set_tag_used, set_tag_amortized, tri_used, tri_amortized;
387 interface->estimated_memory_use( tree_sets,
388 &set_used, &set_amortized,
389 &set_store_used, &set_store_amortized,
390 0, 0, 0, 0,
391 &set_tag_used, &set_tag_amortized );
392 interface->estimated_memory_use( triangles, &tri_used, &tri_amortized );
393
394 int num_tri = 0;
395 interface->get_number_entities_by_type( 0, MBTRI, num_tri );
396
397 tool.stats( root, std::cout );
398
399 unsigned long long real_rss, real_vsize;
400 memory_use( real_vsize, real_rss );
401
402 printf("------------------------------------------------------------------\n");
403 printf("\nmemory: used amortized\n");
404 printf(" ---------- ----------\n");
405 printf("triangles %10s %10s\n",mem_to_string(tri_used).c_str(), mem_to_string(tri_amortized).c_str());
406 printf("sets (total)%10s %10s\n",mem_to_string(set_used).c_str(), mem_to_string(set_amortized).c_str());
407 printf("sets %10s %10s\n",mem_to_string(set_store_used).c_str(), mem_to_string(set_store_amortized).c_str());
408 printf("set tags %10s %10s\n",mem_to_string(set_tag_used).c_str(), mem_to_string(set_tag_amortized).c_str());
409 printf("total real %10s %10s\n",mem_to_string(real_rss).c_str(), mem_to_string(real_vsize).c_str());
410 printf("------------------------------------------------------------------\n");
411 }
412
413
hash_handle(EntityHandle handle)414 static int hash_handle( EntityHandle handle )
415 {
416 EntityID h = ID_FROM_HANDLE(handle);
417 return (int)((h * 13 + 7) % MAX_TAG_VALUE) + 1;
418 }
419
420 class TriTagger : public OrientedBoxTreeTool::Op
421 {
422 private:
423 Interface* mMB;
424 Tag mTag;
425 std::vector<EntityHandle> mHandles;
426 std::vector<int> mTagData;
427 public:
TriTagger(Tag tag,Interface * moab)428 TriTagger( Tag tag, Interface* moab )
429 : mMB(moab), mTag(tag) {}
430
visit(EntityHandle,int,bool & descent)431 ErrorCode visit( EntityHandle, int, bool& descent )
432 { descent = true; return MB_SUCCESS; }
433
leaf(EntityHandle node)434 ErrorCode leaf( EntityHandle node ) {
435 mHandles.clear();
436 mMB->get_entities_by_handle( node, mHandles );
437 mTagData.clear();
438 mTagData.resize( mHandles.size(), hash_handle( node ) );
439 mMB->tag_set_data( mTag, &mHandles[0], mHandles.size(), &mTagData[0] );
440 return MB_SUCCESS;
441 }
442 };
443
444
tag_triangles(Interface * moab)445 void tag_triangles( Interface* moab )
446 {
447 EntityHandle root;
448 ErrorCode rval = get_root( moab, root );
449 if (MB_SUCCESS != rval) {
450 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
451 exit(5);
452 }
453
454 Tag tag;
455 int zero = 0;
456 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE|MB_TAG_CREAT, &zero );
457 TriTagger op( tag, moab );
458
459 OrientedBoxTreeTool tool(moab);
460 rval = tool.preorder_traverse( root, op );
461 if (MB_SUCCESS != rval) {
462 std::cerr << "Internal error tagging triangles" << std::endl;
463 exit(5);
464 }
465 }
466
467
468 class VtxTagger : public OrientedBoxTreeTool::Op
469 {
470 private:
471 Interface* mMB;
472 Tag mTag;
473 std::vector<EntityHandle> mHandles;
474 std::vector<EntityHandle> mConn;
475 std::vector<int> mTagData;
476 public:
VtxTagger(Tag tag,Interface * moab)477 VtxTagger( Tag tag, Interface* moab )
478 : mMB(moab), mTag(tag) {}
479
visit(EntityHandle,int,bool & descent)480 ErrorCode visit( EntityHandle, int, bool& descent )
481 { descent = true; return MB_SUCCESS; }
482
leaf(EntityHandle node)483 ErrorCode leaf( EntityHandle node ) {
484 mHandles.clear();
485 mMB->get_entities_by_handle( node, mHandles );
486 mConn.clear();
487 mMB->get_connectivity( &mHandles[0], mHandles.size(), mConn );
488 mTagData.clear();
489 mTagData.resize( mConn.size(), hash_handle( node ) );
490 mMB->tag_set_data( mTag, &mConn[0], mConn.size(), &mTagData[0] );
491 return MB_SUCCESS;
492 }
493 };
494
495
tag_vertices(Interface * moab)496 void tag_vertices( Interface* moab )
497 {
498 EntityHandle root;
499 ErrorCode rval = get_root( moab, root );
500 if (MB_SUCCESS != rval) {
501 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
502 exit(5);
503 }
504
505 Tag tag;
506 int zero = 0;
507 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE|MB_TAG_CREAT, &zero );
508 VtxTagger op( tag, moab );
509
510 OrientedBoxTreeTool tool(moab);
511 rval = tool.preorder_traverse( root, op );
512 if (MB_SUCCESS != rval) {
513 std::cerr << "Internal error tagging vertices" << std::endl;
514 exit(5);
515 }
516 }
517
518
519 class LeafHexer : public OrientedBoxTreeTool::Op
520 {
521 private:
522 OrientedBoxTreeTool* mTool;
523 Interface* mOut;
524 Tag mTag;
525 std::vector<EntityHandle> mHandles;
526 std::vector<EntityHandle> mConn;
527 std::vector<int> mTagData;
528 public:
LeafHexer(OrientedBoxTreeTool * tool,Interface * mb2,Tag tag)529 LeafHexer( OrientedBoxTreeTool* tool, Interface* mb2, Tag tag )
530 : mTool(tool), mOut( mb2 ), mTag(tag) {}
531
visit(EntityHandle,int,bool & descent)532 ErrorCode visit( EntityHandle, int, bool& descent )
533 { descent = true; return MB_SUCCESS; }
534
leaf(EntityHandle node)535 ErrorCode leaf( EntityHandle node ) {
536 OrientedBox box;
537 ErrorCode rval = mTool->box( node, box );
538 if (MB_SUCCESS !=rval) return rval;
539 EntityHandle h;
540 rval = box.make_hex( h, mOut );
541 if (MB_SUCCESS !=rval) return rval;
542 int i = hash_handle( node );
543 return mOut->tag_set_data( mTag, &h, 1, &i );
544 }
545 };
546
write_tree_blocks(Interface * interface,const char * file)547 void write_tree_blocks( Interface* interface, const char* file )
548 {
549 EntityHandle root;
550 ErrorCode rval = get_root( interface, root );
551 if (MB_SUCCESS != rval) {
552 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
553 exit(5);
554 }
555
556 Core moab2;
557 Tag tag;
558 int zero = 0;
559 moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE|MB_TAG_CREAT, &zero );
560
561 OrientedBoxTreeTool tool(interface);
562 LeafHexer op( &tool, &moab2, tag );
563 rval = tool.preorder_traverse( root, op );
564 if (MB_SUCCESS != rval) {
565 std::cerr << "Internal error: failed to construct leaf hexes" << std::endl;
566 exit(5);
567 }
568
569 moab2.write_mesh( file );
570 }
571
572