kd_tree_tool.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/CartVect.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/GeomUtil.hpp"
#include "Internals.hpp"
#include "moab/Range.hpp"
#include "moab/FileOptions.hpp"
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <limits>
#include <sstream>
#include <cstdlib>
#include <ctime>
#if !defined( _MSC_VER ) && !defined( __MINGW32__ )
#include <unistd.h>
#include <sys/stat.h>
#include <fcntl.h>
#endif
 
using namespace moab;
 
const int MAX_TAG_VALUE = 32; // maximum value to use when tagging entities with tree cell number
const char* TAG_NAME = "TREE_CELL";
 
std::string clock_to_string( clock_t t );
std::string mem_to_string( unsigned long mem );
 
void build_tree( Interface* interface, const Range& elems, FileOptions& opts );
void build_grid( Interface* iface, const double dims[3] );
void delete_existing_tree( Interface* interface );
void print_stats( Interface* interface );
void tag_elements( Interface* interface );
void tag_vertices( Interface* interface );
void write_tree_blocks( Interface* interface, const char* file );
 
static void usage( bool err = true )
{
 std::ostream& s = err ? std::cerr : std::cout;
 s << "kd_tree_tool [-d <int>] [-2|-3] [-n <int>] [-u|-p|-m|-v] [-N <int>] [-s|-S] <input file> "
 "<output file>"
 << std::endl
 << "kd_tree_tool [-d <int>] -G <dims> [-s|-S] <output file>" << std::endl
 << "kd_tree_tool [-h]" << std::endl;
 if( !err )
 {
 s << "Tool to build adaptive kd-Tree" << std::endl;
 s << " -d <int> Specify maximum depth for tree. Default: 30" << std::endl
 << " -n <int> Specify maximum entities per leaf. Default: 6" << std::endl
 << " -2 Build tree from surface elements. Default: yes" << std::endl
 << " -3 Build tree from volume elements. Default: yes" << std::endl
 << " -u Use 'SUBDIVISION' scheme for tree construction" << std::endl
 << " -p Use 'SUBDIVISION_SNAP' tree construction algorithm." << std::endl
 << " -m Use 'VERTEX_MEDIAN' tree construction algorithm." << std::endl
 << " -v Use 'VERTEX_SAMPLE' tree construction algorithm." << std::endl
 << " -N <int> Specify candidate split planes per axis. Default: 3" << std::endl
 << " -t Tag elements will tree cell number." << std::endl
 << " -T Write tree boxes to file." << std::endl
 << " -G <dims> Generate grid - no input elements. Dims must be " << std::endl
 << " HxWxD or H." << std::endl
 << " -s Use range-based sets for tree nodes" << std::endl
 << " -S Use vector-based sets for tree nodes" << std::endl
 << std::endl;
 }
 exit( err );
}
 
/*
#if !defined(_MSC_VER) && !defined(__MINGW32__)
static void memory_use( unsigned long& vsize, unsigned long& rss )
{
 char buffer[512];
 int filp = open( "/proc/self/stat", O_RDONLY );
 ssize_t r = read( filp, buffer, sizeof(buffer)-1 );
 close( filp );
 if (r < 0) r = 0;
 vsize = rss = 0;
 buffer[r] = '\0';
 sscanf( buffer, "%*d %*s %*c " // pid command state
 "%*d %*d " // ppid pgrp
 "%*d %*d %*d " // session tty_nr tpgid
 "%*u " // flags
 "%*u %*u %*u %*u " // minflt cminflt majflt cmajflt
 "%*u %*u %*d %*d " // utime stime cutime cstime
 "%*d %*d %*d " // priority nice (unused)
 "%*d %*u " // itrealval starttime
 "%lu %lu", &vsize, &rss );
 rss *= getpagesize();
}
#else
static void memory_use( unsigned long& vsize, unsigned long& rss )
 { vsize = rss = 0; }
#endif
*/
 
