MOAB: Mesh Oriented datABase  (version 5.3.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             default:
00229                 usage();
00230         }
00231     }
00232 
00233     // this test relies on not cleaning up trees
00234     options << "CLEAN_UP=false;";
00235 
00236     if( make_grid != !output_file ) usage();
00237 
00238     // default to both
00239     if( !surf_elems && !vol_elems ) surf_elems = vol_elems = true;
00240 
00241     ErrorCode rval;
00242     Core moab_core;
00243     Interface* interface = &moab_core;
00244     FileOptions opts( options.str().c_str() );
00245 
00246     if( make_grid )
00247     {
00248         load_time   = 0;
00249         output_file = input_file;
00250         input_file  = 0;
00251         build_time  = clock();
00252         build_grid( interface, dims );
00253         build_time = clock() - build_time;
00254     }
00255     else
00256     {
00257         load_time = clock();
00258         rval      = interface->load_mesh( input_file );
00259         if( MB_SUCCESS != rval )
00260         {
00261             std::cerr << "Error reading file: " << input_file << std::endl;
00262             exit( 2 );
00263         }
00264         load_time = clock() - load_time;
00265 
00266         delete_existing_tree( interface );
00267 
00268         std::cout << "Building tree..." << std::endl;
00269         build_time = clock();
00270         Range elems;
00271         if( !surf_elems ) { interface->get_entities_by_dimension( 0, 3, elems ); }
00272         else
00273         {
00274             interface->get_entities_by_dimension( 0, 2, elems );
00275             if( vol_elems )
00276             {
00277                 Range tmp;
00278                 interface->get_entities_by_dimension( 0, 3, tmp );
00279                 elems.merge( tmp );
00280             }
00281         }
00282 
00283         build_tree( interface, elems, opts );
00284         build_time = clock() - build_time;
00285     }
00286 
00287     std::cout << "Calculating stats..." << std::endl;
00288     AdaptiveKDTree tool( interface );
00289     Range range;
00290     tool.find_all_trees( range );
00291     if( range.size() != 1 )
00292     {
00293         if( range.empty() )
00294             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00295         else
00296             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00297         exit( 5 );
00298     }
00299     tool.print();
00300 
00301     stat_time = clock() - build_time;
00302 
00303     if( tag_elems )
00304     {
00305         std::cout << "Tagging tree..." << std::endl;
00306         tag_elements( interface );
00307         tag_vertices( interface );
00308     }
00309     tag_time = clock() - stat_time;
00310 
00311     std::cout << "Writing file... ";
00312     std::cout.flush();
00313     rval = interface->write_mesh( output_file );
00314     if( MB_SUCCESS != rval )
00315     {
00316         std::cerr << "Error writing file: " << output_file << std::endl;
00317         exit( 3 );
00318     }
00319     write_time = clock() - tag_time;
00320     std::cout << "Wrote " << output_file << std::endl;
00321 
00322     if( tree_file )
00323     {
00324         std::cout << "Writing tree block rep...";
00325         std::cout.flush();
00326         write_tree_blocks( interface, tree_file );
00327         std::cout << "Wrote " << tree_file << std::endl;
00328     }
00329     block_time = clock() - write_time;
00330 
00331     std::cout << "Times:  "
00332               << "    Load"
00333               << "   Build"
00334               << "   Stats"
00335               << "   Write";
00336     if( tag_elems ) std::cout << "Tag Sets";
00337     if( tree_file ) std::cout << "Block   ";
00338     std::cout << std::endl;
00339 
00340     std::cout << "        " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 )
00341               << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 )
00342               << clock_to_string( write_time );
00343     if( tag_elems ) std::cout << std::setw( 8 ) << clock_to_string( tag_time );
00344     if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time );
00345     std::cout << std::endl;
00346 
00347     return 0;
00348 }
00349 
00350 void delete_existing_tree( Interface* interface )
00351 {
00352     Range trees;
00353     AdaptiveKDTree tool( interface );
00354 
00355     tool.find_all_trees( trees );
00356     if( !trees.empty() ) std::cout << "Destroying existing tree(s) contained in file" << std::endl;
00357     for( Range::iterator i = trees.begin(); i != trees.end(); ++i )
00358         tool.reset_tree();
00359 
00360     trees.clear();
00361     tool.find_all_trees( trees );
00362     if( !trees.empty() )
00363     {
00364         std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
00365         exit( 5 );
00366     }
00367 }
00368 
00369 void build_tree( Interface* interface, const Range& elems, FileOptions& opts )
00370 {
00371     EntityHandle root = 0;
00372 
00373     if( elems.empty() )
00374     {
00375         std::cerr << "No elements from which to build tree." << std::endl;
00376         exit( 4 );
00377     }
00378 
00379     AdaptiveKDTree tool( interface, elems, &root, &opts );
00380 }
00381 
00382 void build_grid( Interface* interface, const double dims[3] )
00383 {
00384     ErrorCode rval;
00385     EntityHandle root = 0;
00386     AdaptiveKDTree tool( interface );
00387     AdaptiveKDTreeIter iter;
00388     AdaptiveKDTree::Plane plane;
00389 
00390     double min[3] = { -0.5 * dims[0], -0.5 * dims[1], -0.5 * dims[2] };
00391     double max[3] = { 0.5 * dims[0], 0.5 * dims[1], 0.5 * dims[2] };
00392     rval          = tool.create_root( min, max, root );
00393     if( MB_SUCCESS != rval || !root )
00394     {
00395         std::cerr << "Failed to create tree root." << std::endl;
00396         exit( 4 );
00397     }
00398 
00399     rval = tool.get_tree_iterator( root, iter );
00400     if( MB_SUCCESS != rval ) { std::cerr << "Failed to get tree iterator." << std::endl; }
00401 
00402     do
00403     {
00404         while( iter.depth() < tool.get_max_depth() )
00405         {
00406             plane.norm  = iter.depth() % 3;
00407             plane.coord = 0.5 * ( iter.box_min()[plane.norm] + iter.box_max()[plane.norm] );
00408             rval        = tool.split_leaf( iter, plane );
00409             if( MB_SUCCESS != rval )
00410             {
00411                 std::cerr << "Failed to split tree leaf at depth " << iter.depth() << std::endl;
00412                 exit( 4 );
00413             }
00414         }
00415     } while( ( rval = iter.step() ) == MB_SUCCESS );
00416 
00417     if( rval != MB_ENTITY_NOT_FOUND )
00418     {
00419         std::cerr << "Error stepping KDTree iterator." << std::endl;
00420         exit( 4 );
00421     }
00422 }
00423 
00424 std::string clock_to_string( clock_t t )
00425 {
00426     char unit[5] = "s";
00427     char buffer[256];
00428     double dt = t / (double)CLOCKS_PER_SEC;
00429     // if (dt > 300) {
00430     //  dt /= 60;
00431     //  strcpy( unit, "min" );
00432     //}
00433     // if (dt > 300) {
00434     //  dt /= 60;
00435     //  strcpy( unit, "hr" );
00436     //}
00437     // if (dt > 100) {
00438     //  dt /= 24;
00439     //  strcpy( unit, "days" );
00440     //}
00441     sprintf( buffer, "%0.2f%s", dt, unit );
00442     return buffer;
00443 }
00444 
00445 static int hash_handle( EntityHandle handle )
00446 {
00447     EntityID h = ID_FROM_HANDLE( handle );
00448     return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1;
00449 }
00450 
00451 void tag_elements( Interface* moab )
00452 {
00453     EntityHandle root;
00454     Range range;
00455     AdaptiveKDTree tool( moab );
00456 
00457     tool.find_all_trees( range );
00458     if( range.size() != 1 )
00459     {
00460         if( range.empty() )
00461             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00462         else
00463             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00464         exit( 5 );
00465     }
00466 
00467     root = *range.begin();
00468     range.clear();
00469 
00470     Tag tag;
00471     int zero = 0;
00472     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00473 
00474     AdaptiveKDTreeIter iter;
00475     tool.get_tree_iterator( root, iter );
00476 
00477     std::vector< int > tag_vals;
00478     do
00479     {
00480         range.clear();
00481         moab->get_entities_by_handle( iter.handle(), range );
00482         tag_vals.clear();
00483         tag_vals.resize( range.size(), hash_handle( iter.handle() ) );
00484         moab->tag_set_data( tag, range, &tag_vals[0] );
00485     } while( MB_SUCCESS == iter.step() );
00486 }
00487 
00488 void tag_vertices( Interface* moab )
00489 {
00490     EntityHandle root;
00491     Range range;
00492     AdaptiveKDTree tool( moab );
00493 
00494     tool.find_all_trees( range );
00495     if( range.size() != 1 )
00496     {
00497         if( range.empty() )
00498             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00499         else
00500             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00501         exit( 5 );
00502     }
00503 
00504     root = *range.begin();
00505     range.clear();
00506 
00507     Tag tag;
00508     int zero = 0;
00509     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00510 
00511     AdaptiveKDTreeIter iter;
00512     tool.get_tree_iterator( root, iter );
00513 
00514     do
00515     {
00516         range.clear();
00517         moab->get_entities_by_handle( iter.handle(), range );
00518 
00519         int tag_val = hash_handle( iter.handle() );
00520         Range verts;
00521         moab->get_adjacencies( range, 0, false, verts, Interface::UNION );
00522         for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
00523         {
00524             CartVect coords;
00525             moab->get_coords( &*i, 1, coords.array() );
00526             if( GeomUtil::box_point_overlap( CartVect( iter.box_min() ), CartVect( iter.box_max() ), coords, 1e-7 ) )
00527                 moab->tag_set_data( tag, &*i, 1, &tag_val );
00528         }
00529     } while( MB_SUCCESS == iter.step() );
00530 }
00531 
00532 void write_tree_blocks( Interface* interface, const char* file )
00533 {
00534     EntityHandle root;
00535     Range range;
00536     AdaptiveKDTree tool( interface );
00537 
00538     tool.find_all_trees( range );
00539     if( range.size() != 1 )
00540     {
00541         if( range.empty() )
00542             std::cerr << "Internal error: Failed to retreive tree." << std::endl;
00543         else
00544             std::cerr << "Internal error: Multiple tree roots." << std::endl;
00545         exit( 5 );
00546     }
00547 
00548     root = *range.begin();
00549     range.clear();
00550 
00551     Core moab2;
00552     Tag tag;
00553     int zero = 0;
00554     moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00555 
00556     AdaptiveKDTreeIter iter;
00557     tool.get_tree_iterator( root, iter );
00558 
00559     do
00560     {
00561         double x1         = iter.box_min()[0];
00562         double y1         = iter.box_min()[1];
00563         double z1         = iter.box_min()[2];
00564         double x2         = iter.box_max()[0];
00565         double y2         = iter.box_max()[1];
00566         double z2         = iter.box_max()[2];
00567         double coords[24] = { x1, y1, z1, x2, y1, z1, x2, y2, z1, x1, y2, z1,
00568                               x1, y1, z2, x2, y1, z2, x2, y2, z2, x1, y2, z2 };
00569         EntityHandle verts[8], elem;
00570         for( int i = 0; i < 8; ++i )
00571             moab2.create_vertex( coords + 3 * i, verts[i] );
00572         moab2.create_element( MBHEX, verts, 8, elem );
00573         int tag_val = hash_handle( iter.handle() );
00574         moab2.tag_set_data( tag, &elem, 1, &tag_val );
00575     } while( MB_SUCCESS == iter.step() );
00576 
00577     moab2.write_mesh( file );
00578 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines