MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /** 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in this software. 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 2.1 of the License, or (at your option) any later version. 00013 * 00014 */ 00015 00016 /** 00017 * \class ReadGmsh 00018 * \brief Gmsh (http://www.geuz.org/gmsh) file reader 00019 * 00020 * See: http://geuz.org/gmsh/doc/texinfo/gmsh.html#MSH-ASCII-file-format 00021 * 00022 * \author Jason Kraftcheck 00023 */ 00024 00025 #include "ReadGmsh.hpp" 00026 #include "FileTokenizer.hpp" // for file tokenizer 00027 #include "Internals.hpp" 00028 #include "moab/Interface.hpp" 00029 #include "moab/ReadUtilIface.hpp" 00030 #include "moab/Range.hpp" 00031 #include "MBTagConventions.hpp" 00032 #include "MBParallelConventions.h" 00033 #include "moab/CN.hpp" 00034 #include "GmshUtil.hpp" 00035 00036 #include <cerrno> 00037 #include <cstring> 00038 #include <map> 00039 #include <set> 00040 00041 namespace moab 00042 { 00043 00044 ReaderIface* ReadGmsh::factory( Interface* iface ) 00045 { 00046 return new ReadGmsh( iface ); 00047 } 00048 00049 ReadGmsh::ReadGmsh( Interface* impl ) : mdbImpl( impl ), globalId( 0 ) 00050 { 00051 mdbImpl->query_interface( readMeshIface ); 00052 } 00053 00054 ReadGmsh::~ReadGmsh() 00055 { 00056 if( readMeshIface ) 00057 { 00058 mdbImpl->release_interface( readMeshIface ); 00059 readMeshIface = 0; 00060 } 00061 } 00062 00063 ErrorCode ReadGmsh::read_tag_values( const char* /* file_name */, const char* /* tag_name */, 00064 const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */, 00065 const SubsetList* /* subset_list */ ) 00066 { 00067 return MB_NOT_IMPLEMENTED; 00068 } 00069 00070 ErrorCode ReadGmsh::load_file( const char* filename, const EntityHandle*, const FileOptions&, 00071 const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag ) 00072 { 00073 int num_material_sets = 0; 00074 const int* material_set_list = 0; 00075 00076 if( subset_list ) 00077 { 00078 if( subset_list->tag_list_length > 1 && !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) ) 00079 { 00080 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "GMsh supports subset read only by material ID" ); 00081 } 00082 material_set_list = subset_list->tag_list[0].tag_values; 00083 num_material_sets = subset_list->tag_list[0].num_tag_values; 00084 } 00085 00086 geomSets.clear(); 00087 globalId = mdbImpl->globalId_tag(); 00088 00089 // Create set for more convenient check for material set ids 00090 std::set< int > blocks; 00091 for( const int* mat_set_end = material_set_list + num_material_sets; material_set_list != mat_set_end; 00092 ++material_set_list ) 00093 blocks.insert( *material_set_list ); 00094 00095 // Map of ID->handle for nodes 00096 std::map< long, EntityHandle > node_id_map; 00097 int data_size = 8; 00098 00099 // Open file and hand off pointer to tokenizer 00100 FILE* file_ptr = fopen( filename, "r" ); 00101 if( !file_ptr ) { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) ); } 00102 FileTokenizer tokens( file_ptr, readMeshIface ); 00103 00104 // Determine file format version 00105 const char* const start_tokens[] = { "$NOD", "$MeshFormat", 0 }; 00106 int format_version = tokens.match_token( start_tokens ); 00107 if( !format_version ) return MB_FILE_DOES_NOT_EXIST; 00108 00109 // If version 2.0, read additional header info 00110 if( 2 == format_version ) 00111 { 00112 double version; 00113 if( !tokens.get_doubles( 1, &version ) ) return MB_FILE_WRITE_ERROR; 00114 00115 if( version != 2.0 && version != 2.1 && version != 2.2 ) 00116 { 00117 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": unknown format version: " << version ); 00118 return MB_FILE_DOES_NOT_EXIST; 00119 } 00120 00121 int file_format; 00122 if( !tokens.get_integers( 1, &file_format ) || !tokens.get_integers( 1, &data_size ) || 00123 !tokens.match_token( "$EndMeshFormat" ) ) 00124 return MB_FILE_WRITE_ERROR; 00125 // If physical entities in the gmsh file -> discard this 00126 const char* const phys_tokens[] = { "$Nodes", "$PhysicalNames", 0 }; 00127 int hasPhys = tokens.match_token( phys_tokens ); 00128 00129 if( hasPhys == 2 ) 00130 { 00131 long num_phys; 00132 if( !tokens.get_long_ints( 1, &num_phys ) ) return MB_FILE_WRITE_ERROR; 00133 for( long loop_phys = 0; loop_phys < num_phys; loop_phys++ ) 00134 { 00135 long physDim; 00136 long physGroupNum; 00137 // char const * physName; 00138 if( !tokens.get_long_ints( 1, &physDim ) ) return MB_FILE_WRITE_ERROR; 00139 if( !tokens.get_long_ints( 1, &physGroupNum ) ) return MB_FILE_WRITE_ERROR; 00140 const char* ptc = tokens.get_string(); 00141 if( !ptc ) return MB_FILE_WRITE_ERROR; 00142 // try to get to the end of the line, without reporting errors 00143 // really, we need to skip this 00144 while( !tokens.get_newline( false ) ) 00145 ptc = tokens.get_string(); 00146 } 00147 if( !tokens.match_token( "$EndPhysicalNames" ) || !tokens.match_token( "$Nodes" ) ) 00148 return MB_FILE_WRITE_ERROR; 00149 } 00150 } 00151 00152 // Read number of nodes 00153 long num_nodes; 00154 if( !tokens.get_long_ints( 1, &num_nodes ) ) return MB_FILE_WRITE_ERROR; 00155 00156 // Allocate nodes 00157 std::vector< double* > coord_arrays; 00158 EntityHandle handle = 0; 00159 ErrorCode result = readMeshIface->get_node_coords( 3, num_nodes, MB_START_ID, handle, coord_arrays ); 00160 if( MB_SUCCESS != result ) return result; 00161 00162 // Read nodes 00163 double *x = coord_arrays[0], *y = coord_arrays[1], *z = coord_arrays[2]; 00164 for( long i = 0; i < num_nodes; ++i, ++handle ) 00165 { 00166 long id; 00167 if( !tokens.get_long_ints( 1, &id ) || !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) || 00168 !tokens.get_doubles( 1, z++ ) ) 00169 return MB_FILE_WRITE_ERROR; 00170 00171 if( !node_id_map.insert( std::pair< long, EntityHandle >( id, handle ) ).second ) 00172 { 00173 MB_SET_ERR( MB_FILE_WRITE_ERROR, "Duplicate node ID at line " << tokens.line_number() ); 00174 } 00175 } 00176 00177 // Create reverse map from handle to id 00178 std::vector< int > ids( num_nodes ); 00179 std::vector< int >::iterator id_iter = ids.begin(); 00180 std::vector< EntityHandle > handles( num_nodes ); 00181 std::vector< EntityHandle >::iterator h_iter = handles.begin(); 00182 for( std::map< long, EntityHandle >::iterator i = node_id_map.begin(); i != node_id_map.end(); 00183 ++i, ++id_iter, ++h_iter ) 00184 { 00185 *id_iter = i->first; 00186 *h_iter = i->second; 00187 } 00188 // Store IDs in tags 00189 result = mdbImpl->tag_set_data( globalId, &handles[0], num_nodes, &ids[0] ); 00190 if( MB_SUCCESS != result ) return result; 00191 if( file_id_tag ) 00192 { 00193 result = mdbImpl->tag_set_data( *file_id_tag, &handles[0], num_nodes, &ids[0] ); 00194 if( MB_SUCCESS != result ) return result; 00195 } 00196 ids.clear(); 00197 handles.clear(); 00198 00199 // Get tokens signifying end of node data and start of elements 00200 if( !tokens.match_token( format_version == 1 ? "$ENDNOD" : "$EndNodes" ) || 00201 !tokens.match_token( format_version == 1 ? "$ELM" : "$Elements" ) ) 00202 return MB_FILE_WRITE_ERROR; 00203 00204 // Get element count 00205 long num_elem; 00206 if( !tokens.get_long_ints( 1, &num_elem ) ) return MB_FILE_WRITE_ERROR; 00207 00208 // Lists of data accumulated for elements 00209 std::vector< EntityHandle > connectivity; 00210 std::vector< int > mat_set_list, geom_set_list, part_set_list, id_list; 00211 // Temporary, per-element data 00212 std::vector< int > int_data( 5 ), tag_data( 2 ); 00213 std::vector< long > tmp_conn; 00214 int curr_elem_type = -1; 00215 for( long i = 0; i < num_elem; ++i ) 00216 { 00217 // Read element description 00218 // File format 1.0 00219 if( 1 == format_version ) 00220 { 00221 if( !tokens.get_integers( 5, &int_data[0] ) ) return MB_FILE_WRITE_ERROR; 00222 tag_data[0] = int_data[2]; 00223 tag_data[1] = int_data[3]; 00224 if( (unsigned)tag_data[1] < GmshUtil::numGmshElemType && 00225 GmshUtil::gmshElemTypes[tag_data[1]].num_nodes != (unsigned)int_data[4] ) 00226 { 00227 MB_SET_ERR( MB_FILE_WRITE_ERROR, 00228 "Invalid node count for element type at line " << tokens.line_number() ); 00229 } 00230 } 00231 // File format 2.0 00232 else 00233 { 00234 if( !tokens.get_integers( 3, &int_data[0] ) ) return MB_FILE_WRITE_ERROR; 00235 tag_data.resize( int_data[2] ); 00236 if( !tokens.get_integers( tag_data.size(), &tag_data[0] ) ) return MB_FILE_WRITE_ERROR; 00237 } 00238 00239 // If a list of material sets was specified in the 00240 // argument list, skip any elements for which the 00241 // material set is not specified or is not in the 00242 // passed list. 00243 if( !blocks.empty() && ( tag_data.empty() || blocks.find( tag_data[0] ) != blocks.end() ) ) continue; 00244 00245 // If the next element is not the same type as the last one, 00246 // create a sequence for the block of elements we've read 00247 // to this point (all of the same type), and clear accumulated 00248 // data. 00249 if( int_data[1] != curr_elem_type ) 00250 { 00251 if( !id_list.empty() ) 00252 { // First iteration 00253 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list, 00254 part_set_list, connectivity, file_id_tag ); 00255 if( MB_SUCCESS != result ) return result; 00256 } 00257 00258 id_list.clear(); 00259 mat_set_list.clear(); 00260 geom_set_list.clear(); 00261 part_set_list.clear(); 00262 connectivity.clear(); 00263 curr_elem_type = int_data[1]; 00264 if( (unsigned)curr_elem_type >= GmshUtil::numGmshElemType || 00265 GmshUtil::gmshElemTypes[curr_elem_type].mb_type == MBMAXTYPE ) 00266 { 00267 MB_SET_ERR( MB_FILE_WRITE_ERROR, 00268 "Unsupported element type " << curr_elem_type << " at line " << tokens.line_number() ); 00269 } 00270 tmp_conn.resize( GmshUtil::gmshElemTypes[curr_elem_type].num_nodes ); 00271 } 00272 00273 // Store data from element description 00274 id_list.push_back( int_data[0] ); 00275 if( tag_data.size() > 3 ) 00276 part_set_list.push_back( tag_data[3] ); // it must be new format for gmsh, >= 2.5 00277 // it could have negative partition ids, for ghost elements 00278 else if( tag_data.size() > 2 ) 00279 part_set_list.push_back( tag_data[2] ); // old format, partition id saved in 3rd tag field 00280 else 00281 part_set_list.push_back( 0 ); 00282 geom_set_list.push_back( tag_data.size() > 1 ? tag_data[1] : 0 ); 00283 mat_set_list.push_back( tag_data.size() > 0 ? tag_data[0] : 0 ); 00284 00285 // Get element connectivity 00286 if( !tokens.get_long_ints( tmp_conn.size(), &tmp_conn[0] ) ) return MB_FILE_WRITE_ERROR; 00287 00288 // Convert connectivity from IDs to handles 00289 for( unsigned j = 0; j < tmp_conn.size(); ++j ) 00290 { 00291 std::map< long, EntityHandle >::iterator k = node_id_map.find( tmp_conn[j] ); 00292 if( k == node_id_map.end() ) 00293 { 00294 MB_SET_ERR( MB_FILE_WRITE_ERROR, "Invalid node ID at line " << tokens.line_number() ); 00295 } 00296 connectivity.push_back( k->second ); 00297 } 00298 } // for (num_nodes) 00299 00300 // Create entity sequence for last element(s). 00301 if( !id_list.empty() ) 00302 { 00303 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list, 00304 part_set_list, connectivity, file_id_tag ); 00305 if( MB_SUCCESS != result ) return result; 00306 } 00307 00308 // Construct parent-child relations for geometric sets. 00309 // Note: At the time this comment was written, the following 00310 // function was not implemented. 00311 result = create_geometric_topology(); 00312 geomSets.clear(); 00313 return result; 00314 } 00315 00316 //! Create an element sequence 00317 ErrorCode ReadGmsh::create_elements( const GmshElemType& type, const std::vector< int >& elem_ids, 00318 const std::vector< int >& matl_ids, const std::vector< int >& geom_ids, 00319 const std::vector< int >& prtn_ids, 00320 const std::vector< EntityHandle >& connectivity, const Tag* file_id_tag ) 00321 { 00322 ErrorCode result; 00323 00324 // Make sure input is consistent 00325 const unsigned long num_elem = elem_ids.size(); 00326 const int node_per_elem = type.num_nodes; 00327 if( matl_ids.size() != num_elem || geom_ids.size() != num_elem || prtn_ids.size() != num_elem || 00328 connectivity.size() != num_elem * node_per_elem ) 00329 return MB_FAILURE; 00330 00331 // Create the element sequence 00332 // for points, simply gather the connectivities and create the materials 00333 if( type.mb_type == MBVERTEX ) 00334 { 00335 Range elements; 00336 elements.insert< std::vector< EntityHandle > >( connectivity.begin(), connectivity.end() ); 00337 result = create_sets( type.mb_type, elements, matl_ids, 0 ); 00338 if( MB_SUCCESS != result ) return result; 00339 00340 return MB_SUCCESS; 00341 } 00342 EntityHandle handle = 0; 00343 EntityHandle* conn_array; 00344 result = 00345 readMeshIface->get_element_connect( num_elem, node_per_elem, type.mb_type, MB_START_ID, handle, conn_array ); 00346 if( MB_SUCCESS != result ) return result; 00347 00348 // Copy passed element connectivity into entity sequence data. 00349 if( type.node_order ) 00350 { 00351 for( unsigned long i = 0; i < num_elem; ++i ) 00352 for( int j = 0; j < node_per_elem; ++j ) 00353 conn_array[i * node_per_elem + type.node_order[j]] = connectivity[i * node_per_elem + j]; 00354 } 00355 else 00356 { 00357 memcpy( conn_array, &connectivity[0], connectivity.size() * sizeof( EntityHandle ) ); 00358 } 00359 00360 // Notify MOAB of the new elements 00361 result = readMeshIface->update_adjacencies( handle, num_elem, node_per_elem, conn_array ); 00362 if( MB_SUCCESS != result ) return result; 00363 00364 // Store element IDs 00365 Range elements( handle, handle + num_elem - 1 ); 00366 result = mdbImpl->tag_set_data( globalId, elements, &elem_ids[0] ); 00367 if( MB_SUCCESS != result ) return result; 00368 if( file_id_tag ) 00369 { 00370 result = mdbImpl->tag_set_data( *file_id_tag, elements, &elem_ids[0] ); 00371 if( MB_SUCCESS != result ) return result; 00372 } 00373 00374 // Add elements to material sets 00375 result = create_sets( type.mb_type, elements, matl_ids, 0 ); 00376 if( MB_SUCCESS != result ) return result; 00377 // Add elements to geometric sets 00378 result = create_sets( type.mb_type, elements, geom_ids, 1 ); 00379 if( MB_SUCCESS != result ) return result; 00380 // Add elements to parallel partitions 00381 result = create_sets( type.mb_type, elements, prtn_ids, 2 ); 00382 if( MB_SUCCESS != result ) return result; 00383 00384 return MB_SUCCESS; 00385 } 00386 00387 //! Add elements to sets as dictated by grouping ID in file. 00388 ErrorCode ReadGmsh::create_sets( EntityType type, const Range& elements, const std::vector< int >& set_ids, 00389 int set_type ) 00390 { 00391 ErrorCode result; 00392 00393 // Get a unique list of set IDs 00394 std::set< int > ids; 00395 for( std::vector< int >::const_iterator i = set_ids.begin(); i != set_ids.end(); ++i ) 00396 ids.insert( *i ); 00397 00398 // No Sets? 00399 if( ids.empty() || ( ids.size() == 1 && *ids.begin() == 0 ) ) return MB_SUCCESS; // no sets (all ids are zero) 00400 00401 // Get/create tag handles 00402 int num_tags; 00403 Tag tag_handles[2]; 00404 int tag_val; 00405 const void* tag_values[2] = { &tag_val, NULL }; 00406 00407 switch( set_type ) 00408 { 00409 default: 00410 return MB_FAILURE; 00411 case 0: 00412 case 2: { 00413 const char* name = set_type ? PARALLEL_PARTITION_TAG_NAME : MATERIAL_SET_TAG_NAME; 00414 result = mdbImpl->tag_get_handle( name, 1, MB_TYPE_INTEGER, tag_handles[0], MB_TAG_SPARSE | MB_TAG_CREAT ); 00415 if( MB_SUCCESS != result ) return result; 00416 num_tags = 1; 00417 break; 00418 } 00419 case 1: { 00420 result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handles[1], 00421 MB_TAG_SPARSE | MB_TAG_CREAT ); 00422 if( MB_SUCCESS != result ) return result; 00423 tag_values[1] = NULL; 00424 tag_handles[0] = globalId; 00425 num_tags = 2; 00426 break; 00427 } 00428 } // switch 00429 00430 // For each unique set ID... 00431 for( std::set< int >::iterator i = ids.begin(); i != ids.end(); ++i ) 00432 { 00433 // Skip "null" set ID 00434 if( *i == 0 ) continue; 00435 00436 // Get all entities with the current set ID 00437 Range entities, sets; 00438 std::vector< int >::const_iterator j = set_ids.begin(); 00439 for( Range::iterator k = elements.begin(); k != elements.end(); ++j, ++k ) 00440 if( *i == *j ) entities.insert( *k ); 00441 00442 // Get set by ID 00443 // Cppcheck warning (false positive): variable tag_val is assigned a value that is never 00444 // used 00445 tag_val = *i; 00446 result = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, tag_handles, tag_values, num_tags, sets ); 00447 if( MB_SUCCESS != result && MB_ENTITY_NOT_FOUND != result ) return result; 00448 00449 // Don't use existing geometry sets (from some other file) 00450 if( 1 == set_type ) // Geometry 00451 sets = intersect( sets, geomSets ); 00452 00453 // Get set handle 00454 EntityHandle set; 00455 // If no sets with ID, create one 00456 if( sets.empty() ) 00457 { 00458 result = mdbImpl->create_meshset( MESHSET_SET, set ); 00459 if( MB_SUCCESS != result ) return result; 00460 00461 result = mdbImpl->tag_set_data( tag_handles[0], &set, 1, &*i ); 00462 if( MB_SUCCESS != result ) return result; 00463 00464 if( 1 == set_type ) 00465 { // Geometry 00466 int dim = CN::Dimension( type ); 00467 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim ); 00468 if( MB_SUCCESS != result ) return result; 00469 geomSets.insert( set ); 00470 } 00471 } 00472 else 00473 { 00474 set = *sets.begin(); 00475 if( 1 == set_type ) 00476 { // Geometry 00477 int dim = CN::Dimension( type ); 00478 // Get dimension of set 00479 int dim2; 00480 result = mdbImpl->tag_get_data( tag_handles[1], &set, 1, &dim2 ); 00481 if( MB_SUCCESS != result ) return result; 00482 // If we're putting geometry of a higher dimension into the 00483 // set, increase the dimension of the set. 00484 if( dim > dim2 ) 00485 { 00486 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim ); 00487 if( MB_SUCCESS != result ) return result; 00488 } 00489 } 00490 } 00491 00492 // Put the mesh entities into the set 00493 result = mdbImpl->add_entities( set, entities ); 00494 if( MB_SUCCESS != result ) return result; 00495 } // for (ids) 00496 00497 return MB_SUCCESS; 00498 } 00499 00500 //! NOT IMPLEMENTED 00501 //! Reconstruct parent-child relations for geometry sets from 00502 //! mesh connectivity. 00503 ErrorCode ReadGmsh::create_geometric_topology() 00504 { 00505 if( geomSets.empty() ) return MB_SUCCESS; 00506 00507 // Not implemented yet 00508 geomSets.clear(); 00509 return MB_SUCCESS; 00510 } 00511 00512 } // namespace moab