RequestData

import numpy as np
from vtk.util.numpy_support import numpy_to_vtk
from vtk.util.numpy_support import vtk_to_numpy

# Parallel information, useful for debugging only
controller = vtk.vtkMultiProcessController.GetGlobalController()
rank = controller.GetLocalProcessId()

# Read the binary cartesian file
def readcart(fname, dims):
    dims = np.array(dims)  # Easier to use np array to get count of items
    with open(fname, 'rb') as f:
        # Fortran records have INTEGER*4 headers which are useless here
        f.seek(4)

        data = np.fromfile(f, dtype=np.float32, count=np.prod(dims), sep="")
        # Order="F" is the key here
        # Bare in mind that if I need to save in F order I have to ravel
        # first
        return data.reshape(dims, order="F")


xname = "/home/mfontana/Tesis/Codigo/Modes-Display/odir/field_br.R+1+1+0.out"
yname = "/home/mfontana/Tesis/Codigo/Modes-Display/odir/field_bt.R+1+1+0.out"
zname = "/home/mfontana/Tesis/Codigo/Modes-Display/odir/field_bp.R+1+1+0.out"
dims = (dim_x, dim_y, dim_z)

# Construct grid
output.SetDimensions(dims)
output.SetXCoordinates(numpy_to_vtk(linspace(-1, 1, dims[0]), deep=True,
                       array_type=vtk.VTK_FLOAT))
output.SetYCoordinates(numpy_to_vtk(linspace(-1, 1, dims[1]), deep=True,
                       array_type=vtk.VTK_FLOAT))
output.SetZCoordinates(numpy_to_vtk(linspace(-1, 1, dims[2]), deep=True,
                       array_type=vtk.VTK_FLOAT))

# Get data
# Reshaping is because I need two axes to hstack and get (n, 3). If I ravel
# I get (3n,) with hstack
# Must use "F" order again as reshaping only changes access order, not
# memory location (implying that every reshaping must use "F" order)
bx = readcart(xname, dims).reshape(-1, 1, order="F")
by = readcart(yname, dims).reshape(-1, 1, order="F")
bz = readcart(zname, dims).reshape(-1, 1, order="F")

b = np.hstack((bx, by, bz))

# Convert to VTK and give names
array = numpy_to_vtk(b, deep=True, array_type=vtk.VTK_FLOAT)
array.SetName("B")
array.SetComponentName(0, "r")
array.SetComponentName(1, "t")
array.SetComponentName(2, "p")

# Associate data to grid
output.GetPointData().AddArray(array)

output.SetExtent(0, dims[0]-1 , 0, dims[1]-1 , 0, dims[2]-1)