Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadTetGen.cpp
Go to the documentation of this file.
1 #include "ReadTetGen.hpp"
2 #include "moab/Interface.hpp"
3 #include "moab/Range.hpp"
4 #include "moab/ReadUtilIface.hpp"
5 #include "moab/FileOptions.hpp"
6 #include "MBTagConventions.hpp"
7 #include <iostream>
8 #include <fstream>
9 #include <sstream>
10 #include <cctype>
11 #include <map>
12 
13 namespace moab
14 {
15 
17 {
18  return new ReadTetGen( moab );
19 }
20 
21 ReadTetGen::ReadTetGen( Interface* moab ) : mbIface( moab ), readTool( 0 )
22 {
23  moab->query_interface( readTool );
24 }
25 
27 {
29 }
30 
31 ErrorCode ReadTetGen::open_file( const std::string& filename,
32  const std::string& basename,
33  const std::string& suffix,
34  const char* exp_suffix,
35  const char* opt_name,
36  const FileOptions& opts,
37  std::ifstream& file_stream,
38  bool file_required )
39 {
40  std::string real_file_name;
41  ErrorCode rval = opts.get_option( opt_name, real_file_name );
42  if( MB_ENTITY_NOT_FOUND == rval || real_file_name.empty() )
43  {
44  if( MB_SUCCESS == rval ) file_required = true;
45  if( suffix == exp_suffix )
46  {
47  real_file_name = filename;
48  }
49  else
50  {
51  real_file_name = basename;
52  real_file_name += ".";
53  real_file_name += exp_suffix;
54  }
55  }
56 
57  if( !real_file_name.empty() ) file_stream.open( real_file_name.c_str(), std::ios::in );
58  if( file_required && !file_stream.is_open() )
59  {
60  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, real_file_name << ": cannot read file" );
61  }
62 
63  return MB_SUCCESS;
64 }
65 
66 ErrorCode ReadTetGen::read_tag_values( const char* /* file_name */,
67  const char* /* tag_name */,
68  const FileOptions& /* opts */,
69  std::vector< int >& /* tag_values_out */,
70  const SubsetList* /* subset_list */ )
71 {
72  return MB_NOT_IMPLEMENTED;
73 }
74 
75 ErrorCode ReadTetGen::load_file( const char* file_name_c,
76  const EntityHandle* /* file_set */,
77  const FileOptions& opts,
78  const ReaderIface::SubsetList* subset_list,
79  const Tag* file_id_tag )
80 {
81  std::ifstream node_file, ele_file, face_file, edge_file;
82  ErrorCode rval;
83 
84  if( subset_list )
85  {
86  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for TetGen" );
87  }
88 
89  std::string suffix, base, filename( file_name_c );
90  size_t dot_idx = filename.find_last_of( '.' );
91  if( dot_idx == std::string::npos )
92  {
93  base = filename;
94  }
95  else
96  {
97  suffix = filename.substr( dot_idx + 1 );
98  for( size_t i = 0; i < suffix.length(); ++i )
99  suffix[i] = (char)tolower( suffix[i] );
100  if( suffix == "node" || suffix == "ele" || suffix == "face" || suffix == "edge" )
101  {
102  base = filename.substr( 0, dot_idx );
103  }
104  else
105  {
106  base = filename;
107  suffix.clear();
108  }
109  }
110 
111  rval = open_file( filename, base, suffix, "node", "NODE_FILE", opts, node_file, true );
112  if( MB_SUCCESS != rval ) return rval;
113  rval = open_file( filename, base, suffix, "ele", "ELE_FILE", opts, ele_file );
114  if( MB_SUCCESS != rval ) return rval;
115  rval = open_file( filename, base, suffix, "face", "FACE_FILE", opts, face_file );
116  if( MB_SUCCESS != rval ) return rval;
117  rval = open_file( filename, base, suffix, "edge", "EDGE_FILE", opts, edge_file );
118  if( MB_SUCCESS != rval ) return rval;
119 
120  std::vector< Tag > attr_tags[4];
121  std::vector< int > attr_idx[4];
122  const char* option_names[4] = { "NODE_ATTR_LIST", "EDGE_ATTR_LIST", "TRI_ATTR_LIST", "TET_ATTR_LIST" };
123  const char* group_names[4] = { 0, "CURVE_ID", "SURFACE_ID", "VOLUME_ID" };
124  for( int i = 0; i < 4; ++i )
125  {
126  std::string opt_str;
127  rval = opts.get_str_option( option_names[i], opt_str );
128  if( MB_SUCCESS != rval ) continue;
129  rval = parse_attr_list( opt_str, attr_tags[i], attr_idx[i], group_names[i] );
130  if( MB_SUCCESS != rval )
131  {
132  MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, option_names[i] << ": invalid option value" );
133  }
134  }
135 
136  Range tets, tris, edges;
137  std::vector< EntityHandle > nodes;
138  rval = read_node_file( node_file, &attr_tags[0][0], &attr_idx[0][0], attr_tags[0].size(), nodes );
139  if( MB_SUCCESS == rval && ele_file.is_open() ) rval = read_elem_file( MBTET, ele_file, nodes, tets );
140  if( MB_SUCCESS == rval && face_file.is_open() ) rval = read_elem_file( MBTRI, face_file, nodes, tris );
141  if( MB_SUCCESS == rval && edge_file.is_open() ) rval = read_elem_file( MBEDGE, edge_file, nodes, edges );
142 
143  if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, &nodes[0], nodes.size() );
144  if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, edges );
145  if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tris );
146  if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tets );
147 
148  return rval;
149 }
150 
151 ErrorCode ReadTetGen::parse_attr_list( const std::string& option_str,
152  std::vector< Tag >& tag_list,
153  std::vector< int >& index_list,
154  const char* group_designator )
155 {
156  std::vector< std::string > name_list;
157  size_t prev_pos = 0;
158  while( prev_pos != std::string::npos )
159  {
160  size_t pos = option_str.find_first_of( ',', prev_pos );
161  name_list.push_back( option_str.substr( prev_pos, pos ) );
162  prev_pos = pos + 1;
163  }
164 
165  index_list.resize( name_list.size() );
166  std::map< std::string, int > name_count;
167  for( size_t i = 0; i < name_list.size(); ++i )
168  index_list[i] = name_count[name_list[i]]++;
169 
170  for( size_t i = 0; i < name_list.size(); ++i )
171  {
172  if( group_designator && name_list[i] == group_designator )
173  {
174  tag_list[i] = 0;
175  index_list[i] = -1;
176  }
177  else if( name_list.empty() )
178  {
179  tag_list[i] = 0;
180  index_list[i] = 0;
181  }
182  else
183  {
184  ErrorCode rval = mbIface->tag_get_handle( name_list[i].c_str(), name_count[name_list[i]], MB_TYPE_DOUBLE,
185  tag_list[i], MB_TAG_DENSE | MB_TAG_CREAT );
186  if( MB_SUCCESS != rval ) return rval;
187  }
188  }
189 
190  return MB_SUCCESS;
191 }
192 
193 ErrorCode ReadTetGen::read_line( std::istream& file, std::string& line, int& lineno )
194 {
195  // Loop until we find a non-empty line
196  do
197  {
198  // Read a line
199  line.clear();
200  if( !getline( file, line ) ) return MB_FILE_WRITE_ERROR;
201  ++lineno;
202  // Strip comments from line
203  size_t pos = line.find_first_of( '#' );
204  if( pos != std::string::npos ) line = line.substr( 0, pos );
205  // Strip leading whitespace from line
206  for( pos = 0; pos < line.length() && isspace( line[pos] ); ++pos )
207  ;
208  if( pos == line.length() )
209  line.clear();
210  else if( pos != 0 )
211  line = line.substr( pos );
212  } while( line.empty() );
213 
214  return MB_SUCCESS;
215 }
216 
217 ErrorCode ReadTetGen::read_line( std::istream& file, double* values_out, int num_values, int& lineno )
218 {
219  // Get a line of text
220  std::string line;
221  ErrorCode rval = read_line( file, line, lineno );
222  if( MB_SUCCESS != rval ) return rval;
223 
224  // Tokenize line as doubles
225  std::stringstream str( line );
226  for( int i = 0; i < num_values; ++i )
227  {
228  double v;
229  if( !( str >> v ) )
230  {
231  MB_SET_ERR( MB_FAILURE, "Error reading node data at line " << lineno );
232  }
233  values_out[i] = v;
234  }
235 
236  // Check that we're at the end of the line
237  int junk;
238  if( ( str >> junk ) || !str.eof() )
239  {
240  MB_SET_ERR( MB_FAILURE, "Unexpected trailing data for line " << lineno << " of node data" );
241  }
242 
243  return MB_SUCCESS;
244 }
245 
247  const Tag* attr_tag_list,
248  const int* attr_tag_index,
249  int attr_tag_list_len,
250  std::vector< EntityHandle >& nodes )
251 {
252  int lineno = 0;
253  ErrorCode rval;
254 
255  double header_vals[4];
256  rval = read_line( file, header_vals, 4, lineno );
257  if( MB_SUCCESS != rval ) return rval;
258 
259  const int num_vtx = (int)header_vals[0];
260  const int dim = (int)header_vals[1];
261  const int num_attr = (int)header_vals[2];
262  const int bdry_flag = (int)header_vals[3];
263  if( num_vtx < 1 || dim < 2 || dim > 3 || num_attr < 0 || bdry_flag < 0 || bdry_flag > 1 )
264  {
265  MB_SET_ERR( MB_FAILURE, "Invalid header line for node data" );
266  }
267  if( attr_tag_list_len > num_attr ) attr_tag_list_len = num_attr;
268 
269  // Allocate space for tag data
270  std::map< Tag, int > tag_size;
271  std::map< Tag, std::vector< double > > tag_data;
272  for( int i = 0; i < attr_tag_list_len; ++i )
273  {
274  if( !attr_tag_list[i] || attr_tag_index[i] < 0 ) continue;
275  std::vector< double >& data = tag_data[attr_tag_list[i]];
276  // Increase tag size by one value per vertex for each time
277  // we encounter it in the list.
278  data.resize( data.size() + num_vtx );
279  ++tag_size[attr_tag_list[i]];
280  }
281  std::vector< double* > attr_data( attr_tag_list_len );
282  std::vector< int > attr_size( attr_tag_list_len );
283  for( int i = 0; i < attr_tag_list_len; ++i )
284  {
285  if( !attr_tag_list[i] || attr_tag_index[i] < 0 )
286  {
287  attr_data[i] = 0;
288  attr_size[i] = 0;
289  }
290  else
291  {
292  attr_data[i] = &( tag_data[attr_tag_list[i]] )[0];
293  attr_size[i] = tag_size[attr_tag_list[i]];
294  }
295  }
296 
297  // Allocate vertices
298  std::vector< double* > coords;
299  EntityHandle start_handle;
300  rval = readTool->get_node_coords( dim, num_vtx, 1, start_handle, coords );
301  if( MB_SUCCESS != rval ) return rval;
302 
303  // Read data line for each node
304  nodes.reserve( num_vtx );
305  std::vector< double > data( 1 + dim + num_attr + bdry_flag );
306  std::vector< int > ids( num_vtx );
307  for( int i = 0; i < num_vtx; ++i )
308  {
309  rval = read_line( file, &data[0], data.size(), lineno );
310  if( MB_SUCCESS != rval ) return rval;
311 
312  // Get ID
313  ids[i] = (int)data[0];
314  if( ids[i] >= (int)nodes.size() ) nodes.resize( ids[i] + 1 );
315  nodes[ids[i]] = start_handle + i;
316 
317  // Get coordinates
318  // Cppcheck warning (false positive): variable coords is assigned a value that is never used
319  for( int j = 0; j < dim; ++j )
320  coords[j][i] = data[j + 1];
321 
322  // Get attribute data
323  for( int j = 0; j < attr_tag_list_len; ++j )
324  if( attr_data[j] ) attr_data[j][i * attr_size[j] + attr_tag_index[j]] = data[j + 1 + dim];
325 
326  // Discard boundary bit
327  }
328 
329  // Store tag data
330  Range node_range;
331  node_range.insert( start_handle, start_handle + num_vtx - 1 );
332  for( std::map< Tag, std::vector< double > >::iterator i = tag_data.begin(); i != tag_data.end(); ++i )
333  {
334  rval = mbIface->tag_set_data( i->first, node_range, &i->second[0] );
335  if( MB_SUCCESS != rval ) return rval;
336  }
337  Tag idtag = mbIface->globalId_tag();
338 
339  rval = mbIface->tag_set_data( idtag, node_range, &ids[0] );
340  if( MB_SUCCESS != rval ) return rval;
341 
342  return MB_SUCCESS;
343 }
344 
346  std::istream& file,
347  const std::vector< EntityHandle >& nodes,
348  Range& elems )
349 {
350  int lineno = 0;
351  ErrorCode rval;
352 
353  int node_per_elem, have_group_id, dim;
354  double header_vals[3];
355  switch( type )
356  {
357  case MBTET:
358  rval = read_line( file, header_vals, 3, lineno );
359  node_per_elem = (int)header_vals[1];
360  have_group_id = (int)header_vals[2];
361  dim = 3;
362  break;
363  case MBTRI:
364  rval = read_line( file, header_vals, 2, lineno );
365  node_per_elem = 3;
366  have_group_id = (int)header_vals[1];
367  dim = 2;
368  break;
369  case MBEDGE:
370  rval = read_line( file, header_vals, 1, lineno );
371  node_per_elem = 2;
372  have_group_id = 0;
373  dim = 1;
374  break;
375  default:
376  rval = MB_FAILURE;
377  break;
378  }
379  if( MB_SUCCESS != rval ) return rval;
380  const int num_elem = (int)header_vals[0];
381  if( num_elem < 1 || node_per_elem < 2 || have_group_id < 0 || have_group_id > 1 )
382  {
383  MB_SET_ERR( MB_FAILURE, "Invalid header line for element data" );
384  }
385 
386  // Create group map
387  std::map< double, EntityHandle > groups;
388  Tag dim_tag, id_tag;
389  id_tag = mbIface->globalId_tag();
390 
391  const int negone = -1;
393  &negone );
394  if( MB_SUCCESS != rval ) return rval;
395 
396  // Allocate elements
397  EntityHandle start_handle, *conn_array;
398  rval = readTool->get_element_connect( num_elem, node_per_elem, type, 1, start_handle, conn_array );
399  if( MB_SUCCESS != rval ) return rval;
400  elems.insert( start_handle, start_handle + num_elem - 1 );
401 
402  // Read data line for each node
403  std::vector< double > data( 1 + node_per_elem + have_group_id );
404  std::vector< int > ids( num_elem );
405  for( int i = 0; i < num_elem; ++i )
406  {
407  rval = read_line( file, &data[0], data.size(), lineno );
408  if( MB_SUCCESS != rval ) return rval;
409 
410  // Get ID
411  ids[i] = (int)data[0];
412 
413  // Get connectivity
414  for( int j = 0; j < node_per_elem; ++j )
415  conn_array[node_per_elem * i + j] = nodes[(int)data[j + 1]];
416 
417  // Grouping
418  if( have_group_id && 0.0 != data[node_per_elem + 1] )
419  {
420  double id = data[node_per_elem + 1];
421  EntityHandle grp = groups[id];
422  if( 0 == grp )
423  {
424  rval = mbIface->create_meshset( MESHSET_SET, grp );
425  if( MB_SUCCESS != rval ) return rval;
426  elems.insert( grp );
427  rval = mbIface->tag_set_data( dim_tag, &grp, 1, &dim );
428  if( MB_SUCCESS != rval ) return rval;
429  int iid = (int)id;
430  rval = mbIface->tag_set_data( id_tag, &grp, 1, &iid );
431  if( MB_SUCCESS != rval ) return rval;
432  groups[id] = grp;
433  }
434  EntityHandle handle = start_handle + i;
435  rval = mbIface->add_entities( grp, &handle, 1 );
436  if( MB_SUCCESS != rval ) return rval;
437  }
438  }
439 
440  // Store id data
441  Range elems2;
442  elems2.insert( start_handle, start_handle + num_elem - 1 );
443  rval = mbIface->tag_set_data( id_tag, elems2, &ids[0] );
444  if( MB_SUCCESS != rval ) return rval;
445 
446  return MB_SUCCESS;
447 }
448 
449 } // namespace moab