VTK : how to read grid cells' lenth, width and height?

Question

I have a huge grid in *.pvd format. I would like to ensure some cells size specification have been respected when building said grid. To do so, I should get a cell data array with (dx,dy,dz)

I first tried to check this in Paraview with very little success. Then I resolved to export the mesh in various format (vtk, vtu, ex2) and import things into python using the vtk module, as in the code below. Unfortunately, the size of the mesh forbids it and I get various error messages stating "Unable to allocate n cells of size x".

import vtk
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName("my_mesh.vtu")
reader.Update()

Finally, in Paraview there is a python-shell that allows me to open the grid file in either pvd or vtk format:

>>> from paraview.simple import *
>>> my_vtk = OpenDataFile("my_mesh.vtk")
>>> print dir(my_vtk)

Despite my browsing the methods and attribute of this reader object, I remain clueless about where to fetch any geometry information on the grid. I also browsed through the simple module documentation and I can't really wrap my head around it.

So how can one retrieve information regarding the geometry of cells from a paraview.servermanager.LegacyVTKReader object?

Any clue about how to achieve this with through the paraview GUI, or any kludge to load the vtk object into python vtk despite the memory issue is also very welcome. Sorry for such a hazy question, but I don't really know where to get started...

Unknown Unknown · Accepted Answer · 2017-12-08T23:56:04

You can use GetClientSideObject() (see here) to get a VTK object in the Paraview Python shell. After that you can use all the regular VTK Python functions. For example, you can write the following in the Paraview Python shell

>>> from paraview.simple import *
>>> currentSelection = GetActiveSource()
>>> readerObj = currentSelection.GetClientSideObject()
>>> unstructgrid = readerObj.GetOutput()
>>> firstCell = unstructgrid.GetCell(0)
>>> cellPoints = firstCell.GetPoints()

Alternatively, you can use Programmable Filter in ParaView. This allows access to full VTK python module and even NumPy or other modules. You can enter following script in the script window of the programmable filter:

import vtk as v
import numpy as np

inp = self.GetUnstructuredGridInput()
cells = inp.GetCells()
cells.InitTraversal()
cellPtIds = v.vtkIdList()
lenArr = v.vtkDoubleArray()
lenArr.SetNumberOfComponents(3)
lenArr.SetName('CellSize')
while cells.GetNextCell( cellPtIds ):
    pts = []
    for i in range( cellPtIds.GetNumberOfIds() ):
        ptCoords = inp.GetPoint( cellPtIds.GetId(i) )
        pts.append( ptCoords )    
    pts = np.array( pts )
    dx = np.max(pts[:,0]) - np.min(pts[:,0])
    dy = np.max(pts[:,1]) - np.min(pts[:,1])
    dz = np.max(pts[:,2]) - np.min(pts[:,2])
    lenArr.InsertNextTuple3(dx, dy, dz)
out = self.GetUnstructuredGridOutput()
out.ShallowCopy( inp )
out.GetCellData().AddArray( lenArr )

In Paraview when you select the 'ProgrammableFilter1' icon in your pipeline, a new cell data array will be available to you from the drop-down as shown in the screenshot below. You can modify the script above to save the data to file to analyze externally.

VTK : how to read grid cells' lenth, width and height?

2 Answers