static int parseint( int& i, int argc, char* argv[] )
{
 char* end;
 ++i;
 if( i == argc )
 {
 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 int result = strtol( argv[i], &end, 0 );
 if( result < 0 || *end )
 {
 std::cerr << "Expected positive integer following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 return result;
}
 
static void parsedims( int& i, int argc, char* argv[], double dims[3] )
{
 char* end;
 ++i;
 if( i == argc )
 {
 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 dims[0] = strtod( argv[i], &end );
 if( *end == '\0' )
 {
 dims[2] = dims[1] = dims[0];
 return;
 }
 else if( *end != 'x' && *end != 'X' )
 goto error;
 
 ++end;
 dims[1] = strtod( end, &end );
 if( *end == '\0' )
 {
 dims[2] = dims[1];
 return;
 }
 else if( *end != 'x' && *end != 'X' )
 goto error;
 
 ++end;
 dims[2] = strtod( end, &end );
 if( *end != '\0' ) goto error;
 
 return;
error:
 std::cerr << "Invaild dimension specification." << std::endl;
 usage();
}
 
int main( int argc, char* argv[] )
{
 bool surf_elems = false, vol_elems = false;
 const char* input_file = 0;
 const char* output_file = 0;
 const char* tree_file = 0;
 std::ostringstream options;
 bool tag_elems = false;
 clock_t load_time, build_time, stat_time, tag_time, write_time, block_time;
 bool make_grid = false;
 double dims[3];
 
 for( int i = 1; i < argc; ++i )
 {
 if( argv[i][0] != '-' )
 {
 if( !input_file )
 input_file = argv[i];
 else if( !output_file )
 output_file = argv[i];
 else
 usage();
 continue;
 }
 
 if( !argv[i][1] || argv[i][2] ) usage();
 
 switch( argv[i][1] )
 {
 case '2':
 surf_elems = true;
 break;
 case '3':
 vol_elems = true;
 break;
 case 'S':
 options << "MESHSET_FLAGS=" << MESHSET_ORDERED << ";";
 break;
 case 's':
 break;
 case 'd':
 options << "MAX_DEPTH=" << parseint( i, argc, argv ) << ";";
 break;
 case 'n':
 options << "MAX_PER_LEAF=" << parseint( i, argc, argv ) << ";";
 break;
 case 'u':
 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::SUBDIVISION << ";";
 break;
 case 'p':
 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::SUBDIVISION_SNAP << ";";
 break;
 case 'm':
 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::VERTEX_MEDIAN << ";";
 break;
 case 'v':
 options << "CANDIDATE_PLANE_SET=" << AdaptiveKDTree::VERTEX_SAMPLE << ";";
 break;
 case 'N':
 options << "SPLITS_PER_DIR=" << parseint( i, argc, argv ) << ";";
 break;
 case 't':
 tag_elems = true;
 break;
 case 'T':
 tree_file = argv[++i];
 break;
 case 'G':
 make_grid = true;
 parsedims( i, argc, argv, dims );
 break;
 case 'h':
 usage( false );
 break;
 default:
 usage();
 }
 }
 
 // this test relies on not cleaning up trees
 options << "CLEAN_UP=false;";
 
 if( make_grid != !output_file ) usage();
 
 // default to both
 if( !surf_elems && !vol_elems ) surf_elems = vol_elems = true;
 
 ErrorCode rval;
 Core moab_core;
 Interface* interface = &moab_core;
 FileOptions opts( options.str().c_str() );
 
 if( make_grid )
 {
 load_time = 0;
 output_file = input_file;
 input_file = 0;
 build_time = clock();
 build_grid( interface, dims );
 build_time = clock() - build_time;
 }
 else
 {
 load_time = clock();
 rval = interface->load_mesh( input_file );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Error reading file: " << input_file << std::endl;
 exit( 2 );
 }
 load_time = clock() - load_time;
 
 delete_existing_tree( interface );
 
 std::cout << "Building tree..." << std::endl;
 build_time = clock();
 Range elems;
 if( !surf_elems )
 {
 interface->get_entities_by_dimension( 0, 3, elems );
 }
 else
 {
 interface->get_entities_by_dimension( 0, 2, elems );
 if( vol_elems )
 {
 Range tmp;
 interface->get_entities_by_dimension( 0, 3, tmp );
 elems.merge( tmp );
 }
 }
 
 build_tree( interface, elems, opts );
 build_time = clock() - build_time;
 }
 
 std::cout << "Calculating stats..." << std::endl;
 AdaptiveKDTree tool( interface );
 Range range;
 tool.find_all_trees( range );
 if( range.size() != 1 )
 {
 if( range.empty() )
 std::cerr << "Internal error: Failed to retreive tree." << std::endl;
 else
 std::cerr << "Internal error: Multiple tree roots." << std::endl;
 exit( 5 );
 }
 tool.print();
 
 stat_time = clock() - build_time;
 
 if( tag_elems )
 {
 std::cout << "Tagging tree..." << std::endl;
 tag_elements( interface );
 tag_vertices( interface );
 }
 tag_time = clock() - stat_time;
 
 std::cout << "Writing file... ";
 std::cout.flush();
 rval = interface->write_mesh( output_file );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Error writing file: " << output_file << std::endl;
 exit( 3 );
 }
 write_time = clock() - tag_time;
 std::cout << "Wrote " << output_file << std::endl;
 
 if( tree_file )
 {
 std::cout << "Writing tree block rep...";
 std::cout.flush();
 write_tree_blocks( interface, tree_file );
 std::cout << "Wrote " << tree_file << std::endl;
 }
 block_time = clock() - write_time;
 
 std::cout << "Times: "
 << " Load"
 << " Build"
 << " Stats"
 << " Write";
 if( tag_elems ) std::cout << "Tag Sets";
 if( tree_file ) std::cout << "Block ";
 std::cout << std::endl;
 
 std::cout << " " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 )
 << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 )
 << clock_to_string( write_time );
 if( tag_elems ) std::cout << std::setw( 8 ) << clock_to_string( tag_time );
 if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time );
 std::cout << std::endl;
 
 return 0;
}
 
