Извлечение идентификаторов треугольников из файла VTU - PullRequest
Новая
       76

Извлечение идентификаторов треугольников из файла VTU

0 голосов
/ 30 октября 2019

У меня есть файл pvtu со связанными файлами vtu, из которого я хочу отобразить некоторые данные. Если я загружаю pvtu в Paraview (5.6+), я получаю следующее изображение, когда выбираю Сплошной цвет (белый) и Поверхность с краями: enter image description here Сетка явно анизотропная вблизи верхней границы спочти сплющенные треугольники;это ожидаемое поведение.

Если я сейчас загружу ту же pvtu в Python и выведу сетку следующим образом, 

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

Я получу это: enter image description here enter image description here, где вы можете легко увидеть, что анизотропия отсутствует. Поэтому мой наивный вопрос: как мне воспроизвести сетку, отображаемую в Paraview с Python? Однако, возможно, есть более точный вопрос. Я полностью осознаю, что библиотека триангуляции matplotlib принимает треугольники в качестве аргумента, но я не могу найти команду, чтобы извлечь их из pvtu. Так что, возможно, лучший вопрос - как получить треугольники из файла pvtu?

Любая помощь приветствуется.

Ответы [ 2 ]

1 голос
/ 09 ноября 2019

Ваша проблема в том, что вы не используете triangles параметр matplotlib.tri. Фактически, связность мешей, которые присутствуют в ParaView, теряется, если вы не укажете это в matplotlib. Фактически, вы предоставляете 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)
0 голосов
/ 10 ноября 2019

Исходя из ответа Alone Programmer, следующий код позволил мне получить ту же сетку, что и 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()

Это отображает

enter image description here

...