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

[David Thompson]

Hi Nico,

Adding this line to _read_exodusii_mesh() makes things work for me:

     reader.SetGlobalResultArrayStatus( 'mu', 1 )

The Exodus format calls field data "global" results (since they are  
not defined over points or cells) and they are not read by default.

	David

On Mar 13, 2012, at 8:32 , Nico Schlömer wrote:

>> 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
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>