ReadVtk.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

/**
 * \class ReadVtk
 * \brief VTK reader from Mesquite
 * \author Jason Kraftcheck
 */

#include 
#include 
#include 

#include "ReadVtk.hpp"
#include "moab/Range.hpp"
#include "Internals.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/FileOptions.hpp"
#include "FileTokenizer.hpp"
#include "moab/VtkUtil.hpp"
#include "MBTagConventions.hpp"

// #define MB_VTK_MATERIAL_SETS
#ifdef MB_VTK_MATERIAL_SETS
#include 
#endif

namespace moab
{

#ifdef MB_VTK_MATERIAL_SETS

class Hash
{
  public:
    unsigned long value;

    Hash()
    {
        this->value = 5381L;
    }

    Hash( const unsigned char* bytes, size_t len )
    {  // djb2, a hash by dan bernstein presented on comp.lang.c for hashing strings.
        this->value = 5381L;
        for( ; len; --len, ++bytes )
            this->value = this->value * 33 + ( *bytes );
    }

    Hash( bool duh )
    {
        this->value = duh;  // Hashing a single bit with a long is stupid but easy.
    }

    Hash( const Hash& src )
    {
        this->value = src.value;
    }

    Hash& operator=( const Hash& src )
    {
        this->value = src.value;
        return *this;
    }

    bool operator<( const Hash& other ) const
    {
        return this->value < other.value;
    }
};

// Pass this class opaque data + a handle and it adds the handle to a set
// whose opaque data has the same hash. If no set exists for the hash,
// it creates one. Each set is tagged with the opaque data.
// When entities with different opaque data have the same hash, they
// will be placed into the same set.
// There will be no collisions when the opaque data is shorter than an
// unsigned long, and this is really the only case we need to support.
// The rest is bonus. Hash does quite well with strings, even those
// with identical prefixes.
class Modulator : public std::map< Hash, EntityHandle >
{
  public:
    Modulator( Interface* iface ) : mesh( iface ) {}

    ErrorCode initialize( std::string tag_name, DataType mb_type, size_t sz, size_t per_elem )
    {
        std::vector< unsigned char > default_val;
        default_val.resize( sz * per_elem );
        ErrorCode rval = this->mesh->tag_get_handle( tag_name.c_str(), per_elem, mb_type, this->tag,
                                                     MB_TAG_SPARSE | MB_TAG_BYTES | MB_TAG_CREAT, &default_val[0] );
        return rval;
    }

    void add_entity( EntityHandle ent, const unsigned char* bytes, size_t len )
    {
        Hash h( bytes, len );
        EntityHandle mset = this->congruence_class( h, bytes );
        ErrorCode rval;
        rval = this->mesh->add_entities( mset, &ent, 1 );MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
    }

    void add_entities( Range& range, const unsigned char* bytes, size_t bytes_per_ent )
    {
        for( Range::iterator it = range.begin(); it != range.end(); ++it, bytes += bytes_per_ent )
        {
            Hash h( bytes, bytes_per_ent );
            EntityHandle mset = this->congruence_class( h, bytes );
            ErrorCode rval;
            rval = this->mesh->add_entities( mset, &*it, 1 );MB_CHK_SET_ERR_RET( rval, "Failed to add entities to mesh" );
        }
    }

    EntityHandle congruence_class( const Hash& h, const void* tag_data )
    {
        std::map< Hash, EntityHandle >::iterator it = this->find( h );
        if( it == this->end() )
        {
            EntityHandle mset;
            Range preexist;
            ErrorCode rval;
            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 );
            if( preexist.size() )
            {
                mset = *preexist.begin();
            }
            else
            {
                rval = this->mesh->create_meshset( MESHSET_SET, mset );MB_CHK_SET_ERR_RET_VAL( rval, "Failed to create mesh set", (EntityHandle)0 );
                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 );
            }
            ( *this )[h] = mset;
            return mset;
        }
        return it->second;
    }

    Interface* mesh;
    Tag tag;
};
#endif  // MB_VTK_MATERIAL_SETS

ReaderIface* ReadVtk::factory( Interface* iface )
{
    return new ReadVtk( iface );
}

ReadVtk::ReadVtk( Interface* impl ) : mdbImpl( impl ), mPartitionTagName( MATERIAL_SET_TAG_NAME )
{
    mdbImpl->query_interface( readMeshIface );
}

ReadVtk::~ReadVtk()
{
    if( readMeshIface )
    {
        mdbImpl->release_interface( readMeshIface );
        readMeshIface = 0;
    }
}

const char* const vtk_type_names[] = {
    "bit",  "char",          "unsigned_char", "short",  "unsigned_short", "int", "unsigned_int",
    "long", "unsigned_long", "float",         "double", "vtkIdType",      0 };

ErrorCode ReadVtk::read_tag_values( const char* /* file_name */,
                                    const char* /* tag_name */,
                                    const FileOptions& /* opts */,
                                    std::vector< int >& /* tag_values_out */,
                                    const SubsetList* /* subset_list */ )
{
    return MB_NOT_IMPLEMENTED;
}

