Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadVtk.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 ReadVtk
18  * \brief VTK reader from Mesquite
19  * \author Jason Kraftcheck
20  */
21 
22 #include <cassert>
23 #include <cstdlib>
24 #include <cstring>
25 
26 #include "ReadVtk.hpp"
27 #include "moab/Range.hpp"
28 #include "Internals.hpp"
29 #include "moab/Interface.hpp"
30 #include "moab/ReadUtilIface.hpp"
31 #include "moab/FileOptions.hpp"
32 #include "FileTokenizer.hpp"
33 #include "moab/VtkUtil.hpp"
34 #include "MBTagConventions.hpp"
35 
36 // #define MB_VTK_MATERIAL_SETS
37 #ifdef MB_VTK_MATERIAL_SETS
38 #include <map>
39 #endif
40 
41 namespace moab
42 {
43 
44 #ifdef MB_VTK_MATERIAL_SETS
45 
46 class Hash
47 {
48  public:
49  unsigned long value;
50 
51  Hash()
52  {
53  this->value = 5381L;
54  }
55 
56  Hash( const unsigned char* bytes, size_t len )
57  { // djb2, a hash by dan bernstein presented on comp.lang.c for hashing strings.
58  this->value = 5381L;
59  for( ; len; --len, ++bytes )
60  this->value = this->value * 33 + ( *bytes );
61  }
62 
63  Hash( bool duh )
64  {
65  this->value = duh; // Hashing a single bit with a long is stupid but easy.
66  }
67 
68  Hash( const Hash& src )
69  {
70  this->value = src.value;
71  }
72 
73  Hash& operator=( const Hash& src )
74  {
75  this->value = src.value;
76  return *this;
77  }
78 
79  bool operator<( const Hash& other ) const
80  {
81  return this->value < other.value;
82  }
83 };
84 
85 // Pass this class opaque data + a handle and it adds the handle to a set
86 // whose opaque data has the same hash. If no set exists for the hash,
87 // it creates one. Each set is tagged with the opaque data.
88 // When entities with different opaque data have the same hash, they
89 // will be placed into the same set.
90 // There will be no collisions when the opaque data is shorter than an
91 // unsigned long, and this is really the only case we need to support.
92 // The rest is bonus. Hash does quite well with strings, even those
93 // with identical prefixes.
94 class Modulator : public std::map< Hash, EntityHandle >
95 {
96  public:
97  Modulator( Interface* iface ) : mesh( iface ) {}
98 
99  ErrorCode initialize( std::string tag_name, DataType mb_type, size_t sz, size_t per_elem )
100  {
101  std::vector< unsigned char > default_val;
102  default_val.resize( sz * per_elem );
103  ErrorCode rval = this->mesh->tag_get_handle( tag_name.c_str(), per_elem, mb_type, this->tag,
104  MB_TAG_SPARSE | MB_TAG_BYTES | MB_TAG_CREAT, &default_val[0] );
105  return rval;
106  }
107 
108  void add_entity( EntityHandle ent, const unsigned char* bytes, size_t len )
109  {
110  Hash h( bytes, len );
111  EntityHandle mset = this->congruence_class( h, bytes );
112  ErrorCode rval;
113  rval = this->mesh->add_entities( mset, &ent, 1 );MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
114  }
115 
116  void add_entities( Range& range, const unsigned char* bytes, size_t bytes_per_ent )
117  {
118  for( Range::iterator it = range.begin(); it != range.end(); ++it, bytes += bytes_per_ent )
119  {
120  Hash h( bytes, bytes_per_ent );
121  EntityHandle mset = this->congruence_class( h, bytes );
122  ErrorCode rval;
123  rval = this->mesh->add_entities( mset, &*it, 1 );MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
124  }
125  }
126 
127  EntityHandle congruence_class( const Hash& h, const void* tag_data )
128  {
129  std::map< Hash, EntityHandle >::iterator it = this->find( h );
130  if( it == this->end() )
131  {
132  EntityHandle mset;
133  Range preexist;
134  ErrorCode rval;
135  rval = this->mesh->get_entities_by_type_and_tag( 0, MBENTITYSET, &this->tag, &tag_data, 1, preexist );MB_CHK_SET_ERR_RET_VAL( rval, "Failed to get entities by type and tag", (EntityHandle)0 );
136  if( preexist.size() )
137  {
138  mset = *preexist.begin();
139  }
140  else
141  {
142  rval = this->mesh->create_meshset( MESHSET_SET, mset );MB_CHK_SET_ERR_RET_VAL( rval, "Failed to create mesh set", (EntityHandle)0 );
143  rval = this->mesh->tag_set_data( this->tag, &mset, 1, tag_data );MB_CHK_SET_ERR_RET_VAL( rval, "Failed to set tag data", (EntityHandle)0 );
144  }
145  ( *this )[h] = mset;
146  return mset;
147  }
148  return it->second;
149  }
150 
151  Interface* mesh;
152  Tag tag;
153 };
154 #endif // MB_VTK_MATERIAL_SETS
155 
156 ReaderIface* ReadVtk::factory( Interface* iface )
157 {
158  return new ReadVtk( iface );
159 }
160 
161 ReadVtk::ReadVtk( Interface* impl ) : mdbImpl( impl ), mPartitionTagName( MATERIAL_SET_TAG_NAME )
162 {
163  mdbImpl->query_interface( readMeshIface );
164 }
165 
166 ReadVtk::~ReadVtk()
167 {
168  if( readMeshIface )
169  {
170  mdbImpl->release_interface( readMeshIface );
171  readMeshIface = 0;
172  }
173 }
174 
175 const char* const vtk_type_names[] = {
176  "bit", "char", "unsigned_char", "short", "unsigned_short", "int", "unsigned_int",
177  "long", "unsigned_long", "float", "double", "vtkIdType", 0 };
178 
179 ErrorCode ReadVtk::read_tag_values( const char* /* file_name */,
180  const char* /* tag_name */,
181  const FileOptions& /* opts */,
182  std::vector< int >& /* tag_values_out */,
183  const SubsetList* /* subset_list */ )
184 {
185  return MB_NOT_IMPLEMENTED;
186 }
187 
188 ErrorCode ReadVtk::load_file( const char* filename,
189  const EntityHandle* /* file_set */,
190  const FileOptions& opts,
191  const ReaderIface::SubsetList* subset_list,
192  const Tag* file_id_tag )
193 {
194  ErrorCode result;
195 
196  int major, minor;
197  char vendor_string[257];
198  std::vector< Range > element_list;
199  Range vertices;
200 
201  if( subset_list )
202  {
203  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for VTK" );
204  }
205 
206  // Does the caller want a field to be used for partitioning the entities?
207  // If not, we'll assume any scalar integer field named MATERIAL_SET specifies partitions.
208  std::string partition_tag_name;
209  result = opts.get_option( "PARTITION", partition_tag_name );
210  if( result == MB_SUCCESS ) mPartitionTagName = partition_tag_name;
211 
212  FILE* file = fopen( filename, "r" );
213  if( !file ) return MB_FILE_DOES_NOT_EXIST;
214 
215  // Read file header
216 
217  if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
218  {
219  fclose( file );
220  return MB_FAILURE;
221  }
222 
223  if( !strchr( vendor_string, '\n' ) || 2 != sscanf( vendor_string, "# vtk DataFile Version %d.%d", &major, &minor ) )
224  {
225  fclose( file );
226  return MB_FAILURE;
227  }
228 
229  if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
230  {
231  fclose( file );
232  return MB_FAILURE;
233  }
234 
235  // VTK spec says this should not exceed 256 chars.
236  if( !strchr( vendor_string, '\n' ) )
237  {
238  fclose( file );
239  MB_SET_ERR( MB_FAILURE, "Vendor string (line 2) exceeds 256 characters" );
240  }
241 
242  // Check file type
243 
244  FileTokenizer tokens( file, readMeshIface );
245  const char* const file_type_names[] = { "ASCII", "BINARY", 0 };
246  int filetype = tokens.match_token( file_type_names );
247  switch( filetype )
248  {
249  case 2: // BINARY
250  MB_SET_ERR( MB_FAILURE, "Cannot read BINARY VTK files" );
251  default: // ERROR
252  return MB_FAILURE;
253  case 1: // ASCII
254  break;
255  }
256 
257  // Read the mesh
258  if( !tokens.match_token( "DATASET" ) ) return MB_FAILURE;
259  result = vtk_read_dataset( tokens, vertices, element_list );
260  if( MB_SUCCESS != result ) return result;
261 
262  if( file_id_tag )
263  {
264  result = store_file_ids( *file_id_tag, vertices, element_list );
265  if( MB_SUCCESS != result ) return result;
266  }
267 
268  // Count the number of elements read
269  long elem_count = 0;
270  for( std::vector< Range >::iterator it = element_list.begin(); it != element_list.end(); ++it )
271  elem_count += it->size();
272 
273  // Read attribute data until end of file.
274  const char* const block_type_names[] = { "POINT_DATA", "CELL_DATA", 0 };
275  std::vector< Range > vertex_list( 1 );
276  vertex_list[0] = vertices;
277  int blocktype = 0;
278  while( !tokens.eof() )
279  {
280  // Get POINT_DATA or CELL_DATA
281  int new_block_type = tokens.match_token( block_type_names, false );
282  if( tokens.eof() ) break;
283 
284  if( !new_block_type )
285  {
286  // If next token was neither POINT_DATA nor CELL_DATA,
287  // then there's another attribute under the current one.
288  if( blocktype )
289  tokens.unget_token();
290  else
291  break;
292  }
293  else
294  {
295  blocktype = new_block_type;
296  long count;
297  if( !tokens.get_long_ints( 1, &count ) ) return MB_FAILURE;
298 
299  if( blocktype == 1 && (unsigned long)count != vertices.size() )
300  {
301  MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of vertices at line " << tokens.line_number() );
302  }
303  else if( blocktype == 2 && count != elem_count )
304  {
305  MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of elements at line " << tokens.line_number() );
306  }
307  }
308 
309  if( blocktype == 1 )
310  result = vtk_read_attrib_data( tokens, vertex_list );
311  else
312  result = vtk_read_attrib_data( tokens, element_list );
313 
314  if( MB_SUCCESS != result ) return result;
315  }
316 
317  return MB_SUCCESS;
318 }
319 
320 ErrorCode ReadVtk::allocate_vertices( long num_verts,
321  EntityHandle& start_handle_out,
322  double*& x_coord_array_out,
323  double*& y_coord_array_out,
324  double*& z_coord_array_out )
325 {
326  ErrorCode result;
327 
328  // Create vertices
329  std::vector< double* > arrays;
330  start_handle_out = 0;
331  result = readMeshIface->get_node_coords( 3, num_verts, MB_START_ID, start_handle_out, arrays );
332  if( MB_SUCCESS != result ) return result;
333 
334  x_coord_array_out = arrays[0];
335  y_coord_array_out = arrays[1];
336  z_coord_array_out = arrays[2];
337 
338  return MB_SUCCESS;
339 }
340 
341 ErrorCode ReadVtk::read_vertices( FileTokenizer& tokens, long num_verts, EntityHandle& start_handle_out )
342 {
343  ErrorCode result;
344  double *x, *y, *z;
345 
346  result = allocate_vertices( num_verts, start_handle_out, x, y, z );
347  if( MB_SUCCESS != result ) return result;
348 
349  // Read vertex coordinates
350  for( long vtx = 0; vtx < num_verts; ++vtx )
351  {
352  if( !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) || !tokens.get_doubles( 1, z++ ) )
353  return MB_FAILURE;
354  }
355 
356  return MB_SUCCESS;
357 }
358 
359 ErrorCode ReadVtk::allocate_elements( long num_elements,
360  int vert_per_element,
361  EntityType type,
362  EntityHandle& start_handle_out,
363  EntityHandle*& conn_array_out,
364  std::vector< Range >& append_to_this )
365 {
366  ErrorCode result;
367 
368  start_handle_out = 0;
369  result = readMeshIface->get_element_connect( num_elements, vert_per_element, type, MB_START_ID, start_handle_out,
370  conn_array_out );
371  if( MB_SUCCESS != result ) return result;
372 
373  Range range( start_handle_out, start_handle_out + num_elements - 1 );
374  append_to_this.push_back( range );
375  return MB_SUCCESS;
376 }
377 
378 ErrorCode ReadVtk::vtk_read_dataset( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& element_list )
379 {
380  const char* const data_type_names[] = {
381  "STRUCTURED_POINTS", "STRUCTURED_GRID", "UNSTRUCTURED_GRID", "POLYDATA", "RECTILINEAR_GRID", "FIELD", 0 };
382  int datatype = tokens.match_token( data_type_names );
383  switch( datatype )
384  {
385  case 1:
386  return vtk_read_structured_points( tokens, vertex_list, element_list );
387  case 2:
388  return vtk_read_structured_grid( tokens, vertex_list, element_list );
389  case 3:
390  return vtk_read_unstructured_grid( tokens, vertex_list, element_list );
391  case 4:
392  return vtk_read_polydata( tokens, vertex_list, element_list );
393  case 5:
394  return vtk_read_rectilinear_grid( tokens, vertex_list, element_list );
395  case 6:
396  return vtk_read_field( tokens );
397  default:
398  return MB_FAILURE;
399  }
400 }
401 
402 ErrorCode ReadVtk::vtk_read_structured_points( FileTokenizer& tokens,
403  Range& vertex_list,
404  std::vector< Range >& elem_list )
405 {
406  long i, j, k;
407  long dims[3];
408  double origin[3], space[3];
409  ErrorCode result;
410 
411  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
412  return MB_FAILURE;
413 
414  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
415  {
416  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
417  }
418 
419  if( !tokens.match_token( "ORIGIN" ) || !tokens.get_doubles( 3, origin ) || !tokens.get_newline() )
420  return MB_FAILURE;
421 
422  const char* const spacing_names[] = { "SPACING", "ASPECT_RATIO", 0 };
423  if( !tokens.match_token( spacing_names ) || !tokens.get_doubles( 3, space ) || !tokens.get_newline() )
424  return MB_FAILURE;
425 
426  // Create vertices
427  double *x, *y, *z;
428  EntityHandle start_handle = 0;
429  long num_verts = dims[0] * dims[1] * dims[2];
430  result = allocate_vertices( num_verts, start_handle, x, y, z );
431  if( MB_SUCCESS != result ) return result;
432  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
433 
434  for( k = 0; k < dims[2]; ++k )
435  for( j = 0; j < dims[1]; ++j )
436  for( i = 0; i < dims[0]; ++i )
437  {
438  *x = origin[0] + i * space[0];
439  ++x;
440  *y = origin[1] + j * space[1];
441  ++y;
442  *z = origin[2] + k * space[2];
443  ++z;
444  }
445 
446  return vtk_create_structured_elems( dims, start_handle, elem_list );
447 }
448 
449 ErrorCode ReadVtk::vtk_read_structured_grid( FileTokenizer& tokens,
450  Range& vertex_list,
451  std::vector< Range >& elem_list )
452 {
453  long num_verts, dims[3];
454  ErrorCode result;
455 
456  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
457  return MB_FAILURE;
458 
459  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
460  {
461  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
462  }
463 
464  if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
465  !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
466  return MB_FAILURE;
467 
468  if( num_verts != ( dims[0] * dims[1] * dims[2] ) )
469  {
470  MB_SET_ERR( MB_FAILURE, "Point count not consistent with dimensions at line " << tokens.line_number() );
471  }
472 
473  // Create and read vertices
474  EntityHandle start_handle = 0;
475  result = read_vertices( tokens, num_verts, start_handle );
476  if( MB_SUCCESS != result ) return result;
477  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
478 
479  return vtk_create_structured_elems( dims, start_handle, elem_list );
480 }
481 
482 ErrorCode ReadVtk::vtk_read_rectilinear_grid( FileTokenizer& tokens,
483  Range& vertex_list,
484  std::vector< Range >& elem_list )
485 {
486  int i, j, k;
487  long dims[3];
488  const char* labels[] = { "X_COORDINATES", "Y_COORDINATES", "Z_COORDINATES" };
489  std::vector< double > coords[3];
490  ErrorCode result;
491 
492  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
493  return MB_FAILURE;
494 
495  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
496  {
497  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
498  }
499 
500  for( i = 0; i < 3; i++ )
501  {
502  long count;
503  if( !tokens.match_token( labels[i] ) || !tokens.get_long_ints( 1, &count ) ||
504  !tokens.match_token( vtk_type_names ) )
505  return MB_FAILURE;
506 
507  if( count != dims[i] )
508  {
509  MB_SET_ERR( MB_FAILURE, "Coordinate count inconsistent with dimensions at line " << tokens.line_number() );
510  }
511 
512  coords[i].resize( count );
513  if( !tokens.get_doubles( count, &coords[i][0] ) ) return MB_FAILURE;
514  }
515 
516  // Create vertices
517  double *x, *y, *z;
518  EntityHandle start_handle = 0;
519  long num_verts = dims[0] * dims[1] * dims[2];
520  result = allocate_vertices( num_verts, start_handle, x, y, z );
521  if( MB_SUCCESS != result ) return result;
522  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
523 
524  // Calculate vertex coordinates
525  for( k = 0; k < dims[2]; ++k )
526  for( j = 0; j < dims[1]; ++j )
527  for( i = 0; i < dims[0]; ++i )
528  {
529  *x = coords[0][i];
530  ++x;
531  *y = coords[1][j];
532  ++y;
533  *z = coords[2][k];
534  ++z;
535  }
536 
537  return vtk_create_structured_elems( dims, start_handle, elem_list );
538 }
539 
540 ErrorCode ReadVtk::vtk_read_polydata( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& elem_list )
541 {
542  ErrorCode result;
543  long num_verts;
544  const char* const poly_data_names[] = { "VERTICES", "LINES", "POLYGONS", "TRIANGLE_STRIPS", 0 };
545 
546  if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
547  !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
548  return MB_FAILURE;
549 
550  if( num_verts < 1 )
551  {
552  MB_SET_ERR( MB_FAILURE, "Invalid point count at line " << tokens.line_number() );
553  }
554 
555  // Create vertices and read coordinates
556  EntityHandle start_handle = 0;
557  result = read_vertices( tokens, num_verts, start_handle );
558  if( MB_SUCCESS != result ) return result;
559  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
560 
561  int poly_type = tokens.match_token( poly_data_names );
562  switch( poly_type )
563  {
564  case 0:
565  result = MB_FAILURE;
566  break;
567  case 1:
568  MB_SET_ERR( MB_FAILURE, "Vertex element type at line " << tokens.line_number() );
569  break;
570  case 2:
571  MB_SET_ERR( MB_FAILURE, "Unsupported type: polylines at line " << tokens.line_number() );
572  result = MB_FAILURE;
573  break;
574  case 3:
575  result = vtk_read_polygons( tokens, start_handle, elem_list );
576  break;
577  case 4:
578  MB_SET_ERR( MB_FAILURE, "Unsupported type: triangle strips at line " << tokens.line_number() );
579  result = MB_FAILURE;
580  break;
581  }
582 
583  return result;
584 }
585 
586 ErrorCode ReadVtk::vtk_read_polygons( FileTokenizer& tokens, EntityHandle first_vtx, std::vector< Range >& elem_list )
587 {
588  ErrorCode result;
589  long size[2];
590 
591  if( !tokens.get_long_ints( 2, size ) || !tokens.get_newline() ) return MB_FAILURE;
592 
593  const Range empty;
594  std::vector< EntityHandle > conn_hdl;
595  std::vector< long > conn_idx;
596  EntityHandle first = 0, prev = 0, handle;
597  for( int i = 0; i < size[0]; ++i )
598  {
599  long count;
600  if( !tokens.get_long_ints( 1, &count ) ) return MB_FAILURE;
601  conn_idx.resize( count );
602  conn_hdl.resize( count );
603  if( !tokens.get_long_ints( count, &conn_idx[0] ) ) return MB_FAILURE;
604 
605  for( long j = 0; j < count; ++j )
606  conn_hdl[j] = first_vtx + conn_idx[j];
607 
608  result = mdbImpl->create_element( MBPOLYGON, &conn_hdl[0], count, handle );
609  if( MB_SUCCESS != result ) return result;
610 
611  if( prev + 1 != handle )
612  {
613  if( first )
614  { // True except for first iteration (first == 0)
615  if( elem_list.empty() || first < elem_list.back().front() ) // Only need new range if order would get
616  // mixed up, or we just began inserting
617  elem_list.push_back( empty );
618  elem_list.back().insert( first, prev );
619  }
620  first = handle;
621  }
622  prev = handle;
623  }
624  if( first )
625  { // True unless no elements (size[0] == 0)
626  if( elem_list.empty() || first < elem_list.back().front() ) // Only need new range if order would get mixed
627  // up, or we just began inserting
628  elem_list.push_back( empty );
629  elem_list.back().insert( first, prev );
630  }
631 
632  return MB_SUCCESS;
633 }
634 
635 ErrorCode ReadVtk::vtk_read_unstructured_grid( FileTokenizer& tokens,
636  Range& vertex_list,
637  std::vector< Range >& elem_list )
638 {
639  ErrorCode result;
640  long i, num_verts, num_elems[2];
641  EntityHandle tmp_conn_list[27];
642 
643  // Poorly formatted VTK legacy format document seems to
644  // lead many to think that a FIELD block can occur within
645  // an UNSTRUCTURED_GRID dataset rather than as its own data
646  // set. So allow for field data between other blocks of
647  // data.
648 
649  const char* pts_str[] = { "FIELD", "POINTS", 0 };
650  while( 1 == ( i = tokens.match_token( pts_str ) ) )
651  {
652  result = vtk_read_field( tokens );
653  if( MB_SUCCESS != result ) return result;
654  }
655  if( i != 2 ) return MB_FAILURE;
656 
657  if( !tokens.get_long_ints( 1, &num_verts ) || !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
658  return MB_FAILURE;
659 
660  if( num_verts < 1 )
661  {
662  MB_SET_ERR( MB_FAILURE, "Invalid point count at line " << tokens.line_number() );
663  }
664 
665  // Create vertices and read coordinates
666  EntityHandle first_vertex = 0;
667  result = read_vertices( tokens, num_verts, first_vertex );
668  if( MB_SUCCESS != result ) return result;
669  vertex_list.insert( first_vertex, first_vertex + num_verts - 1 );
670 
671  const char* cell_str[] = { "FIELD", "CELLS", 0 };
672  while( 1 == ( i = tokens.match_token( cell_str ) ) )
673  {
674  result = vtk_read_field( tokens );
675  if( MB_SUCCESS != result ) return result;
676  }
677  if( i != 2 ) return MB_FAILURE;
678 
679  if( !tokens.get_long_ints( 2, num_elems ) || !tokens.get_newline() ) return MB_FAILURE;
680 
681  // Read element connectivity for all elements
682  std::vector< long > connectivity( num_elems[1] );
683  if( !tokens.get_long_ints( num_elems[1], &connectivity[0] ) ) return MB_FAILURE;
684 
685  if( !tokens.match_token( "CELL_TYPES" ) || !tokens.get_long_ints( 1, &num_elems[1] ) || !tokens.get_newline() )
686  return MB_FAILURE;
687 
688  // Read element types
689  std::vector< long > types( num_elems[0] );
690  if( !tokens.get_long_ints( num_elems[0], &types[0] ) ) return MB_FAILURE;
691 
692  // Create elements in blocks of the same type
693  // It is important to preserve the order in
694  // which the elements were read for later reading
695  // attribute data.
696  long id = 0;
697  std::vector< long >::iterator conn_iter = connectivity.begin();
698  while( id < num_elems[0] )
699  {
700  unsigned vtk_type = types[id];
701  if( vtk_type >= VtkUtil::numVtkElemType ) return MB_FAILURE;
702 
703  EntityType type = VtkUtil::vtkElemTypes[vtk_type].mb_type;
704  if( type == MBMAXTYPE )
705  {
706  MB_SET_ERR( MB_FAILURE, "Unsupported VTK element type: " << VtkUtil::vtkElemTypes[vtk_type].name << " ("
707  << vtk_type << ")" );
708  }
709 
710  int num_vtx = *conn_iter;
711  if( type != MBPOLYGON && type != MBPOLYHEDRON && num_vtx != (int)VtkUtil::vtkElemTypes[vtk_type].num_nodes )
712  {
713  MB_SET_ERR( MB_FAILURE, "Cell " << id << " is of type '" << VtkUtil::vtkElemTypes[vtk_type].name
714  << "' but has " << num_vtx << " vertices" );
715  }
716 
717  // Find any subsequent elements of the same type
718  // if polyhedra, need to look at the number of faces to put in the same range
719  std::vector< long >::iterator conn_iter2 = conn_iter + num_vtx + 1;
720  long end_id = id + 1;
721  if( MBPOLYHEDRON != type )
722  {
723  while( end_id < num_elems[0] && (unsigned)types[end_id] == vtk_type && *conn_iter2 == num_vtx )
724  {
725  ++end_id;
726  conn_iter2 += num_vtx + 1;
727  }
728  }
729  else
730  {
731  // advance only if next is polyhedron too, and if number of faces is the same
732  int num_faces = conn_iter[1];
733  while( end_id < num_elems[0] && (unsigned)types[end_id] == vtk_type && conn_iter2[1] == num_faces )
734  {
735  ++end_id;
736  conn_iter2 += num_vtx + 1;
737  }
738  // num_vtx becomes in this case num_faces
739  num_vtx = num_faces; // for polyhedra, this is what we want
740  // trigger vertex adjacency call
741  Range firstFaces;
742  mdbImpl->get_adjacencies( &first_vertex, 1, 2, false, firstFaces );
743  }
744 
745  // Allocate element block
746  long num_elem = end_id - id;
747  EntityHandle start_handle = 0;
748  EntityHandle* conn_array;
749 
750  // if type is MBVERTEX, skip, we will not create elements with one vertex
751  if( MBVERTEX == type )
752  {
753  id += num_elem;
754  conn_iter += 2 * num_elem; // skip 2 * number of 1-vertex elements
755  continue;
756  }
757 
758  result = allocate_elements( num_elem, num_vtx, type, start_handle, conn_array, elem_list );
759  if( MB_SUCCESS != result ) return result;
760 
761  EntityHandle* conn_sav = conn_array;
762 
763  // Store element connectivity
764  if( type != MBPOLYHEDRON )
765  {
766  for( ; id < end_id; ++id )
767  {
768  if( conn_iter == connectivity.end() )
769  {
770  MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at cell " << id );
771  }
772  // Make sure connectivity length is correct.
773  if( *conn_iter != num_vtx )
774  {
775  MB_SET_ERR( MB_FAILURE, "Cell " << id << " is of type '" << VtkUtil::vtkElemTypes[vtk_type].name
776  << "' but has " << num_vtx << " vertices" );
777  }
778  ++conn_iter;
779 
780  for( i = 0; i < num_vtx; ++i, ++conn_iter )
781  {
782  if( conn_iter == connectivity.end() )
783  {
784  MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at cell " << id );
785  }
786 
787  conn_array[i] = *conn_iter + first_vertex;
788  }
789 
790  const unsigned* order = VtkUtil::vtkElemTypes[vtk_type].node_order;
791  if( order )
792  {
793  assert( num_vtx * sizeof( EntityHandle ) <= sizeof( tmp_conn_list ) );
794  memcpy( tmp_conn_list, conn_array, num_vtx * sizeof( EntityHandle ) );
795  for( int j = 0; j < num_vtx; ++j )
796  conn_array[order[j]] = tmp_conn_list[j];
797  }
798 
799  conn_array += num_vtx;
800  }
801  }
802  else // type == MBPOLYHEDRON
803  {
804  // in some cases, we may need to create new elements; will it screw the tags?
805  // not if the file was not from moab
806  ErrorCode rv = MB_SUCCESS;
807  for( ; id < end_id; ++id )
808  {
809  if( conn_iter == connectivity.end() )
810  {
811  MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at polyhedra cell " << id );
812  }
813  ++conn_iter;
814  // iterator is now at number of faces
815  // we should check it is indeed num_vtx
816  int num_faces = *conn_iter;
817  if( num_faces != num_vtx ) MB_SET_ERR( MB_FAILURE, "Connectivity data wrong at polyhedra cell " << id );
818 
819  EntityHandle connec[20]; // we bet we will have only 20 vertices at most, in a
820  // face in a polyhedra
821  for( int j = 0; j < num_faces; j++ )
822  {
823  conn_iter++;
824  int numverticesInFace = (int)*conn_iter;
825  if( numverticesInFace > 20 )
826  MB_SET_ERR( MB_FAILURE,
827  "too many vertices in face index " << j << " for polyhedra cell " << id );
828  // need to find the face, but first fill with vertices
829  for( int k = 0; k < numverticesInFace; k++ )
830  {
831  connec[k] = first_vertex + *( ++conn_iter ); //
832  }
833  Range adjFaces;
834  // find a face with these vertices; if not, we need to create one, on the fly :(
835  rv = mdbImpl->get_adjacencies( connec, numverticesInFace, 2, false, adjFaces );MB_CHK_ERR( rv );
836  if( adjFaces.size() >= 1 )
837  {
838  conn_array[j] = adjFaces[0]; // get the first face found
839  }
840  else
841  {
842  // create the face; tri, quad or polygon
843  EntityType etype = MBTRI;
844  if( 4 == numverticesInFace ) etype = MBQUAD;
845  if( 4 < numverticesInFace ) etype = MBPOLYGON;
846 
847  rv = mdbImpl->create_element( etype, connec, numverticesInFace, conn_array[j] );MB_CHK_ERR( rv );
848  }
849  }
850 
851  conn_array += num_vtx; // advance for next polyhedra
852  conn_iter++; // advance to the next field
853  }
854  }
855 
856  // Notify MOAB of the new elements
857  result = readMeshIface->update_adjacencies( start_handle, num_elem, num_vtx, conn_sav );
858  if( MB_SUCCESS != result ) return result;
859  }
860 
861  return MB_SUCCESS;
862 }
863 
864 ErrorCode ReadVtk::vtk_create_structured_elems( const long* dims,
865  EntityHandle first_vtx,
866  std::vector< Range >& elem_list )
867 {
868  ErrorCode result;
869  // int non_zero[3] = {0, 0, 0}; // True if dim > 0 for x, y, z respectively
870  long elem_dim = 0; // Element dimension (2->quad, 3->hex)
871  long num_elems = 1; // Total number of elements
872  long vert_per_elem; // Element connectivity length
873  long edims[3] = { 1, 1, 1 }; // Number of elements in each grid direction
874 
875  // Populate above data
876  for( int d = 0; d < 3; d++ )
877  {
878  if( dims[d] > 1 )
879  {
880  // non_zero[elem_dim] = d;
881  ++elem_dim;
882  edims[d] = dims[d] - 1;
883  num_elems *= edims[d];
884  }
885  }
886  vert_per_elem = 1 << elem_dim;
887 
888  // Get element type from element dimension
889  EntityType type;
890  switch( elem_dim )
891  {
892  case 1:
893  type = MBEDGE;
894  break;
895  case 2:
896  type = MBQUAD;
897  break;
898  case 3:
899  type = MBHEX;
900  break;
901  default:
902  MB_SET_ERR( MB_FAILURE, "Invalid dimension for structured elements: " << elem_dim );
903  }
904 
905  // Allocate storage for elements
906  EntityHandle start_handle = 0;
907  EntityHandle* conn_array;
908  result = allocate_elements( num_elems, vert_per_elem, type, start_handle, conn_array, elem_list );
909  if( MB_SUCCESS != result ) return MB_FAILURE;
910 
911  EntityHandle* conn_sav = conn_array;
912 
913  // Offsets of element vertices in grid relative to corner closest to origin
914  long k = dims[0] * dims[1];
915  const long corners[8] = { 0, 1, 1 + dims[0], dims[0], k, k + 1, k + 1 + dims[0], k + dims[0] };
916 
917  // Populate element list
918  for( long z = 0; z < edims[2]; ++z )
919  for( long y = 0; y < edims[1]; ++y )
920  for( long x = 0; x < edims[0]; ++x )
921  {
922  const long index = x + y * dims[0] + z * ( dims[0] * dims[1] );
923  for( long j = 0; j < vert_per_elem; ++j, ++conn_array )
924  *conn_array = index + corners[j] + first_vtx;
925  }
926 
927  // Notify MOAB of the new elements
928  result = readMeshIface->update_adjacencies( start_handle, num_elems, vert_per_elem, conn_sav );
929  if( MB_SUCCESS != result ) return result;
930 
931  return MB_SUCCESS;
932 }
933 
934 ErrorCode ReadVtk::vtk_read_field( FileTokenizer& tokens )
935 {
936  // This is not supported yet.
937  // Parse the data but throw it away because
938  // Mesquite has no internal representation for it.
939 
940  // Could save this in tags, but the only useful thing that
941  // could be done with the data is to write it back out
942  // with the modified mesh. As there's no way to save the
943  // type of a tag in Mesquite, it cannot be written back
944  // out correctly either.
945  // FIXME: Don't know what to do with this data.
946  // For now, read it and throw it out.
947 
948  long num_arrays;
949  if( !tokens.get_string() || // Name
950  !tokens.get_long_ints( 1, &num_arrays ) )
951  return MB_FAILURE;
952 
953  for( long i = 0; i < num_arrays; ++i )
954  {
955  /*const char* name =*/tokens.get_string();
956 
957  long dims[2];
958  if( !tokens.get_long_ints( 2, dims ) || !tokens.match_token( vtk_type_names ) ) return MB_FAILURE;
959 
960  long num_vals = dims[0] * dims[1];
961 
962  for( long j = 0; j < num_vals; j++ )
963  {
964  double junk;
965  if( !tokens.get_doubles( 1, &junk ) ) return MB_FAILURE;
966  }
967  }
968 
969  return MB_SUCCESS;
970 }
971 
972 ErrorCode ReadVtk::vtk_read_attrib_data( FileTokenizer& tokens, std::vector< Range >& entities )
973 {
974  const char* const type_names[] = { "SCALARS", "COLOR_SCALARS", "VECTORS", "NORMALS", "TEXTURE_COORDINATES",
975  "TENSORS", "FIELD", 0 };
976 
977  int type = tokens.match_token( type_names );
978  const char* tmp_name = tokens.get_string();
979  if( !type || !tmp_name ) return MB_FAILURE;
980 
981  std::string name_alloc( tmp_name );
982  const char* name = name_alloc.c_str();
983  switch( type )
984  {
985  case 1:
986  return vtk_read_scalar_attrib( tokens, entities, name );
987  case 2:
988  return vtk_read_color_attrib( tokens, entities, name );
989  case 3:
990  return vtk_read_vector_attrib( tokens, entities, name );
991  case 4:
992  return vtk_read_vector_attrib( tokens, entities, name );
993  case 5:
994  return vtk_read_texture_attrib( tokens, entities, name );
995  case 6:
996  return vtk_read_tensor_attrib( tokens, entities, name );
997  case 7:
998  return vtk_read_field_attrib( tokens, entities, name );
999  }
1000 
1001  return MB_FAILURE;
1002 }
1003 
1004 ErrorCode ReadVtk::vtk_read_tag_data( FileTokenizer& tokens,
1005  int type,
1006  size_t per_elem,
1007  std::vector< Range >& entities,
1008  const char* name )
1009 {
1010  ErrorCode result;
1011  DataType mb_type;
1012  if( type == 1 )
1013  {
1014  mb_type = MB_TYPE_BIT;
1015  }
1016  else if( type >= 2 && type <= 9 )
1017  {
1018  mb_type = MB_TYPE_INTEGER;
1019  }
1020  else if( type == 10 || type == 11 )
1021  {
1022  mb_type = MB_TYPE_DOUBLE;
1023  }
1024  else if( type == 12 )
1025  {
1026  mb_type = MB_TYPE_INTEGER;
1027  }
1028  else
1029  return MB_FAILURE;
1030 
1031 #ifdef MB_VTK_MATERIAL_SETS
1032  size_t size;
1033  if( type == 1 )
1034  {
1035  size = sizeof( bool );
1036  }
1037  else if( type >= 2 && type <= 9 )
1038  {
1039  size = sizeof( int );
1040  }
1041  else if( type == 10 || type == 11 )
1042  {
1043  size = sizeof( double );
1044  }
1045  else /* (type == 12) */
1046  {
1047  size = 4; // Could be 4 or 8, but we don't know. Hope it's 4 because MOAB doesn't support
1048  // 64-bit ints.
1049  }
1050  Modulator materialMap( this->mdbImpl );
1051  result = materialMap.initialize( this->mPartitionTagName, mb_type, size, per_elem );MB_CHK_SET_ERR( result, "MaterialMap tag (" << this->mPartitionTagName << ") creation failed." );
1052  bool isMaterial = size * per_elem <= 4 && // Must have int-sized values (ParallelComm requires it)
1053  !this->mPartitionTagName.empty() && // Must have a non-empty field name...
1054  !strcmp( name, this->mPartitionTagName.c_str() ); // ... that matches our spec.
1055 #endif // MB_VTK_MATERIAL_SETS
1056 
1057  // Get/create tag
1058  Tag handle;
1059  result = mdbImpl->tag_get_handle( name, per_elem, mb_type, handle, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( result, "Tag name conflict for attribute \"" << name << "\" at line " << tokens.line_number() );
1060 
1061  std::vector< Range >::iterator iter;
1062 
1063  if( type == 1 )
1064  {
1065  for( iter = entities.begin(); iter != entities.end(); ++iter )
1066  {
1067  bool* data = new bool[iter->size() * per_elem];
1068  if( !tokens.get_booleans( per_elem * iter->size(), data ) )
1069  {
1070  delete[] data;
1071  return MB_FAILURE;
1072  }
1073 
1074  bool* data_iter = data;
1075  Range::iterator ent_iter = iter->begin();
1076  for( ; ent_iter != iter->end(); ++ent_iter )
1077  {
1078  unsigned char bits = 0;
1079  for( unsigned j = 0; j < per_elem; ++j, ++data_iter )
1080  bits |= (unsigned char)( *data_iter << j );
1081 #ifdef MB_VTK_MATERIAL_SETS
1082  if( isMaterial ) materialMap.add_entity( *ent_iter, &bits, 1 );
1083 #endif // MB_VTK_MATERIAL_SETS
1084  result = mdbImpl->tag_set_data( handle, &*ent_iter, 1, &bits );
1085  if( MB_SUCCESS != result )
1086  {
1087  delete[] data;
1088  return result;
1089  }
1090  }
1091  delete[] data;
1092  }
1093  }
1094  else if( ( type >= 2 && type <= 9 ) || type == 12 )
1095  {
1096  std::vector< int > data;
1097  for( iter = entities.begin(); iter != entities.end(); ++iter )
1098  {
1099  data.resize( iter->size() * per_elem );
1100  if( !tokens.get_integers( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
1101 #ifdef MB_VTK_MATERIAL_SETS
1102  if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
1103 #endif // MB_VTK_MATERIAL_SETS
1104  result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
1105  if( MB_SUCCESS != result ) return result;
1106  }
1107  }
1108  else if( type == 10 || type == 11 )
1109  {
1110  std::vector< double > data;
1111  for( iter = entities.begin(); iter != entities.end(); ++iter )
1112  {
1113  data.resize( iter->size() * per_elem );
1114  if( !tokens.get_doubles( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
1115 #ifdef MB_VTK_MATERIAL_SETS
1116  if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
1117 #endif // MB_VTK_MATERIAL_SETS
1118  result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
1119  if( MB_SUCCESS != result ) return result;
1120  }
1121  }
1122 
1123  return MB_SUCCESS;
1124 }
1125 
1126 ErrorCode ReadVtk::vtk_read_scalar_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1127 {
1128  int type = tokens.match_token( vtk_type_names );
1129  if( !type ) return MB_FAILURE;
1130 
1131  long size;
1132  const char* tok = tokens.get_string();
1133  if( !tok ) return MB_FAILURE;
1134 
1135  const char* end = 0;
1136  size = strtol( tok, (char**)&end, 0 );
1137  if( *end )
1138  {
1139  size = 1;
1140  tokens.unget_token();
1141  }
1142 
1143  // VTK spec says cannot be greater than 4--do we care?
1144  if( size < 1 )
1145  { //|| size > 4)
1146  MB_SET_ERR( MB_FAILURE, "Scalar count out of range [1,4] at line " << tokens.line_number() );
1147  }
1148 
1149  if( !tokens.match_token( "LOOKUP_TABLE" ) || !tokens.match_token( "default" ) ) return MB_FAILURE;
1150 
1151  return vtk_read_tag_data( tokens, type, size, entities, name );
1152 }
1153 
1154 ErrorCode ReadVtk::vtk_read_color_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1155 {
1156  long size;
1157  if( !tokens.get_long_ints( 1, &size ) || size < 1 ) return MB_FAILURE;
1158 
1159  return vtk_read_tag_data( tokens, 10, size, entities, name );
1160 }
1161 
1162 ErrorCode ReadVtk::vtk_read_vector_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1163 {
1164  int type = tokens.match_token( vtk_type_names );
1165  if( !type ) return MB_FAILURE;
1166 
1167  return vtk_read_tag_data( tokens, type, 3, entities, name );
1168 }
1169 
1170 ErrorCode ReadVtk::vtk_read_texture_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1171 {
1172  int type, dim;
1173  if( !tokens.get_integers( 1, &dim ) || !( type = tokens.match_token( vtk_type_names ) ) ) return MB_FAILURE;
1174 
1175  if( dim < 1 || dim > 3 )
1176  {
1177  MB_SET_ERR( MB_FAILURE, "Invalid dimension (" << dim << ") at line " << tokens.line_number() );
1178  }
1179 
1180  return vtk_read_tag_data( tokens, type, dim, entities, name );
1181 }
1182 
1183 ErrorCode ReadVtk::vtk_read_tensor_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1184 {
1185  int type = tokens.match_token( vtk_type_names );
1186  if( !type ) return MB_FAILURE;
1187 
1188  return vtk_read_tag_data( tokens, type, 9, entities, name );
1189 }
1190 
1191 ErrorCode ReadVtk::vtk_read_field_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* )
1192 {
1193  long num_fields;
1194  if( !tokens.get_long_ints( 1, &num_fields ) ) return MB_FAILURE;
1195 
1196  long i;
1197  for( i = 0; i < num_fields; ++i )
1198  {
1199  const char* tok = tokens.get_string();
1200  if( !tok ) return MB_FAILURE;
1201 
1202  std::string name_alloc( tok );
1203 
1204  long num_comp;
1205  if( !tokens.get_long_ints( 1, &num_comp ) ) return MB_FAILURE;
1206 
1207  long num_tuples;
1208  if( !tokens.get_long_ints( 1, &num_tuples ) ) return MB_FAILURE;
1209 
1210  int type = tokens.match_token( vtk_type_names );
1211  if( !type ) return MB_FAILURE;
1212 
1213  ErrorCode result = vtk_read_tag_data( tokens, type, num_comp, entities, name_alloc.c_str() );MB_CHK_SET_ERR( result, "Error reading data for field \"" << name_alloc << "\" (" << num_comp << " components, "
1214  << num_tuples << " tuples, type " << type
1215  << ") at line " << tokens.line_number() );
1216  }
1217 
1218  return MB_SUCCESS;
1219 }
1220 
1221 ErrorCode ReadVtk::store_file_ids( Tag tag, const Range& verts, const std::vector< Range >& elems )
1222 {
1223  ErrorCode rval;
1224 
1225  rval = readMeshIface->assign_ids( tag, verts );
1226  if( MB_SUCCESS != rval ) return rval;
1227 
1228  int id = 0;
1229  for( size_t i = 0; i < elems.size(); ++i )
1230  {
1231  rval = readMeshIface->assign_ids( tag, elems[i], id );
1232  id += elems[i].size();
1233  }
1234 
1235  return MB_SUCCESS;
1236 }
1237 
1238 } // namespace moab