从VTU-fi中提取三角形id

Question

我有一个pvtu文件与相关的vtu文件，我想从中显示一些数据。如果我在Paraview（5.6+）中加载pvtu，当我选择纯色（白色）和带边的曲面时，会得到以下图像： 网格在接近顶部边界处具有明显的各向异性，几乎是扁平的三角形；这是预期的行为。你知道吗

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

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文件中获取三角形？你知道吗

谢谢你的帮助。你知道吗

Answer 1

根据程序员的回答，以下代码允许我实现与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()

这将显示

Answer 2

您的问题是您没有使用triangles的matplotlib.tri选项。实际上，如果不在matplotlib中指定，则ParaView中存在的网格的连接性将丢失。事实上，您给matplotlib一个自由来呈现它想要的任何单元，当您知道三角形网格的连通性时，这显然是不正确的。可以使用以下命令提取三角形网格的连接性：

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)