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

[David Thompson]

Nico,

Also, in _read_field_data(), this

     values = array.GetValue( k )

should be

     values = array.GetValue( 0 )

	David

On Mar 13, 2012, at 13:25 , David Thompson wrote:

> 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
>>
>
>
> _______________________________________________
> 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
>