I have a Python code that reads and analyses a VTU file using the vtk library for python. The file has a vector field stored in it. I would like to calculate the divergence of this vector field. I know Paraview can do this and I know how we can do this using the GUI, However, I would like to know if there is a Python API for Paraview which I can use for this.
I would like to achieve the following, I pass the VTU file and the corresponding vector field to a function, which will return the calculated divergence as an array.
Hi Mathieu, Thanks for the superfast response. I’m afraid I did not put my question in the right words, I’m looking for a possible paraview API that I can integrate into my code. Which will calculate the divergence and returns the value as an array.
Hi, Thank you, Math for the response. I found that the vtkGradientFilter of the VTK library can do what I want to achieve. For anyone looking at this in future,
The VTK library can be installed by pip install vtk or , pip3 install vtk . Inside this library, there is a function called vtkGardientFilter, Which takes a field as input and gives the calculated 9-dimensional gradient as an output array.
If you add the 0, 4 and 8th components in this output array you can get the divergence.
Mathieu, correct me if I made any mistakes.
Hi sorry, I did not know that the vtk libarry had such advanced routines. So I was in vein trying to find a way how I can use the Paraview filters instead. If someone in future made the same mistake as me, I thought this might be helpful for them.
from vtk.numpy_interface import algorithms as algs
div_f = algs.divergence(vec_field)
From the notes “Gradient and algorithms that require access to a mesh work whether that mesh is a uniform grid, a curvilinear grid, or an unstructured grid thanks to VTK’s data model.”
Hi, Thank you for your response. However, This is not working.
import numpy
import vtk
def GetField(vtu_filename, fieldname):
""" This will extract the field with the given name from the specified file """
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(vtu_filename)
reader.Update()
field_array = reader.GetOutput().GetPointData().GetArray(fieldname)
return field_array
def ComputeDivergence(input_array):
gradient_9d = vtk.vtkGradientFilter(input_array)
du_dx, dv_dy, dw_dz = gradient_9d[0], gradient_9d[4], gradient_9d[8]
return du_dx + dv_dy + dw_dz
from vtk.numpy_interface import algorithms as algs
Myfile = "hys_Ldiamond080.vtu"
Magnetization = GetField(Myfile, "Magnetization")
div_f = algs.divergence(Magnetization)
I get the following error.
Traceback (most recent call last):
File "Stack.py", line 21, in <module>
div_f = algs.divergence(Magnetization)
File "/home/rajgourav/.local/lib/python2.7/site-packages/vtk/numpy_interface/algorithms.py", line 125, in new_dsfunc
return dsfunc(array, ds)
File "/home/rajgourav/.local/lib/python2.7/site-packages/vtk/numpy_interface/internal_algorithms.py", line 247, in divergence
if not dataset : dataset = narray.DataSet
AttributeError: 'vtkCommonCorePython.vtkDoubleArray' object has no attribute 'DataSet'
You will need a bit more reading of the documentation page. You need to wrap your VTKArray into an object suitable for numpy processing. Thus, the following code should work for your case:
from vtk.numpy_interface import dataset_adapter as dsa
obj = dsa.WrapDataObject(reader.GetOutput())
Magnetization = obj.PointData[“Magnetization”] #we now have a numpy array
isinstance(Magnetization, np.ndarray)
my_div = algs.divergence(Magnetization)
Get back the VTK array equivalent of this numpy array
Thanks, the documentation was really cryptic and hard to understand. Your answer is really helpful, I have tried this and this is working. However, I still have a small technical doubt. In order to calculate the Divergence, the function should know all the mesh information. In this case, how does its NumPy function getting this information? Since we are only passing the vector field as a Numpy array? or does this also contain the point data?
Functions like this [gradient] require access to the dataset containing the array and the associated mesh. This is one of the reasons why we use a subclass of ndarray in dataset_adapter: Each array points to the dataset containing it. Functions such as gradient use the mesh and the array together.