ReadGmsh.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
/**
 * \class ReadGmsh
 * \brief Gmsh (http://www.geuz.org/gmsh) file reader
 *
 * See: http://geuz.org/gmsh/doc/texinfo/gmsh.html#MSH-ASCII-file-format
 *
 * \author Jason Kraftcheck
 */
 
#include "ReadGmsh.hpp"
#include "FileTokenizer.hpp" // for file tokenizer
#include "Internals.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/CN.hpp"
#include "GmshUtil.hpp"
 
#include <cerrno>
#include <cstring>
#include <map>
#include <set>
 
namespace moab
{
 
ReaderIface* ReadGmsh::factory( Interface* iface )
{
 return new ReadGmsh( iface );
}
 
ReadGmsh::ReadGmsh( Interface* impl ) : mdbImpl( impl ), globalId( 0 )
{
 mdbImpl->query_interface( readMeshIface );
}
 
ReadGmsh::~ReadGmsh()
{
 if( readMeshIface )
 {
 mdbImpl->release_interface( readMeshIface );
 readMeshIface = 0;
 }
}
 
ErrorCode ReadGmsh::read_tag_values( const char* /* file_name */,
 const char* /* tag_name */,
 const FileOptions& /* opts */,
 std::vector< int >& /* tag_values_out */,
 const SubsetList* /* subset_list */ )
{
 return MB_NOT_IMPLEMENTED;
}
 
ErrorCode ReadGmsh::load_file( const char* filename,
 const EntityHandle*,
 const FileOptions&,
 const ReaderIface::SubsetList* subset_list,
 const Tag* file_id_tag )
{
 int num_material_sets = 0;
 const int* material_set_list = 0;
 
 if( subset_list )
 {
 if( subset_list->tag_list_length > 1 && !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
 {
 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "GMsh supports subset read only by material ID" );
 }
 material_set_list = subset_list->tag_list[0].tag_values;
 num_material_sets = subset_list->tag_list[0].num_tag_values;
 }
 
 geomSets.clear();
 globalId = mdbImpl->globalId_tag();
 
 // Create set for more convenient check for material set ids
 std::set< int > blocks;
 for( const int* mat_set_end = material_set_list + num_material_sets; material_set_list != mat_set_end;
 ++material_set_list )
 blocks.insert( *material_set_list );
 
 // Map of ID->handle for nodes
 std::map< long, EntityHandle > node_id_map;
 int data_size = 8;
 
 // Open file and hand off pointer to tokenizer
 FILE* file_ptr = fopen( filename, "r" );
 if( !file_ptr )
 {
 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) );
 }
 FileTokenizer tokens( file_ptr, readMeshIface );
 
 // Determine file format version
 const char* const start_tokens[] = { "$NOD", "$MeshFormat", 0 };
 int format_version = tokens.match_token( start_tokens );
 if( !format_version ) return MB_FILE_DOES_NOT_EXIST;
 
 // If version 2.0, read additional header info
 if( 2 == format_version )
 {
 double version;
 if( !tokens.get_doubles( 1, &version ) ) return MB_FILE_WRITE_ERROR;
 
 if( version != 2.0 && version != 2.1 && version != 2.2 )
 {
 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": unknown format version: " << version );
 return MB_FILE_DOES_NOT_EXIST;
 }
 
 int file_format;
 if( !tokens.get_integers( 1, &file_format ) || !tokens.get_integers( 1, &data_size ) ||
 !tokens.match_token( "$EndMeshFormat" ) )
 return MB_FILE_WRITE_ERROR;
 // If physical entities in the gmsh file -> discard this
 const char* const phys_tokens[] = { "$Nodes", "$PhysicalNames", 0 };
 int hasPhys = tokens.match_token( phys_tokens );
 
 if( hasPhys == 2 )
 {
 long num_phys;
 if( !tokens.get_long_ints( 1, &num_phys ) ) return MB_FILE_WRITE_ERROR;
 for( long loop_phys = 0; loop_phys < num_phys; loop_phys++ )
 {
 long physDim;
 long physGroupNum;
 // char const * physName;
 if( !tokens.get_long_ints( 1, &physDim ) ) return MB_FILE_WRITE_ERROR;
 if( !tokens.get_long_ints( 1, &physGroupNum ) ) return MB_FILE_WRITE_ERROR;
 const char* ptc = tokens.get_string();
 if( !ptc ) return MB_FILE_WRITE_ERROR;
 // try to get to the end of the line, without reporting errors
 // really, we need to skip this
 while( !tokens.get_newline( false ) )
 ptc = tokens.get_string();
 }
 if( !tokens.match_token( "$EndPhysicalNames" ) || !tokens.match_token( "$Nodes" ) )
 return MB_FILE_WRITE_ERROR;
 }
 }
 
 // Read number of nodes
 long num_nodes;
 if( !tokens.get_long_ints( 1, &num_nodes ) ) return MB_FILE_WRITE_ERROR;
 
 // Allocate nodes
 std::vector< double* > coord_arrays;
 EntityHandle handle = 0;
 ErrorCode result = readMeshIface->get_node_coords( 3, num_nodes, MB_START_ID, handle, coord_arrays );
 if( MB_SUCCESS != result ) return result;
 
 // Read nodes
 double *x = coord_arrays[0], *y = coord_arrays[1], *z = coord_arrays[2];
 for( long i = 0; i < num_nodes; ++i, ++handle )
 {
 long id;
 if( !tokens.get_long_ints( 1, &id ) || !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) ||
 !tokens.get_doubles( 1, z++ ) )
 return MB_FILE_WRITE_ERROR;
 
 if( !node_id_map.insert( std::pair< long, EntityHandle >( id, handle ) ).second )
 {
 MB_SET_ERR( MB_FILE_WRITE_ERROR, "Duplicate node ID at line " << tokens.line_number() );
 }
 }
 
 // Create reverse map from handle to id
 std::vector< int > ids( num_nodes );
 std::vector< int >::iterator id_iter = ids.begin();
 std::vector< EntityHandle > handles( num_nodes );
 std::vector< EntityHandle >::iterator h_iter = handles.begin();
 for( std::map< long, EntityHandle >::iterator i = node_id_map.begin(); i != node_id_map.end();
 ++i, ++id_iter, ++h_iter )
 {
 *id_iter = i->first;
 *h_iter = i->second;
 }
 // Store IDs in tags
 result = mdbImpl->tag_set_data( globalId, &handles[0], num_nodes, &ids[0] );
 if( MB_SUCCESS != result ) return result;
 if( file_id_tag )
 {
 result = mdbImpl->tag_set_data( *file_id_tag, &handles[0], num_nodes, &ids[0] );
 if( MB_SUCCESS != result ) return result;
 }
 ids.clear();
 handles.clear();
 
 // Get tokens signifying end of node data and start of elements
 if( !tokens.match_token( format_version == 1 ? "$ENDNOD" : "$EndNodes" ) ||
 !tokens.match_token( format_version == 1 ? "$ELM" : "$Elements" ) )
 return MB_FILE_WRITE_ERROR;
 
 // Get element count
 long num_elem;
 if( !tokens.get_long_ints( 1, &num_elem ) ) return MB_FILE_WRITE_ERROR;
 
 // Lists of data accumulated for elements
 std::vector< EntityHandle > connectivity;
 std::vector< int > mat_set_list, geom_set_list, part_set_list, id_list;
 // Temporary, per-element data
 std::vector< int > int_data( 5 ), tag_data( 2 );
 std::vector< long > tmp_conn;
 int curr_elem_type = -1;
 for( long i = 0; i < num_elem; ++i )
 {
 // Read element description
 // File format 1.0
 if( 1 == format_version )
 {
 if( !tokens.get_integers( 5, &int_data[0] ) ) return MB_FILE_WRITE_ERROR;
 tag_data[0] = int_data[2];
 tag_data[1] = int_data[3];
 if( (unsigned)tag_data[1] < GmshUtil::numGmshElemType &&
 GmshUtil::gmshElemTypes[tag_data[1]].num_nodes != (unsigned)int_data[4] )
 {
 MB_SET_ERR( MB_FILE_WRITE_ERROR,
 "Invalid node count for element type at line " << tokens.line_number() );
 }
 }
 // File format 2.0
 else
 {
 if( !tokens.get_integers( 3, &int_data[0] ) ) return MB_FILE_WRITE_ERROR;
 tag_data.resize( int_data[2] );
 if( !tokens.get_integers( tag_data.size(), &tag_data[0] ) ) return MB_FILE_WRITE_ERROR;
 }
 
 // If a list of material sets was specified in the
 // argument list, skip any elements for which the
 // material set is not specified or is not in the
 // passed list.
 if( !blocks.empty() && ( tag_data.empty() || blocks.find( tag_data[0] ) != blocks.end() ) ) continue;
 
 // If the next element is not the same type as the last one,
 // create a sequence for the block of elements we've read
 // to this point (all of the same type), and clear accumulated
 // data.
 if( int_data[1] != curr_elem_type )
 {
 if( !id_list.empty() )
 { // First iteration
 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list,
 part_set_list, connectivity, file_id_tag );
 if( MB_SUCCESS != result ) return result;
 }
 
 id_list.clear();
 mat_set_list.clear();
 geom_set_list.clear();
 part_set_list.clear();
 connectivity.clear();
 curr_elem_type = int_data[1];
 if( (unsigned)curr_elem_type >= GmshUtil::numGmshElemType ||
 GmshUtil::gmshElemTypes[curr_elem_type].mb_type == MBMAXTYPE )
 {
 MB_SET_ERR( MB_FILE_WRITE_ERROR,
 "Unsupported element type " << curr_elem_type << " at line " << tokens.line_number() );
 }
 tmp_conn.resize( GmshUtil::gmshElemTypes[curr_elem_type].num_nodes );
 }
 
 // Store data from element description
 id_list.push_back( int_data[0] );
 if( tag_data.size() > 3 )
 part_set_list.push_back( tag_data[3] ); // it must be new format for gmsh, >= 2.5
 // it could have negative partition ids, for ghost elements
 else if( tag_data.size() > 2 )
 part_set_list.push_back( tag_data[2] ); // old format, partition id saved in 3rd tag field
 else
 part_set_list.push_back( 0 );
 geom_set_list.push_back( tag_data.size() > 1 ? tag_data[1] : 0 );
 mat_set_list.push_back( tag_data.size() > 0 ? tag_data[0] : 0 );
 
 // Get element connectivity
 if( !tokens.get_long_ints( tmp_conn.size(), &tmp_conn[0] ) ) return MB_FILE_WRITE_ERROR;
 
 // Convert connectivity from IDs to handles
 for( unsigned j = 0; j < tmp_conn.size(); ++j )
 {
 std::map< long, EntityHandle >::iterator k = node_id_map.find( tmp_conn[j] );
 if( k == node_id_map.end() )
 {
 MB_SET_ERR( MB_FILE_WRITE_ERROR, "Invalid node ID at line " << tokens.line_number() );
 }
 connectivity.push_back( k->second );
 }
 } // for (num_nodes)
 
 // Create entity sequence for last element(s).
 if( !id_list.empty() )
 {
 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list,
 part_set_list, connectivity, file_id_tag );
 if( MB_SUCCESS != result ) return result;
 }
 
 // Construct parent-child relations for geometric sets.
 // Note: At the time this comment was written, the following
 // function was not implemented.
 result = create_geometric_topology();
 geomSets.clear();
 return result;
}
 
