[Paraview] visualising a numpy array

[Juha Jäykkä]

> How are you getting the PETSc._DMDA_Vec_array into ParaView? The
> Programmable Filter is also available in ParaView...

I was going to use 

my_numpy_data = petsc_dmda_array[:]
vtkobject = paraview.numpy_support.numpy_to_vtk(my_numpy_data.copy())

That works fine, but then I got stuck with how to stick this into paraview's 
pipeline.

> You'll need to set up a vtkImageData as the output of a Programmable
> Source/Filter at the very least with the proper dimensions, then use
> the vtk.numpy_interface.dataset_adapter tools to adapt your numpy
> array as a VTK data array that you add to the vtkImageData's point
> data.

I would have thought I need a vtkStructuredGrid or vtkRectilinearGrid, but I 
don't seem able to figure out how to use them.

With vSG I can

pts = vtk.vtkPoints()
pdo = self.GetOutput()
newArray = vtk.vtkDoubleArray()
# this is how far I've got before...
newArray.SetName("testdata")
newArray.SetNumberOfComponents(1)
pdo.GetPointData().AddArray(newArray)

# stuff our numbers into pdo.GetPointData().GetArray(0).InsertNextValue(float) 
one by one,
# followed by pts.InsertNextPoint(x,y,z), so need to get xyz as well
xmin,xmax,Nx,ymin,ymax,Ny,zmin,zmax,Nz = -1,1,10,-1,2,20,-1,3,30
loci = 
numpy.mgrid[xmin:xmax:1j*Nx,ymin:ymax:1j*Ny,zmin:zmax:1j*Nz].reshape(3,-1)
values = numpy.random.random(loci[0].shape)
for ind in range(loci[0].shape[0]):
    pdo.GetPointData().GetArray(0).InsertNextValue(values[ind])
    pts.InsertNextPoint(loci[0,ind],loci[1,ind],loci[2,ind])
pdo.SetPoints(pts)
exts = [xmin,xmax,ymin,ymax,zmin,zmax]
pdo.SetExtent(exts)

which seems to do the right thing except I can not make paraview plot anything 
using this. In the ProgrammableSource's Information tab I do see the right 
bounds and extents (but I am a bit confused about what the "dimension" means 
in extents: it says 3, 4, 5 for x, y, z, which makes me think it has used unit 
lattice spacing even though I would have thought it uses whatever it finds in 
pts. The "Data Arrays" box has "testdata" in it and has the correct value 
range. Also the number of points is correct (6000). I don't know where it gets 
its 24 cells from, though. (24 = 2*3*4, so I am a bit concerned it gets this 
implicitly from bounds again with unit spacing).

But I don't understand what I'm doing wrong as I'm basically modifying 
examples from the wiki.

Using vtkImageData I can only get a Number of Points: 60 (and no plot) and 
using vtkRectilinearGrid I get the same. I can fix the 60 -> 6000 by fixing 
the extents, but that still leaves the bounds and no plot on screen.

> Parallel is tougher... You can do the full data set in serial and
> distribute using the D3 filter, but obviously that has some downsides.

Unacceptable ones: I cannot even imagine what that bottleneck would do to the 
code.

> I believe it is possible to distribute the dataset using MPI in a
> Programmable Source/Filter, but have never done so myself.

Hmm... so where does the ProgrammableSource run? I was planning on running the 
whole code in parallel, but didn't really think about how to do that yet. Now 
that I start thinking about it, Catalyst is probably the way to do it, isn't 
it?

Cheers,
Juha

> On Tue, Apr 26, 2016 at 8:49 AM, Juha Jäykkä <juha.jaykka at gmail.com> wrote:
> > Dear list,
> > 
> > What's the current best practice to visualise a (3D) numpy array with
> > paraview? In particular, I'm wondering if
> > 
> > vtkdata = numpy_support.numpy_to_vtk(num_array=NumPy_data.ravel(),
> > deep=True, array_type=vtk.VTK_FLOAT)
> > 
> > can somehow be directly visualised? In parallel.
> > 
> > I'm aware of ProgrammableSource() and how to do it using that, but only
> > from a file! What I'd very much like to do is avoid the file (if I cannot
> > avoid the file, I'll just do xdmf and hdf5.)
> > 
> > Just in case it matters, the numpy array is actually not really a numpy
> > array, but petsc4py.PETSc._DMDA_Vec_array which of course can be "cast"
> > as numpy.ndarray.
> > 
> > Best regards,
> > Juha
> > 
> > _______________________________________________
> > Powered by www.kitware.com
> > 
> > Visit other Kitware open-source projects at
> > http://www.kitware.com/opensource/opensource.html
> > 
> > Please keep messages on-topic and check the ParaView Wiki at:
> > http://paraview.org/Wiki/ParaView
> > 
> > Search the list archives at: http://markmail.org/search/?q=ParaView
> > 
> > Follow this link to subscribe/unsubscribe:
> > http://public.kitware.com/mailman/listinfo/paraview