VTK メモ

2018年11月18日

はじめに

データ可視化ライブラリ VTK を Python で使うメモ。

環境

  • Ubuntu 16.04 LTS/macOS Mojave (10.14.1)
  • Anaconda3

VTK は pip か conda でインストールできる。

円筒の表示

VTK では、source あるいは reader から mapper -> actor -> renderer -> interactor とデータをつないでいく。

以下の例では、円筒を作ってウインドウに表示している。

cylinder.py

import vtk
from vtk.util.colors import tomato

# source
cylinder = vtk.vtkCylinderSource()
cylinder.SetResolution(20)

# mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(cylinder.GetOutputPort())

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(tomato)
actor.RotateX(30.)
actor.RotateY(-45.)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)

ren.ResetCamera()
ren.GetActiveCamera().Zoom(1.5)

ren_win.Render()
inter.Initialize()
inter.Start()

マウスイベント

マウスイベントを処理する。以下のコードは、円筒を表示するコードに MouseInteractorStyle を定義して、interactor に渡している。

iren.SetInteractorStyle(MouseInteractorStyle())

処理は、マウスの右クリック、左クリックでメッセージを表示し、それぞれのもとの処理を呼び出しているだけだが、若干挙動がかわる。

mouse.py

import vtk
from vtk.util.colors import tomato

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)
        self.AddObserver("RightButtonPressEvent", self.rightButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        print("LEFT")
        self.OnLeftButtonDown()

    def rightButtonPressEvent(self, obj, event):
        print("RIGHT")
        self.OnRightButtonDown()

# source
cylinder = vtk.vtkCylinderSource()
cylinder.SetResolution(20)

# mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(cylinder.GetOutputPort())

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(tomato)
actor.RotateX(30.)
actor.RotateY(-45.)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

ren.ResetCamera()
ren.GetActiveCamera().Zoom(1.5)

ren_win.Render()
inter.Initialize()
inter.Start()

キー入力

key.py

import vtk

class KeyInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("KeyPressEvent", self.keyPressEvent)

    def keyPressEvent(self, obj, event):
        key = self.GetInteractor().GetKeySym()
        print(key)
        self.OnKeyPress()

# renderer
ren = vtk.vtkRenderer()
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(KeyInteractorStyle())

ren_win.Render()
inter.Initialize()
inter.Start()

次のようにしてもよい。

def keyPressEvent(obj, event):
    key = obj.GetKeySym()
    print(key)

inter.AddObserver("KeyPressEvent", keyPressEvent)

STL の読み書き

球を作り、STL で書き出して、それを読み込む。

stl_writer.py

import vtk

filename = "sphere.stl"

sphere = vtk.vtkSphereSource()
sphere.SetThetaResolution(50)
sphere.SetPhiResolution(50)

writer = vtk.vtkSTLWriter()
writer.SetFileName(filename)
writer.SetInputConnection(sphere.GetOutputPort())
writer.Write()

# reader
reader = vtk.vtkSTLReader()
reader.SetFileName(filename)

# mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(reader.GetOutputPort())

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)

ren_win.Render()
inter.Initialize()
inter.Start()

VTK の読み込み

VTK ファイル (Unstructure Grid) を読み込む。

scalar.vtk

# vtk DataFile Version 2.0
scalar
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
0 0 0
1 0 0
1 1 0
0 1 0
CELLS 1 5
4 0 1 2 3
CELL_TYPES 1
9
POINT_DATA 4
SCALARS point_scalars float
LOOKUP_TABLE default
1
2
3
4
CELL_DATA 1
SCALARS cell_scalars float
LOOKUP_TABLE default
5

vtk_reader.py

import vtk

filename = "scalar.vtk"

# reader
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()

output = reader.GetOutput()
scalar_range = output.GetScalarRange()
print(scalar_range)

# mapper
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(output)
mapper.SetScalarRange(scalar_range)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)

ren_win.Render()
inter.Initialize()
inter.Start()

画像の書き出し

ウインドウの表示を画像ファイルとして書き出す。

image_writer.py

import vtk
from vtk.util.colors import tomato

filename = "cylinder.png"

# source
cylinder = vtk.vtkCylinderSource()
cylinder.SetResolution(20)

# mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(cylinder.GetOutputPort())

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(tomato)
actor.RotateX(30.)
actor.RotateY(-45.)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(0.1, 0.2, 0.4)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

ren.ResetCamera()
ren.GetActiveCamera().Zoom(1.5)
ren_win.Render()

# window to image filter
imageFilter = vtk.vtkWindowToImageFilter()
imageFilter.SetInput(ren_win)
imageFilter.SetInputBufferTypeToRGB()

# image writer
writer = vtk.vtkPNGWriter()
writer.SetFileName(filename)
writer.SetInputConnection(imageFilter.GetOutputPort())
writer.Write()

オブジェクトの見た目の設定

actor からプロパティ (vtkProperty) を受け取れるので、それで設定する。

prop = actor.GetProperty()

色の設定

prop.SetColor(1, 0, 0)

環境光の設定

prop.SetAmbient(0.5)

質感の設定

prop.SetDiffuse(1)
prop.SetSpecular(0.5)
prop.SetSpecularPower(30)

透明度の設定 (0 で完全に透明)

prop.SetOpacity(1)

点の設定

prop.VertexVisibilityOn()
prop.SetVertexColor(0.5, 1, 0.5)
prop.SetPointSize(5)

辺の設定

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

オブジェクトの点表示

prop.SetColor(0, 0, 0)
prop.SetPointSize(3)
prop.SetRepresentationToPoints()

オブジェクトのワイヤーフレーム表示

prop.SetColor(0, 0, 0)
prop.SetLineWidth(1)
prop.SetRepresentationToWireframe()

オブジェクトの面表示

prop.SetRepresentationToSurface()

背景

背景色は renderer で設定する。

ren.SetBackground(0.1, 0.2, 0.4)

グラデーションにするには、次のようにする。

ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

2 のほうが上になる。

グリッドデータを構成する

grid.py

import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

grid = vtk.vtkUnstructuredGrid()
points = vtk.vtkPoints()
points.InsertNextPoint(0, 0, 0)
points.InsertNextPoint(1, 0, 0)
points.InsertNextPoint(1, 1, 0)
points.InsertNextPoint(0, 1, 0)
points.InsertNextPoint(0, 0, 1)
points.InsertNextPoint(1, 0, 1)
points.InsertNextPoint(1, 1, 1)
points.InsertNextPoint(0, 1, 1)

points.InsertNextPoint(1, 0, 0)
points.InsertNextPoint(2, 0, 0)
points.InsertNextPoint(2, 1, 0)
points.InsertNextPoint(1, 1, 0)
points.InsertNextPoint(1, 0, 1)
points.InsertNextPoint(2, 0, 1)
points.InsertNextPoint(2, 1, 1)
points.InsertNextPoint(1, 1, 1)
grid.SetPoints(points)

cells = vtk.vtkCellArray()
hexa = vtk.vtkHexahedron()
hexa.GetPointIds().SetId(0, 0)
hexa.GetPointIds().SetId(1, 1)
hexa.GetPointIds().SetId(2, 2)
hexa.GetPointIds().SetId(3, 3)
hexa.GetPointIds().SetId(4, 4)
hexa.GetPointIds().SetId(5, 5)
hexa.GetPointIds().SetId(6, 6)
hexa.GetPointIds().SetId(7, 7)
cells.InsertNextCell(hexa)

hexa = vtk.vtkHexahedron()
hexa.GetPointIds().SetId(0, 8)
hexa.GetPointIds().SetId(1, 9)
hexa.GetPointIds().SetId(2, 10)
hexa.GetPointIds().SetId(3, 11)
hexa.GetPointIds().SetId(4, 12)
hexa.GetPointIds().SetId(5, 13)
hexa.GetPointIds().SetId(6, 14)
hexa.GetPointIds().SetId(7, 15)
cells.InsertNextCell(hexa)
grid.SetCells(vtk.VTK_HEXAHEDRON, cells)

# filter
id_filter = vtk.vtkIdFilter()
id_filter.PointIdsOff()
id_filter.CellIdsOn()
id_filter.SetInputData(grid)

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

