我有一个带有关联vtu文件的pvtu文件,我要从中显示一些数据。如果在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文件中获得三角形?
感谢任何帮助。
您的问题是您没有使用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)
基于独立程序员的回答,以下代码使我可以实现与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()
此显示