Embedded Files
３次元可視化ツールとして今の所 無料のParaview を使っている。
僕は Fortranで計算したデータは そのままバイナリ出力して保存しているので、Paraviewで読み込むためにはこれらを VTK（Visualization Toolkit）ファイルに変換する必要がある。現在、この作業はFortranでやっているが、近い将来 Paraviewを動かすのも含めて全てPython に移行する予定。
以下のサイトが参考になる。
このページは、FortranでVTKファイルの出力する方法と Paraviewの基本操作を自分用メモとしてまとめたもの。随時更新中。
３次元デカルト座標データ
３次元デカルト座標データ
デカルト座標の場合は単純に上から
ヘッダー（ファイルの説明 / フォーマットはASCII ）
データタイプは STRUCTURED_GRID
x, y, z の次元数（Nx, Ny, Nz）
データの総個数（Nx*Ny*Nz）
座標データ（上のデータ数行）
出力する変数（POINT_DATA）の個数
出力する変数のタイプと名前
上の座標データに対応する変数データ
と並べてやれば良い。
直方体の表面２次元データを出力
直方体の表面２次元データを出力
例として、３次元対流計算結果の３次元ムービーを作りたい場合を考える。
可視化のためだけなら正直３次元データのまま vtkファイルに出力する必要はない。
見える面だけを２次元データとして保存するのが（容量節約になるし Paraviewもサクサク動くので）良い。
以下の例で使用している対流計算データはここからダウンロード可能。
まず鉛直速度 vz(N_x,N_y,N_z) の z=z(N_z) における断面（top view）を出力する場合は以下のようになる。
同様に、x=0 における断面（side view）を vtkファイルに出力するには次のようにする。
これらの vtkファイルは Paraviewで読み込むことができて、適当に色や視点角度を調整すれば下のような図を簡単に作ることができる。
さらに、図の一部領域を切り抜いたように表示させたい場合を考える（下図）。これは、消したい領域に適当な値（例えば不自然なほど大きな値）を与えておいて、これを Paraview の Threshold 機能を使って不可視化することによって実現できる。
ボリュームレンダリング用の３次元データを出力
ボリュームレンダリング用の３次元データを出力
正直上のような図なら matplotlib で頑張って作れないこともない。
Paraview の恩恵を最大限感じられるのは、やはりボリュームレンダリングや磁力線（流線）を可視化したい時だろう。
ここでは、例として上述の３次元対流計算結果の鉛直速度をボリュームレンダリング表示してみる。
Paraview で読み込んだら表示タイプを Volume にして、Volume Rendering / Select Mapper を Projected tetra から Resampled To Image に変更する。
Contour, Clip, Slice など使える機能が豊富。カラーマップの transparacy を変えれば、以下のような図も作れる。
３次元球座標データ
３次元球座標データ
球座標データは UNSTRUCTURED_GRID なので、データ形状を陽的に記述する必要がある。
この場合のVTKファイルは上から順に
ヘッダー（ファイルの説明 / フォーマットはASCII / データタイプは STRUCTURED_GRID）
座標点データ（POINTS）
上の座標点からなるセルの構成要素（CELLS）/ 右回りか左回りでちゃんと一周するように記述
セルの形状（CELL_TYPES）/ セルタイプ番号については公式ドキュメント参照
上の座標データに対応する変数データ（POINT_DATA）
と並べる。
表面２次元データの出力
表面２次元データの出力
球面上定義された２次元データ Q(θ, φ ) で全球を完全に覆うために、
極 (θ=0, π ) 周りは、３角形（セルタイプ番号 : ５）
それ以外は、４辺形（セルタイプ番号 : ９）
を利用する。
ここでは、例として全球回転対流の数値計算データを３次元表示してみる（用いた計算データはここからダウンロード可能）。
まず、準備として スタッガード格子上定義されていた変数を Q[Nr,Nq,Np] ⇒Qex[Nr,0:Nq,0:Np] と拡張する。
経度方向には周期条件を課し、緯度方向には極 (θ=0, π ) を端点に置いた。
VTKファイルは、座標データを４辺形（cells : 4 / cell_type : 9）と３角形（cells : 4 / cell_type : 9）に分けて以下のように記述する。
赤道面２次元データの出力
赤道面２次元データの出力
赤道面は、４辺形のみで表現できるので簡単。VTKファイルの出力ルーチンは以下の通り。
子午面２次元データの出力
子午面２次元データの出力
子午面も４辺形のみで表現できるが、極を跨いで経度座標インデックスが変化することに注意。VTKファイルの出力ルーチンは以下の通り。
これらのVTKファイルを Paraview で読み込むことで（そして少し工夫することで）、次のような図を作ることが可能。
磁力線の３次元可視化
磁力線の３次元可視化
Paraview の最大の魅力は、磁力線を綺麗に可視化できることである（特に我々のように MHD計算をしている人にとっては）。
３次元データの出力
３次元データの出力
ここでは、例として全球対流ダイナモ計算で得られた磁場データを VTK出力し、磁力線として描画してみる。
すでに、球座標からデカルト座標への変換は済んでるものとする。
この場合のVTKファイルは上から順に
ヘッダー（ファイルの説明 / フォーマットはASCII / データタイプは STRUCTURED_GRID）
座標点データ（POINTS）
上の座標データに対応する変数データPOINT_DATA）の個数
ベクトルデータ（VECTORS）/ 同じ行に x, y, z 成分並べる
スカラー（SCALARS）/ 磁力線に色を付ける指標として使う
と並べる。
下の例では、磁場の経度（トロイダル）成分と磁場の強さ（絶対値）をそれぞれスカラー場 B_tor と Bamp として保存している。
Paraview で磁力線の描画方法
Paraview で磁力線の描画方法
手順は以下の通り。
VTKファイルを読み込んだら、変数をベクトルデータ（mag_vec）に指定した上で、Stream Tracer を選択する
Seeds / Seed Type を High Resolution Line Source から Point Source に変更
Spere Parameters の Show Sphere のチェックを外し、Center の座標を原点（0, 0, 0）に合わせ、Radius を１とする
Number of Points は デフォルトは 100 だが、必要に応じて上げる（磁力線の本数が増えるが、その分動作は重くなる）
Coloring で用いるスカラー場として B_tor を選択
Stream Tracer を Apply したら、Filters / Alphabetical / Tube を選択
Scalars を B_tor から Bamp に変更。ただし、Coloring には B_tor を指定
Tube Radius を調整、Vary Radius を Off から By Scalar に変更。Radius Factor も必要に応じて上げる
Styling の Opacity を調整し、Apply する。カラーマップや視点などは適当に調整する
こうして出来たのが下図である。赤 / 青がそれぞれトロイダル磁場の正負、磁力線の太さが磁場の強さを表している。
Google Sites
Report abuse