void delete_existing_tree( Interface* interface )
{
 Range trees;
 AdaptiveKDTree tool( interface );
 
 tool.find_all_trees( trees );
 if( !trees.empty() ) std::cout << "Destroying existing tree(s) contained in file" << std::endl;
 for( Range::iterator i = trees.begin(); i != trees.end(); ++i )
 tool.reset_tree();
 
 trees.clear();
 tool.find_all_trees( trees );
 if( !trees.empty() )
 {
 std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
 exit( 5 );
 }
}
 
void build_tree( Interface* interface, const Range& elems, FileOptions& opts )
{
 EntityHandle root = 0;
 
 if( elems.empty() )
 {
 std::cerr << "No elements from which to build tree." << std::endl;
 exit( 4 );
 }
 
 AdaptiveKDTree tool( interface, elems, &root, &opts );
}
 
void build_grid( Interface* interface, const double dims[3] )
{
 ErrorCode rval;
 EntityHandle root = 0;
 AdaptiveKDTree tool( interface );
 AdaptiveKDTreeIter iter;
 AdaptiveKDTree::Plane plane;
 
 double min[3] = { -0.5 * dims[0], -0.5 * dims[1], -0.5 * dims[2] };
 double max[3] = { 0.5 * dims[0], 0.5 * dims[1], 0.5 * dims[2] };
 rval = tool.create_root( min, max, root );
 if( MB_SUCCESS != rval || !root )
 {
 std::cerr << "Failed to create tree root." << std::endl;
 exit( 4 );
 }
 
 rval = tool.get_tree_iterator( root, iter );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Failed to get tree iterator." << std::endl;
 }
 
 do
 {
 while( iter.depth() < tool.get_max_depth() )
 {
 plane.norm = iter.depth() % 3;
 plane.coord = 0.5 * ( iter.box_min()[plane.norm] + iter.box_max()[plane.norm] );
 rval = tool.split_leaf( iter, plane );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Failed to split tree leaf at depth " << iter.depth() << std::endl;
 exit( 4 );
 }
 }
 } while( ( rval = iter.step() ) == MB_SUCCESS );
 
 if( rval != MB_ENTITY_NOT_FOUND )
 {
 std::cerr << "Error stepping KDTree iterator." << std::endl;
 exit( 4 );
 }
}
 
std::string clock_to_string( clock_t t )
{
 char unit[5] = "s";
 char buffer[256];
 double dt = t / (double)CLOCKS_PER_SEC;
 // if (dt > 300) {
 // dt /= 60;
 // strcpy( unit, "min" );
 //}
 // if (dt > 300) {
 // dt /= 60;
 // strcpy( unit, "hr" );
 //}
 // if (dt > 100) {
 // dt /= 24;
 // strcpy( unit, "days" );
 //}
 sprintf( buffer, "%0.2f%s", dt, unit );
 return buffer;
}
 
