从 VTU 文件中提取三角形 ID [英] Extract Triangles IDs from VTU file

问题描述

我有一个带有相关vtu文件的pvtu文件，我想从中显示一些数据.如果在Paraview(5.6+)中加载pvtu，则在选择纯色(白色)"和带有边缘的曲面"时会得到以下图像:网格显然是各向异性的，靠近顶部边界，三角形几乎是扁平的.这是预期的行为.

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.

如果我现在在Python中加载相同的pvtu并以以下方式显示网格，

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')

我明白了:您可以很容易地看到各向异性不存在.因此，我天真的问题是:如何使用 Python 重现 Paraview 中显示的网格?但是，可能存在一个更准确的问题.我完全知道 matplotlib 的三角剖分库接受三角形作为参数，但我无法找到从 pvtu 中提取它们的命令.因此，也许更好的问题是如何从pvtu文件中获得三角形?

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?

感谢任何帮助.

推荐答案

你的问题是你没有使用 matplotlib.tri 的 triangles 选项.事实上，当您没有在 matplotlib 中指定它时，存在于 ParaView 中的网格的连接会丢失.实际上，您赋予了matplotlib一个自由呈现其所需单元格的自由，当您知道三角形网格的连通性时，这显然是不正确的.您可以使用以下命令提取三角形网格的连通性:

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)

这篇关于从 VTU 文件中提取三角形 ID的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！