Plot VTK with matplotlib

plot_vtk.py

#!/usr/bin/env python

import os
import numpy as np
import vtk
import matplotlib.pyplot as plt


def load_velocity(filename):
    if not os.path.exists(filename):
        return None
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.ReadAllVectorsOn()
    reader.Update()

    data = reader.GetOutput()

    # Extracting triangulation information
    triangles = data.GetPolys().GetData()
    points = data.GetPoints()

    # Mapping data: cell -> point
    mapper = vtk.vtkCellDataToPointData()
    mapper.AddInputData(data)
    mapper.Update()
    mapped_data = mapper.GetOutput()

    # Extracting interpolate point data
    udata = mapped_data.GetPointData().GetArray(0)

    ntri = triangles.GetNumberOfTuples()/4
    npts = points.GetNumberOfPoints()
    nvls = udata.GetNumberOfTuples()

    tri = np.zeros((ntri, 3))
    x = np.zeros(npts)
    y = np.zeros(npts)
    ux = np.zeros(nvls)
    uy = np.zeros(nvls)

    for i in xrange(0, ntri):
        tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
        tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
        tri[i, 2] = triangles.GetTuple(4*i + 3)[0]

    for i in xrange(npts):
        pt = points.GetPoint(i)
        x[i] = pt[0]
        y[i] = pt[1]

    for i in xrange(0, nvls):
        U = udata.GetTuple(i)
        ux[i] = U[0]
        uy[i] = U[1]

    return (x, y, tri, ux, uy)

plt.clf()
x, y, tri, ux, uy = load_velocity('U_test-plane.vtk')
plt.tricontour(x, y, tri, ux, 16, linestyles='-',
               colors='black', linewidths=0.5)
plt.tricontourf(x, y, tri, ux, 16)
plt.rc('text', usetex=True)
plt.xlim([0, 0.1])
plt.ylim([0, 0.1])
plt.gca().set_aspect('equal')
plt.gca().tick_params(direction='out', which='both')
plt.minorticks_on()
plt.gca().set_xticklabels([])
plt.gca().set_yticklabels([])
plt.title('$\mathsf{Cavity\ tutorial,\ u_x}$')
plt.savefig("cavity-ux.png", dpi=300, bbox_inches='tight')
plt.savefig("cavity-ux.pdf", bbox_inches='tight')

Very nice! Thanks for sharing!

Hello, Thank you for this helpful script. I am trying to use this script for OpenFOAM vtk files. I get the error reads as:
TypeError: 'float' object cannot be interpreted as an integer on the "tri = N.zeros((ntri, 3))" .
Does it something to do with the format setting in surface I defined in my OpenFOAM case?

In Python 3 division is floating point by default, no this line

ntri = triangles.GetNumberOfTuples()/4

makes ntri float, yet VTK requires ntri to be int. So, you just use integer division (ntri = triangles.GetNumberOfTuples()//4) to make VTK happy.

Thank you for your help! I have a more broad question on the realm of VTK files. I can see in my .vtk file for velocity on a "yz" plane that the coordinates are saved under there columns and N rows (N as the number of data point ?) as POLYDATA and the U values as POINT_DATA containing the u1,u2,u3 component in the same file. I want to be able to have analysis over this data, including the calculation of volume flux (U * area) for a certain part of the data (around a certain area in which velocity satisfies any criterion, e.g., contour-like ideas). What are the utilities that I need to this end?

It depends. If you have one file, then you can use ParaView and Programmable filter there. If there are 20 files, then you use Python with VTK (though, I knew students, who were ready to post-process manually 100 files). Post-processing pipe should be something like:

• Point data to Cell Data (for velocity).
• Calculator (for velocity magnitude).
• Threshold (for selecting cells for further post-processing).
• Calculate scalar product of cell area vector and velocity to calculate flux.