static int hash_handle( EntityHandle handle )
{
 EntityID h = ID_FROM_HANDLE( handle );
 return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1;
}
 
void tag_elements( Interface* moab )
{
 EntityHandle root;
 Range range;
 AdaptiveKDTree tool( moab );
 
 tool.find_all_trees( range );
 if( range.size() != 1 )
 {
 if( range.empty() )
 std::cerr << "Internal error: Failed to retreive tree." << std::endl;
 else
 std::cerr << "Internal error: Multiple tree roots." << std::endl;
 exit( 5 );
 }
 
 root = *range.begin();
 range.clear();
 
 Tag tag;
 int zero = 0;
 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 
 AdaptiveKDTreeIter iter;
 tool.get_tree_iterator( root, iter );
 
 std::vector< int > tag_vals;
 do
 {
 range.clear();
 moab->get_entities_by_handle( iter.handle(), range );
 tag_vals.clear();
 tag_vals.resize( range.size(), hash_handle( iter.handle() ) );
 moab->tag_set_data( tag, range, &tag_vals[0] );
 } while( MB_SUCCESS == iter.step() );
}
 
void tag_vertices( Interface* moab )
{
 EntityHandle root;
 Range range;
 AdaptiveKDTree tool( moab );
 
 tool.find_all_trees( range );
 if( range.size() != 1 )
 {
 if( range.empty() )
 std::cerr << "Internal error: Failed to retreive tree." << std::endl;
 else
 std::cerr << "Internal error: Multiple tree roots." << std::endl;
 exit( 5 );
 }
 
 root = *range.begin();
 range.clear();
 
 Tag tag;
 int zero = 0;
 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 
 AdaptiveKDTreeIter iter;
 tool.get_tree_iterator( root, iter );
 
 do
 {
 range.clear();
 moab->get_entities_by_handle( iter.handle(), range );
 
 int tag_val = hash_handle( iter.handle() );
 Range verts;
 moab->get_adjacencies( range, 0, false, verts, Interface::UNION );
 for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
 {
 CartVect coords;
 moab->get_coords( &*i, 1, coords.array() );
 if( GeomUtil::box_point_overlap( CartVect( iter.box_min() ), CartVect( iter.box_max() ), coords, 1e-7 ) )
 moab->tag_set_data( tag, &*i, 1, &tag_val );
 }
 } while( MB_SUCCESS == iter.step() );
}
 
void write_tree_blocks( Interface* interface, const char* file )
{
 EntityHandle root;
 Range range;
 AdaptiveKDTree tool( interface );
 
 tool.find_all_trees( range );
 if( range.size() != 1 )
 {
 if( range.empty() )
 std::cerr << "Internal error: Failed to retreive tree." << std::endl;
 else
 std::cerr << "Internal error: Multiple tree roots." << std::endl;
 exit( 5 );
 }
 
 root = *range.begin();
 range.clear();
 
 Core moab2;
 Tag tag;
 int zero = 0;
 moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 
 AdaptiveKDTreeIter iter;
 tool.get_tree_iterator( root, iter );
 
 do
 {
 double x1 = iter.box_min()[0];
 double y1 = iter.box_min()[1];
 double z1 = iter.box_min()[2];
 double x2 = iter.box_max()[0];
 double y2 = iter.box_max()[1];
 double z2 = iter.box_max()[2];
 double coords[24] = { x1, y1, z1, x2, y1, z1, x2, y2, z1, x1, y2, z1,
 x1, y1, z2, x2, y1, z2, x2, y2, z2, x1, y2, z2 };
 EntityHandle verts[8], elem;
 for( int i = 0; i < 8; ++i )
 moab2.create_vertex( coords + 3 * i, verts[i] );
 moab2.create_element( MBHEX, verts, 8, elem );
 int tag_val = hash_handle( iter.handle() );
 moab2.tag_set_data( tag, &elem, 1, &tag_val );
 } while( MB_SUCCESS == iter.step() );
 
 moab2.write_mesh( file );
}