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 Source

Programmable Source は、新しいデータのソースとなったり、ファイルからデータを読み込んだりするために使用されます。これはPythonで書かれています。後で、このデータを修正するためにフィルタが使用されます。Programmable Source の例としては、.csvファイルを paraview に読み込むことがあります。

• data.csv という名前の .csv ファイルを作成します。このファイルに以下を記述します。

  x coord, y coord, z coord, scalar
  0,0,0,0
  1,0,0,1
  0,1,0,2
  1,1,0,2
  -0.5,-0.5,1,4
  -0.5,-0.5,1,5
  -0.5,-0.5,1,6
  -0.5,-0.5,1,7
  

• Sources → Programmable Source を選択します。

• Output Dataset Type を vtkTable に変更します。

• スクリプトウィンドウに以下を入力します。

  import numpy as np
  data = np.genfromtxt("c:/.../../data.csv", dtype=None, names=True,
  delimiter=',', autostrip=True)
  for name in data.dtype.names:
     array = data[name]
     output.RowData.append(array,name)
  

• Filters → Alphabetical → Table to Points filter を選択します。

• X Column == x_coord.

• Y Column == y_coord.

• Z Column == z_coord.

• Apply を実行します。

• 3Dビューを選択し、表示状態をオンにします。

ParaViewユーザーガイドに掲載されているプログラマブルソースの例をいくつか紹介します：

• CSVファイルの読み込み

• CSVファイルシリーズの読み込み

• パーティクルによるCSVファイルの読み込み

• 2値化された2次元画像の読み込み

• ヘリックスソース

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

13.4. 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 に選択します。

  • 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 に関連するブログ記事から興味深いものを選んでみました。

• https://blog.kitware.com/improved-vtk-numpy-integration/

• https://blog.kitware.com/improved-vtk-numpy-integration-part-2/

• https://blog.kitware.com/improved-vtk-numpy-integration-part-3/

• https://blog.kitware.com/mpi4py-and-vtk/

• https://blog.kitware.com/improved-vtk-numpy-integration-part-4/

• https://blog.kitware.com/improved-vtk-numpy-integration-part-5/