vtk形式によるブロック構造格子AMR(Adaptive Mesh Refinement)データの記述方法とParaViewによる可視化方法

本記事の目的

vtk形式のうち、vtkHierarchicalBoxDataSetを用いたAMR構造格子データの記述とParaViewを用いた可視化について解説します。XMLベースのvtk形式については別記事にまとめましたので、vtk形式の基礎的な部分については省略します。

vtkHierarchicalBoxDataSet(.vthb)について

vtkHierarchicalBoxDataSetは、Serial Structuredなvtk形式の種類の1つであり、階層構造を持った3次元直交格子を記述/可視化するために使われます。難点は公式のドキュメント上には存在していない点、かつ仕様がFixされていないためかversionによって大きく記述形式が変わっている上、明確にそれを説明している資料が見当たらない点です。（もし存在していたら教えてください）

本稿執筆時点ではversion="1.1"が最新のようなので、このバージョンを対象にまとめます。

公式のサンプルデータ

VTKData/AMR/HierarchicalBoxDataSet.v1.1.vthb

https://www.vtk.org/Wiki/VTK_Datasets で公開されています。
git clone git://vtk.org/VTKData.gitすることでclone可能です。

基本形式

基本的な形式はいつも通り(VTKFile要素)ですが、versionを必ず指定する必要があります。可視化ソフトのバージョンによって対応しているversionが違うっぽいのが最悪です。（ParaView 5.2.0でv1.0のvthbを読み込んだら落ちました)

basic.vtk

<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1" >
  <vtkHierarchicalBoxDataSet ...>
    ...
  </vtkHierarchicalBoxDataSet>
</VTKFile>

vtkHierarchicalBoxDataSet要素の中身

vtkHierarchicalBoxDataSet要素はアトリビュートとして以下の2つを持てるようです。(それ以外は不明)

origin

起点の座標

grid_description

軸の記述順番。通常はXYZで良さそうです

また、子としてBlock要素を持ちます。Block要素が階層構造の各レベルに対応します。
Block要素はlevelアトリビュートとspacingアトリビュートを持ちます。

ここまでをまとめると、vtkHierarchicalBoxDataSet要素はこんな形の構造になります。

basic2.vthb

<vtkHierarchicalBoxDataSet origin="-1.75 -1.25 0" grid_description="XYZ">
  <Block level="0" spacing="0.833333 0.833333 0.666667">
    ...  
  </Block>
  <Block level="1" spacing="0.416667 0.416667 0.333333">
    ...
  </Block>
</vtkHierarchicalBoxDataSet>

version1.1ではspacingをFloatで記述できるので整数比(1:2, 1:4等)でない場合にも対応したようですが、
数値計算で使う場合には整数比以外必要なさそうなので、ちょっと面倒に感じます。

注意点として、他のVTKのデータ形式ではアトリビュートが"Origin"、"Spacing"なのですが、
vtkHierarchicalBoxDataSetについては"origin", "spacing"です。間違えると正しく読み込まれません。最悪です。

Block要素の中身

Block要素の中には各グリッドに対応するDataSet要素を記述します。
DataSet要素はindexアトリビュートとamr_boxアトリビュートとfileアトリビュートを持ちます。
それぞれ、下記のような役割があります。

index

同Block内でのindex

amr_box

「親グリッドのどの位置を抜き取る」かどうかを表す。データの`WholeExtent`とは別になります。詳しくは後述します。

file

各Gridの情報を記述しているvtkImageData(.vti形式)ファイルへのPath

fileでImageDataを指定するため、基本的に.vthb形式の中には実データを持ちません。
そのため、Pが頭についていませんがParallelな形式の1つのようなものと考えていいのかもしれません。

DataSet要素を記述した各Blockは下記のような形になります。

root.vthb

<Block level="0" spacing="0.833333 0.833333 0.666667">
  <DataSet index="0" amr_box="0 3 0 2 0 2" file="HierarchicalBoxDataset.v1.1/HierarchicalBoxDataset.v1.1_0.vti"/>    
