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