Vtk tools

Analysing fluidity output in NumPy/SciPy using vtktools

The Python module Vtktools allows you to analyse the output data of fluidity (vtu files) such as scalar fields, vector fields and nodal positions, into NumPy arrays using python. SciPy and NumPy form a very powerfull Python library that can be used for all kinds of mathematical analysis. For those of you who are used to do their analysis in Matlab: SciPy/NumPy provides (almost) all of the functionality of Matlab (arrays, matrices and all possible operations and function for them) and combines it with the superior language of Python. If you don't know any Python yet, the basic use of array operations in NumPy is actually quite similar to Matlab.

Using vtktools

Initially, it will probably be easier to use the interactive interface ipython. The following script imports the vtktools module into python

import sys
sys.path.append('<<fluidity_source_path>>/python')
import vtktools

Then reads in a vtu file. The information from the vtu file is then contained in the vtu object "ug".

ug=vtktools.vtu('run_123.vtu')
#to extract from all your files use this inside a loop:
# ug=vtktools.vtu('run_123_'+str(i)+'.vtu')

For a list of all contained scalar and vector fields:

ug.GetFieldNames()

At this stage, scalar and vector quantities can be extracted. Pressure is extracted and put into a NumPy array called p, and velocity is put into the array uvw.

p=ug.GetScalarField('Pressure')
uvw=ug.GetVectorField('Velocity')

It is possible to add a field to the vtu object. It is also very simple to add new fields to the vtu object. If we want for instance the square of the pressure as an extra field we do:

psq=p**2
ug.AddScalarField('PressureSquared', psq)

Finally if we want to write the changed vtu object which also contains the new field PressureSquared. The new vtu file can then be opened in Mayavi or Paraview, and the new field will be available for visual analysis.

ug.Write()

By not specifying a filename, the original vtu file will be overwritten, but since we have only added a field nothing is lost.

For an overview of all possibilities of vtktools.vtu: help('vtktools.vtu')

Example usage of NumPy and vtktools

The following subtracts a linear temperature profile from the linear stratified flow past a Gaussian seamount problem:

 python
 import sys
 sys.path.append('<<fluidity/tools directory path>>')
 import vtktools
 ug=vtktools.vtu('<<vtu file>>')
 ug.GetFieldNames()
 p=ug.GetScalarField('Pressure')
 uvw=ug.GetVectorField('Velocity')
 temp=ug.GetScalarField('Temperature')
 pos=ug.GetLocations()
 z=pos[:,2]
 temp2 = temp  + 5.55555555555*z
 ug.AddScalarField('<<name of new field>>', temp2)
 ug.Write('<<new vtu file name>>')

Vtktools output an array with the values of all the fields at each point of the mesh that was used in the simulation. However, this mesh could be unstructured or adapted: in this case, the points would not be evenly distributed in space. For some diagnostics it may be desirable to have arrays giving the values of the fields at points that are evenly distributed in space. For a two-dimensional problem, this can be done with the following script:

import sys
sys.path.append(’<<fluidity_source_path>>/tools’)
import vtktools
u=vtktools.vtu("run_123.vtu") 
Xlist = arange(0.0,4.001,0.01)# x co-ordinates of the desired array shape
Ylist = arange(0.0,2.001,0.01)# y co-ordinates of the desired array shape
[X,Y] = meshgrid(Xlist,Ylist) 
Z = 0.0*X # This is 2d so z is an array of zeros.
X = reshape(X,(size(X),))
Y = reshape(Y,(size(Y),))
Z = reshape(Z,(size(Z),))
zip(X,Y,Z)
pts = zip(X,Y,Z)
pts = vtktools.arr(pts)
# create arrays of velocity and temperature values at the desired points
vel = u.ProbeData(pts, ’Velocity’)
temperature_structured = u.ProbeData(pts, ’Temperature’)

The script reads again the data in a vtu file using vtktools' method "vtu". Then, it sets the desired x and y co-ordinates by putting them into an array, and uses the function "meshgrid" to create a mesh of size determined by the properties of that array. This mesh will then have its points evenly distributed in space, as required. Once the z coordinate is set to zero to deal with the 2-dimensional problem, the individual 1-dimensional arrays corresponding to the three Cartesian coordinates are zipped together. A 3-dimensional array in space, called "pts", is then created. Now, field data can be probed with the "ProbeData" method, and assigned to the evenly distributed points in space of coordinates stored in the "pts" array.

Python tips

The following can be used to ensure that filenames are treated in the correct order, i.e. nicely sorted so filename_100.vtu comes after filename_33.vtu for example, taken from http://www.codinghorror.com/blog/archives/001018.html.

 #!/usr/bin/env python
 import sys
 import re 
 def sort_nicely( l ): 
   convert = lambda text: int(text) if text.isdigit() else text 
   alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] 
   l.sort( key=alphanum_key ) 
 filelist = sys.argv[1:]
 sort_nicely(filelist)

This allows you to write a python script that takes *.vtu as command line input for example and outputs whatever it is you are doing/calculating in correct time order.