//! Create an element sequence
ErrorCode ReadGmsh::create_elements( const GmshElemType& type,
 const std::vector< int >& elem_ids,
 const std::vector< int >& matl_ids,
 const std::vector< int >& geom_ids,
 const std::vector< int >& prtn_ids,
 const std::vector< EntityHandle >& connectivity,
 const Tag* file_id_tag )
{
 ErrorCode result;
 
 // Make sure input is consistent
 const unsigned long num_elem = elem_ids.size();
 const int node_per_elem = type.num_nodes;
 if( matl_ids.size() != num_elem || geom_ids.size() != num_elem || prtn_ids.size() != num_elem ||
 connectivity.size() != num_elem * node_per_elem )
 return MB_FAILURE;
 
 // Create the element sequence
 // for points, simply gather the connectivities and create the materials
 if( type.mb_type == MBVERTEX )
 {
 Range elements;
 elements.insert< std::vector< EntityHandle > >( connectivity.begin(), connectivity.end() );
 result = create_sets( type.mb_type, elements, matl_ids, 0 );
 if( MB_SUCCESS != result ) return result;
 
 return MB_SUCCESS;
 }
 EntityHandle handle = 0;
 EntityHandle* conn_array;
 result =
 readMeshIface->get_element_connect( num_elem, node_per_elem, type.mb_type, MB_START_ID, handle, conn_array );
 if( MB_SUCCESS != result ) return result;
 
 // Copy passed element connectivity into entity sequence data.
 if( type.node_order )
 {
 for( unsigned long i = 0; i < num_elem; ++i )
 for( int j = 0; j < node_per_elem; ++j )
 conn_array[i * node_per_elem + type.node_order[j]] = connectivity[i * node_per_elem + j];
 }
 else
 {
 memcpy( conn_array, &connectivity[0], connectivity.size() * sizeof( EntityHandle ) );
 }
 
 // Notify MOAB of the new elements
 result = readMeshIface->update_adjacencies( handle, num_elem, node_per_elem, conn_array );
 if( MB_SUCCESS != result ) return result;
 
 // Store element IDs
 Range elements( handle, handle + num_elem - 1 );
 result = mdbImpl->tag_set_data( globalId, elements, &elem_ids[0] );
 if( MB_SUCCESS != result ) return result;
 if( file_id_tag )
 {
 result = mdbImpl->tag_set_data( *file_id_tag, elements, &elem_ids[0] );
 if( MB_SUCCESS != result ) return result;
 }
 
 // Add elements to material sets
 result = create_sets( type.mb_type, elements, matl_ids, 0 );
 if( MB_SUCCESS != result ) return result;
 // Add elements to geometric sets
 result = create_sets( type.mb_type, elements, geom_ids, 1 );
 if( MB_SUCCESS != result ) return result;
 // Add elements to parallel partitions
 result = create_sets( type.mb_type, elements, prtn_ids, 2 );
 if( MB_SUCCESS != result ) return result;
 
 return MB_SUCCESS;
}
 