# mapper
mapper = vtk.vtkDataSetMapper()
mapper.SetInputConnection(id_filter.GetOutputPort())
mapper.SelectColorArray(id_filter.GetIdsArrayName())
mapper.SetScalarRange(id_filter.GetOutput().GetScalarRange())
mapper.SetLookupTable(lut)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.VertexVisibilityOn()
prop.SetVertexColor(0.5, 1, 0.5)
prop.SetPointSize(5)

prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

ren_win.Render()
inter.Initialize()
inter.Start()

lookup table は、スカラー値に対応させる色のテーブルである。デフォルトでは小さい値が青で大きい値が赤になっているので、逆にしている。

lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

これは mapper に設定される。

mapper.SetLookupTable(lut)

セルデータを構築する

cell_data.py

import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

grid = vtk.vtkUnstructuredGrid()
points = vtk.vtkPoints()
points.InsertNextPoint(0, 0, 0)
points.InsertNextPoint(1, 0, 0)
points.InsertNextPoint(1, 1, 0)
points.InsertNextPoint(0, 1, 0)
points.InsertNextPoint(0, 0, 1)
points.InsertNextPoint(1, 0, 1)
points.InsertNextPoint(1, 1, 1)
points.InsertNextPoint(0, 1, 1)

points.InsertNextPoint(1, 0, 0)
points.InsertNextPoint(2, 0, 0)
points.InsertNextPoint(2, 1, 0)
points.InsertNextPoint(1, 1, 0)
points.InsertNextPoint(1, 0, 1)
points.InsertNextPoint(2, 0, 1)
points.InsertNextPoint(2, 1, 1)
points.InsertNextPoint(1, 1, 1)
grid.SetPoints(points)

cells = vtk.vtkCellArray()
hexa = vtk.vtkHexahedron()
hexa.GetPointIds().SetId(0, 0)
hexa.GetPointIds().SetId(1, 1)
hexa.GetPointIds().SetId(2, 2)
hexa.GetPointIds().SetId(3, 3)
hexa.GetPointIds().SetId(4, 4)
hexa.GetPointIds().SetId(5, 5)
hexa.GetPointIds().SetId(6, 6)
hexa.GetPointIds().SetId(7, 7)
cells.InsertNextCell(hexa)

hexa = vtk.vtkHexahedron()
hexa.GetPointIds().SetId(0, 8)
hexa.GetPointIds().SetId(1, 9)
hexa.GetPointIds().SetId(2, 10)
hexa.GetPointIds().SetId(3, 11)
hexa.GetPointIds().SetId(4, 12)
hexa.GetPointIds().SetId(5, 13)
hexa.GetPointIds().SetId(6, 14)
hexa.GetPointIds().SetId(7, 15)
cells.InsertNextCell(hexa)
grid.SetCells(vtk.VTK_HEXAHEDRON, cells)

array = vtk.vtkDoubleArray()
array.SetName("mydata")
array.SetNumberOfComponents(1)
array.InsertNextValue(2)
array.InsertNextValue(1)
grid.GetCellData().AddArray(array)

array = grid.GetCellData().GetArray("mydata")
print(array.GetValue(0), array.GetValue(1))
print(array.GetRange())

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.SetNumberOfColors(2)
lut.Build()

# mapper
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(grid)
mapper.SetScalarModeToUseCellFieldData()
mapper.SelectColorArray("mydata")
mapper.SetScalarRange(array.GetRange())
mapper.SetLookupTable(lut)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.VertexVisibilityOn()
prop.SetVertexColor(0.5, 1, 0.5)
prop.SetPointSize(5)

prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

# scalar bar
scalar_bar = vtk.vtkScalarBarActor()
scalar_bar.SetLookupTable(lut)

scalar_bar_widget = vtk.vtkScalarBarWidget()
scalar_bar_widget.SetInteractor(inter)
scalar_bar_widget.SetScalarBarActor(scalar_bar)
scalar_bar_widget.On()

ren_win.Render()
inter.Initialize()
inter.Start()

カラーバーは widget の形で追加しており、モデルとは独立して存在する。ドラッグで位置を変えることもできる。

座標軸の表示

座標軸はカラーバー同様 widget として追加する。

axes = vtk.vtkAxesActor()
axes_widget = vtk.vtkOrientationMarkerWidget()
axes_widget.SetOutlineColor(1, 1, 1)
axes_widget.SetOrientationMarker(axes)
axes_widget.SetInteractor(inter)
axes_widget.SetViewport(0., 0., 0.4, 0.4)
axes_widget.EnabledOn()
axes_widget.InteractiveOff()

