MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "ReadTetGen.hpp" 00002 #include "moab/Interface.hpp" 00003 #include "moab/Range.hpp" 00004 #include "moab/ReadUtilIface.hpp" 00005 #include "moab/FileOptions.hpp" 00006 #include "MBTagConventions.hpp" 00007 #include <iostream> 00008 #include <fstream> 00009 #include <sstream> 00010 #include <ctype.h> 00011 #include <map> 00012 00013 namespace moab 00014 { 00015 00016 ReaderIface* ReadTetGen::factory( Interface* moab ) 00017 { 00018 return new ReadTetGen( moab ); 00019 } 00020 00021 ReadTetGen::ReadTetGen( Interface* moab ) : mbIface( moab ), readTool( 0 ) 00022 { 00023 moab->query_interface( readTool ); 00024 } 00025 00026 ReadTetGen::~ReadTetGen() 00027 { 00028 if( mbIface && readTool ) mbIface->release_interface( readTool ); 00029 } 00030 00031 ErrorCode ReadTetGen::open_file( const std::string& filename, const std::string& basename, const std::string& suffix, 00032 const char* exp_suffix, const char* opt_name, const FileOptions& opts, 00033 std::ifstream& file_stream, bool file_required ) 00034 { 00035 std::string real_file_name; 00036 ErrorCode rval = opts.get_option( opt_name, real_file_name ); 00037 if( MB_ENTITY_NOT_FOUND == rval || real_file_name.empty() ) 00038 { 00039 if( MB_SUCCESS == rval ) file_required = true; 00040 if( suffix == exp_suffix ) { real_file_name = filename; } 00041 else 00042 { 00043 real_file_name = basename; 00044 real_file_name += "."; 00045 real_file_name += exp_suffix; 00046 } 00047 } 00048 00049 if( !real_file_name.empty() ) file_stream.open( real_file_name.c_str(), std::ios::in ); 00050 if( file_required && !file_stream.is_open() ) 00051 { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, real_file_name << ": cannot read file" ); } 00052 00053 return MB_SUCCESS; 00054 } 00055 00056 ErrorCode ReadTetGen::read_tag_values( const char* /* file_name */, const char* /* tag_name */, 00057 const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */, 00058 const SubsetList* /* subset_list */ ) 00059 { 00060 return MB_NOT_IMPLEMENTED; 00061 } 00062 00063 ErrorCode ReadTetGen::load_file( const char* file_name_c, const EntityHandle* /* file_set */, const FileOptions& opts, 00064 const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag ) 00065 { 00066 std::ifstream node_file, ele_file, face_file, edge_file; 00067 ErrorCode rval; 00068 00069 if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for TetGen" ); } 00070 00071 std::string suffix, base, filename( file_name_c ); 00072 size_t dot_idx = filename.find_last_of( '.' ); 00073 if( dot_idx == std::string::npos ) { base = filename; } 00074 else 00075 { 00076 suffix = filename.substr( dot_idx + 1 ); 00077 for( size_t i = 0; i < suffix.length(); ++i ) 00078 suffix[i] = (char)tolower( suffix[i] ); 00079 if( suffix == "node" || suffix == "ele" || suffix == "face" || suffix == "edge" ) 00080 { base = filename.substr( 0, dot_idx ); } 00081 else 00082 { 00083 base = filename; 00084 suffix.clear(); 00085 } 00086 } 00087 00088 rval = open_file( filename, base, suffix, "node", "NODE_FILE", opts, node_file, true ); 00089 if( MB_SUCCESS != rval ) return rval; 00090 rval = open_file( filename, base, suffix, "ele", "ELE_FILE", opts, ele_file ); 00091 if( MB_SUCCESS != rval ) return rval; 00092 rval = open_file( filename, base, suffix, "face", "FACE_FILE", opts, face_file ); 00093 if( MB_SUCCESS != rval ) return rval; 00094 rval = open_file( filename, base, suffix, "edge", "EDGE_FILE", opts, edge_file ); 00095 if( MB_SUCCESS != rval ) return rval; 00096 00097 std::vector< Tag > attr_tags[4]; 00098 std::vector< int > attr_idx[4]; 00099 const char* option_names[4] = { "NODE_ATTR_LIST", "EDGE_ATTR_LIST", "TRI_ATTR_LIST", "TET_ATTR_LIST" }; 00100 const char* group_names[4] = { 0, "CURVE_ID", "SURFACE_ID", "VOLUME_ID" }; 00101 for( int i = 0; i < 4; ++i ) 00102 { 00103 std::string opt_str; 00104 rval = opts.get_str_option( option_names[i], opt_str ); 00105 if( MB_SUCCESS != rval ) continue; 00106 rval = parse_attr_list( opt_str, attr_tags[i], attr_idx[i], group_names[i] ); 00107 if( MB_SUCCESS != rval ) { MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, option_names[i] << ": invalid option value" ); } 00108 } 00109 00110 Range tets, tris, edges; 00111 std::vector< EntityHandle > nodes; 00112 rval = read_node_file( node_file, &attr_tags[0][0], &attr_idx[0][0], attr_tags[0].size(), nodes ); 00113 if( MB_SUCCESS == rval && ele_file.is_open() ) rval = read_elem_file( MBTET, ele_file, nodes, tets ); 00114 if( MB_SUCCESS == rval && face_file.is_open() ) rval = read_elem_file( MBTRI, face_file, nodes, tris ); 00115 if( MB_SUCCESS == rval && edge_file.is_open() ) rval = read_elem_file( MBEDGE, edge_file, nodes, edges ); 00116 00117 if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, &nodes[0], nodes.size() ); 00118 if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, edges ); 00119 if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tris ); 00120 if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tets ); 00121 00122 return rval; 00123 } 00124 00125 ErrorCode ReadTetGen::parse_attr_list( const std::string& option_str, std::vector< Tag >& tag_list, 00126 std::vector< int >& index_list, const char* group_designator ) 00127 { 00128 std::vector< std::string > name_list; 00129 size_t prev_pos = 0; 00130 while( prev_pos != std::string::npos ) 00131 { 00132 size_t pos = option_str.find_first_of( ',', prev_pos ); 00133 name_list.push_back( option_str.substr( prev_pos, pos ) ); 00134 prev_pos = pos + 1; 00135 } 00136 00137 index_list.resize( name_list.size() ); 00138 std::map< std::string, int > name_count; 00139 for( size_t i = 0; i < name_list.size(); ++i ) 00140 index_list[i] = name_count[name_list[i]]++; 00141 00142 for( size_t i = 0; i < name_list.size(); ++i ) 00143 { 00144 if( group_designator && name_list[i] == group_designator ) 00145 { 00146 tag_list[i] = 0; 00147 index_list[i] = -1; 00148 } 00149 else if( name_list.empty() ) 00150 { 00151 tag_list[i] = 0; 00152 index_list[i] = 0; 00153 } 00154 else 00155 { 00156 ErrorCode rval = mbIface->tag_get_handle( name_list[i].c_str(), name_count[name_list[i]], MB_TYPE_DOUBLE, 00157 tag_list[i], MB_TAG_DENSE | MB_TAG_CREAT ); 00158 if( MB_SUCCESS != rval ) return rval; 00159 } 00160 } 00161 00162 return MB_SUCCESS; 00163 } 00164 00165 ErrorCode ReadTetGen::read_line( std::istream& file, std::string& line, int& lineno ) 00166 { 00167 // Loop until we find a non-empty line 00168 do 00169 { 00170 // Read a line 00171 line.clear(); 00172 if( !getline( file, line ) ) return MB_FILE_WRITE_ERROR; 00173 ++lineno; 00174 // Strip comments from line 00175 size_t pos = line.find_first_of( '#' ); 00176 if( pos != std::string::npos ) line = line.substr( 0, pos ); 00177 // Strip leading whitespace from line 00178 for( pos = 0; pos < line.length() && isspace( line[pos] ); ++pos ) 00179 ; 00180 if( pos == line.length() ) 00181 line.clear(); 00182 else if( pos != 0 ) 00183 line = line.substr( pos ); 00184 } while( line.empty() ); 00185 00186 return MB_SUCCESS; 00187 } 00188 00189 ErrorCode ReadTetGen::read_line( std::istream& file, double* values_out, int num_values, int& lineno ) 00190 { 00191 // Get a line of text 00192 std::string line; 00193 ErrorCode rval = read_line( file, line, lineno ); 00194 if( MB_SUCCESS != rval ) return rval; 00195 00196 // Tokenize line as doubles 00197 std::stringstream str( line ); 00198 for( int i = 0; i < num_values; ++i ) 00199 { 00200 double v; 00201 if( !( str >> v ) ) { MB_SET_ERR( MB_FAILURE, "Error reading node data at line " << lineno ); } 00202 values_out[i] = v; 00203 } 00204 00205 // Check that we're at the end of the line 00206 int junk; 00207 if( ( str >> junk ) || !str.eof() ) 00208 { MB_SET_ERR( MB_FAILURE, "Unexpected trailing data for line " << lineno << " of node data" ); } 00209 00210 return MB_SUCCESS; 00211 } 00212 00213 ErrorCode ReadTetGen::read_node_file( std::istream& file, const Tag* attr_tag_list, const int* attr_tag_index, 00214 int attr_tag_list_len, std::vector< EntityHandle >& nodes ) 00215 { 00216 int lineno = 0; 00217 ErrorCode rval; 00218 00219 double header_vals[4]; 00220 rval = read_line( file, header_vals, 4, lineno ); 00221 if( MB_SUCCESS != rval ) return rval; 00222 00223 const int num_vtx = (int)header_vals[0]; 00224 const int dim = (int)header_vals[1]; 00225 const int num_attr = (int)header_vals[2]; 00226 const int bdry_flag = (int)header_vals[3]; 00227 if( num_vtx < 1 || dim < 2 || dim > 3 || num_attr < 0 || bdry_flag < 0 || bdry_flag > 1 ) 00228 { MB_SET_ERR( MB_FAILURE, "Invalid header line for node data" ); } 00229 if( attr_tag_list_len > num_attr ) attr_tag_list_len = num_attr; 00230 00231 // Allocate space for tag data 00232 std::map< Tag, int > tag_size; 00233 std::map< Tag, std::vector< double > > tag_data; 00234 for( int i = 0; i < attr_tag_list_len; ++i ) 00235 { 00236 if( !attr_tag_list[i] || attr_tag_index[i] < 0 ) continue; 00237 std::vector< double >& data = tag_data[attr_tag_list[i]]; 00238 // Increase tag size by one value per vertex for each time 00239 // we encounter it in the list. 00240 data.resize( data.size() + num_vtx ); 00241 ++tag_size[attr_tag_list[i]]; 00242 } 00243 std::vector< double* > attr_data( attr_tag_list_len ); 00244 std::vector< int > attr_size( attr_tag_list_len ); 00245 for( int i = 0; i < attr_tag_list_len; ++i ) 00246 { 00247 if( !attr_tag_list[i] || attr_tag_index[i] < 0 ) 00248 { 00249 attr_data[i] = 0; 00250 attr_size[i] = 0; 00251 } 00252 else 00253 { 00254 attr_data[i] = &( tag_data[attr_tag_list[i]] )[0]; 00255 attr_size[i] = tag_size[attr_tag_list[i]]; 00256 } 00257 } 00258 00259 // Allocate vertices 00260 std::vector< double* > coords; 00261 EntityHandle start_handle; 00262 rval = readTool->get_node_coords( dim, num_vtx, 1, start_handle, coords ); 00263 if( MB_SUCCESS != rval ) return rval; 00264 00265 // Read data line for each node 00266 nodes.reserve( num_vtx ); 00267 std::vector< double > data( 1 + dim + num_attr + bdry_flag ); 00268 std::vector< int > ids( num_vtx ); 00269 for( int i = 0; i < num_vtx; ++i ) 00270 { 00271 rval = read_line( file, &data[0], data.size(), lineno ); 00272 if( MB_SUCCESS != rval ) return rval; 00273 00274 // Get ID 00275 ids[i] = (int)data[0]; 00276 if( ids[i] >= (int)nodes.size() ) nodes.resize( ids[i] + 1 ); 00277 nodes[ids[i]] = start_handle + i; 00278 00279 // Get coordinates 00280 // Cppcheck warning (false positive): variable coords is assigned a value that is never used 00281 for( int j = 0; j < dim; ++j ) 00282 coords[j][i] = data[j + 1]; 00283 00284 // Get attribute data 00285 for( int j = 0; j < attr_tag_list_len; ++j ) 00286 if( attr_data[j] ) attr_data[j][i * attr_size[j] + attr_tag_index[j]] = data[j + 1 + dim]; 00287 00288 // Discard boundary bit 00289 } 00290 00291 // Store tag data 00292 Range node_range; 00293 node_range.insert( start_handle, start_handle + num_vtx - 1 ); 00294 for( std::map< Tag, std::vector< double > >::iterator i = tag_data.begin(); i != tag_data.end(); ++i ) 00295 { 00296 rval = mbIface->tag_set_data( i->first, node_range, &i->second[0] ); 00297 if( MB_SUCCESS != rval ) return rval; 00298 } 00299 Tag idtag = mbIface->globalId_tag(); 00300 00301 rval = mbIface->tag_set_data( idtag, node_range, &ids[0] ); 00302 if( MB_SUCCESS != rval ) return rval; 00303 00304 return MB_SUCCESS; 00305 } 00306 00307 ErrorCode ReadTetGen::read_elem_file( EntityType type, std::istream& file, const std::vector< EntityHandle >& nodes, 00308 Range& elems ) 00309 { 00310 int lineno = 0; 00311 ErrorCode rval; 00312 00313 int node_per_elem, have_group_id, dim; 00314 double header_vals[3]; 00315 switch( type ) 00316 { 00317 case MBTET: 00318 rval = read_line( file, header_vals, 3, lineno ); 00319 node_per_elem = (int)header_vals[1]; 00320 have_group_id = (int)header_vals[2]; 00321 dim = 3; 00322 break; 00323 case MBTRI: 00324 rval = read_line( file, header_vals, 2, lineno ); 00325 node_per_elem = 3; 00326 have_group_id = (int)header_vals[1]; 00327 dim = 2; 00328 break; 00329 case MBEDGE: 00330 rval = read_line( file, header_vals, 1, lineno ); 00331 node_per_elem = 2; 00332 have_group_id = 0; 00333 dim = 1; 00334 break; 00335 default: 00336 rval = MB_FAILURE; 00337 break; 00338 } 00339 if( MB_SUCCESS != rval ) return rval; 00340 const int num_elem = (int)header_vals[0]; 00341 if( num_elem < 1 || node_per_elem < 2 || have_group_id < 0 || have_group_id > 1 ) 00342 { MB_SET_ERR( MB_FAILURE, "Invalid header line for element data" ); } 00343 00344 // Create group map 00345 std::map< double, EntityHandle > groups; 00346 Tag dim_tag, id_tag; 00347 id_tag = mbIface->globalId_tag(); 00348 00349 const int negone = -1; 00350 rval = mbIface->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag, MB_TAG_SPARSE | MB_TAG_CREAT, 00351 &negone ); 00352 if( MB_SUCCESS != rval ) return rval; 00353 00354 // Allocate elements 00355 EntityHandle start_handle, *conn_array; 00356 rval = readTool->get_element_connect( num_elem, node_per_elem, type, 1, start_handle, conn_array ); 00357 if( MB_SUCCESS != rval ) return rval; 00358 elems.insert( start_handle, start_handle + num_elem - 1 ); 00359 00360 // Read data line for each node 00361 std::vector< double > data( 1 + node_per_elem + have_group_id ); 00362 std::vector< int > ids( num_elem ); 00363 for( int i = 0; i < num_elem; ++i ) 00364 { 00365 rval = read_line( file, &data[0], data.size(), lineno ); 00366 if( MB_SUCCESS != rval ) return rval; 00367 00368 // Get ID 00369 ids[i] = (int)data[0]; 00370 00371 // Get connectivity 00372 for( int j = 0; j < node_per_elem; ++j ) 00373 conn_array[node_per_elem * i + j] = nodes[(int)data[j + 1]]; 00374 00375 // Grouping 00376 if( have_group_id && 0.0 != data[node_per_elem + 1] ) 00377 { 00378 double id = data[node_per_elem + 1]; 00379 EntityHandle grp = groups[id]; 00380 if( 0 == grp ) 00381 { 00382 rval = mbIface->create_meshset( MESHSET_SET, grp ); 00383 if( MB_SUCCESS != rval ) return rval; 00384 elems.insert( grp ); 00385 rval = mbIface->tag_set_data( dim_tag, &grp, 1, &dim ); 00386 if( MB_SUCCESS != rval ) return rval; 00387 int iid = (int)id; 00388 rval = mbIface->tag_set_data( id_tag, &grp, 1, &iid ); 00389 if( MB_SUCCESS != rval ) return rval; 00390 groups[id] = grp; 00391 } 00392 EntityHandle handle = start_handle + i; 00393 rval = mbIface->add_entities( grp, &handle, 1 ); 00394 if( MB_SUCCESS != rval ) return rval; 00395 } 00396 } 00397 00398 // Store id data 00399 Range elems2; 00400 elems2.insert( start_handle, start_handle + num_elem - 1 ); 00401 rval = mbIface->tag_set_data( id_tag, elems2, &ids[0] ); 00402 if( MB_SUCCESS != rval ) return rval; 00403 00404 return MB_SUCCESS; 00405 } 00406 00407 } // namespace moab