//! Add elements to sets as dictated by grouping ID in file.
ErrorCode ReadGmsh::create_sets( EntityType type,
 const Range& elements,
 const std::vector< int >& set_ids,
 int set_type )
{
 ErrorCode result;
 
 // Get a unique list of set IDs
 std::set< int > ids;
 for( std::vector< int >::const_iterator i = set_ids.begin(); i != set_ids.end(); ++i )
 ids.insert( *i );
 
 // No Sets?
 if( ids.empty() || ( ids.size() == 1 && *ids.begin() == 0 ) ) return MB_SUCCESS; // no sets (all ids are zero)
 
 // Get/create tag handles
 int num_tags;
 Tag tag_handles[2];
 int tag_val;
 const void* tag_values[2] = { &tag_val, NULL };
 
 switch( set_type )
 {
 default:
 return MB_FAILURE;
 case 0:
 case 2: {
 const char* name = set_type ? PARALLEL_PARTITION_TAG_NAME : MATERIAL_SET_TAG_NAME;
 result = mdbImpl->tag_get_handle( name, 1, MB_TYPE_INTEGER, tag_handles[0], MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 num_tags = 1;
 break;
 }
 case 1: {
 result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handles[1],
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 tag_values[1] = NULL;
 tag_handles[0] = globalId;
 num_tags = 2;
 break;
 }
 } // switch
 
 // For each unique set ID...
 for( std::set< int >::iterator i = ids.begin(); i != ids.end(); ++i )
 {
 // Skip "null" set ID
 if( *i == 0 ) continue;
 
 // Get all entities with the current set ID
 Range entities, sets;
 std::vector< int >::const_iterator j = set_ids.begin();
 for( Range::iterator k = elements.begin(); k != elements.end(); ++j, ++k )
 if( *i == *j ) entities.insert( *k );
 
 // Get set by ID
 // Cppcheck warning (false positive): variable tag_val is assigned a value that is never
 // used
 tag_val = *i;
 result = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, tag_handles, tag_values, num_tags, sets );
 if( MB_SUCCESS != result && MB_ENTITY_NOT_FOUND != result ) return result;
 
 // Don't use existing geometry sets (from some other file)
 if( 1 == set_type ) // Geometry
 sets = intersect( sets, geomSets );
 
 // Get set handle
 EntityHandle set;
 // If no sets with ID, create one
 if( sets.empty() )
 {
 result = mdbImpl->create_meshset( MESHSET_SET, set );
 if( MB_SUCCESS != result ) return result;
 
 result = mdbImpl->tag_set_data( tag_handles[0], &set, 1, &*i );
 if( MB_SUCCESS != result ) return result;
 
 if( 1 == set_type )
 { // Geometry
 int dim = CN::Dimension( type );
 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim );
 if( MB_SUCCESS != result ) return result;
 geomSets.insert( set );
 }
 }
 else
 {
 set = *sets.begin();
 if( 1 == set_type )
 { // Geometry
 int dim = CN::Dimension( type );
 // Get dimension of set
 int dim2;
 result = mdbImpl->tag_get_data( tag_handles[1], &set, 1, &dim2 );
 if( MB_SUCCESS != result ) return result;
 // If we're putting geometry of a higher dimension into the
 // set, increase the dimension of the set.
 if( dim > dim2 )
 {
 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim );
 if( MB_SUCCESS != result ) return result;
 }
 }
 }
 
 // Put the mesh entities into the set
 result = mdbImpl->add_entities( set, entities );
 if( MB_SUCCESS != result ) return result;
 } // for (ids)
 
 return MB_SUCCESS;
}
 
//! NOT IMPLEMENTED
//! Reconstruct parent-child relations for geometry sets from
//! mesh connectivity.
ErrorCode ReadGmsh::create_geometric_topology()
{
 if( geomSets.empty() ) return MB_SUCCESS;
 
 // Not implemented yet
 geomSets.clear();
 return MB_SUCCESS;
}
 
} // namespace moab