Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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