Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadGmsh.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 /**
17  * \class ReadGmsh
18  * \brief Gmsh (http://www.geuz.org/gmsh) file reader
19  *
20  * See: http://geuz.org/gmsh/doc/texinfo/gmsh.html#MSH-ASCII-file-format
21  *
22  * \author Jason Kraftcheck
23  */
24 
25 #include "ReadGmsh.hpp"
26 #include "FileTokenizer.hpp" // for file tokenizer
27 #include "Internals.hpp"
28 #include "moab/Interface.hpp"
29 #include "moab/ReadUtilIface.hpp"
30 #include "moab/Range.hpp"
31 #include "MBTagConventions.hpp"
32 #include "MBParallelConventions.h"
33 #include "moab/CN.hpp"
34 #include "GmshUtil.hpp"
35 
36 #include <cerrno>
37 #include <cstring>
38 #include <map>
39 #include <set>
40 
41 namespace moab
42 {
43 
45 {
46  return new ReadGmsh( iface );
47 }
48 
49 ReadGmsh::ReadGmsh( Interface* impl ) : mdbImpl( impl ), globalId( 0 )
50 {
52 }
53 
55 {
56  if( readMeshIface )
57  {
59  readMeshIface = 0;
60  }
61 }
62 
63 ErrorCode ReadGmsh::read_tag_values( const char* /* file_name */,
64  const char* /* tag_name */,
65  const FileOptions& /* opts */,
66  std::vector< int >& /* tag_values_out */,
67  const SubsetList* /* subset_list */ )
68 {
69  return MB_NOT_IMPLEMENTED;
70 }
71 
72 ErrorCode ReadGmsh::load_file( const char* filename,
73  const EntityHandle*,
74  const FileOptions&,
75  const ReaderIface::SubsetList* subset_list,
76  const Tag* file_id_tag )
77 {
78  int num_material_sets = 0;
79  const int* material_set_list = 0;
80 
81  if( subset_list )
82  {
83  if( subset_list->tag_list_length > 1 && !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
84  {
85  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "GMsh supports subset read only by material ID" );
86  }
87  material_set_list = subset_list->tag_list[0].tag_values;
88  num_material_sets = subset_list->tag_list[0].num_tag_values;
89  }
90 
91  geomSets.clear();
93 
94  // Create set for more convenient check for material set ids
95  std::set< int > blocks;
96  for( const int* mat_set_end = material_set_list + num_material_sets; material_set_list != mat_set_end;
97  ++material_set_list )
98  blocks.insert( *material_set_list );
99 
100  // Map of ID->handle for nodes
101  std::map< long, EntityHandle > node_id_map;
102  int data_size = 8;
103 
104  // Open file and hand off pointer to tokenizer
105  FILE* file_ptr = fopen( filename, "r" );
106  if( !file_ptr )
107  {
108  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) );
109  }
110  FileTokenizer tokens( file_ptr, readMeshIface );
111 
112  // Determine file format version
113  const char* const start_tokens[] = { "$NOD", "$MeshFormat", 0 };
114  int format_version = tokens.match_token( start_tokens );
115  if( !format_version ) return MB_FILE_DOES_NOT_EXIST;
116 
117  // If version 2.0, read additional header info
118  if( 2 == format_version )
119  {
120  double version;
121  if( !tokens.get_doubles( 1, &version ) ) return MB_FILE_WRITE_ERROR;
122 
123  if( version != 2.0 && version != 2.1 && version != 2.2 )
124  {
125  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": unknown format version: " << version );
126  return MB_FILE_DOES_NOT_EXIST;
127  }
128 
129  int file_format;
130  if( !tokens.get_integers( 1, &file_format ) || !tokens.get_integers( 1, &data_size ) ||
131  !tokens.match_token( "$EndMeshFormat" ) )
132  return MB_FILE_WRITE_ERROR;
133  // If physical entities in the gmsh file -> discard this
134  const char* const phys_tokens[] = { "$Nodes", "$PhysicalNames", 0 };
135  int hasPhys = tokens.match_token( phys_tokens );
136 
137  if( hasPhys == 2 )
138  {
139  long num_phys;
140  if( !tokens.get_long_ints( 1, &num_phys ) ) return MB_FILE_WRITE_ERROR;
141  for( long loop_phys = 0; loop_phys < num_phys; loop_phys++ )
142  {
143  long physDim;
144  long physGroupNum;
145  // char const * physName;
146  if( !tokens.get_long_ints( 1, &physDim ) ) return MB_FILE_WRITE_ERROR;
147  if( !tokens.get_long_ints( 1, &physGroupNum ) ) return MB_FILE_WRITE_ERROR;
148  const char* ptc = tokens.get_string();
149  if( !ptc ) return MB_FILE_WRITE_ERROR;
150  // try to get to the end of the line, without reporting errors
151  // really, we need to skip this
152  while( !tokens.get_newline( false ) )
153  ptc = tokens.get_string();
154  }
155  if( !tokens.match_token( "$EndPhysicalNames" ) || !tokens.match_token( "$Nodes" ) )
156  return MB_FILE_WRITE_ERROR;
157  }
158  }
159 
160  // Read number of nodes
161  long num_nodes;
162  if( !tokens.get_long_ints( 1, &num_nodes ) ) return MB_FILE_WRITE_ERROR;
163 
164  // Allocate nodes
165  std::vector< double* > coord_arrays;
166  EntityHandle handle = 0;
167  ErrorCode result = readMeshIface->get_node_coords( 3, num_nodes, MB_START_ID, handle, coord_arrays );
168  if( MB_SUCCESS != result ) return result;
169 
170  // Read nodes
171  double *x = coord_arrays[0], *y = coord_arrays[1], *z = coord_arrays[2];
172  for( long i = 0; i < num_nodes; ++i, ++handle )
173  {
174  long id;
175  if( !tokens.get_long_ints( 1, &id ) || !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) ||
176  !tokens.get_doubles( 1, z++ ) )
177  return MB_FILE_WRITE_ERROR;
178 
179  if( !node_id_map.insert( std::pair< long, EntityHandle >( id, handle ) ).second )
180  {
181  MB_SET_ERR( MB_FILE_WRITE_ERROR, "Duplicate node ID at line " << tokens.line_number() );
182  }
183  }
184 
185  // Create reverse map from handle to id
186  std::vector< int > ids( num_nodes );
187  std::vector< int >::iterator id_iter = ids.begin();
188  std::vector< EntityHandle > handles( num_nodes );
189  std::vector< EntityHandle >::iterator h_iter = handles.begin();
190  for( std::map< long, EntityHandle >::iterator i = node_id_map.begin(); i != node_id_map.end();
191  ++i, ++id_iter, ++h_iter )
192  {
193  *id_iter = i->first;
194  *h_iter = i->second;
195  }
196  // Store IDs in tags
197  result = mdbImpl->tag_set_data( globalId, &handles[0], num_nodes, &ids[0] );
198  if( MB_SUCCESS != result ) return result;
199  if( file_id_tag )
200  {
201  result = mdbImpl->tag_set_data( *file_id_tag, &handles[0], num_nodes, &ids[0] );
202  if( MB_SUCCESS != result ) return result;
203  }
204  ids.clear();
205  handles.clear();
206 
207  // Get tokens signifying end of node data and start of elements
208  if( !tokens.match_token( format_version == 1 ? "$ENDNOD" : "$EndNodes" ) ||
209  !tokens.match_token( format_version == 1 ? "$ELM" : "$Elements" ) )
210  return MB_FILE_WRITE_ERROR;
211 
212  // Get element count
213  long num_elem;
214  if( !tokens.get_long_ints( 1, &num_elem ) ) return MB_FILE_WRITE_ERROR;
215 
216  // Lists of data accumulated for elements
217  std::vector< EntityHandle > connectivity;
218  std::vector< int > mat_set_list, geom_set_list, part_set_list, id_list;
219  // Temporary, per-element data
220  std::vector< int > int_data( 5 ), tag_data( 2 );
221  std::vector< long > tmp_conn;
222  int curr_elem_type = -1;
223  for( long i = 0; i < num_elem; ++i )
224  {
225  // Read element description
226  // File format 1.0
227  if( 1 == format_version )
228  {
229  if( !tokens.get_integers( 5, &int_data[0] ) ) return MB_FILE_WRITE_ERROR;
230  tag_data[0] = int_data[2];
231  tag_data[1] = int_data[3];
232  if( (unsigned)tag_data[1] < GmshUtil::numGmshElemType &&
233  GmshUtil::gmshElemTypes[tag_data[1]].num_nodes != (unsigned)int_data[4] )
234  {
236  "Invalid node count for element type at line " << tokens.line_number() );
237  }
238  }
239  // File format 2.0
240  else
241  {
242  if( !tokens.get_integers( 3, &int_data[0] ) ) return MB_FILE_WRITE_ERROR;
243  tag_data.resize( int_data[2] );
244  if( !tokens.get_integers( tag_data.size(), &tag_data[0] ) ) return MB_FILE_WRITE_ERROR;
245  }
246 
247  // If a list of material sets was specified in the
248  // argument list, skip any elements for which the
249  // material set is not specified or is not in the
250  // passed list.
251  if( !blocks.empty() && ( tag_data.empty() || blocks.find( tag_data[0] ) != blocks.end() ) ) continue;
252 
253  // If the next element is not the same type as the last one,
254  // create a sequence for the block of elements we've read
255  // to this point (all of the same type), and clear accumulated
256  // data.
257  if( int_data[1] != curr_elem_type )
258  {
259  if( !id_list.empty() )
260  { // First iteration
261  result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list,
262  part_set_list, connectivity, file_id_tag );
263  if( MB_SUCCESS != result ) return result;
264  }
265 
266  id_list.clear();
267  mat_set_list.clear();
268  geom_set_list.clear();
269  part_set_list.clear();
270  connectivity.clear();
271  curr_elem_type = int_data[1];
272  if( (unsigned)curr_elem_type >= GmshUtil::numGmshElemType ||
273  GmshUtil::gmshElemTypes[curr_elem_type].mb_type == MBMAXTYPE )
274  {
276  "Unsupported element type " << curr_elem_type << " at line " << tokens.line_number() );
277  }
278  tmp_conn.resize( GmshUtil::gmshElemTypes[curr_elem_type].num_nodes );
279  }
280 
281  // Store data from element description
282  id_list.push_back( int_data[0] );
283  if( tag_data.size() > 3 )
284  part_set_list.push_back( tag_data[3] ); // it must be new format for gmsh, >= 2.5
285  // it could have negative partition ids, for ghost elements
286  else if( tag_data.size() > 2 )
287  part_set_list.push_back( tag_data[2] ); // old format, partition id saved in 3rd tag field
288  else
289  part_set_list.push_back( 0 );
290  geom_set_list.push_back( tag_data.size() > 1 ? tag_data[1] : 0 );
291  mat_set_list.push_back( tag_data.size() > 0 ? tag_data[0] : 0 );
292 
293  // Get element connectivity
294  if( !tokens.get_long_ints( tmp_conn.size(), &tmp_conn[0] ) ) return MB_FILE_WRITE_ERROR;
295 
296  // Convert connectivity from IDs to handles
297  for( unsigned j = 0; j < tmp_conn.size(); ++j )
298  {
299  std::map< long, EntityHandle >::iterator k = node_id_map.find( tmp_conn[j] );
300  if( k == node_id_map.end() )
301  {
302  MB_SET_ERR( MB_FILE_WRITE_ERROR, "Invalid node ID at line " << tokens.line_number() );
303  }
304  connectivity.push_back( k->second );
305  }
306  } // for (num_nodes)
307 
308  // Create entity sequence for last element(s).
309  if( !id_list.empty() )
310  {
311  result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], id_list, mat_set_list, geom_set_list,
312  part_set_list, connectivity, file_id_tag );
313  if( MB_SUCCESS != result ) return result;
314  }
315 
316  // Construct parent-child relations for geometric sets.
317  // Note: At the time this comment was written, the following
318  // function was not implemented.
319  result = create_geometric_topology();
320  geomSets.clear();
321  return result;
322 }
323 
324 //! Create an element sequence
326  const std::vector< int >& elem_ids,
327  const std::vector< int >& matl_ids,
328  const std::vector< int >& geom_ids,
329  const std::vector< int >& prtn_ids,
330  const std::vector< EntityHandle >& connectivity,
331  const Tag* file_id_tag )
332 {
333  ErrorCode result;
334 
335  // Make sure input is consistent
336  const unsigned long num_elem = elem_ids.size();
337  const int node_per_elem = type.num_nodes;
338  if( matl_ids.size() != num_elem || geom_ids.size() != num_elem || prtn_ids.size() != num_elem ||
339  connectivity.size() != num_elem * node_per_elem )
340  return MB_FAILURE;
341 
342  // Create the element sequence
343  // for points, simply gather the connectivities and create the materials
344  if( type.mb_type == MBVERTEX )
345  {
346  Range elements;
347  elements.insert< std::vector< EntityHandle > >( connectivity.begin(), connectivity.end() );
348  result = create_sets( type.mb_type, elements, matl_ids, 0 );
349  if( MB_SUCCESS != result ) return result;
350 
351  return MB_SUCCESS;
352  }
353  EntityHandle handle = 0;
354  EntityHandle* conn_array;
355  result =
356  readMeshIface->get_element_connect( num_elem, node_per_elem, type.mb_type, MB_START_ID, handle, conn_array );
357  if( MB_SUCCESS != result ) return result;
358 
359  // Copy passed element connectivity into entity sequence data.
360  if( type.node_order )
361  {
362  for( unsigned long i = 0; i < num_elem; ++i )
363  for( int j = 0; j < node_per_elem; ++j )
364  conn_array[i * node_per_elem + type.node_order[j]] = connectivity[i * node_per_elem + j];
365  }
366  else
367  {
368  memcpy( conn_array, &connectivity[0], connectivity.size() * sizeof( EntityHandle ) );
369  }
370 
371  // Notify MOAB of the new elements
372  result = readMeshIface->update_adjacencies( handle, num_elem, node_per_elem, conn_array );
373  if( MB_SUCCESS != result ) return result;
374 
375  // Store element IDs
376  Range elements( handle, handle + num_elem - 1 );
377  result = mdbImpl->tag_set_data( globalId, elements, &elem_ids[0] );
378  if( MB_SUCCESS != result ) return result;
379  if( file_id_tag )
380  {
381  result = mdbImpl->tag_set_data( *file_id_tag, elements, &elem_ids[0] );
382  if( MB_SUCCESS != result ) return result;
383  }
384 
385  // Add elements to material sets
386  result = create_sets( type.mb_type, elements, matl_ids, 0 );
387  if( MB_SUCCESS != result ) return result;
388  // Add elements to geometric sets
389  result = create_sets( type.mb_type, elements, geom_ids, 1 );
390  if( MB_SUCCESS != result ) return result;
391  // Add elements to parallel partitions
392  result = create_sets( type.mb_type, elements, prtn_ids, 2 );
393  if( MB_SUCCESS != result ) return result;
394 
395  return MB_SUCCESS;
396 }
397 
398 //! Add elements to sets as dictated by grouping ID in file.
400  const Range& elements,
401  const std::vector< int >& set_ids,
402  int set_type )
403 {
404  ErrorCode result;
405 
406  // Get a unique list of set IDs
407  std::set< int > ids;
408  for( std::vector< int >::const_iterator i = set_ids.begin(); i != set_ids.end(); ++i )
409  ids.insert( *i );
410 
411  // No Sets?
412  if( ids.empty() || ( ids.size() == 1 && *ids.begin() == 0 ) ) return MB_SUCCESS; // no sets (all ids are zero)
413 
414  // Get/create tag handles
415  int num_tags;
416  Tag tag_handles[2];
417  int tag_val;
418  const void* tag_values[2] = { &tag_val, NULL };
419 
420  switch( set_type )
421  {
422  default:
423  return MB_FAILURE;
424  case 0:
425  case 2: {
426  const char* name = set_type ? PARALLEL_PARTITION_TAG_NAME : MATERIAL_SET_TAG_NAME;
427  result = mdbImpl->tag_get_handle( name, 1, MB_TYPE_INTEGER, tag_handles[0], MB_TAG_SPARSE | MB_TAG_CREAT );
428  if( MB_SUCCESS != result ) return result;
429  num_tags = 1;
430  break;
431  }
432  case 1: {
433  result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handles[1],
435  if( MB_SUCCESS != result ) return result;
436  tag_values[1] = NULL;
437  tag_handles[0] = globalId;
438  num_tags = 2;
439  break;
440  }
441  } // switch
442 
443  // For each unique set ID...
444  for( std::set< int >::iterator i = ids.begin(); i != ids.end(); ++i )
445  {
446  // Skip "null" set ID
447  if( *i == 0 ) continue;
448 
449  // Get all entities with the current set ID
450  Range entities, sets;
451  std::vector< int >::const_iterator j = set_ids.begin();
452  for( Range::iterator k = elements.begin(); k != elements.end(); ++j, ++k )
453  if( *i == *j ) entities.insert( *k );
454 
455  // Get set by ID
456  // Cppcheck warning (false positive): variable tag_val is assigned a value that is never
457  // used
458  tag_val = *i;
459  result = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, tag_handles, tag_values, num_tags, sets );
460  if( MB_SUCCESS != result && MB_ENTITY_NOT_FOUND != result ) return result;
461 
462  // Don't use existing geometry sets (from some other file)
463  if( 1 == set_type ) // Geometry
464  sets = intersect( sets, geomSets );
465 
466  // Get set handle
467  EntityHandle set;
468  // If no sets with ID, create one
469  if( sets.empty() )
470  {
471  result = mdbImpl->create_meshset( MESHSET_SET, set );
472  if( MB_SUCCESS != result ) return result;
473 
474  result = mdbImpl->tag_set_data( tag_handles[0], &set, 1, &*i );
475  if( MB_SUCCESS != result ) return result;
476 
477  if( 1 == set_type )
478  { // Geometry
479  int dim = CN::Dimension( type );
480  result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim );
481  if( MB_SUCCESS != result ) return result;
482  geomSets.insert( set );
483  }
484  }
485  else
486  {
487  set = *sets.begin();
488  if( 1 == set_type )
489  { // Geometry
490  int dim = CN::Dimension( type );
491  // Get dimension of set
492  int dim2;
493  result = mdbImpl->tag_get_data( tag_handles[1], &set, 1, &dim2 );
494  if( MB_SUCCESS != result ) return result;
495  // If we're putting geometry of a higher dimension into the
496  // set, increase the dimension of the set.
497  if( dim > dim2 )
498  {
499  result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim );
500  if( MB_SUCCESS != result ) return result;
501  }
502  }
503  }
504 
505  // Put the mesh entities into the set
506  result = mdbImpl->add_entities( set, entities );
507  if( MB_SUCCESS != result ) return result;
508  } // for (ids)
509 
510  return MB_SUCCESS;
511 }
512 
513 //! NOT IMPLEMENTED
514 //! Reconstruct parent-child relations for geometry sets from
515 //! mesh connectivity.
517 {
518  if( geomSets.empty() ) return MB_SUCCESS;
519 
520  // Not implemented yet
521  geomSets.clear();
522  return MB_SUCCESS;
523 }
524 
525 } // namespace moab