 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /* Simple example of use of moab::AdaptiveKDTree class.
2
3 Given a hexahedral mesh, find the hexahedron containing each
4 input position.
5 */
6
7 #include "moab/Core.hpp"
8 #include "moab/AdaptiveKDTree.hpp"
9 #include "moab/Range.hpp"
10 #include "moab/GeomUtil.hpp"
11
12 #include <iostream>
13 #include <string>
14
15 const double EPSILON = 1e-6; // tolerance to use in intersection checks
16
17 // Help with error handling. Given ErrorCode, print
18 // corresponding string and any available message.
19 void print_error( moab::Interface& mb, moab::ErrorCode err )
20 {
21 std::string message;
22 std::string code;
23 if( moab::MB_SUCCESS != mb.get_last_error( message ) ) message.clear();
24 code = mb.get_error_string( err );
25 std::cerr << "Error: " << code << std::endl;
26 if( !message.empty() ) std::cerr << " " << message << std::endl;
27 }
28
29 // Print diagnostic info for unexpected failures.
30 #define CHKERR( err ) \
31 do \
32 { \
33 if( moab::MB_SUCCESS != ( err ) ) \
34 { \
35 print_error( mb, ( err ) ); \
36 std::cerr << "Unexpected failure at: " << __FILE__ << ":" << __LINE__ << std::endl; \
37 return 2; \
38 } \
39 } while( false )
40
41 // Given an entity set and a point, find the hex contained in the
42 // entity set which in turn contains the specified point. Returns
43 // 0 if point is not in any hexahedron.
44 moab::EntityHandle hex_containing_point( moab::Interface& mb, moab::EntityHandle set, const double point[3] );
45
46 // Print hex containing point.
47 void print_hex( moab::Interface& mb, moab::EntityHandle hex );
48
49 int main()
50 {
51 // Ask user for file to read
52 std::string filename;
53 std::cout << "Hex mesh file name: ";
54 std::cin >> filename;
55
56 // Read file into MOAB instance
57 moab::ErrorCode rval;
58 moab::Core moab;
59 moab::Interface& mb = moab;
60 rval = mb.load_file( filename.c_str() );
61 if( moab::MB_SUCCESS != rval )
62 {
63 print_error( mb, rval );
64 std::cerr << filename << ": file load failed" << std::endl;
65 return 1;
66 }
67
68 // Get all hex elemeents
69 moab::Range elems;
70 rval = mb.get_entities_by_type( 0, moab::MBHEX, elems );CHKERR( rval );
71 if( elems.empty() )
72 {
73 std::cerr << filename << ": file containd no hexahedra" << std::endl;
74 return 1;
75 }
76
77 // Build a kD-tree from hex elements
78 moab::EntityHandle tree_root;
79 moab::AdaptiveKDTree tool( &mb );
80 rval = tool.build_tree( elems, tree_root );CHKERR( rval );
81
82 // Loop forever (or until EOF), asking user for a point
83 // to query and printing the hex element containing that
84 // point.
85 for( ;; )
86 {
87 double point[3];
88 std::cout << "Point coordinates: ";
89 if( !( std::cin >> point[0] >> point[1] >> point[2] ) ) break;
90
91 moab::EntityHandle leaf;
92 rval = tool.leaf_containing_point( tree_root, point, leaf );CHKERR( rval );
93 moab::EntityHandle hex = hex_containing_point( mb, leaf, point );
94 if( 0 == hex )
95 std::cout << "Point is not contained in any hexahedron." << std::endl;
96 else
97 print_hex( mb, hex );
98 }
99
100 return 0;
101 }
102
103 moab::EntityHandle hex_containing_point( moab::Interface& mb, moab::EntityHandle set, const double point[3] )
104 {
105 moab::ErrorCode rval;
106 moab::CartVect pt( point ); // input location
107 moab::CartVect coords[8]; // coordinates of corners of hexahedron
108 const moab::EntityHandle* conn; // hex connectivity
109 int conn_len;
110
111 // Get hexes in leaf
112 std::vector< moab::EntityHandle > hexes;
113 rval = mb.get_entities_by_type( set, moab::MBHEX, hexes );CHKERR( rval );
114
115 // Check which hex the point is in
116 std::vector< moab::EntityHandle >::const_iterator i;
117 for( i = hexes.begin(); i != hexes.end(); ++i )
118 {
119 rval = mb.get_connectivity( *i, conn, conn_len );CHKERR( rval );
120 rval = mb.get_coords( conn, 8, &coords[0][0] );CHKERR( rval );
121 if( moab::GeomUtil::point_in_trilinear_hex( coords, pt, EPSILON ) ) return *i;
122 }
123
124 // Return 0 if no hex contains point.
125 return 0;
126 }
127
128 void print_hex( moab::Interface& mb, moab::EntityHandle hex )
129 {
130 // Get MOAB's internal ID for hex element
131 int id = mb.id_from_handle( hex );
132
133 // Get vertex handles for hex corners
134 const moab::EntityHandle* conn; // hex connectivity
135 int conn_len;
136 mb.get_connectivity( hex, conn, conn_len );
137
138 // Get coordinates of vertices
139 double coords[3 * 8];
140 mb.get_coords( conn, 8, coords );
141
142 // Print
143 std::cout << " Point is in hex " << id << " with corners: " << std::endl;
144 for( int i = 0; i < 8; ++i )
145 {
146 std::cout << " (" << coords[3 * i] << ", " << coords[3 * i + 1] << ", " << coords[3 * i + 2] << ")"
147 << std::endl;
148 }
149 }
1.9.1.