MOAB: Mesh Oriented datABase  (version 5.3.0)
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         {
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines