ここでは，多数の等分割の直方体の格子から構成されるデータの時間発展を記述する例を示します．詳細はXCodeでVTK楽園を参照
概要
//ヘッダー{{{
#include <boost/filesystem.hpp> // フォルダーの操作を行うので, boostを使います．
#include <iostream> // この辺は，おなじみですよね♪
#include <vector> // もちろんSTLを使っても大丈夫です. ただしMicrosoft-CLIを使って良いか？不明です．
#include <cmath>
#include <sstream>
#include <vtkVersion.h> // VTKのヘッダーです
#include <vtkSmartPointer.h>
#include <vtkProperty.h>
#include <vtkImageData.h> // 等分割直方格子の領域を使います
#include <vtkMultiBlockDataSet.h> // １ケースに多数の領域が含まれる場合，これが必要
#include <vtkXMLMultiBlockDataWriter.h>// 多数の領域を一発で保存する
#include <vtkDoubleArray.h> // セルや格子点に取り付ける物理量を定義します
#include <vtkPointData.h> // 格子点のデータ
#include <vtkCellData.h> // セルのデータ
#include <vtkCell.h> // セル
//}}}
int main() // 複雑形状を構成 {{{
{
+--- 7 行:時間ステップ----------------------------------------------------------------------------------------------------
+--- 12 行:フォルダーを作成------------------------------------------------------------------------------------------------
// 時間ステップだけ作業する
for(unsigned int nt=1;nt<=ntMax;nt++)
{
// この時間のデータを保存するMultiBlockDataSetを定義
vtkSmartPointer<vtkMultiBlockDataSet> multi=vtkSmartPointer<vtkMultiBlockDataSet>::New();
// 領域達を定義 {{{
for(int idom=0;idom<領域総数;idom++)
{
+---- 22 行:ImageData形状を定義--------------------------------------------------------------------------------------------
+---- 3 行:MultiBlockDataSetに追加----------------------------------------------------------------------------------------
+---- 35 行:ImageDataに物理量を定義 ---------------------------------------------------------------------------------------
std::cout << "作製:" << idom << std::endl;
}//}}}
+--- 16 行:ファイルに保存 -------------------------------------------------------------------------------------------------
}
return EXIT_SUCCESS;
}//}}}
時間ステップの設定およびフォルダー作成
// 時間ステップ{{{
unsigned int ntMax=100;
boost::filesystem::path folderName("myDSMC");
std::string stemName="time"; // ファイル名は folderName/stemName時間ステップ(N桁の数字)extName
std::string extName=".vtm"; // 複数領域では，拡張子は .vtm です
int fileN=4; // まあテストなので, 4桁あれば十分
//}}}
// フォルダーを作成(既存のフォルダーは黙って削除){{{
try
{
if (boost::filesystem::exists(folderName)) boost::filesystem::remove_all(folderName);
boost::filesystem::create_directory(folderName);
}
catch (std::exception& e)
{
std::cerr << folderName << "フォルダがエラー:" << e.what() << std::endl;
return 1;
}//}}}
ImageData形状を定義
// ImageData形状を定義{{{
typedef vtkSmartPointer<vtkImageData> ImageDataPtr; //行が長くなるのでtypedef
ImageDataPtr image=ImageDataPtr::New(); //上のtypedefを用いてImageDataを作成
image->SetSpacing(0.05,0.05,0.05); // もちろんX-Y-Z方向の格子の間隔
int imageExtent[6]={0,0,0,0,0,0}; // あとで各領域の原点や格子点数を設定します
double imageOrigin[3]={0.,0.,0.};
switch (idom)
{
case 0: // 領域0のパラメータを設定
imageExtent[1]=20; imageExtent[3]=10; imageExtent[5]=10;
imageOrigin[0]=0; imageOrigin[1]=0; imageOrigin[2]=0;
break;
;;
case 1: // 領域1のパラメータを設定
imageExtent[1]=10; imageExtent[3]=40; imageExtent[5]=40;
imageOrigin[0]=1; imageOrigin[1]=-0.5; imageOrigin[2]=-0.5;
break;
;;
// その他の領域も必要なら書く
}
image->SetExtent(imageExtent);
image->SetOrigin(imageOrigin);
//}}}
MultiBlockDataSetに追加
// MultiBlockDataSetに追加{{{
multi->SetBlock(idom,image);
//}}}
ImageDataに物理量を定義
// ImageDataに物理量を定義 {{{ // ここでは，セルに物理量を準備してみましょう．格子点の場合 Cell をPointに書き換えます
vtkDoubleArray* pressure = vtkDoubleArray::New(); // 圧力の配列を作成
vtkDoubleArray* velocity = vtkDoubleArray::New(); // 流速の配列を作成
pressure->SetNumberOfComponents(0); // 圧力はスケーラなので成分は無いっす
velocity->SetNumberOfComponents(3); // 流速はベクタで3成分すね
pressure->SetName("p"); // 名前を付けておきませう
velocity->SetName("v");
velocity->SetComponentName(0,"v_1");
velocity->SetComponentName(1,"v_2");
velocity->SetComponentName(2,"v_3");
pressure->SetNumberOfTuples(image->GetNumberOfCells()); // 圧力の要素数は, セルの総数
velocity->SetNumberOfTuples(image->GetNumberOfCells()); // 流速の要素数は, セルの総数
image->GetCellData()->AddArray(pressure); //領域のセルデータに，さっき作った圧力を追加
image->GetCellData()->AddArray(velocity); // さっき作った流速を追加
// まあ，テストのために値を定義しておきましょう
double* weight = new double [image->GetMaxCellSize()];
for(vtkIdType cellId=0; cellId < image->GetNumberOfCells(); cellId++)
{ // セルの中心座標を取得します. 格子点の場合には image->GetPoint(pointId,Xcoord)で取得できますが, セルは面倒
vtkCell* cell=image->GetCell(cellId);
double Pcoord[3]; // セルは体積要素ですので，内部にいちいち曲線座標が存在します
int subId= cell->GetParametricCenter(Pcoord); //曲線座標表示で, セルの中心を求めます．
// セルは複合セルであったりするので, subIdも取得します(普通はsubId=0)
double Xcoord[3];
cell->EvaluateLocation(subId,Pcoord,Xcoord,weight); // ようやく, セル中心のXcoordが取得できました♪
double T=pow(Xcoord[0],2)+pow(Xcoord[1],2)+pow(Xcoord[2],2); // てきとーです．サンプルですので
double p=exp(-0.1*T);
pressure->SetValue(cellId,p); // まあ，要するにこうやって値を設定するわけ
double v1=0.1;
double v2=p*Xcoord[1];
double v3=-0.1;
velocity->SetComponent(cellId,0,v1);
velocity->SetComponent(cellId,1,v2);
velocity->SetComponent(cellId,2,v3);
}
//}}}
std::cout << "作製:" << idom << std::endl;
}//}}}
ファイルに保存
// ファイルに保存 {{{
std::stringstream sst;
sst<< stemName << std::setfill('0') << std::setw(fileN) << nt << extName;
boost::filesystem::path fileName(sst.str());
fileName = folderName / fileName;
// 保存を実行
vtkSmartPointer<vtkXMLMultiBlockDataWriter> writer = vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
writer->SetFileName(fileName.c_str()); // ここでファイル名を設定
#if VTK_MAJOR_VERSION <= 5 // 2012年現在，たまたま移行期なので・・・今はまだVersion5です
writer->SetInputConnection(multi->GetProducerPort());
#else
writer->SetInputData(multi);
#endif
writer->Write(); // これで大量のファイルが生成されます
//}}}
コンパイル
上記プログラムは, 次のようにコンパイルできます]
ic++ --boost --vtk myProgram.cxx -o myExecutable