13. Python & Batch: Python Calculator, Programmable Source, Filter

13.1. はじめに

ParaView には、基礎となるVTKライブラリと同様にpythonの数学関数にアクセスするための3つのフィルタがあります。これらは、Python Calculator、Programmable Source、Filterです。このチュートリアルでは、これらについて説明します。

13.2. Python Calculator

Python Calculatorは、Pythonで利用可能な計算を適用することができます。これらは、セルの体積を取得し、セルの面積を取得し、クロス積、ドット積、カールなどを取得するような関数が含まれています。計算式はすべて1行に収まらなければなりません。Python Calculator を例として見てみましょう。

新しいポイント変数を作成し、5をロードすることができます。
- can.ex2を開きます。
- すべての変数をオンにします。
- Apply を実行します。
- Filters → Alphabetical → Python Calculator.
- Expressionを 5 に変更します。
- Array Association を Point Data に設定します。
- Array Name を Calculated Variable に変更します。
- Apply を実行します。
- Calculated Variable で色付けします。
変位量の2倍に相当する変数を作成しましょう。
- Expressionを DISPL*2 に変更します。
- Apply を実行します。
- 最後のタイムステップへ移動します。
- データ範囲に再スケールします。
- 注：変位データにアクセスする、より完全な方法。フィルタへの最初の入力から明示的に DISPL を引き出しています。最初の入力は[0]、2番目は[1]、などであることを忘れないでください。
- Expressionを (inputs[0].PointData['DISPL']) * 2 に変更します。
- Apply を実行します。
- データ範囲に再スケールします。
ここでは、ベクトルとグローバル変数の掛け算の方法を説明します。DISPLにタイムステップ(TMSTEP)を掛けてみます。
- Expressionを inputs[0].FieldData['NSTEPS'][time_index] * DISPL に変更します。
- Apply を実行します。
- データ範囲に再スケールします。
- 要約すると、関数が変数の入力を必要とする場合、上記の文字列を使用します。関数が入力メッシュを必要とする場合は、inputs[]を使用します。
Did you know?

Pythonでは、配列には角括弧を、数式の一部をグループ化するには括弧を使用します。また、文字列（変数名）を指定する場合は、シングルクォーテーションまたはダブルクォーテーションを使用します。
セルの体積を格納する変数を作成しましょう。
- Expressionを volume(inputs[0]) に変更します。
- Apply を実行します。
- データ範囲に再スケールします。
多数の関数を1つの式にまとめることができます。例えば、カールの発散のsinを取るというのはナンセンスな表現です。以下はその式です。
- Expressionを sin(divergence(curl('ACCL'))) に変更します。Apply -> Rescale to Data Rangeを実行します。

Python Calculatorから利用できる関数が興味深いです。

area(dataset)
aspect(dataset)
cos(array)
cross(X,Y) where X and Y are two 3D vector arrays
curl(array)
divergence(array)
dot(a1,a2)
eigenvalue and eigenvector(array)
gradient(array)
max(array)
mean(array)
min(array)
norm(array)
sin(array)
strain(array)
volume(array)
vorticity(array)

完全なリストはユーザーズガイドの 5.9.3 章で見ることができます。

13.3. Programmable Filter

Programmable Filterは、Pythonを使ってパイプラインのデータを変更するために使用されます。変数を2で割るProgrammable Filterの例。

ACCL を2で割ります。
- can.ex2を開きます。
- すべての変数をオンにします。
- Apply を実行します。
- Filters → Alphabetical → Programmable Filter を選択します。
- 出力データを残します。
- Typeを Same as Input に設定します。
- スクリプトウィンドウに以下を入力します。
```
input0 = inputs[0]
dataArray = input0.PointData["ACCL"] / 2.0
output.PointData.append(dataArray, "ACCL_half")
```
- Coloring を ACCL_half で設定する。
2つのデータセットから互いに引き算をします。disk_out_ref.ex2の2つのインスタンスを2つのデータセットとして使用し、AsH3 から GaMe3 を引き算する。
- disk_out_ref.ex2を開きます。
- すべての変数をオンにします。
- Apply を実行します。
- disk_out_ref.ex2を開きます。
- すべての変数をオンにします。
- Apply を実行します。
- 両方のデータセットをハイライト表示します。マウスでキーを使って、複数の入力を選択します。
- Filters → Alphabetical → Programmable Filter を選択します。
- 出力データを残します。
- Typeを Same as Input に設定します。
- スクリプトウィンドウに以下を入力します。
```
v_0 = inputs[0].PointData['AsH3']
v_1 = inputs[1].PointData['GaMe3']
output.PointData.append(v_1 - v_0, 'difference')
```
- Coloring を difference に設定します。

テンソルを作成する

from paraview.vtk.numpy_interface import dataset_adapter as dsa
import numpy
def make_tensor(xx,yy,zz, xy, yz, xz):

      t = numpy.vstack([xx,yy,zz,xy, yz,
      xz]).transpose().view(dsa.VTKArray)
      t.DataSet = xx.DataSet
      t.Association = xx.Association
      return t

xx = inputs[0].PointData["sigma_xx"]
yy = inputs[0].PointData["sigma_yy"]
zz = inputs[0].PointData["sigma_zz"]
xy = inputs[0].PointData["sigma_xy"]
yz = inputs[0].PointData["sigma_yz"]
xz = inputs[0].PointData["sigma_xz"]
output.PointData.append(make_tensor(xx,yy,zz,xy,yz,xz), "tensor")

