从VTU文件中提取三角形的ID

1 人关注

我有一个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文件中获得三角形？
希望得到任何帮助。


         python


         matplotlib

vtk


         paraview


        2
        
        个回答


          
           
           
            Mithridates the Great
           
          
          
           发布于
           
           2019-11-10


          已采纳


         0
         
         人赞同


          
           你的问题是，你没有使用
           
            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)


           
            
             非常感谢你的回答，我确实在期待一个提供三角形的方法，但我不知道如何用vtk获得它们。然而，请注意，你的代码对我来说并不能正常工作。我不得不稍作修改以使其工作，请看我下面的回答。


          
           
            根据Alone程序员的回答，下面的代码让我实现了与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()