Extract Triangles IDs from VTU file

Question

I have a pvtu file with associated vtu files from which I want to display some data. If I load the pvtu in Paraview (5.6+), I get the following image when I choose Solid Color (White) and Surface With Edges: The mesh is clearly anisotropic close to the top boundary with almost flattened triangles; this is the expected behaviour.

If I now load the same pvtu in Python and display the mesh in the following manner,

import numpy
import matplotlib.pyplot as plt
import vtk
gridreader = vtk.vtkXMLPUnstructuredGridReader()
gridreader.SetFileName('whatever.pvtu')
gridreader.Update()
vtkOut = gridreader.GetOutput()
vtkData = vtkOut.GetPoints().GetData()
coords = numpy.array([vtkData.GetTuple3(x)
                      for x in range(vtkData.GetNumberOfTuples())])
plt.triplot(coords[:, 0], coords[:, 1])
plt.gcf().set_size_inches(16, 8)
plt.gca().set_aspect('equal')
plt.savefig('meshPython1.png', bbox_inches='tight')
plt.gca().set_xlim((5e5, 3e6))
plt.gca().set_ylim((6e5, 1e6))
plt.savefig('meshPython2.png', bbox_inches='tight')

I get that: where you can readily see that the anisotropy is not present. Therefore, my naive question is: how do I reproduce the mesh displayed in Paraview with Python? However, there is probably a more accurate question. I am fully aware that the triangulation library of matplotlib accepts triangles as an argument, but I am unable to find a command to extract them from the pvtu. So maybe a better question would be how to obtain the triangles from a pvtu file?

Any help appreciated.

score 4 · Accepted Answer · answered Nov 09 '19 at 16:47

Your problem is that you don't use triangles option of matplotlib.tri. In fact, connectivity of meshes that are present in the ParaView is lost when you don't specify it in matplotlib. In fact, you give matplotlib a freedom to present cells as whatever it wants, which is not correct obviously when you know the connectivity of your triangular meshes. You can extract the connectivity of triangular meshes by using this command:

cell_connecitivty_matrix = []

for i in range(vtOut.GetNumberOfCells()):
 assert vtkOut.GetCell(i).GetNumberOfPoints() == 3
 cell_connecitivty_matrix.append(vtkOut.GetCell(i).GetPointIds())

cell_connecitivty_matrix = np.array(cell_connecitivty_matrix, dtype=np.float).reshape((vtOut.GetNumberOfCells(),3))

#plot triangles with their connectivity matrix

plt.triplot(coords[:, 0], coords[:, 1], triangles=cell_connectivity_matrix)

Thanks a lot for the answer, I was indeed expecting a way to provide the triangles, but I had no idea of how to obtain them using vtk. However, please note that your code does not work as is for me. I had to slightly modify it to make it work, please see my answer below. — Patol75, Nov 10 '19 at 09:44

score 1 · Answer 2 · answered Nov 10 '19 at 09:42

Based on Alone Programmer's answer, the following code allowed me to achieve the same mesh as Paraview:

import numpy
import matplotlib.pyplot as plt
import vtk
gridreader = vtk.vtkXMLPUnstructuredGridReader()
gridreader.SetFileName('whatever.pvtu')
gridreader.Update()
vtkOut = gridreader.GetOutput()
vtkData = vtkOut.GetPoints().GetData()
coords = numpy.array([vtkData.GetTuple3(x)
                      for x in range(vtkData.GetNumberOfTuples())])
cell_connectivity_matrix = []
for i in range(vtkOut.GetNumberOfCells()):
    assert vtkOut.GetCell(i).GetNumberOfPoints() == 3
    cell_connectivity_matrix.append(
        [vtkOut.GetCell(i).GetPointIds().GetId(j)
         for j in range(vtkOut.GetCell(i).GetPointIds().GetNumberOfIds())])
cell_connectivity_matrix = numpy.array(cell_connectivity_matrix,
                                       dtype=numpy.float)
plt.triplot(coords[:, 0], coords[:, 1], triangles=cell_connectivity_matrix)
plt.gcf().set_size_inches(16, 8)
plt.gca().set_aspect('equal')
plt.show()

This displays

Extract Triangles IDs from VTU file

2 Answers2

Linked