Feature Request - Support for POLYDATA in vtk format

[zhaofenqiang]

Hi, thank you for the project for reading UNSTRUCTURED_GRID data in vtk format. I am just wondering could you help extend the code for reading the POLYDATA in vtk format too? Thank you!

The error occured when I try to read POLYDATA is shown below.
Traceback (most recent call last):

  File "<ipython-input-1-a89949ddc2bd>", line 3, in <module>
    mesh = meshio.read('smooth100.vtk')  # optionally specify file_format

  File "/home/fenqiang/anaconda3/envs/latest/lib/python3.6/site-packages/meshio-2.1.0-py3.6.egg/meshio/helpers.py", line 162, in read
    return format_to_reader[file_format].read(filename)

  File "/home/fenqiang/anaconda3/envs/latest/lib/python3.6/site-packages/meshio-2.1.0-py3.6.egg/meshio/vtk_io.py", line 112, in read
    out = read_buffer(f)

  File "/home/fenqiang/anaconda3/envs/latest/lib/python3.6/site-packages/meshio-2.1.0-py3.6.egg/meshio/vtk_io.py", line 161, in read_buffer
    dataset_type

AssertionError: Only VTK UNSTRUCTURED_GRID supported (not POLYDATA).

[kayarre]

if I have a surface mesh defined by [nx3] points and [nx3] triangles how can I export it to the vtk xml polydata format (".vtp" extension) from meshio? is this related to the above?

[keurfonluu]

Not sure whether it can be done with meshio. You could give pyvista a try but it requires the library vtk to be installed.

If you have two arrays points and cells of shape (n, 3), you can try:

import numpy
import pyvista

poly = pyvista.PolyData(points, numpy.column_stack((numpy.full(len(cells), 3), cells)))
poly.save("filename.vtp")

Did not try it myself, though.

[kayarre]

It’s really easy if vtk is installed. Mainly trying to keep it simple. Once you import in paraview one can extract the surface which gives a polydata.

-Kurt

On Thu, Mar 12, 2020 at 12:53 PM Keurfon Luu [email protected] wrote:

Not sure whether it can be done with meshio. You could give pyvista a try but it requires the library vtk to be installed.

If you have two arrays points and cells of shape (n, 3), you can try:

import numpyimport pyvista

poly = pyvista.PolyData(points, numpy.column_stack((numpy.full(len(cells), 3), cells))) poly.save("filename.vtp")

Did not try it myself, though.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/nschloe/meshio/issues/527#issuecomment-598334665, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACE44VJFNDQTKBLPSJ22ZRTRHEOQVANCNFSM4JPVN54Q .

-- Kurt Sansom

[tomfrankkirk]

I've also been coming up against this problem. I've previously used PyVista for this, but my understanding now is that PyVista uses meshio to perform VTK file IO? If so, it seems a bit counter-intuitive PyVista can read triangle PolyData but meshio cannot (though, it can read the same triangle surfaces if they are marked as UnstructuredGrid, not PolyData)

[keurfonluu]

I've also been coming up against this problem. I've previously used PyVista for this, but my understanding now is that PyVista uses meshio to perform VTK file IO?

As far as I know, pyvista is still using VTK readers for every VTK files. It only switches to meshio if the format is not natively supported by vtk (e.g. Dolfin XML, Gmsh...).

[T4mmi]

If I may, what would be the choking point in supporting PolyData (*.vtp) files in meshio ?

is there any technical reason this isn't included ? (i mean except for the time needed to implement it)

[T4mmi]

I started to work on this, but i'm not sure what is the "correct" meshio way to handle polygon cells...

From meshio types they should be "polygon" but since the length of vertices might vary it would mean handle them not as numpy array but as a list (or at least as an array of obj, not numerical type). Looking at the VTU reader I saw that meshio clustered polygons by length to force np.array but here it's a bit problematic since the cell_data ordering has no direct link with cells (based on vtkPolyData description)...

So I'm looking for advice on the best option between the headache with "polygonX" and simple "polygon" with dtype=np.object

[edit] same goes for lines and triangle stripes in fact ...

[nschloe]

Looking at the VTU reader I saw that meshio clustered polygons by length

The order is kept as it is in the file, it's just that blocks of the same length are packed together in one CellBlock.

[T4mmi]

I might not have read it correctly but it seems to me that there is a reodering grouping of all cells by length here: https://github.com/nschloe/meshio/blob/b1753bdf450413b4590ecaced6dcf4ce65ffc7d2/meshio/vtu/_vtu.py#L94 unless i'm doing it wrong ... here is an example:

import vtk 
import meshio

def foo(file="polygons.vtu"):

    # Setup four points
    points = vtk.vtkPoints()
    points.InsertNextPoint(0.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 1.0, 0.0)
    points.InsertNextPoint(0.0, 1.0, 0.0)

    # Create the polygons
    quad = vtk.vtkPolygon()
    quad.GetPointIds().SetNumberOfIds(4)  # make a quad
    quad.GetPointIds().SetId(0, 0)
    quad.GetPointIds().SetId(1, 1)
    quad.GetPointIds().SetId(2, 2)
    quad.GetPointIds().SetId(3, 3)
    tri = vtk.vtkPolygon()
    tri.GetPointIds().SetNumberOfIds(3)  # make a triangle
    tri.GetPointIds().SetId(0, 0)
    tri.GetPointIds().SetId(1, 1)
    tri.GetPointIds().SetId(2, 2)

    # Add the polygon to a list of polygons
    polygons = vtk.vtkCellArray()
    polygons.InsertNextCell(quad)
    polygons.InsertNextCell(tri)
    polygons.InsertNextCell(quad)

    # Create a VTU
    vtu = vtk.vtkUnstructuredGrid()
    vtu.SetPoints(points)
    vtu.SetCells(7, polygons)

    # Write VTU
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetInputData(vtu)
    writer.SetFileName(file)
    writer.Write()

    mesh = meshio.read(file)
    mesh.cells



if __name__ == "__main__":
    foo()

which produces :

[<meshio CellBlock, type: polygon3, num cells: 1>,
 <meshio CellBlock, type: polygon4, num cells: 2>]

here the polygons should be [polygon4, polygon3, polygon4] ...

[nschloe]

Good catch! That has to go away. I'll check it out now.

[T4mmi]

ok then i'll stay tuned ;)

[T4mmi]

But in the VTU it's no big deal since the cell_data is also re-ordered ... the problem in vtp is that is has to reorder cell_array but only in the slice of polygons (and after that there will be the same problem for lines ... and triangle strips)

[nschloe]

But in the VTU it's no big deal since the cell_data is also re-ordered .

meshio promises not to mess with the data, and that includes the order of the cells. Definitely a bug.

[nschloe]

Alright, fixed now. (At least in VTU.)

[kayarre]

any updates on supporting the vtp format with meshio. vtk is a heavyweight library and it would be nice in meshio or maybe derived from meshio? I am giving it a go at an implementation.