2つのタイムステップを互いに引き算する
- can.ex2の読み込み
- Filters → Alphabetical → Force Time を選択します。
- タイムステップには、変位をゼロにしたいものを設定します。
- このフィルタの出力は、選択したタイムステップで "凍結された" can.ex2 データセットです。これは、時間を進めても変化しない。
- パイプラインブラウザで、can.ex2 (CTRL キーを使用) と ForceTime1 ソースの両方を選択します。元のcan.ex2ソースはタイムステップの変更に伴って更新されますが、ForceTime1 ソースは変更されないことに注意してください。
- Filters → Alphabetical → Python Calculator に選択します。
- Substract the displacement in the "frozen" dataset from the current timestep in can.ex2.
- Expressionを inputs[0].PointData['DISPL'] - inputs[1].PointData['DISPL'] にセットします。
- Python Calculator への入力の順序はうまく定義されていないので、結果の符号を正しくするためにinputs[0]とinputs[1]のインデックスを入れ替える必要があるかもしれませんが、これでうまくいくはずです。
- (必要であれば、) 最後に、Plot Selection over Timeフィルタを追加します。このフィルタは、すべてのタイムステップにわたって実行され、現在のタイムステップのデータから、時間強制フィルタで生成された "凍結" タイムステップのデータを引き、その結果をグラフにプロットします。最初のタイムステップの値は0であるべきです。

固有ベクトルを読み、固有値を計算し、3変数に入れる

# ParaView Programmable Filter script.  ParaView 5.0.1.
#
# This code will read in an eigenvector, calculate an
#    eigenvalue, and place it into three variables.
#
# Be sure to correct the input variables below.  Also, note that
#   the code uses ZX, not XZ.
#
# This code only works with any multiblock or vtk
#   datasets (including ones with only one block - i.e.,
#   Exodus datasets as input).
#
# Usage:  Run the Programmable Filter.
#   Cut and paste this file in the section named "Script".
#   Leave "Output Data Set Type" as "Same as Input".
#   Click Apply button
#
# Written by Jeff Mauldin and Alan Scott
#

import numpy as np

def process_composite_dataset(input0):
  # Pick up input arrays
  xxar = input0.CellData["EPSXX"]
  xyar = input0.CellData["EPSXY"]
  zxar = input0.CellData["EPSZX"]
  yyar = input0.CellData["EPSYY"]
  yzar = input0.CellData["EPSYZ"]
  zzar = input0.CellData["EPSZZ"]

  #print `xxar`
  #print len(xxar.Arrays)

  # Set output arrays to same type as input array.
  # Do a multiply to make sure we don't just have a
  # pointer to the original.
  outarray0 = xxar*0.5
  outarray1 = xxar*0.5
  outarray2 = xxar*0.5


  # Run a for loop over all blocks
  numsubarrays = len(xxar.Arrays)
  for ii in range(0, numsubarrays):
    # pick up input arrays for each block.
    xxarsub = xxar.Arrays[ii]
    xyarsub = xyar.Arrays[ii]
    zxarsub = zxar.Arrays[ii]
    yyarsub = yyar.Arrays[ii]
    yzarsub = yzar.Arrays[ii]
    zzarsub = zzar.Arrays[ii]

    #print `xxarsub`

    # Transpose and calculate the principle strain.
    strain = np.transpose(
        np.array(
          [ [xxarsub, xyarsub, zxarsub],
            [xyarsub, yyarsub, yzarsub],
            [zxarsub, yzarsub, zzarsub] ] ),
              (2,0,1))

    principal_strain = np.linalg.eigvalsh(strain)

    # Move principle strain to temp output arrays for this block
    outarray0.Arrays[ii] = principal_strain[:,0]
    outarray1.Arrays[ii] = principal_strain[:,1]
    outarray2.Arrays[ii] = principal_strain[:,2]

  #ps0 = principal_strain[:,0]
  #print "ps0 len: " + str(len(ps0))

  # Finally, move the temp arrays to output arrays
  output.CellData.append(outarray0, "principal_strain_0")
  output.CellData.append(outarray1, "principal_strain_1")
  output.CellData.append(outarray2, "principal_strain_2")


def process_unstructured_dataset(input0):
  # Pick up input arrays
  xxar = input0.CellData["EPSXX"]
  xyar = input0.CellData["EPSXY"]
  zxar = input0.CellData["EPSZX"]
  yyar = input0.CellData["EPSYY"]
  yzar = input0.CellData["EPSYZ"]
  zzar = input0.CellData["EPSZZ"]

  #print `xxar`
  #print len(xxar.Arrays)

  # Transpose and calculate the principle strain.
  strain = np.transpose(
        np.array(
          [ [xxar, xyar, zxar],
            [xyar, yyar, yzar],
            [zxar, yzar, zzar] ] ),
              (2,0,1))

  principal_strain = np.linalg.eigvalsh(strain)

  #ps0 = principal_strain[:,0]
  #print "ps0 len: " + str(len(ps0))

  # Finally, move the temp arrays to output arrays
  output.CellData.append(principal_strain[:,0],
     "principal_strain_0")
  output.CellData.append(principal_strain[:,1],
     "principal_strain_1")
  output.CellData.append(principal_strain[:,2],
     "principal_strain_2")

input0 = inputs[0]

if input0.IsA("vtkCompositeDataSet"):
  process_composite_dataset(input0)
elif input0.IsA("vtkUnstructuredGrid"):
  process_unstructured_dataset(input0)
else:
  print "Bad dataset type for this script"

重要： Programmable Source と Programmable Filter はサーバーレベルで動作するため、paraview.simple はロードまたは使用できません。

Programmable Filter に関するモードの詳細は、リファレンスマニュアルの 5 章に記載されています。

Programmable Filter に関連するブログ記事から興味深いものを選んでみました。