paraview からスカラやベクトルのデータを直接読みだすために servermanager を使う

18165 ワード
Python Visualization paraview Mesh CAE Python テキストリンク

始めに

シミュレーションデータの可視化に使われる ParaView で、フィールドのデータを python スクリプトから読みだす方法は初見ではわかりにくく、なかでも servermanager を介して配列に直接アクセスする方法についての情報が見つけにくかったので、servermanager を使って動作するコードを書いてみました。
ここでのポイントは、
* servermanager を呼び出すことで、PointData, CellData にアクセスできるようになる
* paraview.vtk.numpy_interface.dataset_adapter を用いて、numpy でアクセスできる配列が生成できる
という2点です。
全体の骨組みとしては、ParaView の [Tools] - [Start Trace] で記録したコードを下敷きにしました。コードは、Windows10 （64 ビット）上の paraview-5.5.1 で動作を確認しています。

情報源

  • dataset_adapter については、The ParaView Guide の Chapter 14 に詳しく書かれています。
  • servermanager については、cfd-online の投稿 を参考にしました。

デモコードの内容

既存のフィルターのままでは実現できない処理として、以下のような課題を実行します。

CellData としてあるスカラ値が入っているときに、このスカラ値の範囲ごとに振り分けて、該当するセル面積の和を求めてヒストグラムのように表示すること。

python のコードで直接スカラ値を読んで、どの範囲に該当するか判別して、該当する面積をインクリメントしていく処理です。
ParaView での集計結果は csv ファイルに書き出し、ヒストグラムとしての描画は、外部で実施することとします。
ちなみに以下のコードでは、スカラ値そのものも流速分布から作っているので、処理の流れは次のようになっています。

  • 例題データを読み込む
    • スライスする
    • スカラ値として、スライス面に平行な速度の大きさを計算する。（平行な速度成分を選んだのに、特段の意味はありませんが、例示として）
    • ポイントデータをセルデータに変換する

ここまで行っているのは、例示用のデータの作成です。
CellData として、振り分けの基準に使うスカラ値が入っているデータセットがあれば、それを使えば十分です。（そもそもそれが、個人的な本来の目的でした）

  • CellSize フィルタを使って、各セルの面積を求める
  • 結果の配列を python で読み出して集計する
  • numpy の機能を使って、CSV ファイルを保存。
  • 確認のため、注目領域だけが見えるような ParaView 設定を残しておく
  • グラフ描画は外部で行います。

コードの実行方法

ParaView を立ち上げて、[View] - [Python Shell] で シェルを開き、シェルのウィンドウ下側のRun Script ボタンから、デモコードのファイルを選んで実行します。
7行目、8行目に、データファイルと出力ファイルの名前をハードコードしているので、ここを適宜修正してください。また、ヒストグラムの上限 dmaxとビン幅 step も終盤でハードコードしています。

# based on a trace generated using paraview version 5.5.1-4

from paraview.simple import *
import paraview.vtk.numpy_interface.dataset_adapter as dsa
import numpy as np

flowfile = r'C:\Users\Your_datadir\cyl_flow.vtu'
histfile = r'C:\Users\Your_destdir\hist.csv'

paraview.simple._DisableFirstRenderCameraReset()

cyl_flowvtu = XMLUnstructuredGridReader(FileName=[flowfile])
renderView1 = GetActiveViewOrCreate('RenderView')

cyl_flowvtuDisplay = Show(cyl_flowvtu, renderView1)
Hide(cyl_flowvtu, renderView1)

slice1 = Slice(Input=cyl_flowvtu)
slice1.SliceType.Origin = [0.015, 0.0, 0.01]
slice1Display = Show(slice1, renderView1)
slice1Display.Representation = 'Surface'
Hide3DWidgets(proxy=slice1.SliceType)
Hide(slice1, renderView1)

surfaceVectors1 = SurfaceVectors(Input=slice1)
surfaceVectors1Display = Show(surfaceVectors1, renderView1)
surfaceVectors1Display.Representation = 'Surface'
Hide(surfaceVectors1, renderView1)

pointDatatoCellData1 = PointDatatoCellData(Input=surfaceVectors1)
pointDatatoCellData1Display = Show(pointDatatoCellData1, renderView1)
pointDatatoCellData1Display.Representation = 'Surface'
Hide(pointDatatoCellData1, renderView1)

calculator1 = Calculator(Input=pointDatatoCellData1)
calculator1.AttributeType = 'Cell Data'
calculator1.ResultArrayName = 'v_para'
calculator1.Function = 'mag(VELOCITY)'
calculator1Display = Show(calculator1, renderView1)
calculator1Display.Representation = 'Surface'
calculator1Display.SetScalarBarVisibility(renderView1, True)
Hide(calculator1, renderView1)
# ここまでで、スカラ値 v_para が CellData に入る。Spread Sheet View で確認するとわかりやすい

# 各セルの面積を計算させる。
cellSize1 = CellSize(Input=calculator1)
cellSize1Display = Show(cellSize1, renderView1)
cellSize1Display.Representation = 'Surface'
cellSize1Display.SetScalarBarVisibility(renderView1, True)

renderView1.Update()

# ここからが、データを読みだすコード
cellsizeobj = servermanager.Fetch(cellSize1)
m_data = dsa.WrapDataObject(cellsizeobj)

# VTKArray として配列を取り出す (numpy の Array として振舞う)
area = m_data.GetCellData().GetArray('Quality')
v_para = m_data.GetCellData().GetArray('v_para')

# ヒストグラムの上限とビンの幅は、データを見て適切に選ぶ。
dmax, step = 16, 2
div = dmax / step + 1
areabin = np.zeros(div)
binlabel = np.linspace(0., dmax, div)
for v, a in zip(v_para, area):
    binnum = int(np.ceil(v / step))
    areabin[binnum] += a

hist = np.vstack([binlabel, areabin]).T

np.savetxt(histfile, hist, fmt="%f", delimiter=",")

# スカラ値が小さいところは、Opacity を小さくして見えなくする。ここでは、2 以下を隠している
# set scalar coloring
ColorBy(cellSize1Display, ('CELLS', 'v_para'))

v_paraLUT = GetColorTransferFunction('v_para')
v_paraLUT.AutomaticRescaleRangeMode = 'Never'
v_paraLUT.EnableOpacityMapping = 1
v_paraLUT.RescaleTransferFunction(0.0, dmax)
v_paraPWF = GetOpacityTransferFunction('v_para')
v_paraPWF.RescaleTransferFunction(0.0, dmax)

v_paraPWF.Points = [0.0, 0.0, 0.5, 0.0, 2., 0.0, 0.5, 0.0, 8., 1.0, 0.5, 0.0, dmax, 1.0, 0.5, 0.0, 16., 0.0, 0.5, 0.0]
v_paraLUT.RGBPoints = [0.0, 0.231373, 0.298039, 0.752941, 2.0, 0.865003, 0.865003, 0.865003, dmax, 0.705882, 0.0156863, 0.14902]

renderView1.CameraPosition = [0.093222, -0.0269241, 0.0564596112]
renderView1.CameraFocalPoint = [-5.4620, 1.5775, -3.3081]
renderView1.CameraViewUp = [0.16518, 0.96800, 0.188874]

結果

cellSize の結果を ParaView からスプレッドシートとして吐き出して、手元で集計するよりも手軽で時間のかからない処理ができるようになりました。