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