座標軸の大きさは SetViewport() で調整できる。SetViewport(0., 0., 0.4, 0.4) というのは、画面の (0, 0)-(0.4, 0.4) の範囲で表示するという意味である (画面全体は (0, 0)-(1, 1) の範囲)。もし右端に表示したければ、SetViewport(0.6, 0., 1., 0.4) などとすればいい。InteractiveOff() の代わりに InteractiveOn() とすると、座標軸を動かせるようになる。

注意: 関数の中で上記を書く場合には注意が必要。vtkOrientationMarkerWidget のオブジェクトを保持するオブジェクトがないので、関数を抜けた途端このオブジェクトは消滅してしまう。対策としては、おそらくクラスのメンバー関数の中で処理を書くことになるはずなので、クラスのメンバー変数に入れてしまう。

参考: Using vtkOrientationMarkerWidget with QVTKRenderWindowInteractor [PyQt4/PySide] (Stack Overflow)

ラベルの色などを変えたい場合は、次のようにする。

prop = vtk.vtkTextProperty()
prop.SetColor(0, 0, 0)
prop.BoldOn()
axes.GetXAxisCaptionActor2D().SetCaptionTextProperty(prop)
axes.GetYAxisCaptionActor2D().SetCaptionTextProperty(prop)
axes.GetZAxisCaptionActor2D().SetCaptionTextProperty(prop)

OpenFOAM ケースの読み込み

OpenFOAM の計算結果を読み込む。ここでは cavity ケースを用いる。

of_reader.py

import vtk
from math import sqrt

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

case = "cavity"
filename = "{}/system/controlDict".format(case)

# reader
reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()

#n = reader.GetTimeValues().GetNumberOfValues()
#latest_time = reader.GetTimeValues().GetValue(n-1)
latest_time = reader.GetTimeValues().GetRange()[1]

reader.UpdateTimeStep(latest_time)
reader.Update()

block = reader.GetOutput()
grid = block.GetBlock(0)

field = "U"
data_type = "cell"

if field == "p":
    if data_type == "cell":
        p = grid.GetCellData().GetArray("p")
    elif data_type == "point":
        p = grid.GetPointData().GetArray("p")
    array_name = "p"
    array = p
elif field == "U":
    magU = vtk.vtkDoubleArray()
    magU.SetName("magU")
    magU.SetNumberOfComponents(1)

    if data_type == "cell":
        U = grid.GetCellData().GetArray("U")

        for i in range(U.GetNumberOfTuples()):
            Ux = U.GetComponent(i, 0)
            Uy = U.GetComponent(i, 1)
            Uz = U.GetComponent(i, 2)
            magU.InsertNextValue(sqrt(Ux*Ux + Uy*Uy + Uz*Uz))

        grid.GetCellData().AddArray(magU)
    elif data_type == "point":
        U = grid.GetPointData().GetArray("U")

        for i in range(U.GetNumberOfTuples()):
            Ux = U.GetComponent(i, 0)
            Uy = U.GetComponent(i, 1)
            Uz = U.GetComponent(i, 2)
            magU.InsertNextValue(sqrt(Ux*Ux + Uy*Uy + Uz*Uz))

        grid.GetPointData().AddArray(magU)
    else:
        assert(0)

    array_name = "magU"
    array = magU
else:
    assert(0)

# filter
geom_filter = vtk.vtkGeometryFilter()
geom_filter.SetInputData(reader.GetOutput())

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

# mapper
mapper = vtk.vtkCompositePolyDataMapper()
mapper.SetInputConnection(geom_filter.GetOutputPort())
if data_type == "cell":
    mapper.SetScalarModeToUseCellFieldData()
elif data_type == "point":
    mapper.SetScalarModeToUsePointFieldData()
mapper.SelectColorArray(array_name)
mapper.SetScalarRange(array.GetRange())
mapper.SetLookupTable(lut)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

# scalar bar
scalar_bar = vtk.vtkScalarBarActor()
scalar_bar.SetOrientationToVertical()
scalar_bar.SetLookupTable(lut)
prop = vtk.vtkTextProperty()
prop.SetColor(0, 0, 0)
prop.SetFontSize(20)
prop.SetFontFamilyToArial()
prop.ItalicOff()
prop.BoldOff()
scalar_bar.UnconstrainedFontSizeOn()
scalar_bar.SetTitleTextProperty(prop)
scalar_bar.SetLabelTextProperty(prop)
scalar_bar.SetLabelFormat("%5.2f")
scalar_bar.SetTitle(array_name)
ren.AddActor(scalar_bar)

# text
text = vtk.vtkTextActor()
text.SetInput(case)
text.GetTextProperty().SetColor(0, 0, 0)
text.GetTextProperty().SetFontSize(24)
text.SetPosition(300, 30)
ren.AddActor(text)

ren_win.Render()
inter.Initialize()
inter.Start()

ケースは case で指定している。時刻データがいくつかある可能性があるが、最後の時刻を取り出している。filed に "p"、"U" を指定して、圧力か速度を表示する。速度はベクトルなので、そのノルム magU を作っている。data_type で "cell" か "point" を指定してセル値と節点値を切り替えている。

カラーバーは、今回は widget ではなく直接 renderer に渡している。下部に vtkTextActor でケース名を表示している。テキストを widget にすることもできるが、テキストを動かしたい機会はあまりなさそう。

OpenFOAM ケースのブロック情報の取得

OpenFOAM のケースのブロックの情報 (パッチ) を得る。

get_of_patches.py

import sys
import vtk

if len(sys.argv) < 2:
    print("usage: get_of_patches.py case")
    exit(0)

case = sys.argv[1]
filename = "{}/system/controlDict".format(case)

reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()

n = reader.GetNumberOfPatchArrays()
for i in range(n):
    name = reader.GetPatchArrayName(i)
    print(name)

実行例

$ python get_of_patches.py cavity
internalMesh
movingWall
fixedWalls
frontAndBack

OpenFOAM のマルチリージョンケースの読み込み

OpenFOAM のマルチリージョンケースの計算結果を読み込む。ここでは multiRegionHeater ケースを用いる。

of_reader_multi.py

import sys
import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

case = "multiRegionHeater"
filename = "{}/system/controlDict".format(case)

# reader
reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()

latest_time = reader.GetTimeValues().GetRange()[1]
print("Latest time:", latest_time)

reader.UpdateTimeStep(latest_time)
reader.Update()


n = reader.GetNumberOfPatchArrays()
for i in range(n):
    name = reader.GetPatchArrayName(i)
    print(name)

enabled_block = [
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
]

reader.DisableAllPatchArrays()
for b in enabled_block:
    reader.SetPatchArrayStatus(b, 1)
reader.Update()


field = "T"
data_type = "point"

block = reader.GetOutput()

array_name = field
array_min = sys.float_info.max
array_max = sys.float_info.min

itr = block.NewIterator()
itr.InitTraversal()
while not itr.IsDoneWithTraversal():
    grid = itr.GetCurrentDataObject()

    if data_type == "cell":
        T = grid.GetCellData().GetArray(field)
    elif data_type == "point":
        T = grid.GetCellData().GetArray(field)
    else:
        assert(0)

    T_min, T_max = T.GetRange()

    if T_min < array_min:
        array_min = T_min
    if T_max > array_max:
        array_max = T_max

    itr.GoToNextItem()


# filter
geom_filter = vtk.vtkGeometryFilter()
geom_filter.SetInputData(reader.GetOutput())

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

# mapper
mapper = vtk.vtkCompositePolyDataMapper()
mapper.SetInputConnection(geom_filter.GetOutputPort())
if data_type == "cell":
    mapper.SetScalarModeToUseCellFieldData()
elif data_type == "point":
    mapper.SetScalarModeToUsePointFieldData()
mapper.SelectColorArray(array_name)
mapper.SetScalarRange(array_min, array_max)
mapper.SetLookupTable(lut)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

# scalar bar
scalar_bar = vtk.vtkScalarBarActor()
scalar_bar.SetOrientationToVertical()
scalar_bar.SetLookupTable(lut)
prop = vtk.vtkTextProperty()
prop.SetColor(0, 0, 0)
prop.SetFontSize(20)
prop.SetFontFamilyToArial()
prop.ItalicOff()
prop.BoldOff()
scalar_bar.UnconstrainedFontSizeOn()
scalar_bar.SetTitleTextProperty(prop)
scalar_bar.SetLabelTextProperty(prop)
scalar_bar.SetLabelFormat("%5.2f")
scalar_bar.SetTitle(array_name)
ren.AddActor(scalar_bar)

# text
text = vtk.vtkTextActor()
text.SetInput(case)
text.GetTextProperty().SetColor(0, 0, 0)
text.GetTextProperty().SetFontSize(24)
text.SetPosition(200, 30)
ren.AddActor(text)

ren_win.Render()
inter.Initialize()
inter.Start()

以下では、一旦すべてのブロックを無効にして、表示させたいものだけ有効にしている。

enabled_block = [
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
]

reader.DisableAllPatchArrays()
for b in enabled_block:
    reader.SetPatchArrayStatus(b, 1)
reader.Update()

イテレータでブロックの巡回をしているところがある。

itr = block.NewIterator()
itr.InitTraversal()
while not itr.IsDoneWithTraversal():
    grid = itr.GetCurrentDataObject()
    ...
    itr.GoToNextItem()

これは SetPachArrayStatus() で有効にされたブロックだけ巡回する。

ブロックの抽出

上の例では、ブロックの選択を reader の SetPatchArrayStatus() で行なったが、これを vtkExtractBlock で行うこともできる。ParaView で、OpenFOAM リーダーの設定でブロックを選択するか、Extract Block フィルタで選択するかというのと同じである。

vtkExtracBlock でブロックを抽出するには、インデックスで指定するしかないが、わかりにくいので、以下の例では一度パッチ名とインデックスの対応テーブルを作っている。

extract_block.py

import sys
import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

case = "multiRegionHeater"
filename = "{}/system/controlDict".format(case)

# reader
reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()

latest_time = reader.GetTimeValues().GetRange()[1]
print("Latest time:", latest_time)

reader.UpdateTimeStep(latest_time)
reader.Update()


#
# create patch to index table
#
reader.EnableAllPatchArrays()
reader.Update()

block_index = {}
index = 1
for i in range(reader.GetNumberOfPatchArrays()):
    name = reader.GetPatchArrayName(i)
    r = re.match(".*internalMesh$", name)
    if r is not None:
        index += 1
    block_index[name] = index
    if r is not None:
        index += 1
    index += 1

for name in block_index:
    print(name, block_index[name])


enabled_block = [
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
]

extract_block = vtk.vtkExtractBlock()
extract_block.SetInputData(reader.GetOutput())
for b in enabled_block:
    extract_block.AddIndex(block_index[b])
extract_block.Update()


field = "T"
data_type = "cell"

block = extract_block.GetOutput()

array_name = field
array_min = sys.float_info.max
array_max = sys.float_info.min

itr = block.NewIterator()
itr.InitTraversal()
while not itr.IsDoneWithTraversal():
    grid = itr.GetCurrentDataObject()

    if data_type == "cell":
        T = grid.GetCellData().GetArray(field)
    elif data_type == "point":
        T = grid.GetPointData().GetArray(field)
    else:
        assert(0)

    T_min, T_max = T.GetRange()

    if T_min < array_min:
        array_min = T_min
    if T_max > array_max:
        array_max = T_max

    itr.GoToNextItem()


# filter
geom_filter = vtk.vtkGeometryFilter()
geom_filter.SetInputConnection(extract_block.GetOutputPort())

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

# mapper
mapper = vtk.vtkCompositePolyDataMapper()
mapper.SetInputConnection(geom_filter.GetOutputPort())
if data_type == "cell":
    mapper.SetScalarModeToUseCellFieldData()
elif data_type == "point":
    mapper.SetScalarModeToUsePointFieldData()
mapper.SelectColorArray(array_name)
mapper.SetScalarRange(array_min, array_max)
mapper.SetLookupTable(lut)

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

# scalar bar
scalar_bar = vtk.vtkScalarBarActor()
scalar_bar.SetOrientationToVertical()
scalar_bar.SetLookupTable(lut)
prop = vtk.vtkTextProperty()
prop.SetColor(0, 0, 0)
prop.SetFontSize(20)
prop.SetFontFamilyToArial()
prop.ItalicOff()
prop.BoldOff()
scalar_bar.UnconstrainedFontSizeOn()
scalar_bar.SetTitleTextProperty(prop)
scalar_bar.SetLabelTextProperty(prop)
scalar_bar.SetLabelFormat("%5.2f")
scalar_bar.SetTitle(array_name)
ren.AddActor(scalar_bar)

# text
text = vtk.vtkTextActor()
text.SetInput(case)
text.GetTextProperty().SetColor(0, 0, 0)
text.GetTextProperty().SetFontSize(24)
text.SetPosition(200, 30)
ren.AddActor(text)

ren_win.Render()
inter.Initialize()
inter.Start()

実行結果は of_reader_multi.py と同じである。

OpenFOAM のケースを読んだ時、reader.GetOutpu() で得られるのは vtkMultiBlockDataSet のオブジェクトである。これから vtkUnstructuredGrid を取り出して、セル値の表示の設定などを行う。of_reader.py の例では、ただちに GetBlock(0) でグリッドデータを取り出しているが、一般的には、マルチブロックデータは階層構造になっている。vtkExtractBlock ではこの構造に flat index でアクセスできるのだが、どのブロックがどの index なのかが見えない。OpenFOAM ケースの場合は、次のようになっているようである。

Multi-block Dataset 0
  defaultRegion 1
    internalMesh 2
    Patches 3
      maxY 4
      minX 5
      maxX 6
      ...

これは ParaView の Extract Block を見ればわかる。ParaView では vtkCompositeIndex というものを表示できるが、それがここで言っているブロックの index である。ルートのインデックスが 0、次のレペルで vtkMultiBlockDataSet のリージョンブロックがあり、その下に vtkUnstructuredGrid である internalMesh、その後に vtkMultiBlockDataSet の Patches、その下に vtkPolyData のパッチが並んでいる。名前は ParaView が勝手につけいているようである。データセットには名前情報はなさそう。

block_index の作成部分は、ブロックの構造を辿りながら処理をするとしたら、次のように書ける。

patch = []
n = reader.GetNumberOfPatchArrays()
for i in range(n):
    name = reader.GetPatchArrayName(i)
    patch.append(name)

block = reader.GetOutput()
block_index = {}
patchi = 0
index = 1
for i in range(block.GetNumberOfBlocks()):
    subblock = block.GetBlock(i)
    index += 1
    for j in range(subblock.GetNumberOfBlocks()):
        subsubblock = subblock.GetBlock(j)
        if isinstance(subsubblock, vtk.vtkMultiBlockDataSet):
            index += 1
            for k in range(subsubblock.GetNumberOfBlocks()):
                block_index[patch[patchi]] = index
                patchi += 1
                index += 1
        else:
            block_index[patch[patchi]] = index
            patchi += 1
            index += 1

注意: 時刻の設定を reader.UpdateTimeStep() で行うが、これは vtkExtractBlock のほうには反映されない。もし vtkExtractBlock でブロックを抽出したあとに reader.UpdateTimeStep() 呼ぶ場合は、一緒に vtkExtractBlock のほうの UpdateTimeStep() も呼び出す必要がある。

ハイライト表示

注目しているブロックをハイライト表示したい。対象のブロックだけ Extract Block で取り出して、色をつければよい。

highlight.py

import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

case = "multiRegionHeater"
filename = "{}/system/controlDict".format(case)

# reader
reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()


#
# create patch to index table
#
reader.EnableAllPatchArrays()
reader.Update()

block_index = {}
index = 1
for i in range(reader.GetNumberOfPatchArrays()):
    name = reader.GetPatchArrayName(i)
    r = re.match(".*internalMesh$", name)
    if r is not None:
        index += 1
    block_index[name] = index
    if r is not None:
        index += 1
    index += 1

for name in block_index:
    print(name, block_index[name])


enabled_block = [
    "bottomWater/internalMesh",
    "topAir/internalMesh",
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
]

highlight_block = [
    "topAir/minX",
]

extract_block = vtk.vtkExtractBlock()
extract_block.SetInputData(reader.GetOutput())
for b in enabled_block:
    extract_block.AddIndex(block_index[b])
extract_block.Update()

extract_block2 = vtk.vtkExtractBlock()
extract_block2.SetInputData(reader.GetOutput())
for b in highlight_block:
    extract_block2.AddIndex(block_index[b])
extract_block2.Update()


#
# block index
#
block = extract_block.GetOutput()

array_name = "index"
array_min = min(block_index.values())
array_max = max(block_index.values())

itr = block.NewIterator()
itr.InitTraversal()
i = 0
while not itr.IsDoneWithTraversal():
    grid = itr.GetCurrentDataObject()

    index = vtk.vtkIntArray()
    index.SetName(array_name)
    index.SetNumberOfComponents(1)

    for j in range(grid.GetNumberOfCells()):
        index.InsertNextValue(block_index[enabled_block[i]])
    grid.GetCellData().AddArray(index)

    itr.GoToNextItem()
    i += 1


# filter
geom_filter = vtk.vtkGeometryFilter()
geom_filter.SetInputConnection(extract_block.GetOutputPort())

geom_filter2 = vtk.vtkGeometryFilter()
geom_filter2.SetInputConnection(extract_block2.GetOutputPort())

# lookup table
lut = vtk.vtkLookupTable()
lut.SetHueRange(0.667, 0)
lut.Build()

# mapper
mapper = vtk.vtkCompositePolyDataMapper()
mapper.SetInputConnection(geom_filter.GetOutputPort())
mapper.SetScalarModeToUseCellFieldData()
mapper.SelectColorArray(array_name)
mapper.SetScalarRange(array_min, array_max)
mapper.SetLookupTable(lut)

mapper2 = vtk.vtkCompositePolyDataMapper()
mapper2.SetInputConnection(geom_filter2.GetOutputPort())
mapper2.SetScalarModeToUseCellFieldData()

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.BackfaceCullingOn()
prop.SetAmbient(0.5)

prop.EdgeVisibilityOn()
prop.SetEdgeColor(0, 0, 0)
prop.SetLineWidth(2)

prop.SetRepresentationToSurface()

actor2 = vtk.vtkActor()
actor2.SetMapper(mapper2)

prop = actor2.GetProperty()
prop.SetColor(1, 0, 1)
prop.SetOpacity(0.8)

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.AddActor(actor2)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

ren_win.Render()
inter.Initialize()
inter.Start()

見栄えのためにブロックを色分けしている。組み込みの機能でできそうではあるが、やりかたがよくわからなかった。

特徴線

特徴線は vtkFeatureEdges で描画できる。OpenFOAM のマルチリージョンケースで特徴線を描画してみる。

feature_edge.py

import vtk

class MouseInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
    def __init__(self, parent=None):
        self.AddObserver("LeftButtonPressEvent", self.leftButtonPressEvent)

    def leftButtonPressEvent(self, obj, event):
        self.OnLeftButtonDown()

case = "multiRegionHeater"
filename = "{}/system/controlDict".format(case)

# reader
reader = vtk.vtkOpenFOAMReader()
reader.SetFileName(filename)
reader.Update()

enabled_block = [
    "bottomWater/internalMesh",
    "topAir/internalMesh",
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
]

reader.DisableAllPatchArrays()
for b in enabled_block:
    reader.SetPatchArrayStatus(b, 1)
reader.Update()

# filter
geom_filter = vtk.vtkGeometryFilter()
geom_filter.SetInputData(reader.GetOutput())

feature_edge = vtk.vtkFeatureEdges()
feature_edge.SetInputConnection(geom_filter.GetOutputPort())

# mapper
mapper = vtk.vtkCompositePolyDataMapper()
mapper.SetInputConnection(feature_edge.GetOutputPort())
mapper.SetScalarModeToUseCellFieldData()

# actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)

prop = actor.GetProperty()
prop.SetColor(0, 0, 0)
prop.SetLineWidth(2)
prop.SetRepresentationToSurface()

# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.GradientBackgroundOn()
ren.SetBackground(1, 1, 1)
ren.SetBackground2(0.4, 0.55, .75)

ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetSize(640, 480)

# interactor
inter = vtk.vtkRenderWindowInteractor()
inter.SetRenderWindow(ren_win)
inter.SetInteractorStyle(MouseInteractorStyle())

ren_win.Render()
inter.Initialize()
inter.Start()

参考