[vtkusers] VTK: Field data I/O with Exodus files

[Nico Schlömer]

> A non-text attachment was scrubbed...

Oh well, here you go
================== *snip* ==================
# start of file
# -*- coding: utf-8 -*-
# # ===================================================
def _main():
    filetype = 'vtu'
    #filetype = 'exodus'

    # create some points and cells
    points = [[0.0, 0.0, 0.0],
              [1.0, 0.0, 0.0],
              [1.0, 1.0, 0.0],
              [0.0, 1.0, 0.0]]
    cells = [[0,1,3],
             [1,2,3]]

    # create vtk mesh
    vtk_mesh = _generate_vtk_mesh(points, cells)

    # add field data mu=0.5 to the mesh
    from vtk import vtkDoubleArray
    scalars0 = vtkDoubleArray()
    scalars0.SetName('mu')
    scalars0.SetNumberOfComponents( 1 )
    scalars0.InsertNextValue( 0.5 )
    vtk_mesh.GetFieldData().AddArray( scalars0 )

    # instantiate the writer
    if filetype == 'vtu':
        filename = 'test.vtu'
        from vtk import vtkXMLUnstructuredGridWriter
        writer = vtkXMLUnstructuredGridWriter()
    elif filetype == 'exodus':
        filename = 'test.e'
        from vtk import vtkExodusIIWriter
        writer = vtkExodusIIWriter()
        writer.WriteAllTimeStepsOn()
    else:
        raise ValueError('File type error.')

    # write it
    writer.SetFileName( filename )
    writer.SetInput( vtk_mesh )
    writer.Write()

    # aaaand read it again
    if filetype == 'vtu':
        from vtk import vtkXMLUnstructuredGridReader
        reader = vtkXMLUnstructuredGridReader()
        reader.SetFileName(filename)
        reader.Update()
        vtk_mesh = reader.GetOutput()
    elif filetype == 'exodus':
        from vtk import vtkExodusIIReader
        reader = vtkExodusIIReader()
        vtk_mesh = _read_exodusii_mesh(reader, filename)
    else:
        raise ValueError('File type error.')

    # read the field data
    field_data = _read_field_data( vtk_mesh )

    # expected output: {'mu': 0.5}
    print field_data

    return
# ===================================================
def _generate_vtk_mesh(points, cellsNodes):
    from vtk import vtkUnstructuredGrid, VTK_TRIANGLE, VTK_TETRA,
vtkIdList, vtkPoints
    mesh = vtkUnstructuredGrid()

    # set points
    vtk_points = vtkPoints()
    for point in points:
        vtk_points.InsertNextPoint(point[0], point[1], point[2])
    mesh.SetPoints( vtk_points )

    # set cells
    for cellNodes in cellsNodes:
        pts = vtkIdList()
        num_local_nodes = len(cellNodes)
        pts.SetNumberOfIds(num_local_nodes)
        # get the connectivity for this element
        for k, node_index in enumerate(cellNodes):
            pts.InsertId(k, node_index)
        mesh.InsertNextCell(VTK_TRIANGLE, pts)

    return mesh
# ===================================================
def _read_exodusii_mesh( reader, file_name ):
    '''Uses a vtkExodusIIReader to return a vtkUnstructuredGrid.
    '''
    reader.SetFileName( file_name )

    # Fetch metadata.
    reader.UpdateInformation()

    # Make sure the point fields are read during Update().
    for k in xrange( reader.GetNumberOfPointResultArrays() ):
        arr_name = reader.GetPointResultArrayName( k )
        reader.SetPointResultArrayStatus( arr_name, 1 )

    # Read the file.
    reader.Update()
    out = reader.GetOutput()

    # Loop through the blocks and search for a vtkUnstructuredGrid.
    vtk_mesh = []
    for i in xrange( out.GetNumberOfBlocks() ):
        blk = out.GetBlock( i )
        for j in xrange( blk.GetNumberOfBlocks() ):
            sub_block = blk.GetBlock( j )
            if sub_block.IsA( 'vtkUnstructuredGrid' ):
                vtk_mesh.append( sub_block )

    if len(vtk_mesh) == 0:
        raise IOError( 'No \'vtkUnstructuredGrid\' found!' )
    elif len(vtk_mesh) > 1:
        raise IOError( 'More than one \'vtkUnstructuredGrid\' found!' )

    return vtk_mesh[0]
# ===================================================
def _read_field_data( vtk_data ):
    '''Gather field data.
    '''
    vtk_field_data = vtk_data.GetFieldData()
    num_arrays = vtk_field_data.GetNumberOfArrays()

    field_data = {}
    for k in xrange( num_arrays ):
        array  = vtk_field_data.GetArray(k)
        name   = array.GetName()
        num_values = array.GetDataSize()
        if num_values == 1:
            values = array.GetValue( k )
        else:
            values = np.zeros( num_values )
            for i in xrange( num_values ):
                values[i] = array.GetValue(i)
        field_data[ name ] = values

    return field_data
# ===================================================
if __name__ == '__main__':
   _main()
# ===================================================
# eof
================== *snap* ==================

--Nico


On Tue, Mar 13, 2012 at 2:01 PM, Nico Schlömer <nico.schloemer at gmail.com> wrote:
> Hi,
>
> I'm having trouble reading field data from Exodus files using VTK.
>
> When I do field data I/O with VTU files, all seems fine; the attached
> Python code writes mu=0.5 to a VTU file and then reads it out
> correctly.
> When I do the same thing with exodus files (set filetype='exodus' in
> the attached example code), the field data is not retrieved correctly.
> Judging from what I see in ParaView, it seems to be written out
> alright -- at least a field with the name 'mu' is displayed.
>
> Any hints?
>
> If you notice that I'm using something else unlike intended, let me
> know of course.
>
> Cheers,
> Nico