MOAB: Mesh Oriented datABase  (version 5.4.0)
kd_tree_tool.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/CartVect.hpp"
00003 #include "moab/AdaptiveKDTree.hpp"
00004 #include "moab/GeomUtil.hpp"
00005 #include "Internals.hpp"
00006 #include "moab/Range.hpp"
00007 #include "moab/FileOptions.hpp"
00008 #include <iostream>
00009 #include <iomanip>
00010 #include <cstdio>
00011 #include <limits>
00012 #include <sstream>
00013 #include <cstdlib>
00014 #include <ctime>
00015 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00016 #include <unistd.h>
00017 #include <sys/stat.h>
00018 #include <fcntl.h>
00019 #endif
00020 
00021 using namespace moab;
00022 
00023 const int MAX_TAG_VALUE = 32;  // maximum value to use when tagging entities with tree cell number
00024 const char* TAG_NAME    = "TREE_CELL";
00025 
00026 std::string clock_to_string( clock_t t );
00027 std::string mem_to_string( unsigned long mem );
00028 
00029 void build_tree( Interface* interface, const Range& elems, FileOptions& opts );
00030 void build_grid( Interface* iface, const double dims[3] );
00031 void delete_existing_tree( Interface* interface );
00032 void print_stats( Interface* interface );
00033 void tag_elements( Interface* interface );
00034 void tag_vertices( Interface* interface );
00035 void write_tree_blocks( Interface* interface, const char* file );
00036 
00037 static void usage( bool err = true )
00038 {
00039     std::ostream& s = err ? std::cerr : std::cout;
00040     s << "kd_tree_tool [-d <int>] [-2|-3] [-n <int>] [-u|-p|-m|-v] [-N <int>] [-s|-S] <input file> "
00041          "<output file>"
00042       << std::endl
00043       << "kd_tree_tool [-d <int>] -G <dims> [-s|-S] <output file>" << std::endl
00044       << "kd_tree_tool [-h]" << std::endl;
00045     if( !err )
00046     {
00047         s << "Tool to build adaptive kd-Tree" << std::endl;
00048         s << "  -d <int>  Specify maximum depth for tree. Default: 30" << std::endl
00049           << "  -n <int>  Specify maximum entities per leaf. Default: 6" << std::endl
00050           << "  -2        Build tree from surface elements. Default: yes" << std::endl
00051           << "  -3        Build tree from volume elements. Default: yes" << std::endl
00052           << "  -u        Use 'SUBDIVISION' scheme for tree construction" << std::endl
00053           << "  -p        Use 'SUBDIVISION_SNAP' tree construction algorithm." << std::endl
00054           << "  -m        Use 'VERTEX_MEDIAN' tree construction algorithm." << std::endl
00055           << "  -v        Use 'VERTEX_SAMPLE' tree construction algorithm." << std::endl
00056           << "  -N <int>  Specify candidate split planes per axis.  Default: 3" << std::endl
00057           << "  -t        Tag elements will tree cell number." << std::endl
00058           << "  -T        Write tree boxes to file." << std::endl
00059           << "  -G <dims> Generate grid - no input elements.  Dims must be " << std::endl
00060           << "            HxWxD or H." << std::endl
00061           << "  -s        Use range-based sets for tree nodes" << std::endl
00062           << "  -S        Use vector-based sets for tree nodes" << std::endl
00063           << std::endl;
00064     }
00065     exit( err );
00066 }
00067 
00068 /*
00069 #if !defined(_MSC_VER) && !defined(__MINGW32__)
00070 static void memory_use( unsigned long& vsize, unsigned long& rss )
00071 {
00072   char buffer[512];
00073   int filp = open( "/proc/self/stat", O_RDONLY );
00074   ssize_t r = read( filp, buffer, sizeof(buffer)-1 );
00075   close( filp );
00076   if (r < 0) r = 0;
00077   vsize = rss = 0;
00078   buffer[r] = '\0';
00079   sscanf( buffer, "%*d %*s %*c "         // pid command state
00080                   "%*d %*d "             // ppid pgrp
00081                   "%*d %*d %*d "         // session tty_nr tpgid
00082                   "%*u "                // flags
00083                   "%*u %*u %*u %*u " // minflt cminflt majflt cmajflt
00084                   "%*u %*u %*d %*d " // utime stime cutime cstime
00085                   "%*d %*d %*d "      // priority nice (unused)
00086                   "%*d %*u "           // itrealval starttime
00087                   "%lu %lu",             &vsize, &rss );
00088   rss *= getpagesize();
00089 }
00090 #else
00091 static void memory_use( unsigned long& vsize, unsigned long& rss )
00092   { vsize = rss = 0; }
00093 #endif
00094 */
00095 
00096 static int parseint( int& i, int argc, char* argv[] )
00097 {
00098     char* end;
00099     ++i;
00100     if( i == argc )
00101     {
00102         std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
00103         usage();
00104     }
00105 
00106     int result = strtol( argv[i], &end, 0 );
00107     if( result < 0 || *end )
00108     {
00109         std::cerr << "Expected positive integer following '" << argv[i - 1] << "'" << std::endl;
00110         usage();
00111     }
00112 
00113     return result;
00114 }
00115 
00116 static void parsedims( int& i, int argc, char* argv[], double dims[3] )
00117 {
00118     char* end;
00119     ++i;
00120     if( i == argc )
00121     {
00122         std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
00123         usage();
00124     }
00125 
00126     dims[0] = strtod( argv[i], &end );
00127     if( *end == '\0' )
00128     {
00129         dims[2] = dims[1] = dims[0];
00130         return;
00131     }
00132     else if( *end != 'x' && *end != 'X' )
00133         goto error;
00134 
00135     ++end;
00136     dims[1] = strtod( end, &end );
00137     if( *end == '\0' )
00138     {
00139         dims[2] = dims[1];
00140         return;
00141     }
00142     else if( *end != 'x' && *end != 'X' )
00143         goto error;
00144 
00145     ++end;
00146     dims[2] = strtod( end, &end );
00147     if( *end != '\0' ) goto error;
00148 
00149     return;
00150 error:
00151     std::cerr << "Invaild dimension specification." << std::endl;
00152     usage();
00153 }
00154 
00155 int main( int argc, char* argv[] )
00156 {
00157     bool surf_elems = false, vol_elems = false;
00158     const char* input_file  = 0;
00159     const char* output_file = 0;
00160     const char* tree_file   = 0;
00161     std::ostringstream options;
00162     bool tag_elems = false;
00163     clock_t load_time, build_time, stat_time, tag_time, write_time, block_time;
00164     bool make_grid = false;
00165     double dims[3];
00166 
00167     for( int i = 1; i < argc; ++i )
00168     {
00169         if( argv[i][0] != '-' )
00170         {
00171             if( !input_file )
00172                 input_file = argv[i];
00173             else if( !output_file )
00174                 output_file = argv[i];
00175             else
00176                 usage();
00177             continue;
00178         }
00179 
00180         if( !argv[i][1] || argv[i][2] ) usage();
00181 
00182         switch( argv[i][1] )
00183         {
00184             case '2':
00185                 surf_elems = true;
00186                 break;
00187             case '3':
00188                 vol_elems = true;
00189                 break;
00190             case 'S':
00191                 options << "MESHSET_FLAGS=" << MESHSET_ORDERED << ";";
00192                 break;
00193             case 's':
00194                 break;
00195             case 'd':
00196                 options << "MAX_DEPTH=" << parseint( i, argc, argv ) << ";";
00197                 break;
00198             case 'n':
00199                 options << "MAX_PER_LEAF=" << parseint( i, argc, argv ) << ";";
00200                 break;
00201             case 'u':
00202                 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::SUBDIVISION << ";";
00203                 break;
00204             case 'p':
00205                 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::SUBDIVISION_SNAP << ";";
00206                 break;
00207             case 'm':
00208                 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::VERTEX_MEDIAN << ";";
00209                 break;
00210             case 'v':
00211                 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::VERTEX_SAMPLE << ";";
00212                 break;
00213             case 'N':
00214                 options << "SPLITS_PER_DIR=" << parseint( i, argc, argv ) << ";";
00215                 break;
00216             case 't':
00217                 tag_elems = true;
00218                 break;
00219             case 'T':
00220                 tree_file = argv[++i];
00221                 break;
00222             case 'G':
00223                 make_grid = true;
00224                 parsedims( i, argc, argv, dims );
00225                 break;
00226             case 'h':
00227                 usage( false );
00228                 break;
00229             default:
00230                 usage();
00231         }
00232     }
00233 
00234     // this test relies on not cleaning up trees
00235     options << "CLEAN_UP=false;";
00236 
00237     if( make_grid != !output_file ) usage();
00238 
00239     // default to both
00240     if( !surf_elems && !vol_elems ) surf_elems = vol_elems = true;
00241 
00242     ErrorCode rval;
00243     Core moab_core;
00244     Interface* interface = &moab_core;
00245     FileOptions opts( options.str().c_str() );
00246 
00247     if( make_grid )
00248     {
00249         load_time   = 0;
00250         output_file = input_file;
00251         input_file  = 0;
00252         build_time  = clock();
00253         build_grid( interface, dims );
00254         build_time = clock() - build_time;
00255     }
00256     else
00257     {
00258         load_time = clock();
00259         rval      = interface->load_mesh( input_file );
00260         if( MB_SUCCESS != rval )
00261         {
00262             std::cerr << "Error reading file: " << input_file << std::endl;
00263             exit( 2 );
00264         }
00265         load_time = clock() - load_time;
00266 
00267         delete_existing_tree( interface );
00268 
00269         std::cout << "Building tree..." << std::endl;
00270         build_time = clock();
00271         Range elems;
00272         if( !surf_elems )
00273         {
00274             interface->get_entities_by_dimension( 0, 3, elems );
00275         }
00276         else
00277         {
00278             interface->get_entities_by_dimension( 0, 2, elems );
00279             if( vol_elems )
00280             {
00281                 Range tmp;
00282                 interface->get_entities_by_dimension( 0, 3, tmp );
00283                 elems.merge( tmp );
00284             }
00285         }
00286 
00287         build_tree( interface, elems, opts );
00288         build_time = clock() - build_time;
00289     }
00290 
00291     std::cout << "Calculating stats..." << std::endl;
00292     AdaptiveKDTree tool( interface );
00293     Range range;
00294     tool.find_all_trees( range );
00295     if( range.size() != 1 )
00296     {
00297         if( range.empty() )
00298             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00299         else
00300             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00301         exit( 5 );
00302     }
00303     tool.print();
00304 
00305     stat_time = clock() - build_time;
00306 
00307     if( tag_elems )
00308     {
00309         std::cout << "Tagging tree..." << std::endl;
00310         tag_elements( interface );
00311         tag_vertices( interface );
00312     }
00313     tag_time = clock() - stat_time;
00314 
00315     std::cout << "Writing file... ";
00316     std::cout.flush();
00317     rval = interface->write_mesh( output_file );
00318     if( MB_SUCCESS != rval )
00319     {
00320         std::cerr << "Error writing file: " << output_file << std::endl;
00321         exit( 3 );
00322     }
00323     write_time = clock() - tag_time;
00324     std::cout << "Wrote " << output_file << std::endl;
00325 
00326     if( tree_file )
00327     {
00328         std::cout << "Writing tree block rep...";
00329         std::cout.flush();
00330         write_tree_blocks( interface, tree_file );
00331         std::cout << "Wrote " << tree_file << std::endl;
00332     }
00333     block_time = clock() - write_time;
00334 
00335     std::cout << "Times:  "
00336               << "    Load"
00337               << "   Build"
00338               << "   Stats"
00339               << "   Write";
00340     if( tag_elems ) std::cout << "Tag Sets";
00341     if( tree_file ) std::cout << "Block   ";
00342     std::cout << std::endl;
00343 
00344     std::cout << "        " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 )
00345               << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 )
00346               << clock_to_string( write_time );
00347     if( tag_elems ) std::cout << std::setw( 8 ) << clock_to_string( tag_time );
00348     if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time );
00349     std::cout << std::endl;
00350 
00351     return 0;
00352 }
00353 
00354 void delete_existing_tree( Interface* interface )
00355 {
00356     Range trees;
00357     AdaptiveKDTree tool( interface );
00358 
00359     tool.find_all_trees( trees );
00360     if( !trees.empty() ) std::cout << "Destroying existing tree(s) contained in file" << std::endl;
00361     for( Range::iterator i = trees.begin(); i != trees.end(); ++i )
00362         tool.reset_tree();
00363 
00364     trees.clear();
00365     tool.find_all_trees( trees );
00366     if( !trees.empty() )
00367     {
00368         std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
00369         exit( 5 );
00370     }
00371 }
00372 
00373 void build_tree( Interface* interface, const Range& elems, FileOptions& opts )
00374 {
00375     EntityHandle root = 0;
00376 
00377     if( elems.empty() )
00378     {
00379         std::cerr << "No elements from which to build tree." << std::endl;
00380         exit( 4 );
00381     }
00382 
00383     AdaptiveKDTree tool( interface, elems, &root, &opts );
00384 }
00385 
00386 void build_grid( Interface* interface, const double dims[3] )
00387 {
00388     ErrorCode rval;
00389     EntityHandle root = 0;
00390     AdaptiveKDTree tool( interface );
00391     AdaptiveKDTreeIter iter;
00392     AdaptiveKDTree::Plane plane;
00393 
00394     double min[3] = { -0.5 * dims[0], -0.5 * dims[1], -0.5 * dims[2] };
00395     double max[3] = { 0.5 * dims[0], 0.5 * dims[1], 0.5 * dims[2] };
00396     rval          = tool.create_root( min, max, root );
00397     if( MB_SUCCESS != rval || !root )
00398     {
00399         std::cerr << "Failed to create tree root." << std::endl;
00400         exit( 4 );
00401     }
00402 
00403     rval = tool.get_tree_iterator( root, iter );
00404     if( MB_SUCCESS != rval )
00405     {
00406         std::cerr << "Failed to get tree iterator." << std::endl;
00407     }
00408 
00409     do
00410     {
00411         while( iter.depth() < tool.get_max_depth() )
00412         {
00413             plane.norm  = iter.depth() % 3;
00414             plane.coord = 0.5 * ( iter.box_min()[plane.norm] + iter.box_max()[plane.norm] );
00415             rval        = tool.split_leaf( iter, plane );
00416             if( MB_SUCCESS != rval )
00417             {
00418                 std::cerr << "Failed to split tree leaf at depth " << iter.depth() << std::endl;
00419                 exit( 4 );
00420             }
00421         }
00422     } while( ( rval = iter.step() ) == MB_SUCCESS );
00423 
00424     if( rval != MB_ENTITY_NOT_FOUND )
00425     {
00426         std::cerr << "Error stepping KDTree iterator." << std::endl;
00427         exit( 4 );
00428     }
00429 }
00430 
00431 std::string clock_to_string( clock_t t )
00432 {
00433     char unit[5] = "s";
00434     char buffer[256];
00435     double dt = t / (double)CLOCKS_PER_SEC;
00436     // if (dt > 300) {
00437     //  dt /= 60;
00438     //  strcpy( unit, "min" );
00439     //}
00440     // if (dt > 300) {
00441     //  dt /= 60;
00442     //  strcpy( unit, "hr" );
00443     //}
00444     // if (dt > 100) {
00445     //  dt /= 24;
00446     //  strcpy( unit, "days" );
00447     //}
00448     sprintf( buffer, "%0.2f%s", dt, unit );
00449     return buffer;
00450 }
00451 
00452 static int hash_handle( EntityHandle handle )
00453 {
00454     EntityID h = ID_FROM_HANDLE( handle );
00455     return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1;
00456 }
00457 
00458 void tag_elements( Interface* moab )
00459 {
00460     EntityHandle root;
00461     Range range;
00462     AdaptiveKDTree tool( moab );
00463 
00464     tool.find_all_trees( range );
00465     if( range.size() != 1 )
00466     {
00467         if( range.empty() )
00468             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00469         else
00470             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00471         exit( 5 );
00472     }
00473 
00474     root = *range.begin();
00475     range.clear();
00476 
00477     Tag tag;
00478     int zero = 0;
00479     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00480 
00481     AdaptiveKDTreeIter iter;
00482     tool.get_tree_iterator( root, iter );
00483 
00484     std::vector< int > tag_vals;
00485     do
00486     {
00487         range.clear();
00488         moab->get_entities_by_handle( iter.handle(), range );
00489         tag_vals.clear();
00490         tag_vals.resize( range.size(), hash_handle( iter.handle() ) );
00491         moab->tag_set_data( tag, range, &tag_vals[0] );
00492     } while( MB_SUCCESS == iter.step() );
00493 }
00494 
00495 void tag_vertices( Interface* moab )
00496 {
00497     EntityHandle root;
00498     Range range;
00499     AdaptiveKDTree tool( moab );
00500 
00501     tool.find_all_trees( range );
00502     if( range.size() != 1 )
00503     {
00504         if( range.empty() )
00505             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00506         else
00507             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00508         exit( 5 );
00509     }
00510 
00511     root = *range.begin();
00512     range.clear();
00513 
00514     Tag tag;
00515     int zero = 0;
00516     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00517 
00518     AdaptiveKDTreeIter iter;
00519     tool.get_tree_iterator( root, iter );
00520 
00521     do
00522     {
00523         range.clear();
00524         moab->get_entities_by_handle( iter.handle(), range );
00525 
00526         int tag_val = hash_handle( iter.handle() );
00527         Range verts;
00528         moab->get_adjacencies( range, 0, false, verts, Interface::UNION );
00529         for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
00530         {
00531             CartVect coords;
00532             moab->get_coords( &*i, 1, coords.array() );
00533             if( GeomUtil::box_point_overlap( CartVect( iter.box_min() ), CartVect( iter.box_max() ), coords, 1e-7 ) )
00534                 moab->tag_set_data( tag, &*i, 1, &tag_val );
00535         }
00536     } while( MB_SUCCESS == iter.step() );
00537 }
00538 
00539 void write_tree_blocks( Interface* interface, const char* file )
00540 {
00541     EntityHandle root;
00542     Range range;
00543     AdaptiveKDTree tool( interface );
00544 
00545     tool.find_all_trees( range );
00546     if( range.size() != 1 )
00547     {
00548         if( range.empty() )
00549             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00550         else
00551             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00552         exit( 5 );
00553     }
00554 
00555     root = *range.begin();
00556     range.clear();
00557 
00558     Core moab2;
00559     Tag tag;
00560     int zero = 0;
00561     moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00562 
00563     AdaptiveKDTreeIter iter;
00564     tool.get_tree_iterator( root, iter );
00565 
00566     do
00567     {
00568         double x1         = iter.box_min()[0];
00569         double y1         = iter.box_min()[1];
00570         double z1         = iter.box_min()[2];
00571         double x2         = iter.box_max()[0];
00572         double y2         = iter.box_max()[1];
00573         double z2         = iter.box_max()[2];
00574         double coords[24] = { x1, y1, z1, x2, y1, z1, x2, y2, z1, x1, y2, z1,
00575                               x1, y1, z2, x2, y1, z2, x2, y2, z2, x1, y2, z2 };
00576         EntityHandle verts[8], elem;
00577         for( int i = 0; i < 8; ++i )
00578             moab2.create_vertex( coords + 3 * i, verts[i] );
00579         moab2.create_element( MBHEX, verts, 8, elem );
00580         int tag_val = hash_handle( iter.handle() );
00581         moab2.tag_set_data( tag, &elem, 1, &tag_val );
00582     } while( MB_SUCCESS == iter.step() );
00583 
00584     moab2.write_mesh( file );
00585 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines