 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
44 ReaderIface* ReadGmsh::factory( Interface* iface )
45 {
46 return new ReadGmsh( iface );
47 }
48
49 ReadGmsh::ReadGmsh( Interface* impl ) : mdbImpl( impl ), globalId( 0 )
50 {
51 mdbImpl->query_interface( readMeshIface );
52 }
53
54 ReadGmsh::~ReadGmsh()
55 {
56 if( readMeshIface )
57 {
58 mdbImpl->release_interface( readMeshIface );
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();
92 globalId = mdbImpl->globalId_tag();
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 {
235 MB_SET_ERR( MB_FILE_WRITE_ERROR,
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 {
275 MB_SET_ERR( MB_FILE_WRITE_ERROR,
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
325 ErrorCode ReadGmsh::create_elements( const GmshElemType& type,
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.
399 ErrorCode ReadGmsh::create_sets( EntityType type,
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],
434 MB_TAG_SPARSE | MB_TAG_CREAT );
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.
516 ErrorCode ReadGmsh::create_geometric_topology()
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
1.9.1.