I am looking to calculate the velocity gradients (du/dx, du/dy, dv,dx, dv/dy) at each point in an Unstructured Grid using Python. I have arrays for the x and y coordinate points as well as the x and y velocities at each point. However, I am not sure which VTK function to use for an unstructured grid. I have tried using the Numpy gradient function which is unsuccessful. Any help is appreciated!

Youâll want to use the `vtkGradientFilter`

filter:

vtkGradientFilter - A general filter for gradient estimation.

Superclass: vtkDataSetAlgorithm

Estimates the gradient of a field in a data set. The gradient

calculation is dependent on the input dataset type. The created

gradient array is of the same type as the array it is calculated from

(e.g. point data or cell data) but the data type will be either float

or double. At the boundary the gradient is not central differencing.

The output gradient array has 3*number of components of the input

data array. The ordering for the output gradient tuple will be

{du/dx, du/dy, du/dz, dv/dx, dv/dy, dv/dz, dw/dx, dw/dy, dw/dz} for

an input array {u, v, w}.

If in Python, PyVista provides a helper function for this filter.

```
from pyvista import examples
# A vtkStructuredGrid - but could be any mesh type
mesh = examples.download_kitchen()
mesh_g = mesh.compute_gradient(scalars="velocity")
mesh_g["gradient"]
```

`mesh_g["gradient"]`

is an N by 9 NumPy array of the gradients, so we could make a dictionary of numpy arrays of the gradients

```
>>> keys = "du/dx", "du/dy", "du/dz", "dv/dx", "dv/dy", "dv/dz", "dw/dx", "dw/dy", "dw/dz"
>>> gradients = dict(zip(keys, mesh_g["gradient"].T))
>>> gradients
{'du/dx': array([ 0.68209255, 0.14414474, -0.08285546, ..., -0.01418129,
-0.02921473, -0.42019978], dtype=float32),
'du/dy': array([-0.09281642, 0.17003751, 0.16026665, ..., -0.3221813 ,
-0.17364405, 0. ], dtype=float32),
'du/dz': array([ 0.17440438, -0.16565122, 0.01180762, ..., -1.3813452 ,
-0.04220725, 0. ], dtype=float32),
'dv/dx': array([-0.10317238, 0.04822876, 0.05100548, ..., 0. ,
0. , 0. ], dtype=float32),
'dv/dy': array([ 0.63790584, 0.86880773, 0.39256015, ..., -0.7274904 ,
-1.1117566 , -0.9115953 ], dtype=float32),
'dv/dz': array([0.17219502, 0.2421124 , 0.14371887, ..., 0. , 0. ,
0. ], dtype=float32),
'dw/dx': array([-0.2244322 , -0.08625986, -0.02245897, ..., 0. ,
0. , 0. ], dtype=float32),
'dw/dy': array([-0.21746539, -0.5370083 , -0.13856363, ..., 0. ,
0. , 0. ], dtype=float32),
'dw/dz': array([-0.15360962, -0.2934037 , -0.2604678 , ..., 0.8079373 ,
1.043386 , 1.3317952 ], dtype=float32)}
```

Can you post the file you are reading? And Iâll set up the vtk pipeline for you.

In brief, you need to pass the vtk mesh to the filter with the arrays to compute the gradient of as point data

Awesome. So I presume you want to compute the gradient of the âvelocityâ scalar field?

I did this two ways. The first is with PyVista and the second is with pure VTK code. PyVista is just a layer on top of VTK to make things a little easier.

Hope this helps!

```
import pyvista as pv
mesh = pv.read("soln_volume_CFD_00290.vtk")
print(mesh.array_names)
# Compute velocity
mesh["Velocity"] = mesh["Momentum"] / mesh["Density"][:,None]
# Compute gradient of `Velocity` array
mesh_grad = mesh.compute_gradient(scalars="Velocity")
# Save out!
mesh_grad.save("Output.vtk")
```

and with just VTK

