はじめに
何をしたいか
２次元の温度シミュレーション結果のVTKファイルをjupyterで表示したい。
どうやって？
VTKファイルをPythonで読み込み、座標データと温度データを取り出して、Matplotlibで表示する。
環境
OS: Windows 10 home edition
Tool: Anaconda(Python 3.7.4)
シミュレーション結果のVTKファイル
有限要素法フリーソフトウェアであるFreefem++を利用して、熱拡散方程式を解き、その結果のVTKファイルを作成した。熱拡散方程式のコードはexampleのHeat.edpを利用した。Freefem++のインストールと実行方法は割愛する。利用したVTKファイルは、ここにおいておく。
Python VTKライブラリのインストール
conda を使ってインストールする。
conda install vtk
import 関係
操作に必要なNumpy、Matplotlib、VTKのライブラリをimportする。
import numpy as np
import maptlotlib.pyplot as plt
import vtk
from vtk.util import numpy_support
%matplotlib inline
VTKファイルの読み込み
VTKファイルが非構造格子だったため、Unstructuredで読み込む。
#ファイル読み込み
filename = "heat_result.vtk"
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
温度のCell DataのPoint Dataへの変換
今回のVTKファイルの座標データがPoint Dataで、温度データがCell Dataとなっている。座標データと温度データを一対一の関係にしたいので、全部Point Dataに変換する。
#cell data から point data変換
cell2point = vtk.vtkCellDataToPointData()
cell2point.SetInputData(reader.GetOutput())
cell2point.Update()
変換したデータからの座標と温度データの抜き出し
VTKファイルを読み込み、vtk_supportを利用して、VTKデータから座標データと温度データを取り出す。
#座標と温度データのnumpy化
coord = numpy_support.vtk_to_numpy(cell2point.GetOutput().GetPoints().GetData())
x = coord[:,0]
y = coord[:,1]
z = coord[:,2]#使わない。２次元の熱拡散方程式の結果なのですべて0である。
#GetAbstractArray(0)にはLabel、GetAbstractArray(1)に温度データが入っている。
temperature = numpy_support.vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray(1))
Matplotlibでのカラーマップ表示
準備が整ったので、あとはMatplotlibで表示するだけである。非構造データのため、tricontourfを利用してカラーマップを表示する。
##カラーマップ出力
plt.tricontourf(x,y,temperature,levels=15,cmap="jet")
plt.colorbar()
全コード
import numpy as np
import matplotlib.pyplot as plt
import vtk
from vtk.util import numpy_support
%matplotlib inline
#ファイル読み込み
filename = "heat_result.vtk"
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
#cell data から point data変換
cell2point = vtk.vtkCellDataToPointData()
cell2point.SetInputData(reader.GetOutput())
cell2point.Update()
#座標と温度データのNumpy化
coord = numpy_support.vtk_to_numpy(cell2point.GetOutput().GetPoints().GetData())
x = coord[:,0]
y = coord[:,1]
z = coord[:,2]#使わない。２次元の熱拡散方程式の結果なのですべて0である。
#GetAbstractArray(0)にはLabel、GetAbstractArray(1)に温度データが入っている。
temperature = numpy_support.vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray(1))
#カラーマップ出力
plt.tricontourf(x,y,temperature,levels=15,cmap="jet")
plt.colorbar()
参考
VTKクラスは、c++のclassリファレンスが公開されている。Python VTKの構造と名前はc++リファレンスと同じなので、参考にはできる。主に参考にしたの下記のサイトである。