Gradient of Unstructured Grid in Python

Support
1

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!

1 Like
2

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)}
4

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

7

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()
8

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.

9

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

10

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

12

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.

13

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

14

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:

  1. How to generate polygon mesh based on the introduced coordinates.
  2. 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)

15
  1. 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

  1. 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"]
16

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?

17

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?

18

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.