```
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from vtk.util import numpy_support
from vtk.util.numpy_support import numpy_to_vtk
from operator import truediv
# File reading into a mesh data structure
filename = 'soln_volume_CFD_00290.vtk'
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename)
reader.ReadAllScalarsOn()
reader.ReadAllVectorsOn()
reader.ReadAllTensorsOn()
reader.Update()
data = reader.GetOutput() # the mesh
# Fetch arrays
density_vtk_array = data.GetPointData().GetArray(0)
momentum_vtk_array = data.GetPointData().GetArray(1)
density_numpy_array = np.array(density_vtk_array)
momentum_numpy_array = np.array(momentum_vtk_array)
# Make velocity array
res = list(map(truediv, momentum_numpy_array, density_numpy_array))
velocity = np.array(res)
velocity_vtk_arr = numpy_to_vtk(velocity)
velocity_vtk_arr.SetName("Velocity")
# Add velocity array to the mesh
data.GetPointData().AddArray(velocity_vtk_arr)
# Input mesh with velocity to the gradient filter
grad = vtk.vtkGradientFilter()
grad.SetInputData(data)
grad.SetInputScalars(0, "Velocity") # Zero b/c point data
grad.SetResultArrayName("gradient")
grad.Update()
# fetched the new mesh with the gradient calculated
computed = grad.GetOutput()
# Save out
writer = vtk.vtkDataSetWriter()
writer.SetFileName("Output.vtk")
writer.SetInputData(computed)
writer.Update()
writer.Write()
```

Perfect thanks a lot! Iâve visualised it in Paraview and it works.

My final question is how to extract each of the components e.g. du/dx, du/dy etc as I need to manipulate them for a Q-Criterion calculation. They appear in Paraview numbered 0-9 corresponding to each gradient.

If you convert that gradients array to a NumPy array then itâs an N (number of points) by 9 array where each column is a gradient so you can extract each column and do whatever

FYI: if youâre already using ParaView, you can run the gradient filter in ParaView via the GUI menu

How is it possible to adapt this procedure to compute gradient of a scalar field when corresponding cell center coordinates are available in .txt file. In fact, the input txt file contains (x,y,T) as columns, and gradients in output file (dT/dx, dT/dy) are required. I highly appreciate your response.

Youâll need to load the data from the text file (using Pandas for instance), pass the X and Y coordinates along with a constant Z coordinate to a `vtkPolyData`

mesh, add the T values as scalars, and then run the filter

Thank you so much for your response. I followed the recommended procedure. But I got trouble generating mesh by vtkPolyData. It would help a lot if I have your feedback on my two questions:

- How to generate polygon mesh based on the introduced coordinates.
- Is the way I assigned scalar temperature values to the generated mesh correct?

The script along with sample 3D input file are attached.

untitled11.py (1.4 KB)

snp-4264.txt (143.7 KB)

- How to generate polygon mesh based on the introduced coordinates.

Youâll need to triangulate the coordinates you have since you do not have the original cell connectivity

- Is the way I assigned scalar temperature values to the generated mesh correct?

The script along with sample 3D input file are attached.

Iâm not sureâŚ here is a simplified snippet using pandas to read the txt file.

```
import pyvista as pv
import numpy as np
import pandas as pd
df = pd.read_csv("snp-4264.txt", sep=",", skipinitialspace=True)
xyz = df[["x-coordinate", "y-coordinate", "z-coordinate"]].values
velocity = df[["x-velocity", "y-velocity", "z-velocity"]].values
temp = df["total-temperature"].values
vol = df["cell-volume"].values
mesh = pv.PolyData(xyz)
mesh["velocity"] = velocity
mesh["total-temperature"] = temp
mesh["cell-volume"] = vol
mesh_grad = mesh.compute_gradient(scalars='total-temperature')
grad = mesh_grad["gradient"]
```

Thanks, but still the obtained solution is not desirable (zero gradient is calculated). It seems the problem is in mesh generation. In your proposed script, if instead of âpv.PolyData(xyz)â, the following command is used:

point_cloud = pv.PolyData(xyz)

mesh = point_cloud.delaunay_3d(alpha=1)

the calculated gradient field is nonzero, but still it is not correct. How can I use pv.PolyData(xyz) so as not to obtain zero gradient as solution?

Hi,

I have a very similar task and followed the steps above, but I get the error: ââPolyDataâ object has no attribute âcompute_gradientââ. I also get this error when just running the âdownload_kitchenâ example above. Was the gradient function removed?

Hi All,

I am an active researcher in the field of ocean science. This time I stuck at one of problem related to gradient issue. I need to calculate pressure gradient in the x and y direction on a skew grid (unstructured grid). I tried may way from Matlab and python from last more than one months. I read the previous discussion about it but I am not clearly understand how can i implement it on my actual data. Can any one please help me to resolve this issue. If any one need my data I can share also.