Saving data as .csv files

MFiX
#7

I used the same code as you provided above:

import vtk
from vtk.util.numpy_support  import vtk_to_numpy

# File to read
fname = 'DES_FB1_DES_00019.vtp'

# --- read a vtp file ---
points = vtk.vtkXMLPolyDataReader()
points.SetFileName(fname)
points.Update()

# print the arrays out
data = points.GetOutput()
point_data = data.GetPointData()
for i in range(point_data.GetNumberOfArrays()):
    print(i, point_data.GetArrayName(i), 'Components:', point_data.GetArray(i).GetNumberOfComponents())

points = point_data.GetAbstractArray('array name')
array = vtk_to_numpy(points)
#8

you need to change the 'array name' in this line points = point_data.GetAbstractArray('array name') to one of the names that is printed out.

#9

Sorry on my version I changed this, but I pasted the original above.
In my version it reads:

vel = vtk_to_numpy(point_data.GetAbstractArray('Velocity'))
points = vtk_to_numpy(data.GetPoints().GetData())
#10

how big are your files? can you post one?

#11

Hi Justin,

I actually had a similar question. Could you help me with it? Thanks!

The website does not allow me to post vtu file, but I got my vtu using this mfx file. vtu.mfx (13.6 KB)

Then I used the following code to extract data

import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy

reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(“BACKGROUND_IC_0005.vtu”)
reader.Update()
data = reader.GetOutput()

point_data = data.GetPointData()
for i in range(point_data.GetNumberOfArrays()):
print(i, point_data.GetArrayName(i), ‘Components:’, point_data.GetArray(i).GetNumberOfComponents())

But I didn’t see anything printing out. Did I misunderstand somewhere on how to use the code?

Thanks,
Jack

#12

There’s no point data in that file (it contains cell data only). It’s useful to open an interactive Python session for data exploration:

% python
>>> import vtk
>>> reader = vtk.vtkXMLUnstructuredGridReader()
>>> reader.SetFileName("BACKGROUND_IC_0005.vtu")  
>>> reader.Update()
>>> data = reader.GetOutput()
>>> data.GetPointData().GetNumberOfArrays()
0
>>> data.GetCellData().GetNumberOfArrays()
1

Also note that although you cannot post VTU files, you can post zipfiles, so just use ‘zip’ to create a zipfile if you need to upload data files.

– Charles

How to read the vtu files with python?
#13

Thanks a lot Charles! I see where the problem is now. I got the cell array data using triangles=vtk_to_numpy(data.GetCells().GetData()). However, I am very new to vtk and still got really confused by the result. Could you give me some hint on what the result mean? I noticed that the my result of vtk_to_numpy(data.GetCells().GetData()) was a series of numbers with many 4s. Each 4 was followed by 4 integers. What information is it providing? Thank you!

I checked out the kitware tutorials that you sent to me, but I didn’t find how GetCells() or GetPoints() work. What are the points and cells? In my understanding, point data stores information about where the points are. Namely it stores arrays of tuples that represent coordinates of points. Then cell data can be parsed into tuples of integers, and each integer represents the index of a point in the point data. Did I misunderstand how vtk works?

Thanks a lot,
Jack

#14

Jack -

You are calling vtk_to_numpy on the wrong object. You need to call GetCellData, not GetCells.

Try this:

import vtk
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName("BACKGROUND_IC_0005.vtu")
reader.Update()
output = reader.GetOutput()
data = output.GetCellData()
print(data.GetNumberOfArrays())
arr = data.GetArray(0)
print(arr)
print(arr.GetTuple(0))
print(vtk_to_numpy(arr))

these values are just the x-velocity components that you’re writing in the VTK file. If you had selected more variables the tuples would contain multiple values. CORRECTION: each variable is stored in its own array.

– Charles

1 Like
How to plot time average particle velocities along x-direction
How to have access to all cell data?
#15

This is how this histogram viewer works in the gui, if it is a vtu file, use read_griddata, if it is a vtp file, use read_partdata. The functions will return a dict of variables names with some other info including the actual array, data. Then just use the vtk_to_numpy function on that data array.

import vtk

def read_griddata(path):
    ugrid_reader = (
        vtk.vtkXMLPUnstructuredGridReader()
        if path.suffix == "pvtu"
        else vtk.vtkXMLUnstructuredGridReader()
    )

    ugrid_reader.SetFileName(path)
    ugrid_reader.Update()
    data = ugrid_reader.GetOutput()

    cell_data = data.GetCellData()
    array_info = {}
    for i in range(cell_data.GetNumberOfArrays()):
        array = cell_data.GetArray(i)
        n_comp = array.GetNumberOfComponents()
        name = cell_data.GetArrayName(i)
        array_info[name] = {
            'magnitude':  array.GetRange(-1),
            'components': n_comp,
            'range':      [array.GetRange(i) for i in range(n_comp)],
            'data':       cell_data.GetAbstractArray(name),
        }
    return array_info

def read_partdata(path):
    particle_reader = (
        vtk.vtkXMLPPolyDataReader()
        if path.suffix == ".pvtp"
        else vtk.vtkXMLPolyDataReader()
    )
    particle_reader.SetFileName(path)
    particle_reader.Update()

    data = particle_reader.GetOutput()
    point_data = data.GetPointData()
    array_info = {}
    n_tuples = None
    for i in range(point_data.GetNumberOfArrays()):
        array = point_data.GetArray(i)
        n_comp = array.GetNumberOfComponents()
        n_tuples = array.GetNumberOfTuples()

        name = point_data.GetArrayName(i)
        array_info[name] = {
            'number_of_tuples': n_tuples,
            'components':       n_comp,
            'range':            [array.GetRange(i) for i in range(n_comp)],
            'magnitude':        array.GetRange(-1),
            'data':             point_data.GetAbstractArray(name),
        }
    return array_info

if __name__ == "__main__":
    vtp_arrays = read_partdata('/path/to/output.vtp')
1 Like
#16

Yeah that works. I can retrieve the data now. Thanks a lot Charles! One more follow up question though. If I selected multiple variables, how can I figure out which value corresponds to which variable in each tuple?

Thanks,
Jack

#17

I was mistaken when I said VTK would store a tuple. Each variable goes into its own array, and every array has a name.

Use GetArrayName:

>>> print(data.GetNumberOfArrays())
3
>>> print(data.GetArrayName(0))
U_G

Note that this is contained in the code Justin posted above.

Also if you’re trying to learn VTK, at an interactive Python prompt you can type
>>> help(data)

(where data is a VTK data object) you will get a list of all the methods available with a description of each. Or consult a VTK manual.

#18

I see. That works now. Thanks Charles! I was trying out Justin’s code and it gave me an error saying NameError: name ‘GridArray’ is not defined. I guess I was missing some library. But your solution works now.

OK cool I will use the help() function.

Thank you,
Jack

#19

Yes, GridArray and ParticleArray are not part of VTK, they are defined elsewhere. The code Justin posted was incomplete, just intended as a code sample. Sorry for the confusion!

#20

Sorry. You can replace those “data classes” with regular dictionaries. I edited the code above.

2 Likes
#21

Hi, Justin, May I ask a question? I have used the code obtaining the information about “Diameter” , but how can I obtain the coordinate of each particle from vtp files? The coordinate is not a component name.

#22

you could store the particle x y z coordinates as 3 distinct DES_USR_VAR that can thus be saved as vtk along with velocity, diameter etc

#23

Thanks very much. But in this case, I will recalculate all the cases.

#24

The points (particle locations) are stored in the vtp files differently then the variable arrays. You can access the points like this:

import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy

particle_reader= vtk.vtkXMLPolyDataReader()
particle_reader.SetFileName(fname)
particle_reader.Update()
data = particle_reader.GetOutput()

# get the positions!
particle_positions = vtk_to_numpy(data.GetPoints().GetData())

FYI, I learned most of this from reading the VTK docs and the examples.

#25

Thanks very much for your valuable instructions. It really helps me.

#26

@onlyjus How can I get the velocity of the particles? I can get the positions of the particles by the piece of code you provided, but I am not able to get the velocity?