</Block>
<Block level="1" spacing="0.416667 0.416667 0.333333">
　　<DataSet index="0" amr_box="0 2 0 2 0 2" file="HierarchicalBoxDataset.v1.1/HierarchicalBoxDataset.v1.1_1.vti"/>
　　<DataSet index="1" amr_box="3 7 0 2 0 2" file="HierarchicalBoxDataset.v1.1/HierarchicalBoxDataset.v1.1_2.vti"/>
　　<DataSet index="2" amr_box="0 2 3 5 0 2" file="HierarchicalBoxDataset.v1.1/HierarchicalBoxDataset.v1.1_3.vti"/>
</Block>

ParaViewでの可視化

ParaView5.2.0では、vthb形式を読み込むと自動的にAMR Datasetとして読み込まれます。
（InformationにDeprecatedの文字があるのが気になります……）

Sliceを追加すれば普通に切り出してくれるのですが、

このように一番粗いGridのみが描画されてしまう場合は、Sourceの「Default Number Of Levels」を1以上の値に変更するとサブグリッドまで描画されます。

他、AMR用のFilterの使い方については、使い方を習得し次第追記します。

amr_boxアトリビュートの意味

一番小さいグリッドの単位はシンプルなSerial ImageDataとして表せることが分かったので、AMRデータ化する場合に考慮する必要があるのはamr_boxアトリビュートだけでいいことになります。
amr_boxアトリビュートは親グリッド上のどの位置を非表示にするかを指定するために使われるアトリビュートです。
具体例でどのように描画されるかを見ていきます。

サンプルデータは要素数が多すぎるため、ここからは自前の小さいデータを使います。
親グリッドとして以下のデータを用意します。

parent.vti

<?xml version="1.0"?>
<VTKFile type="ImageData">
    <ImageData WholeExtent="0 2 0 2 0 2" Origin="0.0 0.0 0.0" Spacing="1.0 1.0 1.0">
        <Piece Extent="0 2 0 2 0 2">
            <PointData Scalars="test_node_values">
                <DataArray Name="test_node_values" type="Int32" format="ascii">
                    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
                </DataArray>
            </PointData>
        </Piece>
    </ImageData>
</VTKFile>

子グリッドとして以下のデータを用意します。

child.vti

<?xml version="1.0"?>
<VTKFile type="ImageData">
    <ImageData WholeExtent="0 2 0 2 0 2" Origin="0 0 0" Spacing="0.5 0.5 0.5">
        <Piece Extent="0 2 0 2 0 2">
            <PointData Scalars="test_node_values">
                <DataArray Name="test_node_values" type="Int32" format="ascii">
                    -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
                </DataArray>
            </PointData>
        </Piece>
    </ImageData>
</VTKFile>

