Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
KDTree.cpp
Go to the documentation of this file.
00001 /* Simple example of use of moab::AdaptiveKDTree class.
00002 
00003    Given a hexahedral mesh, find the hexahedron containing each
00004    input position.
00005  */
00006 
00007 #include "moab/Core.hpp"
00008 #include "moab/AdaptiveKDTree.hpp"
00009 #include "moab/Range.hpp"
00010 #include "moab/GeomUtil.hpp"
00011 
00012 #include <iostream>
00013 #include <string>
00014 
00015 const double EPSILON = 1e-6;  // tolerance to use in intersection checks
00016 
00017 // Help with error handling.  Given ErrorCode, print
00018 // corresponding string and any available message.
00019 void print_error( moab::Interface& mb, moab::ErrorCode err )
00020 {
00021     std::string message;
00022     std::string code;
00023     if( moab::MB_SUCCESS != mb.get_last_error( message ) ) message.clear();
00024     code = mb.get_error_string( err );
00025     std::cerr << "Error: " << code << std::endl;
00026     if( !message.empty() ) std::cerr << "  " << message << std::endl;
00027 }
00028 
00029 // Print diagnostic info for unexpected failures.
00030 #define CHKERR( err )                                                                           \
00031     do                                                                                          \
00032     {                                                                                           \
00033         if( moab::MB_SUCCESS != ( err ) )                                                       \
00034         {                                                                                       \
00035             print_error( mb, ( err ) );                                                         \
00036             std::cerr << "Unexpected failure at: " << __FILE__ << ":" << __LINE__ << std::endl; \
00037             return 2;                                                                           \
00038         }                                                                                       \
00039     } while( false )
00040 
00041 // Given an entity set and a point, find the hex contained in the
00042 // entity set which in turn contains the specified point.  Returns
00043 // 0 if point is not in any hexahedron.
00044 moab::EntityHandle hex_containing_point( moab::Interface& mb, moab::EntityHandle set, const double point[3] );
00045 
00046 // Print hex containing point.
00047 void print_hex( moab::Interface& mb, moab::EntityHandle hex );
00048 
00049 int main()
00050 {
00051     // Ask user for file to read
00052     std::string filename;
00053     std::cout << "Hex mesh file name: ";
00054     std::cin >> filename;
00055 
00056     // Read file into MOAB instance
00057     moab::ErrorCode rval;
00058     moab::Core moab;
00059     moab::Interface& mb = moab;
00060     rval                = mb.load_file( filename.c_str() );
00061     if( moab::MB_SUCCESS != rval )
00062     {
00063         print_error( mb, rval );
00064         std::cerr << filename << ": file load failed" << std::endl;
00065         return 1;
00066     }
00067 
00068     // Get all hex elemeents
00069     moab::Range elems;
00070     rval = mb.get_entities_by_type( 0, moab::MBHEX, elems );CHKERR( rval );
00071     if( elems.empty() )
00072     {
00073         std::cerr << filename << ": file containd no hexahedra" << std::endl;
00074         return 1;
00075     }
00076 
00077     // Build a kD-tree from hex elements
00078     moab::EntityHandle tree_root;
00079     moab::AdaptiveKDTree tool( &mb );
00080     rval = tool.build_tree( elems, tree_root );CHKERR( rval );
00081 
00082     // Loop forever (or until EOF), asking user for a point
00083     // to query and printing the hex element containing that
00084     // point.
00085     for( ;; )
00086     {
00087         double point[3];
00088         std::cout << "Point coordinates: ";
00089         if( !( std::cin >> point[0] >> point[1] >> point[2] ) ) break;
00090 
00091         moab::EntityHandle leaf;
00092         rval = tool.leaf_containing_point( tree_root, point, leaf );CHKERR( rval );
00093         moab::EntityHandle hex = hex_containing_point( mb, leaf, point );
00094         if( 0 == hex )
00095             std::cout << "Point is not contained in any hexahedron." << std::endl;
00096         else
00097             print_hex( mb, hex );
00098     }
00099 
00100     return 0;
00101 }
00102 
00103 moab::EntityHandle hex_containing_point( moab::Interface& mb, moab::EntityHandle set, const double point[3] )
00104 {
00105     moab::ErrorCode rval;
00106     moab::CartVect pt( point );      // input location
00107     moab::CartVect coords[8];        // coordinates of corners of hexahedron
00108     const moab::EntityHandle* conn;  // hex connectivity
00109     int conn_len;
00110 
00111     // Get hexes in leaf
00112     std::vector< moab::EntityHandle > hexes;
00113     rval = mb.get_entities_by_type( set, moab::MBHEX, hexes );CHKERR( rval );
00114 
00115     // Check which hex the point is in
00116     std::vector< moab::EntityHandle >::const_iterator i;
00117     for( i = hexes.begin(); i != hexes.end(); ++i )
00118     {
00119         rval = mb.get_connectivity( *i, conn, conn_len );CHKERR( rval );
00120         rval = mb.get_coords( conn, 8, &coords[0][0] );CHKERR( rval );
00121         if( moab::GeomUtil::point_in_trilinear_hex( coords, pt, EPSILON ) ) return *i;
00122     }
00123 
00124     // Return 0 if no hex contains point.
00125     return 0;
00126 }
00127 
00128 void print_hex( moab::Interface& mb, moab::EntityHandle hex )
00129 {
00130     // Get MOAB's internal ID for hex element
00131     int id = mb.id_from_handle( hex );
00132 
00133     // Get vertex handles for hex corners
00134     const moab::EntityHandle* conn;  // hex connectivity
00135     int conn_len;
00136     mb.get_connectivity( hex, conn, conn_len );
00137 
00138     // Get coordinates of vertices
00139     double coords[3 * 8];
00140     mb.get_coords( conn, 8, coords );
00141 
00142     // Print
00143     std::cout << " Point is in hex " << id << " with corners: " << std::endl;
00144     for( int i = 0; i < 8; ++i )
00145     {
00146         std::cout << " (" << coords[3 * i] << ", " << coords[3 * i + 1] << ", " << coords[3 * i + 2] << ")"
00147                   << std::endl;
00148     }
00149 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines