I have run a DEM model in mfix and now want to save the particle data for each time step.
I have seen that I can do this through importing the .pvd file into paraview and see the spreadsheet of data.
However, I want to save all the velocity/coordinates of each particle for each time step. I see I have a series of .vtp files containing this information. Do i have to individually put these into paraview and save the data, or is there a way to save all the data at once?
I suggest just using a vtp reader instead going through paraview. Since you already have vtk installed in the MFiX environment, just use Python. To get you started:
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)
Thank you for your response.
Is there anything in the user guide about the vtk/vtp reader, as I could only find information on how to use paraview or how to build postmfix to output data files?
I am new to coding in python, so was hoping that a full set of instructions on how to do this may be available.
Thanks
no we do not have any documentation included with MFIX. VTK is the C++ library behind Paraview. The documentation is here: VTK: VTK 9.0.20201226 Documentation
There are a collection of examples (C++, Python, Java, C#, js) here: https://lorensen.github.io/VTKExamples/site/
I have been trying to use python to get the data stored in the vtp files as you suggested but whenever I use vtk_to_numpy I get the error:
array = vtk_to_numpy(points)
AttributeError: ‘NoneType’ object has no attribute ‘GetDataType’
Could you provide any information on how to fix this?
looks like “points” is None, what is the rest of your code?
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)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.
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())how big are your files? can you post one?
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
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
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
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
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')
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
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.
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
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!
Sorry. You can replace those “data classes” with regular dictionaries. I edited the code above.