numpy arrays to tvtk/vtk vector on nonuniform rectilinear grid

Question

Actually, the ultimate task for me is to plot streamline of nonuniform rectilinear data.

My field data are stored in a giant array, myData, of shape Nx * Ny * Nz * nComponents. My coordinates are stored in 1d arrays of shapes Nx, Ny, and Nz respectively.

I was able to plot the scalar data using the following code:

from tvtk.api import tvtk
r = tvtk.RectilinearGrid()
# getCoordinates is my custom function
# to get 1d arrays of nonuniform grid coordinates
r.x_coordinates = getCoordinates('x')
r.y_coordinates = getCoordinates('y')
r.z_coordinates = getCoordinates('z')
# get a component, in this case, the density, 'rho'
fldIndex = 0
fldName = 'rho'
s = myData[:, :, :, fldIndex]
# set scalar, I have to pass the parameter 'F' to ravel
r.point_data.scalars = s.ravel(order='F')
r.point_data.scalars.name = fldName
r.dimensions = s.shape
from mayavi import mlab
surf = mlab.pipeline.iso_surface(r, opacity=0.1)

Now I would like add a vector, velocity, to r, so that I can use it as input of

vnorm = mlab.pipeline.extract_vector_norm(velocity)
vline = mlab.pipeline.streamline(vnorm) # plus many tuning parameters

Here, the three components of velocity are stored in myData[:,:,1], myData[:,:,2], myData[:,:3]. How can I utilize them? I naively tried the following codes but could not get any streamline

# the following code did not lead to correct streamline
u = myData[:,:,:,1].ravel(order='F')
v = myData[:,:,:,2].ravel(order='F')
w = myData[:,:,:,3].ravel(order='F')
import numpy as np
vec = np.array([u,v,w]).T # now vec is NOT contiguous and
                          # NOT accepted by the internal numpy_to_vtk
                          # conversion
r.point_data.vectors = vec
r.point_data.vectors.name = 'velocity'

PS: Initially I tried to use mayavi's vector_field. However, it depends on vtk.ImageData, i.e., it only supports uniform grid.

Answer 1

The following code works

from tvtk.api import tvtk
r = tvtk.RectilinearGrid()
r.x_coordinates = getCoordinates('x')
r.y_coordinates = getCoordinates('y')
r.z_coordinates = getCoordinates('z')
# vec is a numpy array of the right shape and storage layout
import numpy as np
vec = np.zeros([size, 3],dtype='float64') #FIXME: 
vec[:, 0] = myData[:,:,:,1].ravel(order='F')
vec[:, 1] = myData[:,:,:,2].ravel(order='F')
vec[:, 2] = myData[:,:,:,3].ravel(order='F')
# vec from the following code also works;
# the key is to generate a contiguous array, which was done by 'zeros'
# above, and is done by ascontiguousarray below
# vec = np.ascontiguousarray( np.array([
#    myData[:,:,:,1].ravel(order='F'),
#    myData[:,:,:,2].ravel(order='F'),
#    myData[:,:,:,3].ravel(order='F')]).T )
# use numpy_to_vtk to get the vtk array
r.point_data.vectors = vec
r.point_data.vectors.name = 'velocity'
Nx,Ny,Nz,Ncomp = MyData.shape
r.dimensions = (Nx, Ny, Nz)

# make streamline
from mayavi import mlab
magnitude = mlab.pipeline.extract_vector_norm(r)
field_lines = mlab.pipeline.streamline(magnitude) # streamline options neglected