Hi Juha, This is a bit beyond my current numpy-foo. Hopefully I'll get a chance to look at it in more depth in the next couple months.
If you figure anything else out, I and others on the list would surely appreciate knowing what you find. Thanks, Cory On Fri, Apr 29, 2016 at 2:19 PM, Cory Quammen <[email protected]> wrote: > On Wed, Apr 27, 2016 at 11:39 AM, Juha Jäykkä <[email protected]> wrote: >>> 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. > > If your grid points are regularly spaced in each dimension, which I > think they are in this case, you can use a vtkImageData. See the > ParaViewGuide (http://www.paraview.org/paraview-guide/) pages 33 and > 34 for pictorial representations of where vtkStructuredGrid or > vtkRectilinearGrid are needed. > >> 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) > > Note that the extents are in the integer index space (e.g., i-j-k) of > the grid points rather than physical space, so exts should be > something like > > [0, 9, 0, 19, 0, 29] > > Also note that xmax, ymax, and zmax should be the highest possible > index, not the number of grid positions in that dimension. > >> 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? > > The Programmable Source is a VTK filter that runs on all nodes on > which you are running the ParaView server. Catalyst may indeed be > useful for you. > > Cory > >> Cheers, >> Juha >> >>> On Tue, Apr 26, 2016 at 8:49 AM, Juha Jäykkä <[email protected]> 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 >> > > > > -- > Cory Quammen > R&D Engineer > Kitware, Inc. -- Cory Quammen R&D Engineer Kitware, Inc. _______________________________________________ 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
