Mesh Oriented datABase  (version 5.6.0)
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 );
114  MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
115  }
116 
117  void add_entities( Range& range, const unsigned char* bytes, size_t bytes_per_ent )
118  {
119  for( Range::iterator it = range.begin(); it != range.end(); ++it, bytes += bytes_per_ent )
120  {
121  Hash h( bytes, bytes_per_ent );
122  EntityHandle mset = this->congruence_class( h, bytes );
123  ErrorCode rval;
124  rval = this->mesh->add_entities( mset, &*it, 1 );
125  MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
126  }
127  }
128 
129  EntityHandle congruence_class( const Hash& h, const void* tag_data )
130  {
131  std::map< Hash, EntityHandle >::iterator it = this->find( h );
132  if( it == this->end() )
133  {
134  EntityHandle mset;
135  Range preexist;
136  ErrorCode rval;
137  rval = this->mesh->get_entities_by_type_and_tag( 0, MBENTITYSET, &this->tag, &tag_data, 1, preexist );
138  MB_CHK_SET_ERR_RET_VAL( rval, "Failed to get entities by type and tag", (EntityHandle)0 );
139  if( preexist.size() )
140  {
141  mset = *preexist.begin();
142  }
143  else
144  {
145  rval = this->mesh->create_meshset( MESHSET_SET, mset );
146  MB_CHK_SET_ERR_RET_VAL( rval, "Failed to create mesh set", (EntityHandle)0 );
147  rval = this->mesh->tag_set_data( this->tag, &mset, 1, tag_data );
148  MB_CHK_SET_ERR_RET_VAL( rval, "Failed to set tag data", (EntityHandle)0 );
149  }
150  ( *this )[h] = mset;
151  return mset;
152  }
153  return it->second;
154  }
155 
156  Interface* mesh;
157  Tag tag;
158 };
159 #endif // MB_VTK_MATERIAL_SETS
160 
161 ReaderIface* ReadVtk::factory( Interface* iface )
162 {
163  return new ReadVtk( iface );
164 }
165 
166 ReadVtk::ReadVtk( Interface* impl ) : mdbImpl( impl ), mPartitionTagName( MATERIAL_SET_TAG_NAME )
167 {
168  mdbImpl->query_interface( readMeshIface );
169 }
170 
171 ReadVtk::~ReadVtk()
172 {
173  if( readMeshIface )
174  {
175  mdbImpl->release_interface( readMeshIface );
176  readMeshIface = 0;
177  }
178 }
179 
180 const char* const vtk_type_names[] = {
181  "bit", "char", "unsigned_char", "short", "unsigned_short", "int", "unsigned_int",
182  "long", "unsigned_long", "float", "double", "vtkIdType", 0 };
183 
184 ErrorCode ReadVtk::read_tag_values( const char* /* file_name */,
185  const char* /* tag_name */,
186  const FileOptions& /* opts */,
187  std::vector< int >& /* tag_values_out */,
188  const SubsetList* /* subset_list */ )
189 {
190  return MB_NOT_IMPLEMENTED;
191 }
192 
193 ErrorCode ReadVtk::load_file( const char* filename,
194  const EntityHandle* /* file_set */,
195  const FileOptions& opts,
196  const ReaderIface::SubsetList* subset_list,
197  const Tag* file_id_tag )
198 {
199  ErrorCode result;
200 
201  int major, minor;
202  char vendor_string[257];
203  std::vector< Range > element_list;
204  Range vertices;
205 
206  if( subset_list )
207  {
208  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for VTK" );
209  }
210 
211  // Does the caller want a field to be used for partitioning the entities?
212  // If not, we'll assume any scalar integer field named MATERIAL_SET specifies partitions.
213  std::string partition_tag_name;
214  result = opts.get_option( "PARTITION", partition_tag_name );
215  if( result == MB_SUCCESS ) mPartitionTagName = partition_tag_name;
216 
217  FILE* file = fopen( filename, "r" );
218  if( !file ) return MB_FILE_DOES_NOT_EXIST;
219 
220  // Read file header
221 
222  if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
223  {
224  fclose( file );
225  return MB_FAILURE;
226  }
227 
228  if( !strchr( vendor_string, '\n' ) || 2 != sscanf( vendor_string, "# vtk DataFile Version %d.%d", &major, &minor ) )
229  {
230  fclose( file );
231  return MB_FAILURE;
232  }
233 
234  if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
235  {
236  fclose( file );
237  return MB_FAILURE;
238  }
239 
240  // VTK spec says this should not exceed 256 chars.
241  if( !strchr( vendor_string, '\n' ) )
242  {
243  fclose( file );
244  MB_SET_ERR( MB_FAILURE, "Vendor string (line 2) exceeds 256 characters" );
245  }
246 
247  // Check file type
248 
249  FileTokenizer tokens( file, readMeshIface );
250  const char* const file_type_names[] = { "ASCII", "BINARY", 0 };
251  int filetype = tokens.match_token( file_type_names );
252  switch( filetype )
253  {
254  case 2: // BINARY
255  MB_SET_ERR( MB_FAILURE, "Cannot read BINARY VTK files" );
256  default: // ERROR
257  return MB_FAILURE;
258  case 1: // ASCII
259  break;
260  }
261 
262  // Read the mesh
263  if( !tokens.match_token( "DATASET" ) ) return MB_FAILURE;
264  result = vtk_read_dataset( tokens, vertices, element_list );
265  if( MB_SUCCESS != result ) return result;
266 
267  if( file_id_tag )
268  {
269  result = store_file_ids( *file_id_tag, vertices, element_list );
270  if( MB_SUCCESS != result ) return result;
271  }
272 
273  // Count the number of elements read
274  long elem_count = 0;
275  for( std::vector< Range >::iterator it = element_list.begin(); it != element_list.end(); ++it )
276  elem_count += it->size();
277 
278  // Read attribute data until end of file.
279  const char* const block_type_names[] = { "POINT_DATA", "CELL_DATA", 0 };
280  std::vector< Range > vertex_list( 1 );
281  vertex_list[0] = vertices;
282  int blocktype = 0;
283  while( !tokens.eof() )
284  {
285  // Get POINT_DATA or CELL_DATA
286  int new_block_type = tokens.match_token( block_type_names, false );
287  if( tokens.eof() ) break;
288 
289  if( !new_block_type )
290  {
291  // If next token was neither POINT_DATA nor CELL_DATA,
292  // then there's another attribute under the current one.
293  if( blocktype )
294  tokens.unget_token();
295  else
296  break;
297  }
298  else
299  {
300  blocktype = new_block_type;
301  long count;
302  if( !tokens.get_long_ints( 1, &count ) ) return MB_FAILURE;
303 
304  if( blocktype == 1 && (unsigned long)count != vertices.size() )
305  {
306  MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of vertices at line " << tokens.line_number() );
307  }
308  else if( blocktype == 2 && count != elem_count )
309  {
310  MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of elements at line " << tokens.line_number() );
311  }
312  }
313 
314  if( blocktype == 1 )
315  result = vtk_read_attrib_data( tokens, vertex_list );
316  else
317  result = vtk_read_attrib_data( tokens, element_list );
318 
319  if( MB_SUCCESS != result ) return result;
320  }
321 
322  return MB_SUCCESS;
323 }
324 
325 ErrorCode ReadVtk::allocate_vertices( long num_verts,
326  EntityHandle& start_handle_out,
327  double*& x_coord_array_out,
328  double*& y_coord_array_out,
329  double*& z_coord_array_out )
330 {
331  ErrorCode result;
332 
333  // Create vertices
334  std::vector< double* > arrays;
335  start_handle_out = 0;
336  result = readMeshIface->get_node_coords( 3, num_verts, MB_START_ID, start_handle_out, arrays );
337  if( MB_SUCCESS != result ) return result;
338 
339  x_coord_array_out = arrays[0];
340  y_coord_array_out = arrays[1];
341  z_coord_array_out = arrays[2];
342 
343  return MB_SUCCESS;
344 }
345 
346 ErrorCode ReadVtk::read_vertices( FileTokenizer& tokens, long num_verts, EntityHandle& start_handle_out )
347 {
348  ErrorCode result;
349  double *x, *y, *z;
350 
351  result = allocate_vertices( num_verts, start_handle_out, x, y, z );
352  if( MB_SUCCESS != result ) return result;
353 
354  // Read vertex coordinates
355  for( long vtx = 0; vtx < num_verts; ++vtx )
356  {
357  if( !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) || !tokens.get_doubles( 1, z++ ) )
358  return MB_FAILURE;
359  }
360 
361  return MB_SUCCESS;
362 }
363 
364 ErrorCode ReadVtk::allocate_elements( long num_elements,
365  int vert_per_element,
366  EntityType type,
367  EntityHandle& start_handle_out,
368  EntityHandle*& conn_array_out,
369  std::vector< Range >& append_to_this )
370 {
371  ErrorCode result;
372 
373  start_handle_out = 0;
374  result = readMeshIface->get_element_connect( num_elements, vert_per_element, type, MB_START_ID, start_handle_out,
375  conn_array_out );
376  if( MB_SUCCESS != result ) return result;
377 
378  Range range( start_handle_out, start_handle_out + num_elements - 1 );
379  append_to_this.push_back( range );
380  return MB_SUCCESS;
381 }
382 
383 ErrorCode ReadVtk::vtk_read_dataset( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& element_list )
384 {
385  const char* const data_type_names[] = {
386  "STRUCTURED_POINTS", "STRUCTURED_GRID", "UNSTRUCTURED_GRID", "POLYDATA", "RECTILINEAR_GRID", "FIELD", 0 };
387  int datatype = tokens.match_token( data_type_names );
388  switch( datatype )
389  {
390  case 1:
391  return vtk_read_structured_points( tokens, vertex_list, element_list );
392  case 2:
393  return vtk_read_structured_grid( tokens, vertex_list, element_list );
394  case 3:
395  return vtk_read_unstructured_grid( tokens, vertex_list, element_list );
396  case 4:
397  return vtk_read_polydata( tokens, vertex_list, element_list );
398  case 5:
399  return vtk_read_rectilinear_grid( tokens, vertex_list, element_list );
400  case 6:
401  return vtk_read_field( tokens );
402  default:
403  return MB_FAILURE;
404  }
405 }
406 
407 ErrorCode ReadVtk::vtk_read_structured_points( FileTokenizer& tokens,
408  Range& vertex_list,
409  std::vector< Range >& elem_list )
410 {
411  long i, j, k;
412  long dims[3];
413  double origin[3], space[3];
414  ErrorCode result;
415 
416  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
417  return MB_FAILURE;
418 
419  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
420  {
421  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
422  }
423 
424  if( !tokens.match_token( "ORIGIN" ) || !tokens.get_doubles( 3, origin ) || !tokens.get_newline() )
425  return MB_FAILURE;
426 
427  const char* const spacing_names[] = { "SPACING", "ASPECT_RATIO", 0 };
428  if( !tokens.match_token( spacing_names ) || !tokens.get_doubles( 3, space ) || !tokens.get_newline() )
429  return MB_FAILURE;
430 
431  // Create vertices
432  double *x, *y, *z;
433  EntityHandle start_handle = 0;
434  long num_verts = dims[0] * dims[1] * dims[2];
435  result = allocate_vertices( num_verts, start_handle, x, y, z );
436  if( MB_SUCCESS != result ) return result;
437  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
438 
439  for( k = 0; k < dims[2]; ++k )
440  for( j = 0; j < dims[1]; ++j )
441  for( i = 0; i < dims[0]; ++i )
442  {
443  *x = origin[0] + i * space[0];
444  ++x;
445  *y = origin[1] + j * space[1];
446  ++y;
447  *z = origin[2] + k * space[2];
448  ++z;
449  }
450 
451  return vtk_create_structured_elems( dims, start_handle, elem_list );
452 }
453 
454 ErrorCode ReadVtk::vtk_read_structured_grid( FileTokenizer& tokens,
455  Range& vertex_list,
456  std::vector< Range >& elem_list )
457 {
458  long num_verts, dims[3];
459  ErrorCode result;
460 
461  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
462  return MB_FAILURE;
463 
464  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
465  {
466  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
467  }
468 
469  if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
470  !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
471  return MB_FAILURE;
472 
473  if( num_verts != ( dims[0] * dims[1] * dims[2] ) )
474  {
475  MB_SET_ERR( MB_FAILURE, "Point count not consistent with dimensions at line " << tokens.line_number() );
476  }
477 
478  // Create and read vertices
479  EntityHandle start_handle = 0;
480  result = read_vertices( tokens, num_verts, start_handle );
481  if( MB_SUCCESS != result ) return result;
482  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
483 
484  return vtk_create_structured_elems( dims, start_handle, elem_list );
485 }
486 
487 ErrorCode ReadVtk::vtk_read_rectilinear_grid( FileTokenizer& tokens,
488  Range& vertex_list,
489  std::vector< Range >& elem_list )
490 {
491  int i, j, k;
492  long dims[3];
493  const char* labels[] = { "X_COORDINATES", "Y_COORDINATES", "Z_COORDINATES" };
494  std::vector< double > coords[3];
495  ErrorCode result;
496 
497  if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
498  return MB_FAILURE;
499 
500  if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
501  {
502  MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
503  }
504 
505  for( i = 0; i < 3; i++ )
506  {
507  long count;
508  if( !tokens.match_token( labels[i] ) || !tokens.get_long_ints( 1, &count ) ||
509  !tokens.match_token( vtk_type_names ) )
510  return MB_FAILURE;
511 
512  if( count != dims[i] )
513  {
514  MB_SET_ERR( MB_FAILURE, "Coordinate count inconsistent with dimensions at line " << tokens.line_number() );
515  }
516 
517  coords[i].resize( count );
518  if( !tokens.get_doubles( count, &coords[i][0] ) ) return MB_FAILURE;
519  }
520 
521  // Create vertices
522  double *x, *y, *z;
523  EntityHandle start_handle = 0;
524  long num_verts = dims[0] * dims[1] * dims[2];
525  result = allocate_vertices( num_verts, start_handle, x, y, z );
526  if( MB_SUCCESS != result ) return result;
527  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
528 
529  // Calculate vertex coordinates
530  for( k = 0; k < dims[2]; ++k )
531  for( j = 0; j < dims[1]; ++j )
532  for( i = 0; i < dims[0]; ++i )
533  {
534  *x = coords[0][i];
535  ++x;
536  *y = coords[1][j];
537  ++y;
538  *z = coords[2][k];
539  ++z;
540  }
541 
542  return vtk_create_structured_elems( dims, start_handle, elem_list );
543 }
544 
545 ErrorCode ReadVtk::vtk_read_polydata( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& elem_list )
546 {
547  ErrorCode result;
548  long num_verts;
549  const char* const poly_data_names[] = { "VERTICES", "LINES", "POLYGONS", "TRIANGLE_STRIPS", 0 };
550 
551  if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
552  !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
553  return MB_FAILURE;
554 
555  if( num_verts < 1 )
556  {
557  MB_SET_ERR( MB_FAILURE, "Invalid point count at line " << tokens.line_number() );
558  }
559 
560  // Create vertices and read coordinates
561  EntityHandle start_handle = 0;
562  result = read_vertices( tokens, num_verts, start_handle );
563  if( MB_SUCCESS != result ) return result;
564  vertex_list.insert( start_handle, start_handle + num_verts - 1 );
565 
566  int poly_type = tokens.match_token( poly_data_names );
567  switch( poly_type )
568  {
569  case 0:
570  result = MB_FAILURE;
571  break;
572  case 1:
573  MB_SET_ERR( MB_FAILURE, "Vertex element type at line " << tokens.line_number() );
574  case 2:
575  MB_SET_ERR( MB_FAILURE, "Unsupported type: polylines at line " << tokens.line_number() );
576  case 3:
577  result = vtk_read_polygons( tokens, start_handle, elem_list );
578  break;
579  case 4:
580  MB_SET_ERR( MB_FAILURE, "Unsupported type: triangle strips at line " << tokens.line_number() );
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 );
836  MB_CHK_ERR( rv );
837  if( adjFaces.size() >= 1 )
838  {
839  conn_array[j] = adjFaces[0]; // get the first face found
840  }
841  else
842  {
843  // create the face; tri, quad or polygon
844  EntityType etype = MBTRI;
845  if( 4 == numverticesInFace ) etype = MBQUAD;
846  if( 4 < numverticesInFace ) etype = MBPOLYGON;
847 
848  rv = mdbImpl->create_element( etype, connec, numverticesInFace, conn_array[j] );
849  MB_CHK_ERR( rv );
850  }
851  }
852 
853  conn_array += num_vtx; // advance for next polyhedra
854  conn_iter++; // advance to the next field
855  }
856  }
857 
858  // Notify MOAB of the new elements
859  result = readMeshIface->update_adjacencies( start_handle, num_elem, num_vtx, conn_sav );
860  if( MB_SUCCESS != result ) return result;
861  }
862 
863  return MB_SUCCESS;
864 }
865 
866 ErrorCode ReadVtk::vtk_create_structured_elems( const long* dims,
867  EntityHandle first_vtx,
868  std::vector< Range >& elem_list )
869 {
870  ErrorCode result;
871  // int non_zero[3] = {0, 0, 0}; // True if dim > 0 for x, y, z respectively
872  long elem_dim = 0; // Element dimension (2->quad, 3->hex)
873  long num_elems = 1; // Total number of elements
874  long vert_per_elem; // Element connectivity length
875  long edims[3] = { 1, 1, 1 }; // Number of elements in each grid direction
876 
877  // Populate above data
878  for( int d = 0; d < 3; d++ )
879  {
880  if( dims[d] > 1 )
881  {
882  // non_zero[elem_dim] = d;
883  ++elem_dim;
884  edims[d] = dims[d] - 1;
885  num_elems *= edims[d];
886  }
887  }
888  vert_per_elem = 1 << elem_dim;
889 
890  // Get element type from element dimension
891  EntityType type;
892  switch( elem_dim )
893  {
894  case 1:
895  type = MBEDGE;
896  break;
897  case 2:
898  type = MBQUAD;
899  break;
900  case 3:
901  type = MBHEX;
902  break;
903  default:
904  MB_SET_ERR( MB_FAILURE, "Invalid dimension for structured elements: " << elem_dim );
905  }
906 
907  // Allocate storage for elements
908  EntityHandle start_handle = 0;
909  EntityHandle* conn_array;
910  result = allocate_elements( num_elems, vert_per_elem, type, start_handle, conn_array, elem_list );
911  if( MB_SUCCESS != result ) return MB_FAILURE;
912 
913  EntityHandle* conn_sav = conn_array;
914 
915  // Offsets of element vertices in grid relative to corner closest to origin
916  long k = dims[0] * dims[1];
917  const long corners[8] = { 0, 1, 1 + dims[0], dims[0], k, k + 1, k + 1 + dims[0], k + dims[0] };
918 
919  // Populate element list
920  for( long z = 0; z < edims[2]; ++z )
921  for( long y = 0; y < edims[1]; ++y )
922  for( long x = 0; x < edims[0]; ++x )
923  {
924  const long index = x + y * dims[0] + z * ( dims[0] * dims[1] );
925  for( long j = 0; j < vert_per_elem; ++j, ++conn_array )
926  *conn_array = index + corners[j] + first_vtx;
927  }
928 
929  // Notify MOAB of the new elements
930  result = readMeshIface->update_adjacencies( start_handle, num_elems, vert_per_elem, conn_sav );
931  if( MB_SUCCESS != result ) return result;
932 
933  return MB_SUCCESS;
934 }
935 
936 ErrorCode ReadVtk::vtk_read_field( FileTokenizer& tokens )
937 {
938  // This is not supported yet.
939  // Parse the data but throw it away because
940  // Mesquite has no internal representation for it.
941 
942  // Could save this in tags, but the only useful thing that
943  // could be done with the data is to write it back out
944  // with the modified mesh. As there's no way to save the
945  // type of a tag in Mesquite, it cannot be written back
946  // out correctly either.
947  // FIXME: Don't know what to do with this data.
948  // For now, read it and throw it out.
949 
950  long num_arrays;
951  if( !tokens.get_string() || // Name
952  !tokens.get_long_ints( 1, &num_arrays ) )
953  return MB_FAILURE;
954 
955  for( long i = 0; i < num_arrays; ++i )
956  {
957  /*const char* name =*/tokens.get_string();
958 
959  long dims[2];
960  if( !tokens.get_long_ints( 2, dims ) || !tokens.match_token( vtk_type_names ) ) return MB_FAILURE;
961 
962  long num_vals = dims[0] * dims[1];
963 
964  for( long j = 0; j < num_vals; j++ )
965  {
966  double junk;
967  if( !tokens.get_doubles( 1, &junk ) ) return MB_FAILURE;
968  }
969  }
970 
971  return MB_SUCCESS;
972 }
973 
974 ErrorCode ReadVtk::vtk_read_attrib_data( FileTokenizer& tokens, std::vector< Range >& entities )
975 {
976  const char* const type_names[] = { "SCALARS", "COLOR_SCALARS", "VECTORS", "NORMALS", "TEXTURE_COORDINATES",
977  "TENSORS", "FIELD", 0 };
978 
979  int type = tokens.match_token( type_names );
980  const char* tmp_name = tokens.get_string();
981  if( !type || !tmp_name ) return MB_FAILURE;
982 
983  std::string name_alloc( tmp_name );
984  const char* name = name_alloc.c_str();
985  switch( type )
986  {
987  case 1:
988  return vtk_read_scalar_attrib( tokens, entities, name );
989  case 2:
990  return vtk_read_color_attrib( tokens, entities, name );
991  case 3:
992  return vtk_read_vector_attrib( tokens, entities, name );
993  case 4:
994  return vtk_read_vector_attrib( tokens, entities, name );
995  case 5:
996  return vtk_read_texture_attrib( tokens, entities, name );
997  case 6:
998  return vtk_read_tensor_attrib( tokens, entities, name );
999  case 7:
1000  return vtk_read_field_attrib( tokens, entities, name );
1001  }
1002 
1003  return MB_FAILURE;
1004 }
1005 
1006 ErrorCode ReadVtk::vtk_read_tag_data( FileTokenizer& tokens,
1007  int type,
1008  size_t per_elem,
1009  std::vector< Range >& entities,
1010  const char* name )
1011 {
1012  ErrorCode result;
1013  DataType mb_type;
1014  if( type == 1 )
1015  {
1016  mb_type = MB_TYPE_BIT;
1017  }
1018  else if( type >= 2 && type <= 9 )
1019  {
1020  mb_type = MB_TYPE_INTEGER;
1021  }
1022  else if( type == 10 || type == 11 )
1023  {
1024  mb_type = MB_TYPE_DOUBLE;
1025  }
1026  else if( type == 12 )
1027  {
1028  mb_type = MB_TYPE_INTEGER;
1029  }
1030  else
1031  return MB_FAILURE;
1032 
1033 #ifdef MB_VTK_MATERIAL_SETS
1034  size_t size;
1035  if( type == 1 )
1036  {
1037  size = sizeof( bool );
1038  }
1039  else if( type >= 2 && type <= 9 )
1040  {
1041  size = sizeof( int );
1042  }
1043  else if( type == 10 || type == 11 )
1044  {
1045  size = sizeof( double );
1046  }
1047  else /* (type == 12) */
1048  {
1049  size = 4; // Could be 4 or 8, but we don't know. Hope it's 4 because MOAB doesn't support
1050  // 64-bit ints.
1051  }
1052  Modulator materialMap( this->mdbImpl );
1053  result = materialMap.initialize( this->mPartitionTagName, mb_type, size, per_elem );
1054  MB_CHK_SET_ERR( result, "MaterialMap tag (" << this->mPartitionTagName << ") creation failed." );
1055  bool isMaterial = size * per_elem <= 4 && // Must have int-sized values (ParallelComm requires it)
1056  !this->mPartitionTagName.empty() && // Must have a non-empty field name...
1057  !strcmp( name, this->mPartitionTagName.c_str() ); // ... that matches our spec.
1058 #endif // MB_VTK_MATERIAL_SETS
1059 
1060  // Get/create tag
1061  Tag handle;
1062  result = mdbImpl->tag_get_handle( name, per_elem, mb_type, handle, MB_TAG_DENSE | MB_TAG_CREAT );
1063  MB_CHK_SET_ERR( result, "Tag name conflict for attribute \"" << name << "\" at line " << tokens.line_number() );
1064 
1065  std::vector< Range >::iterator iter;
1066 
1067  if( type == 1 )
1068  {
1069  for( iter = entities.begin(); iter != entities.end(); ++iter )
1070  {
1071  bool* data = new bool[iter->size() * per_elem];
1072  if( !tokens.get_booleans( per_elem * iter->size(), data ) )
1073  {
1074  delete[] data;
1075  return MB_FAILURE;
1076  }
1077 
1078  bool* data_iter = data;
1079  Range::iterator ent_iter = iter->begin();
1080  for( ; ent_iter != iter->end(); ++ent_iter )
1081  {
1082  unsigned char bits = 0;
1083  for( unsigned j = 0; j < per_elem; ++j, ++data_iter )
1084  bits |= (unsigned char)( *data_iter << j );
1085 #ifdef MB_VTK_MATERIAL_SETS
1086  if( isMaterial ) materialMap.add_entity( *ent_iter, &bits, 1 );
1087 #endif // MB_VTK_MATERIAL_SETS
1088  result = mdbImpl->tag_set_data( handle, &*ent_iter, 1, &bits );
1089  if( MB_SUCCESS != result )
1090  {
1091  delete[] data;
1092  return result;
1093  }
1094  }
1095  delete[] data;
1096  }
1097  }
1098  else if( ( type >= 2 && type <= 9 ) || type == 12 )
1099  {
1100  std::vector< int > data;
1101  for( iter = entities.begin(); iter != entities.end(); ++iter )
1102  {
1103  data.resize( iter->size() * per_elem );
1104  if( !tokens.get_integers( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
1105 #ifdef MB_VTK_MATERIAL_SETS
1106  if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
1107 #endif // MB_VTK_MATERIAL_SETS
1108  result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
1109  if( MB_SUCCESS != result ) return result;
1110  }
1111  }
1112  else if( type == 10 || type == 11 )
1113  {
1114  std::vector< double > data;
1115  for( iter = entities.begin(); iter != entities.end(); ++iter )
1116  {
1117  data.resize( iter->size() * per_elem );
1118  if( !tokens.get_doubles( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
1119 #ifdef MB_VTK_MATERIAL_SETS
1120  if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
1121 #endif // MB_VTK_MATERIAL_SETS
1122  result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
1123  if( MB_SUCCESS != result ) return result;
1124  }
1125  }
1126 
1127  return MB_SUCCESS;
1128 }
1129 
1130 ErrorCode ReadVtk::vtk_read_scalar_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1131 {
1132  int type = tokens.match_token( vtk_type_names );
1133  if( !type ) return MB_FAILURE;
1134 
1135  long size;
1136  const char* tok = tokens.get_string();
1137  if( !tok ) return MB_FAILURE;
1138 
1139  const char* end = 0;
1140  size = strtol( tok, (char**)&end, 0 );
1141  if( *end )
1142  {
1143  size = 1;
1144  tokens.unget_token();
1145  }
1146 
1147  // VTK spec says cannot be greater than 4--do we care?
1148  if( size < 1 )
1149  { //|| size > 4)
1150  MB_SET_ERR( MB_FAILURE, "Scalar count out of range [1,4] at line " << tokens.line_number() );
1151  }
1152 
1153  if( !tokens.match_token( "LOOKUP_TABLE" ) || !tokens.match_token( "default" ) ) return MB_FAILURE;
1154 
1155  return vtk_read_tag_data( tokens, type, size, entities, name );
1156 }
1157 
1158 ErrorCode ReadVtk::vtk_read_color_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1159 {
1160  long size;
1161  if( !tokens.get_long_ints( 1, &size ) || size < 1 ) return MB_FAILURE;
1162 
1163  return vtk_read_tag_data( tokens, 10, size, entities, name );
1164 }
1165 
1166 ErrorCode ReadVtk::vtk_read_vector_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1167 {
1168  int type = tokens.match_token( vtk_type_names );
1169  if( !type ) return MB_FAILURE;
1170 
1171  return vtk_read_tag_data( tokens, type, 3, entities, name );
1172 }
1173 
1174 ErrorCode ReadVtk::vtk_read_texture_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1175 {
1176  int type, dim;
1177  if( !tokens.get_integers( 1, &dim ) || !( type = tokens.match_token( vtk_type_names ) ) ) return MB_FAILURE;
1178 
1179  if( dim < 1 || dim > 3 )
1180  {
1181  MB_SET_ERR( MB_FAILURE, "Invalid dimension (" << dim << ") at line " << tokens.line_number() );
1182  }
1183 
1184  return vtk_read_tag_data( tokens, type, dim, entities, name );
1185 }
1186 
1187 ErrorCode ReadVtk::vtk_read_tensor_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
1188 {
1189  int type = tokens.match_token( vtk_type_names );
1190  if( !type ) return MB_FAILURE;
1191 
1192  return vtk_read_tag_data( tokens, type, 9, entities, name );
1193 }
1194 
1195 ErrorCode ReadVtk::vtk_read_field_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* )
1196 {
1197  long num_fields;
1198  if( !tokens.get_long_ints( 1, &num_fields ) ) return MB_FAILURE;
1199 
1200  long i;
1201  for( i = 0; i < num_fields; ++i )
1202  {
1203  const char* tok = tokens.get_string();
1204  if( !tok ) return MB_FAILURE;
1205 
1206  std::string name_alloc( tok );
1207 
1208  long num_comp;
1209  if( !tokens.get_long_ints( 1, &num_comp ) ) return MB_FAILURE;
1210 
1211  long num_tuples;
1212  if( !tokens.get_long_ints( 1, &num_tuples ) ) return MB_FAILURE;
1213 
1214  int type = tokens.match_token( vtk_type_names );
1215  if( !type ) return MB_FAILURE;
1216 
1217  ErrorCode result = vtk_read_tag_data( tokens, type, num_comp, entities, name_alloc.c_str() );
1218  MB_CHK_SET_ERR( result, "Error reading data for field \"" << name_alloc << "\" (" << num_comp << " components, "
1219  << num_tuples << " tuples, type " << type
1220  << ") at line " << tokens.line_number() );
1221  }
1222 
1223  return MB_SUCCESS;
1224 }
1225 
1226 ErrorCode ReadVtk::store_file_ids( Tag tag, const Range& verts, const std::vector< Range >& elems )
1227 {
1228  ErrorCode rval;
1229 
1230  rval = readMeshIface->assign_ids( tag, verts );
1231  if( MB_SUCCESS != rval ) return rval;
1232 
1233  int id = 0;
1234  for( size_t i = 0; i < elems.size(); ++i )
1235  {
1236  rval = readMeshIface->assign_ids( tag, elems[i], id );
1237  id += elems[i].size();
1238  }
1239 
1240  return MB_SUCCESS;
1241 }
1242 
1243 } // namespace moab