(Spacingを半分にしてPointの値が反転している以外は同じ構造です）
vthbファイルを次のように書いてみます。

test1.vthb

<?xml version="1.0"?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1">
    <vtkHierarchicalBoxDataSet origin="0.0 0.0 0.0" grid_description="XYZ">
        <Block level="0" spacing="1.0 1.0 1.0">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="parent.vti">
            </DataSet>
        </Block>
        <Block level="1" spacing="0.5 0.5 0.5">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="child.vti">
            </DataSet>
        </Block>
    </vtkHierarchicalBoxDataSet>
</VTKFile>

このファイルは下図のように描画されます。

これだけではamr_boxアトリビュートの動きがよく分からないので、vthbファイルを少し弄ってみます。
まずamr_box="0 1 0 1 0 1"からamr_box="0 1 0 2 0 1"に変化させた場合。

test1.vthb

<?xml version="1.0"?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1">
    <vtkHierarchicalBoxDataSet origin="0.0 0.0 0.0" grid_description="XYZ">
        <Block level="0" spacing="1.0 1.0 1.0">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="parent.vti">
            </DataSet>
        </Block>
        <Block level="1" spacing="0.5 0.5 0.5">
            <DataSet index="0" amr_box="0 1 0 2 0 1" file="child.vti">
            </DataSet>
        </Block>
    </vtkHierarchicalBoxDataSet>
</VTKFile>

子の格子がy方向に1つ拡大され、対応する親の領域のみ描画されなくなってしまいました。
amr_box="0 1 0 3 0 1"に変えてみます。

ほぼ同じですが、y方向の格子が更に1つ増えて親格子のサイズと一致しました。

ここで親格子のWholeExtentは"0 2 0 2 0 2"なので、1方向のノードとしては[0, 1, 2]の3点になります。
これをRefinementRatio=2で分割した場合、1方向のノードは[0, 0.5, 1, 1.5, 2.0]の5点になります。(親がn点の時は(n-1)*2 + 1点です)
細分化されたノードに0からインデックスをつけると[0, 1, 2, 3, 4]となります。
ここで、ノードではなくてセルのインデックスは2ノードで1ペアなので1つ減って、[0, 1, 2, 3]になります。
amr_boxアトリビュートは、この「現在の格子のセルインデックス」を指定して、その部分を親グリッドで描画しないようにする、という動作を行うアトリビュートです。
具体的にはamr_box="imin imax jmin jmax kmin kmax"と指定した時、

• x方向: iminからimaxのセルまで (= iminからimax+1のノードまで)
• y方向: jminからjmaxのセルまで (= jminからjmax+1のノードまで)
• z方向: kminからkmaxのセルまで (= kminからkmax+1のノードまで)

を非表示にします。
試しにamr_box="0 0 0 0 0 0"としてみると、下図のようになりました。

白い枠線が子グリッド格子に割り当てられた広がりを示しており、各方向の0番目のセルだけが（ノードだと0から1）、実座標で言えば0と0.5の部分が囲まれていることが分かります。

これまで貼り付けた図からも分かる通り、amr_boxアトリビュートによって非表示にされる領域と実際のvtiデータが描画される位置は無関係であることが分かります。amr_box="0 0 0 0 0 0"とした例では、vtiデータのWholeExtentが各方向3グリッド分あるのに対し、amr_boxアトリビュートで割り当てられているのは2グリッド分であるため親グリッドの描画と重なり、正しくない描画が発生しています。

基本的には、子グリッドとして割り当てるvtiデータのWholeExtent="0 xmax 0 ymax 0 zmax"だったなら、amr_box="0 xmax-1 0 ymax-1 0 zmax-1"とすれば要素数としてはぴったり合います。位置を動かしたい場合は以下の2つを行う必要があります。

• vtiデータのoriginを移動（空間中の正しい位置に描画されるようにする）
• amr_boxアトリビュートにオフセットを追加

子グリッドをx方向に移動してみましょう。まずchild.vtiのoriginをorigin="1.0 0.0 0.0"に変更します。

child.vti

<?xml version="1.0"?>
<VTKFile type="ImageData">
    <ImageData WholeExtent="0 2 0 2 0 2" Origin="1.0 0 0" Spacing="0.5 0.5 0.5">
        <Piece Extent="0 2 0 2 0 2">
            <PointData Scalars="test_node_values">
                <DataArray Name="test_node_values" type="Int32" format="ascii">
                    -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
                </DataArray>
            </PointData>
        </Piece>
    </ImageData>
</VTKFile>

次に、vthb側のamr_boxをamr_box="2 3 0 1 0 1"に変更します。

test1.vthb

<?xml version="1.0"?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1">
    <vtkHierarchicalBoxDataSet origin="0.0 0.0 0.0" grid_description="XYZ">
        <Block level="0" spacing="1.0 1.0 1.0">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="parent.vti">
            </DataSet>
        </Block>
        <Block level="1" spacing="0.5 0.5 0.5">
            <DataSet index="0" amr_box="2 3 0 1 0 1" file="child.vti">
            </DataSet>
        </Block>
    </vtkHierarchicalBoxDataSet>
</VTKFile>

このように描画されます。

上手く移動できました。

最後に、2レベル以上のレベルがある時の指定方法について学びましょう。
孫グリッドファイルを用意します。

grandchild.vti

<?xml version="1.0"?>
<VTKFile type="ImageData">
    <ImageData WholeExtent="0 2 0 2 0 2" Origin="1.0 0 0" Spacing="0.25 0.25 0.25">
        <Piece Extent="0 2 0 2 0 2">
            <PointData Scalars="test_node_values">
                <DataArray Name="test_node_values" type="Int32" format="ascii">
                    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
                </DataArray>
            </PointData>
        </Piece>
    </ImageData>
</VTKFile>

vthbを以下のように変更します。

test1.vthb

<?xml version="1.0"?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1">
    <vtkHierarchicalBoxDataSet origin="0.0 0.0 0.0" grid_description="XYZ">
        <Block level="0" spacing="1.0 1.0 1.0">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="parent.vti">
            </DataSet>
        </Block>
        <Block level="1" spacing="0.5 0.5 0.5">
            <DataSet index="0" amr_box="2 3 0 1 0 1" file="child.vti">
            </DataSet>
        </Block>
        <Block level="2" spacing="0.25 0.25 0.25">
            <DataSet index="0" amr_box="4 5 0 1 0 1" file="grandchild.vti">
            </DataSet>
        </Block>
    </vtkHierarchicalBoxDataSet>
</VTKFile>

こうなります。

孫のamr_boxアトリビュートに注目しましょう。現在のvtkHierarchicalBoxDataSetでは階層毎にindexがあるのみで、親子関係を持つことができません。そのため、amr_boxアトリビュートはそのグリッドの親に対してどこにいるのかではなく、空間全体(Level0のグリッド)に対してどこにいるのかを指定することになります。また孫のvtiファイルのOriginも同様で、空間全体のどこから始まるのかを記述する必要があります。

親子関係を持てるなら、孫グリッドのアトリビュートはamr_box="0 1 0 1 0 1"のように単純に親との関係のみで記述できるのですが、現在は無理のようです。

まとめると

• amr_boxアトリビュートは親のグリッドに対するインデックスではなく、Level0の空間に対するインデックス
• 個々のグリッドのextentと起点が全体空間中でどの位置になるのかは、個々のグリッドが知っておく必要がある

となります。

ちなみに、Block要素は各レベル1つという訳ではなく、順不同で記述することができます。

拙作の宣伝: SimpleVTKによるAMRデータの可視化

基本的にAMRのような階層構造グリッドを作る際、全空間中でどの位置になるのかというのは保持しなくても問題がありません。しかし、vtkHierarchicalBoxDataSetとして出力する場合には何らかの方法で全空間中での座標を解決しなければなりません。Block中のDataSet要素は既にindexを持っているので、親グリッドのIndexさえ指定できれば全空間中におけるamr_boxを解決できます。

ですので、拙作のライブラリ SimpleVTK において「親のグリッド上でどこにいるのかというのを指定したら自動的にamr_boxアトリビュートを解決してくれる」機能を実装しました。
hsimyu/SimpleVTK - https://github.com/hsimyu/SimpleVTK

また、セルベースではなくてノードベースでも指定が可能になっています。
2レベルのvthbを出力するサンプルは以下のようになります。

#include <simple_vtk.hpp>
int main() {
    SimpleVTK gen;
    gen.enableExtentManagement();
    gen.changeBaseExtent(0, 2, 0, 2, 0, 2);
    gen.changeBaseOrigin(0.0, 0.0, 0.0);
    gen.changeBaseSpacing(1.0, 1.0, 0.5);
    gen.changeRefinementRatio(2.0);

    gen.beginVTK("vtkHierarchicalBoxDataSet");
    gen.beginContent();
    gen.setGridDescription("XYZ");
        gen.beginBlock();
            gen.beginDataSet(0);
            gen.setAMRBoxNode(0, 2, 0, 2, 0, 2); // can use node index
            gen.setFile("parent.vti");
            gen.endDataSet();
        gen.endBlock();
        gen.beginBlock();
            gen.beginDataSet(0);
            gen.setAMRBoxNodeFromParentIndex(0, 0, 1, 1, 2, 0, 1); // 親グリッド(0-0)の 0から1, 1から2, 0から1に乗っている)
            gen.setFile("child.vti");
            gen.endDataSet();
        gen.endBlock();
        gen.beginBlock();
            gen.beginDataSet(0);
            gen.setAMRBoxNodeFromParentIndex(0, 0, 1, 0, 1, 0, 1); // 親グリッド(1-0)の 0から1, 0から1, 0から1に乗っている)
            gen.setFile("grandchild.vti");
            gen.endDataSet();
        gen.endBlock();
    gen.endContent();
    gen.endVTK();

    gen.generate("test_vthb_amr_support");
}

この出力は以下のようになります。
(parent.vti, child.vti, grandchild.vtiはそれぞれWholeExtent="0 2 0 2 0 2"のImageDataとしています)

<?xml version="1.0" ?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="1.1">
    <vtkHierarchicalBoxDataSet origin="0 0 0" grid_description="XYZ">
        <Block level="0" spacing="1 1 0.5">
            <DataSet index="0" amr_box="0 1 0 1 0 1" file="parent.vti">
            </DataSet>
        </Block>
        <Block level="1" spacing="0.5 0.5 0.25">
            <DataSet index="0" amr_box="0 1 2 3 0 1" file="child.vti">
            </DataSet>
        </Block>
        <Block level="2" spacing="0.25 0.25 0.125">
            <DataSet index="0" amr_box="0 1 4 5 0 1" file="grandchild.vti">
            </DataSet>
        </Block>
    </vtkHierarchicalBoxDataSet>
</VTKFile>

amr_boxアトリビュートのうちy座標が自動的に解決されていることが分かります。また、spacingについても自動計算をするようになっています。
可視化するとこうなります。

便利ですね！！
お疲れ様でした。

(おまけ) 古いvtkHierarchicalBoxDataSet形式(.vtm形式)

2007年とかそのくらいのリンクを遡ると、vtkHierarchicalBoxDataSetは(.vthbではなく).vtmとして扱われていたようです。

• [Paraview] [vtm] file formats: https://public.kitware.com/pipermail/paraview/2007-June/005153.html
• VTK - Users - Vtk and AMR: http://vtk.1045678.n5.nabble.com/Vtk-and-AMR-td1231671.html

1つめのリンク先で言及されているchombo3dデータについて、2011年にVTKDataをCloneして放置しているUserのRepositoryを発見したので漁ってみると、こんな形式だったようです。

oldoldvtm.vtm

<?xml version="1.0"?>
<VTKFile type="vtkHierarchicalBoxDataSet" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">
  <vtkHierarchicalBoxDataSet>
    <RefinementRatio level="0" refinement="2"/>
    <RefinementRatio level="1" refinement="2"/>
    <DataSet group="0" dataset="0" amr_box="0 7 0 7 0 7" file="chombo3d_0.vti"/>
    <DataSet group="1" dataset="0" amr_box="0 15 0 15 0 15" file="chombo3d_1.vti"/>
    <DataSet group="2" dataset="0" amr_box="0 11 8 23 8 23" file="chombo3d_2.vti"/>
    <DataSet group="2" dataset="1" amr_box="8 11 0 7 0 7" file="chombo3d_3.vti"/>
    <DataSet group="2" dataset="2" amr_box="12 19 0 7 0 7" file="chombo3d_4.vti"/>
    <DataSet group="2" dataset="3" amr_box="12 15 8 15 8 15" file="chombo3d_5.vti"/>
    <DataSet group="2" dataset="4" amr_box="12 15 8 15 16 23" file="chombo3d_6.vti"/>
    <DataSet group="2" dataset="5" amr_box="12 15 16 23 8 15" file="chombo3d_7.vti"/>
    <DataSet group="2" dataset="6" amr_box="12 15 16 23 16 23" file="chombo3d_8.vti"/>
    <DataSet group="2" dataset="7" amr_box="16 19 12 15 12 15" file="chombo3d_9.vti"/>
    <DataSet group="2" dataset="8" amr_box="16 19 12 15 16 19" file="chombo3d_10.vti"/>
    <DataSet group="2" dataset="9" amr_box="16 19 16 19 12 15" file="chombo3d_11.vti"/>
    <DataSet group="2" dataset="10" amr_box="16 19 16 19 16 19" file="chombo3d_12.vti"/>
    <DataSet group="2" dataset="11" amr_box="20 23 8 11 8 11" file="chombo3d_13.vti"/>
    <DataSet group="2" dataset="12" amr_box="20 23 12 19 12 19" file="chombo3d_14.vti"/>
    <DataSet group="2" dataset="13" amr_box="24 31 8 11 8 11" file="chombo3d_15.vti"/>
  </vtkHierarchicalBoxDataSet>
</VTKFile>

Block要素でのまとまりがないため、現在のVersionより使いにくそうです。
（その後のVersionの発展については見つけられませんでした)