ErrorCode ReadVtk::load_file( const char* filename,
                              const EntityHandle* /* file_set */,
                              const FileOptions& opts,
                              const ReaderIface::SubsetList* subset_list,
                              const Tag* file_id_tag )
{
    ErrorCode result;

    int major, minor;
    char vendor_string[257];
    std::vector< Range > element_list;
    Range vertices;

    if( subset_list )
    {
        MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for VTK" );
    }

    // Does the caller want a field to be used for partitioning the entities?
    // If not, we'll assume any scalar integer field named MATERIAL_SET specifies partitions.
    std::string partition_tag_name;
    result = opts.get_option( "PARTITION", partition_tag_name );
    if( result == MB_SUCCESS ) mPartitionTagName = partition_tag_name;

    FILE* file = fopen( filename, "r" );
    if( !file ) return MB_FILE_DOES_NOT_EXIST;

    // Read file header

    if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
    {
        fclose( file );
        return MB_FAILURE;
    }

    if( !strchr( vendor_string, '\n' ) || 2 != sscanf( vendor_string, "# vtk DataFile Version %d.%d", &major, &minor ) )
    {
        fclose( file );
        return MB_FAILURE;
    }

    if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
    {
        fclose( file );
        return MB_FAILURE;
    }

    // VTK spec says this should not exceed 256 chars.
    if( !strchr( vendor_string, '\n' ) )
    {
        fclose( file );
        MB_SET_ERR( MB_FAILURE, "Vendor string (line 2) exceeds 256 characters" );
    }

    // Check file type

    FileTokenizer tokens( file, readMeshIface );
    const char* const file_type_names[] = { "ASCII", "BINARY", 0 };
    int filetype                        = tokens.match_token( file_type_names );
    switch( filetype )
    {
        case 2:  // BINARY
            MB_SET_ERR( MB_FAILURE, "Cannot read BINARY VTK files" );
        default:  // ERROR
            return MB_FAILURE;
        case 1:  // ASCII
            break;
    }

    // Read the mesh
    if( !tokens.match_token( "DATASET" ) ) return MB_FAILURE;
    result = vtk_read_dataset( tokens, vertices, element_list );
    if( MB_SUCCESS != result ) return result;

    if( file_id_tag )
    {
        result = store_file_ids( *file_id_tag, vertices, element_list );
        if( MB_SUCCESS != result ) return result;
    }

    // Count the number of elements read
    long elem_count = 0;
    for( std::vector< Range >::iterator it = element_list.begin(); it != element_list.end(); ++it )
        elem_count += it->size();

    // Read attribute data until end of file.
    const char* const block_type_names[] = { "POINT_DATA", "CELL_DATA", 0 };
    std::vector< Range > vertex_list( 1 );
    vertex_list[0] = vertices;
    int blocktype  = 0;
    while( !tokens.eof() )
    {
        // Get POINT_DATA or CELL_DATA
        int new_block_type = tokens.match_token( block_type_names, false );
        if( tokens.eof() ) break;

        if( !new_block_type )
        {
            // If next token was neither POINT_DATA nor CELL_DATA,
            // then there's another attribute under the current one.
            if( blocktype )
                tokens.unget_token();
            else
                break;
        }
        else
        {
            blocktype = new_block_type;
            long count;
            if( !tokens.get_long_ints( 1, &count ) ) return MB_FAILURE;

            if( blocktype == 1 && (unsigned long)count != vertices.size() )
            {
                MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of vertices at line " << tokens.line_number() );
            }
            else if( blocktype == 2 && count != elem_count )
            {
                MB_SET_ERR( MB_FAILURE, "Count inconsistent with number of elements at line " << tokens.line_number() );
            }
        }

        if( blocktype == 1 )
            result = vtk_read_attrib_data( tokens, vertex_list );
        else
            result = vtk_read_attrib_data( tokens, element_list );

        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::allocate_vertices( long num_verts,
                                      EntityHandle& start_handle_out,
                                      double*& x_coord_array_out,
                                      double*& y_coord_array_out,
                                      double*& z_coord_array_out )
{
    ErrorCode result;

    // Create vertices
    std::vector< double* > arrays;
    start_handle_out = 0;
    result           = readMeshIface->get_node_coords( 3, num_verts, MB_START_ID, start_handle_out, arrays );
    if( MB_SUCCESS != result ) return result;

    x_coord_array_out = arrays[0];
    y_coord_array_out = arrays[1];
    z_coord_array_out = arrays[2];

    return MB_SUCCESS;
}

ErrorCode ReadVtk::read_vertices( FileTokenizer& tokens, long num_verts, EntityHandle& start_handle_out )
{
    ErrorCode result;
    double *x, *y, *z;

    result = allocate_vertices( num_verts, start_handle_out, x, y, z );
    if( MB_SUCCESS != result ) return result;

    // Read vertex coordinates
    for( long vtx = 0; vtx < num_verts; ++vtx )
    {
        if( !tokens.get_doubles( 1, x++ ) || !tokens.get_doubles( 1, y++ ) || !tokens.get_doubles( 1, z++ ) )
            return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::allocate_elements( long num_elements,
                                      int vert_per_element,
                                      EntityType type,
                                      EntityHandle& start_handle_out,
                                      EntityHandle*& conn_array_out,
                                      std::vector< Range >& append_to_this )
{
    ErrorCode result;

    start_handle_out = 0;
    result = readMeshIface->get_element_connect( num_elements, vert_per_element, type, MB_START_ID, start_handle_out,
                                                 conn_array_out );
    if( MB_SUCCESS != result ) return result;

    Range range( start_handle_out, start_handle_out + num_elements - 1 );
    append_to_this.push_back( range );
    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_read_dataset( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& element_list )
{
    const char* const data_type_names[] = {
        "STRUCTURED_POINTS", "STRUCTURED_GRID", "UNSTRUCTURED_GRID", "POLYDATA", "RECTILINEAR_GRID", "FIELD", 0 };
    int datatype = tokens.match_token( data_type_names );
    switch( datatype )
    {
        case 1:
            return vtk_read_structured_points( tokens, vertex_list, element_list );
        case 2:
            return vtk_read_structured_grid( tokens, vertex_list, element_list );
        case 3:
            return vtk_read_unstructured_grid( tokens, vertex_list, element_list );
        case 4:
            return vtk_read_polydata( tokens, vertex_list, element_list );
        case 5:
            return vtk_read_rectilinear_grid( tokens, vertex_list, element_list );
        case 6:
            return vtk_read_field( tokens );
        default:
            return MB_FAILURE;
    }
}

ErrorCode ReadVtk::vtk_read_structured_points( FileTokenizer& tokens,
                                               Range& vertex_list,
                                               std::vector< Range >& elem_list )
{
    long i, j, k;
    long dims[3];
    double origin[3], space[3];
    ErrorCode result;

    if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
    }

    if( !tokens.match_token( "ORIGIN" ) || !tokens.get_doubles( 3, origin ) || !tokens.get_newline() )
        return MB_FAILURE;

    const char* const spacing_names[] = { "SPACING", "ASPECT_RATIO", 0 };
    if( !tokens.match_token( spacing_names ) || !tokens.get_doubles( 3, space ) || !tokens.get_newline() )
        return MB_FAILURE;

    // Create vertices
    double *x, *y, *z;
    EntityHandle start_handle = 0;
    long num_verts            = dims[0] * dims[1] * dims[2];
    result                    = allocate_vertices( num_verts, start_handle, x, y, z );
    if( MB_SUCCESS != result ) return result;
    vertex_list.insert( start_handle, start_handle + num_verts - 1 );

    for( k = 0; k < dims[2]; ++k )
        for( j = 0; j < dims[1]; ++j )
            for( i = 0; i < dims[0]; ++i )
            {
                *x = origin[0] + i * space[0];
                ++x;
                *y = origin[1] + j * space[1];
                ++y;
                *z = origin[2] + k * space[2];
                ++z;
            }

    return vtk_create_structured_elems( dims, start_handle, elem_list );
}

ErrorCode ReadVtk::vtk_read_structured_grid( FileTokenizer& tokens,
                                             Range& vertex_list,
                                             std::vector< Range >& elem_list )
{
    long num_verts, dims[3];
    ErrorCode result;

    if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
    }

    if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
        !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( num_verts != ( dims[0] * dims[1] * dims[2] ) )
    {
        MB_SET_ERR( MB_FAILURE, "Point count not consistent with dimensions at line " << tokens.line_number() );
    }

    // Create and read vertices
    EntityHandle start_handle = 0;
    result                    = read_vertices( tokens, num_verts, start_handle );
    if( MB_SUCCESS != result ) return result;
    vertex_list.insert( start_handle, start_handle + num_verts - 1 );

    return vtk_create_structured_elems( dims, start_handle, elem_list );
}

ErrorCode ReadVtk::vtk_read_rectilinear_grid( FileTokenizer& tokens,
                                              Range& vertex_list,
                                              std::vector< Range >& elem_list )
{
    int i, j, k;
    long dims[3];
    const char* labels[] = { "X_COORDINATES", "Y_COORDINATES", "Z_COORDINATES" };
    std::vector< double > coords[3];
    ErrorCode result;

    if( !tokens.match_token( "DIMENSIONS" ) || !tokens.get_long_ints( 3, dims ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid dimension at line " << tokens.line_number() );
    }

    for( i = 0; i < 3; i++ )
    {
        long count;
        if( !tokens.match_token( labels[i] ) || !tokens.get_long_ints( 1, &count ) ||
            !tokens.match_token( vtk_type_names ) )
            return MB_FAILURE;

        if( count != dims[i] )
        {
            MB_SET_ERR( MB_FAILURE, "Coordinate count inconsistent with dimensions at line " << tokens.line_number() );
        }

        coords[i].resize( count );
        if( !tokens.get_doubles( count, &coords[i][0] ) ) return MB_FAILURE;
    }

    // Create vertices
    double *x, *y, *z;
    EntityHandle start_handle = 0;
    long num_verts            = dims[0] * dims[1] * dims[2];
    result                    = allocate_vertices( num_verts, start_handle, x, y, z );
    if( MB_SUCCESS != result ) return result;
    vertex_list.insert( start_handle, start_handle + num_verts - 1 );

    // Calculate vertex coordinates
    for( k = 0; k < dims[2]; ++k )
        for( j = 0; j < dims[1]; ++j )
            for( i = 0; i < dims[0]; ++i )
            {
                *x = coords[0][i];
                ++x;
                *y = coords[1][j];
                ++y;
                *z = coords[2][k];
                ++z;
            }

    return vtk_create_structured_elems( dims, start_handle, elem_list );
}

ErrorCode ReadVtk::vtk_read_polydata( FileTokenizer& tokens, Range& vertex_list, std::vector< Range >& elem_list )
{
    ErrorCode result;
    long num_verts;
    const char* const poly_data_names[] = { "VERTICES", "LINES", "POLYGONS", "TRIANGLE_STRIPS", 0 };

    if( !tokens.match_token( "POINTS" ) || !tokens.get_long_ints( 1, &num_verts ) ||
        !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( num_verts < 1 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid point count at line " << tokens.line_number() );
    }

    // Create vertices and read coordinates
    EntityHandle start_handle = 0;
    result                    = read_vertices( tokens, num_verts, start_handle );
    if( MB_SUCCESS != result ) return result;
    vertex_list.insert( start_handle, start_handle + num_verts - 1 );

    int poly_type = tokens.match_token( poly_data_names );
    switch( poly_type )
    {
        case 0:
            result = MB_FAILURE;
            break;
        case 1:
            MB_SET_ERR( MB_FAILURE, "Vertex element type at line " << tokens.line_number() );
            break;
        case 2:
            MB_SET_ERR( MB_FAILURE, "Unsupported type: polylines at line " << tokens.line_number() );
            result = MB_FAILURE;
            break;
        case 3:
            result = vtk_read_polygons( tokens, start_handle, elem_list );
            break;
        case 4:
            MB_SET_ERR( MB_FAILURE, "Unsupported type: triangle strips at line " << tokens.line_number() );
            result = MB_FAILURE;
            break;
    }

    return result;
}

ErrorCode ReadVtk::vtk_read_polygons( FileTokenizer& tokens, EntityHandle first_vtx, std::vector< Range >& elem_list )
{
    ErrorCode result;
    long size[2];

    if( !tokens.get_long_ints( 2, size ) || !tokens.get_newline() ) return MB_FAILURE;

    const Range empty;
    std::vector< EntityHandle > conn_hdl;
    std::vector< long > conn_idx;
    EntityHandle first = 0, prev = 0, handle;
    for( int i = 0; i < size[0]; ++i )
    {
        long count;
        if( !tokens.get_long_ints( 1, &count ) ) return MB_FAILURE;
        conn_idx.resize( count );
        conn_hdl.resize( count );
        if( !tokens.get_long_ints( count, &conn_idx[0] ) ) return MB_FAILURE;

        for( long j = 0; j < count; ++j )
            conn_hdl[j] = first_vtx + conn_idx[j];

        result = mdbImpl->create_element( MBPOLYGON, &conn_hdl[0], count, handle );
        if( MB_SUCCESS != result ) return result;

        if( prev + 1 != handle )
        {
            if( first )
            {  // True except for first iteration (first == 0)
                if( elem_list.empty() || first < elem_list.back().front() )  // Only need new range if order would get
                                                                             // mixed up, or we just began inserting
                    elem_list.push_back( empty );
                elem_list.back().insert( first, prev );
            }
            first = handle;
        }
        prev = handle;
    }
    if( first )
    {                                                                // True unless no elements (size[0] == 0)
        if( elem_list.empty() || first < elem_list.back().front() )  // Only need new range if order would get mixed
                                                                     // up, or we just began inserting
            elem_list.push_back( empty );
        elem_list.back().insert( first, prev );
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_read_unstructured_grid( FileTokenizer& tokens,
                                               Range& vertex_list,
                                               std::vector< Range >& elem_list )
{
    ErrorCode result;
    long i, num_verts, num_elems[2];
    EntityHandle tmp_conn_list[27];

    // Poorly formatted VTK legacy format document seems to
    // lead many to think that a FIELD block can occur within
    // an UNSTRUCTURED_GRID dataset rather than as its own data
    // set. So allow for field data between other blocks of
    // data.

    const char* pts_str[] = { "FIELD", "POINTS", 0 };
    while( 1 == ( i = tokens.match_token( pts_str ) ) )
    {
        result = vtk_read_field( tokens );
        if( MB_SUCCESS != result ) return result;
    }
    if( i != 2 ) return MB_FAILURE;

    if( !tokens.get_long_ints( 1, &num_verts ) || !tokens.match_token( vtk_type_names ) || !tokens.get_newline() )
        return MB_FAILURE;

    if( num_verts < 1 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid point count at line " << tokens.line_number() );
    }

    // Create vertices and read coordinates
    EntityHandle first_vertex = 0;
    result                    = read_vertices( tokens, num_verts, first_vertex );
    if( MB_SUCCESS != result ) return result;
    vertex_list.insert( first_vertex, first_vertex + num_verts - 1 );

    const char* cell_str[] = { "FIELD", "CELLS", 0 };
    while( 1 == ( i = tokens.match_token( cell_str ) ) )
    {
        result = vtk_read_field( tokens );
        if( MB_SUCCESS != result ) return result;
    }
    if( i != 2 ) return MB_FAILURE;

    if( !tokens.get_long_ints( 2, num_elems ) || !tokens.get_newline() ) return MB_FAILURE;

    // Read element connectivity for all elements
    std::vector< long > connectivity( num_elems[1] );
    if( !tokens.get_long_ints( num_elems[1], &connectivity[0] ) ) return MB_FAILURE;

    if( !tokens.match_token( "CELL_TYPES" ) || !tokens.get_long_ints( 1, &num_elems[1] ) || !tokens.get_newline() )
        return MB_FAILURE;

    // Read element types
    std::vector< long > types( num_elems[0] );
    if( !tokens.get_long_ints( num_elems[0], &types[0] ) ) return MB_FAILURE;

    // Create elements in blocks of the same type
    // It is important to preserve the order in
    // which the elements were read for later reading
    // attribute data.
    long id                                 = 0;
    std::vector< long >::iterator conn_iter = connectivity.begin();
    while( id < num_elems[0] )
    {
        unsigned vtk_type = types[id];
        if( vtk_type >= VtkUtil::numVtkElemType ) return MB_FAILURE;

        EntityType type = VtkUtil::vtkElemTypes[vtk_type].mb_type;
        if( type == MBMAXTYPE )
        {
            MB_SET_ERR( MB_FAILURE, "Unsupported VTK element type: " << VtkUtil::vtkElemTypes[vtk_type].name << " ("
                                                                     << vtk_type << ")" );
        }

        int num_vtx = *conn_iter;
        if( type != MBPOLYGON && type != MBPOLYHEDRON && num_vtx != (int)VtkUtil::vtkElemTypes[vtk_type].num_nodes )
        {
            MB_SET_ERR( MB_FAILURE, "Cell " << id << " is of type '" << VtkUtil::vtkElemTypes[vtk_type].name
                                            << "' but has " << num_vtx << " vertices" );
        }

        // Find any subsequent elements of the same type
        // if polyhedra, need to look at the number of faces to put in the same range
        std::vector< long >::iterator conn_iter2 = conn_iter + num_vtx + 1;
        long end_id                              = id + 1;
        if( MBPOLYHEDRON != type )
        {
            while( end_id < num_elems[0] && (unsigned)types[end_id] == vtk_type && *conn_iter2 == num_vtx )
            {
                ++end_id;
                conn_iter2 += num_vtx + 1;
            }
        }
        else
        {
            // advance only if next is polyhedron too, and if number of faces is the same
            int num_faces = conn_iter[1];
            while( end_id < num_elems[0] && (unsigned)types[end_id] == vtk_type && conn_iter2[1] == num_faces )
            {
                ++end_id;
                conn_iter2 += num_vtx + 1;
            }
            // num_vtx becomes in this case num_faces
            num_vtx = num_faces;  // for polyhedra, this is what we want
            // trigger vertex adjacency call
            Range firstFaces;
            mdbImpl->get_adjacencies( &first_vertex, 1, 2, false, firstFaces );
        }

        // Allocate element block
        long num_elem             = end_id - id;
        EntityHandle start_handle = 0;
        EntityHandle* conn_array;

        // if type is MBVERTEX, skip, we will not create elements with one vertex
        if( MBVERTEX == type )
        {
            id += num_elem;
            conn_iter += 2 * num_elem;  // skip 2 * number of 1-vertex elements
            continue;
        }

        result = allocate_elements( num_elem, num_vtx, type, start_handle, conn_array, elem_list );
        if( MB_SUCCESS != result ) return result;

        EntityHandle* conn_sav = conn_array;

        // Store element connectivity
        if( type != MBPOLYHEDRON )
        {
            for( ; id < end_id; ++id )
            {
                if( conn_iter == connectivity.end() )
                {
                    MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at cell " << id );
                }
                // Make sure connectivity length is correct.
                if( *conn_iter != num_vtx )
                {
                    MB_SET_ERR( MB_FAILURE, "Cell " << id << " is of type '" << VtkUtil::vtkElemTypes[vtk_type].name
                                                    << "' but has " << num_vtx << " vertices" );
                }
                ++conn_iter;

                for( i = 0; i < num_vtx; ++i, ++conn_iter )
                {
                    if( conn_iter == connectivity.end() )
                    {
                        MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at cell " << id );
                    }

                    conn_array[i] = *conn_iter + first_vertex;
                }

                const unsigned* order = VtkUtil::vtkElemTypes[vtk_type].node_order;
                if( order )
                {
                    assert( num_vtx * sizeof( EntityHandle ) <= sizeof( tmp_conn_list ) );
                    memcpy( tmp_conn_list, conn_array, num_vtx * sizeof( EntityHandle ) );
                    for( int j = 0; j < num_vtx; ++j )
                        conn_array[order[j]] = tmp_conn_list[j];
                }

                conn_array += num_vtx;
            }
        }
        else  // type == MBPOLYHEDRON
        {
            // in some cases, we may need to create new elements; will it screw the tags?
            // not if the file was not from moab
            ErrorCode rv = MB_SUCCESS;
            for( ; id < end_id; ++id )
            {
                if( conn_iter == connectivity.end() )
                {
                    MB_SET_ERR( MB_FAILURE, "Connectivity data truncated at polyhedra cell " << id );
                }
                ++conn_iter;
                // iterator is now at number of faces
                // we should check it is indeed num_vtx
                int num_faces = *conn_iter;
                if( num_faces != num_vtx ) MB_SET_ERR( MB_FAILURE, "Connectivity data wrong at polyhedra cell " << id );

                EntityHandle connec[20];  // we bet we will have only 20 vertices at most, in a
                                          // face in a polyhedra
                for( int j = 0; j < num_faces; j++ )
                {
                    conn_iter++;
                    int numverticesInFace = (int)*conn_iter;
                    if( numverticesInFace > 20 )
                        MB_SET_ERR( MB_FAILURE,
                                    "too many vertices in face index " << j << " for polyhedra cell " << id );
                    // need to find the face, but first fill with vertices
                    for( int k = 0; k < numverticesInFace; k++ )
                    {
                        connec[k] = first_vertex + *( ++conn_iter );  //
                    }
                    Range adjFaces;
                    // find a face with these vertices; if not, we need to create one, on the fly :(
                    rv = mdbImpl->get_adjacencies( connec, numverticesInFace, 2, false, adjFaces );MB_CHK_ERR( rv );
                    if( adjFaces.size() >= 1 )
                    {
                        conn_array[j] = adjFaces[0];  // get the first face found
                    }
                    else
                    {
                        // create the face; tri, quad or polygon
                        EntityType etype = MBTRI;
                        if( 4 == numverticesInFace ) etype = MBQUAD;
                        if( 4 < numverticesInFace ) etype = MBPOLYGON;

                        rv = mdbImpl->create_element( etype, connec, numverticesInFace, conn_array[j] );MB_CHK_ERR( rv );
                    }
                }

                conn_array += num_vtx;  // advance for next polyhedra
                conn_iter++;            // advance to the next field
            }
        }

        // Notify MOAB of the new elements
        result = readMeshIface->update_adjacencies( start_handle, num_elem, num_vtx, conn_sav );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_create_structured_elems( const long* dims,
                                                EntityHandle first_vtx,
                                                std::vector< Range >& elem_list )
{
    ErrorCode result;
    // int non_zero[3] = {0, 0, 0}; // True if dim > 0 for x, y, z respectively
    long elem_dim  = 0;           // Element dimension (2->quad, 3->hex)
    long num_elems = 1;           // Total number of elements
    long vert_per_elem;           // Element connectivity length
    long edims[3] = { 1, 1, 1 };  // Number of elements in each grid direction

    // Populate above data
    for( int d = 0; d < 3; d++ )
    {
        if( dims[d] > 1 )
        {
            // non_zero[elem_dim] = d;
            ++elem_dim;
            edims[d] = dims[d] - 1;
            num_elems *= edims[d];
        }
    }
    vert_per_elem = 1 << elem_dim;

    // Get element type from element dimension
    EntityType type;
    switch( elem_dim )
    {
        case 1:
            type = MBEDGE;
            break;
        case 2:
            type = MBQUAD;
            break;
        case 3:
            type = MBHEX;
            break;
        default:
            MB_SET_ERR( MB_FAILURE, "Invalid dimension for structured elements: " << elem_dim );
    }

    // Allocate storage for elements
    EntityHandle start_handle = 0;
    EntityHandle* conn_array;
    result = allocate_elements( num_elems, vert_per_elem, type, start_handle, conn_array, elem_list );
    if( MB_SUCCESS != result ) return MB_FAILURE;

    EntityHandle* conn_sav = conn_array;

    // Offsets of element vertices in grid relative to corner closest to origin
    long k                = dims[0] * dims[1];
    const long corners[8] = { 0, 1, 1 + dims[0], dims[0], k, k + 1, k + 1 + dims[0], k + dims[0] };

    // Populate element list
    for( long z = 0; z < edims[2]; ++z )
        for( long y = 0; y < edims[1]; ++y )
            for( long x = 0; x < edims[0]; ++x )
            {
                const long index = x + y * dims[0] + z * ( dims[0] * dims[1] );
                for( long j = 0; j < vert_per_elem; ++j, ++conn_array )
                    *conn_array = index + corners[j] + first_vtx;
            }

    // Notify MOAB of the new elements
    result = readMeshIface->update_adjacencies( start_handle, num_elems, vert_per_elem, conn_sav );
    if( MB_SUCCESS != result ) return result;

    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_read_field( FileTokenizer& tokens )
{
    // This is not supported yet.
    // Parse the data but throw it away because
    // Mesquite has no internal representation for it.

    // Could save this in tags, but the only useful thing that
    // could be done with the data is to write it back out
    // with the modified mesh. As there's no way to save the
    // type of a tag in Mesquite, it cannot be written back
    // out correctly either.
    // FIXME: Don't know what to do with this data.
    // For now, read it and throw it out.

    long num_arrays;
    if( !tokens.get_string() ||  // Name
        !tokens.get_long_ints( 1, &num_arrays ) )
        return MB_FAILURE;

    for( long i = 0; i < num_arrays; ++i )
    {
        /*const char* name =*/tokens.get_string();

        long dims[2];
        if( !tokens.get_long_ints( 2, dims ) || !tokens.match_token( vtk_type_names ) ) return MB_FAILURE;

        long num_vals = dims[0] * dims[1];

        for( long j = 0; j < num_vals; j++ )
        {
            double junk;
            if( !tokens.get_doubles( 1, &junk ) ) return MB_FAILURE;
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_read_attrib_data( FileTokenizer& tokens, std::vector< Range >& entities )
{
    const char* const type_names[] = { "SCALARS", "COLOR_SCALARS", "VECTORS", "NORMALS", "TEXTURE_COORDINATES",
                                       "TENSORS", "FIELD",         0 };

    int type             = tokens.match_token( type_names );
    const char* tmp_name = tokens.get_string();
    if( !type || !tmp_name ) return MB_FAILURE;

    std::string name_alloc( tmp_name );
    const char* name = name_alloc.c_str();
    switch( type )
    {
        case 1:
            return vtk_read_scalar_attrib( tokens, entities, name );
        case 2:
            return vtk_read_color_attrib( tokens, entities, name );
        case 3:
            return vtk_read_vector_attrib( tokens, entities, name );
        case 4:
            return vtk_read_vector_attrib( tokens, entities, name );
        case 5:
            return vtk_read_texture_attrib( tokens, entities, name );
        case 6:
            return vtk_read_tensor_attrib( tokens, entities, name );
        case 7:
            return vtk_read_field_attrib( tokens, entities, name );
    }

    return MB_FAILURE;
}

ErrorCode ReadVtk::vtk_read_tag_data( FileTokenizer& tokens,
                                      int type,
                                      size_t per_elem,
                                      std::vector< Range >& entities,
                                      const char* name )
{
    ErrorCode result;
    DataType mb_type;
    if( type == 1 )
    {
        mb_type = MB_TYPE_BIT;
    }
    else if( type >= 2 && type <= 9 )
    {
        mb_type = MB_TYPE_INTEGER;
    }
    else if( type == 10 || type == 11 )
    {
        mb_type = MB_TYPE_DOUBLE;
    }
    else if( type == 12 )
    {
        mb_type = MB_TYPE_INTEGER;
    }
    else
        return MB_FAILURE;

#ifdef MB_VTK_MATERIAL_SETS
    size_t size;
    if( type == 1 )
    {
        size = sizeof( bool );
    }
    else if( type >= 2 && type <= 9 )
    {
        size = sizeof( int );
    }
    else if( type == 10 || type == 11 )
    {
        size = sizeof( double );
    }
    else /* (type == 12) */
    {
        size = 4;  // Could be 4 or 8, but we don't know. Hope it's 4 because MOAB doesn't support
                   // 64-bit ints.
    }
    Modulator materialMap( this->mdbImpl );
    result = materialMap.initialize( this->mPartitionTagName, mb_type, size, per_elem );MB_CHK_SET_ERR( result, "MaterialMap tag (" << this->mPartitionTagName << ") creation failed." );
    bool isMaterial = size * per_elem <= 4 &&              // Must have int-sized values (ParallelComm requires it)
                      !this->mPartitionTagName.empty() &&  // Must have a non-empty field name...
                      !strcmp( name, this->mPartitionTagName.c_str() );  // ... that matches our spec.
#endif                                                                   // MB_VTK_MATERIAL_SETS

    // Get/create tag
    Tag handle;
    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() );

    std::vector< Range >::iterator iter;

    if( type == 1 )
    {
        for( iter = entities.begin(); iter != entities.end(); ++iter )
        {
            bool* data = new bool[iter->size() * per_elem];
            if( !tokens.get_booleans( per_elem * iter->size(), data ) )
            {
                delete[] data;
                return MB_FAILURE;
            }

            bool* data_iter          = data;
            Range::iterator ent_iter = iter->begin();
            for( ; ent_iter != iter->end(); ++ent_iter )
            {
                unsigned char bits = 0;
                for( unsigned j = 0; j < per_elem; ++j, ++data_iter )
                    bits |= (unsigned char)( *data_iter << j );
#ifdef MB_VTK_MATERIAL_SETS
                if( isMaterial ) materialMap.add_entity( *ent_iter, &bits, 1 );
#endif  // MB_VTK_MATERIAL_SETS
                result = mdbImpl->tag_set_data( handle, &*ent_iter, 1, &bits );
                if( MB_SUCCESS != result )
                {
                    delete[] data;
                    return result;
                }
            }
            delete[] data;
        }
    }
    else if( ( type >= 2 && type <= 9 ) || type == 12 )
    {
        std::vector< int > data;
        for( iter = entities.begin(); iter != entities.end(); ++iter )
        {
            data.resize( iter->size() * per_elem );
            if( !tokens.get_integers( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
#ifdef MB_VTK_MATERIAL_SETS
            if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
#endif  // MB_VTK_MATERIAL_SETS
            result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
            if( MB_SUCCESS != result ) return result;
        }
    }
    else if( type == 10 || type == 11 )
    {
        std::vector< double > data;
        for( iter = entities.begin(); iter != entities.end(); ++iter )
        {
            data.resize( iter->size() * per_elem );
            if( !tokens.get_doubles( iter->size() * per_elem, &data[0] ) ) return MB_FAILURE;
#ifdef MB_VTK_MATERIAL_SETS
            if( isMaterial ) materialMap.add_entities( *iter, (unsigned char*)&data[0], per_elem * size );
#endif  // MB_VTK_MATERIAL_SETS
            result = mdbImpl->tag_set_data( handle, *iter, &data[0] );
            if( MB_SUCCESS != result ) return result;
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::vtk_read_scalar_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
{
    int type = tokens.match_token( vtk_type_names );
    if( !type ) return MB_FAILURE;

    long size;
    const char* tok = tokens.get_string();
    if( !tok ) return MB_FAILURE;

    const char* end = 0;
    size            = strtol( tok, (char**)&end, 0 );
    if( *end )
    {
        size = 1;
        tokens.unget_token();
    }

    // VTK spec says cannot be greater than 4--do we care?
    if( size < 1 )
    {  //|| size > 4)
        MB_SET_ERR( MB_FAILURE, "Scalar count out of range [1,4] at line " << tokens.line_number() );
    }

    if( !tokens.match_token( "LOOKUP_TABLE" ) || !tokens.match_token( "default" ) ) return MB_FAILURE;

    return vtk_read_tag_data( tokens, type, size, entities, name );
}

ErrorCode ReadVtk::vtk_read_color_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
{
    long size;
    if( !tokens.get_long_ints( 1, &size ) || size < 1 ) return MB_FAILURE;

    return vtk_read_tag_data( tokens, 10, size, entities, name );
}

ErrorCode ReadVtk::vtk_read_vector_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
{
    int type = tokens.match_token( vtk_type_names );
    if( !type ) return MB_FAILURE;

    return vtk_read_tag_data( tokens, type, 3, entities, name );
}

ErrorCode ReadVtk::vtk_read_texture_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
{
    int type, dim;
    if( !tokens.get_integers( 1, &dim ) || !( type = tokens.match_token( vtk_type_names ) ) ) return MB_FAILURE;

    if( dim < 1 || dim > 3 )
    {
        MB_SET_ERR( MB_FAILURE, "Invalid dimension (" << dim << ") at line " << tokens.line_number() );
    }

    return vtk_read_tag_data( tokens, type, dim, entities, name );
}

ErrorCode ReadVtk::vtk_read_tensor_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* name )
{
    int type = tokens.match_token( vtk_type_names );
    if( !type ) return MB_FAILURE;

    return vtk_read_tag_data( tokens, type, 9, entities, name );
}

ErrorCode ReadVtk::vtk_read_field_attrib( FileTokenizer& tokens, std::vector< Range >& entities, const char* )
{
    long num_fields;
    if( !tokens.get_long_ints( 1, &num_fields ) ) return MB_FAILURE;

    long i;
    for( i = 0; i < num_fields; ++i )
    {
        const char* tok = tokens.get_string();
        if( !tok ) return MB_FAILURE;

        std::string name_alloc( tok );

        long num_comp;
        if( !tokens.get_long_ints( 1, &num_comp ) ) return MB_FAILURE;

        long num_tuples;
        if( !tokens.get_long_ints( 1, &num_tuples ) ) return MB_FAILURE;

        int type = tokens.match_token( vtk_type_names );
        if( !type ) return MB_FAILURE;

        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, "
                                                                  << num_tuples << " tuples, type " << type
                                                                  << ") at line " << tokens.line_number() );
    }

    return MB_SUCCESS;
}

ErrorCode ReadVtk::store_file_ids( Tag tag, const Range& verts, const std::vector< Range >& elems )
{
    ErrorCode rval;

    rval = readMeshIface->assign_ids( tag, verts );
    if( MB_SUCCESS != rval ) return rval;

    int id = 0;
    for( size_t i = 0; i < elems.size(); ++i )
    {
        rval = readMeshIface->assign_ids( tag, elems[i], id );
        id += elems[i].size();
    }

    return MB_SUCCESS;
}

}  // namespace moab