Programable Sourceの使い方

可視化ソフトウェアParaViewでは、"Programmable Source"を使用して、データを作成することができます。

ここでは、"Programmable Source"でImageデータ(Uniform Rectilinear Grid)を作成してみます。データの各点には、Riemann Zeta関数の絶対値を与えます。

Programmable Sourceの場所

Programmable Sourceはメニューバーの

 Sources >> Programmable Source

にあります。

programmable source

Script (RequestInformation)の記述

まず、300x500x1ピクセルの2次元画像を作成します。

Propertiesタブにある"Output Data Set Type"に"vtkImageData"を設定し、"Script (RequestInformation)"に以下を記述します。

Everything is expanded.Everything is shortened.
  1
  2
 
 
from paraview import util
util.SetOutputWholeExtent(self, [0, 299, 0, 499, 0, 0])

 

RequestInformation

Scriptの記述

Propertiesタブにある"Script"には以下を記述します。

Everything is expanded.Everything is shortened.
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
import math
 
lowerThreshold = 1.0e-6
upperBound = 1.0e+4
infinity = 300
 
def binom(n,k):
  m = 1
  if n < 2 * k:
    k = n - k
  for i in range(1, k + 1):
    m = m * (n - i + 1) / i
  return m
 
def zeta(s):
 
  outerSum= 0.0
  prev = 1000000000
 
  for m in range(1, infinity+1):
 
    innerSum = 0.0
    for j in range(1, m+1):
      c1 = 1 if (j-1)%2 == 0 else -1
      c2 = binom(m-1, j-1)
      c3 = j**(-s)
      innerSum += c1*c2*c3
 
    try:
      innerSum = innerSum * (2**(-m)) / (1-2**(1-s))
    except:
      return 0+0j
 
    outerSum += innerSum
 
    if abs(prev - innerSum) < lowerThreshold:
      break
 
    if abs(outerSum) > upperBound:
      break
 
    prev = innerSum
 
  return outerSum
 
 
out = self.GetOutput()
out.SetDimensions(300, 500, 1)
out.SetOrigin(0.0001, 0.0001, 0)
out.SetSpacing(0.1, 0.1, 1)
 
zetaArray = vtk.vtkDoubleArray()
zetaArray.SetName('zeta')
zetaArray.SetNumberOfComponents(1)
zetaArray.SetNumberOfTuples(out.GetNumberOfPoints())
 
dims = out.GetDimensions()
 
for z in range(0, dims[2]):
  for y in range(0, dims[1]):
    for x in range(0, dims[0]):
      pntId = out.ComputePointId([x, y, z])
      point = out.GetPoint(pntId)
      s = point[0] + point[1]*1j
      val = abs(zeta(s))
      zetaArray.SetTuple1(pntId, val)
 
out.GetPointData().AddArray(zetaArray)

これは、画像データ上の点でリーマン・ゼータ関数を求めて、その絶対値をその点に設定します。Riemann Zeta関数は、

\[ \zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots \]

で定義される無限級数の複素関数です。

詳しくは、リーマンのゼータ関数で遊び倒そう (Ruby編)を参照のこと。

なお、メッシュの次元(300, 500, 1)、原点(0.0001, 0.0001, 0)、間隔(0.1, 0.1, 1)はコード内の以下で設定しています。

 out.SetDimensions(300, 500, 1)
 out.SetOrigin(0.0001, 0.0001, 0)
 out.SetSpacing(0.1, 0.1, 1)

 

script

上記の設定が完了しましたら、"Apply"します。(Riemann Zeta関数の計算は時間がかかりますので気長に待ちましょう)。

描画の設定

Riemann Zeta関数の計算が終了し、画像データが作成されましたら、表示を変数"zeta"にしコンター図を表示しましょう。

ここで、Riemann Zeta関数は、s=1が極となり、その周辺で絶対値が非常に大きいため、コンター表示する場合は、値の範囲を"Rescale"する必要があります。

また、形状に座標軸(実部(real)と虚部(imaginary))を表示するとよいでしょう。

以下の図のように、値の範囲を[0, 4]にリスケールし、続けて、Axes Gridsを有効にしてから、"Edit"ボタンから軸の名前を(x軸のタイトル->real, y軸のタイトル->imaginary)に変更します。Axes GridsのEditでは、フォントの種類や大きさも変更できますので適宜調節すると良いでしょう。

rescale and enable axis_grid

以上で、Programable Sourceを使用して、画像データを作成し、その各点にRiemann Zeta関数の絶対値を与え、最期に、そのコンターを見ることができました。

riemann-zeta contours

Riemann予想

Riemann予想によるとRiemann Zeta関数の零点

\[ \zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots = 0を満たすsの値 \]

の実部は全て1/2となります。

Slice機能で、実部が1/2の位置で切断して、表示してみると確かにその線上に、零点(青色)が並んでいます。

real part=0.5

また、Warp by Scalarフィルターなどを使うと、零点の位置が見易いかもしれません。

warp by scalar

添付ファイル: fileRequestInformation.png 24件 [詳細] filewarp_by_scalar.png 27件 [詳細] fileRiemann_zero_point.png 30件 [詳細] fileriemann_zeta_contours.png 26件 [詳細] filescript.png 28件 [詳細] filerescale_and_axis_grid.png 24件 [詳細] filemenue_programmable_source.png 29件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2018-05-03 (木) 08:15:48 (77d)