vtkファイルの出力(python) #Python

VisItで3次元可視化をするためにvtkファイルを作る必要があった。
pythonでやる方法を調べていたらpyevtkを使うのが良さそう。

ここらへんを参考にpyevtkをインストール

ここでは、簡単なつまらない例を紹介。
3次元ガウス関数で球をつくる。

import numpy as np
from pyevtk.hl import *

ix = 128
jx = 128
kx = 128

xmin = 0.
xmax = 1.
ymin = 0.
ymax = 1.
zmin = 0.
zmax = 1.

dx = (xmax-xmin)/ix
dy = (ymax-ymin)/jx
dz = (zmax-zmin)/kx

x = xmin + 0.5*dx + np.arange(ix)*dx
y = ymin + 0.5*dy + np.arange(jx)*dy
z = zmin + 0.5*dz + np.arange(kx)*dz

xx,yy,zz = np.meshgrid(x,y,z)

ds = 0.1

x0 = (xmax+xmin)/2.
y0 = (ymax+ymin)/2.
z0 = (zmax+zmin)/2.

amp = 1.0

rr = (xx-x0)**2+(yy-y0)**2+(zz-z0)**2
qq = amp*np.exp(-rr/(2.*ds**2))

imageToVTK('./test',pointData={'quantity':qq})

大事なのは上から2行目と最後の行。

from pyevtk.hl import *

は何も考えずに書いておけばOK。

imageToVTK('./test',pointData={'quantity':qq})

でvtkファイルを作っている。
'./test' 部分は保存するファイルの名前を指定。カレントディレクトリにtest.vtiというファイルが作られる。
'quantity'部分はVisItでファイルを開いたときの変数の名前。好きな変数名でOK。
qq には保存したい変数を指定。

実際にVisItでtest.vtiを読み、可視化したのが下図。

imageToVTKでは座標を指定できないので、座標を指定したい場合はgridToVTKを使って、

gridToVTK('./test',x,y,z,pointData={'quantity':qq})

とすればよい。カレントディレクトリにはtest.vtrができる。

特殊なメッシュの張り方をしている場合などは
https://www.vtk.org/Wiki/VTK/Writing_VTK_files_using_python
を参照。