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