可視化

パラメータにずらずら書くとバグの素なのでクラスに入れといた. 

//  vtk_myIO.hpp    僕のVTKの入出力
#pragma once
#include <fstream>
#include <string>
#include <Eigen/Dense>
class vtk_myIO
{
public:
   std::filesystem::path filename;
   Eigen::Vector3d ox;
   Eigen::Vector3d dx;
   Eigen::Vector3i nx;
   void write_legacy_vti2d(std::function<double(int,int)> data)
   {
       std::filesystem::create_directories(filename.parent_path());
       std::ofstream ofs(filename);
       ofs << "# vtk DataFile Version 3.0\n";
       ofs << "Legacy VTK file generated by write_legacy_vtk2d\n";
       ofs << "ASCII\n";
       ofs << "DATASET STRUCTURED_POINTS\n";
       ofs << "DIMENSIONS " << nx[0] << " " << nx[1] << " 1\n";
       ofs << "ORIGIN " << ox[0] << " " << ox[1] << " 0.0\n";
       ofs << "SPACING " << dx[0] << " " << dx[1] << " 1.0\n";
       ofs << "POINT_DATA " << nx[0] * nx[1] << "\n";
       ofs << "SCALARS scalars double\n";
       ofs << "LOOKUP_TABLE default\n";
       ofs << std::setprecision(15) << std::scientific;
       for (int j = 0; j < nx[1]; ++j)
           for (int i = 0; i < nx[0]; ++i)
               ofs << data(i,j) << "\n";
       ofs.close();
   }
};

いろいろ試す時に面倒なので，潰しが効くようにしといた．ここは，呼び出す時に，なんかi, jから値が出るのを与えれば良いことにしとく．これは入れといた方が安全だということを前に学んだよね．ここは，標準だと7-8桁しか保存せず結局不足することが多いので，全桁書くように設定した．するってえと, main()は

// main.cpp
#include <iostream>
#include "../include/vtk_myIO.hpp"
int main()
{
   vtk_myIO Writer;
   Writer.ox << -3,-3, 0;      //めんどくさがってvtk_write.oxをEigen::Vector3d にしたので，3成分入れてね
   Writer.dx << 0.1, 0.1, 0.;
   Writer.nx << 60, 60, 1;
   std::filesystem::path myHome = getenv("HOME");
   Writer.filename=(myHome/"Data")/"test05";       // あなたが好きなデータフォルダーを設定
   Writer.filename/="LEGACY_2D.vtk";　　　　　　　　　// ファイル名は  ほげほげ.vtk にしてね！
   Writer.write_legacy_vti2d([&Writer](int i, int j) -> double
   {
       double x=Writer.ox[0]+Writer.dx[0]*i;
       double y=Writer.ox[1]+Writer.dx[1]*j;
       return std::exp(-x*x-y*y);
   });
   return 0;
}

これで LEGACY_2D.vtk ファイルが完成した！

追加してみた. ん？というところもあったので，少し猿と対話したので，上のchat結果とはちょっと違うかも

// vtk_myIO.hpp    僕のVTKの入出力
...
#include <vtkSmartPointer.h>
#include <vtkImageData.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkXMLImageDataWriter.h>
...
    void write_xml_vti2d(std::function<double(int, int)> data)
   {
       // VTKのImageDataオブジェクト作成
       auto image = vtkSmartPointer<vtkImageData>::New();
       image->SetDimensions(nx[0], nx[1], nx[2]);
       image->SetOrigin(ox[0], ox[1], ox[2]);
       image->SetSpacing(dx[0], dx[1], dx[2]);
       // スカラー値の配列を作成
       auto ARR= vtkSmartPointer<vtkDoubleArray>::New();
       ARR->SetName("pressure");
       ARR->SetNumberOfComponents(1);
       ARR->SetNumberOfTuples(nx[0] * nx[1]);
       image->GetPointData()->AddArray(ARR);
       for (int j = 0; j < nx[1]; ++j)
           for (int i = 0; i < nx[0]; ++i)
               ARR->SetValue(j * nx[0] + i, data(i, j));
       // スカラー値の配列をアクティブにする
       image->GetPointData()->SetActiveScalars("pressure");
       // XML VTKファイル（.vti）で出力
       auto writer = vtkSmartPointer<vtkXMLImageDataWriter>::New();
       std::filesystem::create_directories(filename.parent_path());
       writer->SetFileName(filename.c_str());
       writer->SetInputData(image);
       writer->SetCompressorTypeToLZMA(); // LZMA圧縮（zlibより高圧縮）
       writer->Write();
   }
...

ここでは，まず①配列を作成する, ②配列に名を与える, ③配列の要素数を指定する, ④その配列を格子点データに追加する，⑤配列の値を設定する（ここで学んだように1次元化されていることに注意）　少し簡単な方法もある（スカラー1個だけは簡単)が，まあ流体計算では，セルデータに圧力だけじゃなくて流速ベクトルとか設定しないとだしね．ここの方法だと，複数のデータを一つのVTKに追加できるし，追加した関数は全てParaViewで一覧に出現する．これは省略しても大して変わらない．

この例では格子点データに追加しましたが，もちろん，セルデータに追加することもできますが・・・少し面倒なので，ここではやめときましょう．

main()はこんな感じ：

#include <iostream>
#include "../include/vtk_myIO.hpp"
int main()
{
   vtk_myIO Writer;
   Writer.ox << -3, -3, 0;
   Writer.dx << 0.01, 0.01, 0.;
   Writer.nx << 600, 600, 1;
   std::filesystem::path myHome = getenv("HOME");
   auto myFolder = (myHome / "Data") / "test05";
   auto FUNC = [&Writer](int i, int j) -> double
   {
       double x = Writer.ox[0] + Writer.dx[0] * i;
       double y = Writer.ox[1] + Writer.dx[1] * j;
       return std::exp(-x * x - y * y);
   };
   Writer.filename = myFolder / "LEGACY_2D.vtk";
   Writer.write_legacy_vti2d(FUNC);
   Writer.filename = myFolder / "XML_2D.vti";
   Writer.write_xml_vti2d(FUNC);
   return 0;
}

LegacyタイプとXMLタイプの2通り使うので，使いまわせるようにFUNC()というのにしておきました．

こっちは，先ほどと同じLegacyVTK(ファイル名が.vtk)

こっちが，XML-VTK(ファイル名が.vti)    ⭐︎ファイル名を正しく設定しておくと, ParaViewで自動で判断してくれるので便利

PVDwriteは, 最初にPVDファイル名を指定しておきます．で, 時刻t_1のデータがFILE_1にあるのであれば, Append(FILE_1, t_1)します．で，どんどんAppend()してから最後に Write() すると, 指定した名前.pvd というファイルを，指定したフォルダーに作成します．

まあ使い勝手から言って, あるデータフォルダに

• 指定した.PVD
• Contents/ フォルダ
• Contents/FILE_1　　　ノーマルなVTKファイル
• Contents/FILE_2　　   ノーマルなVTKファイル
• ...

というのが，使いやすいと思う．指定した.PVDファイルに，時刻の一覧と多応するファイルが書いてあるので，PVDファイルをParaviewに読み込ませれば動作する．

// main.cpp
#include <iostream>
#include "../include/vtk_myIO.hpp"
#include "../include/pvdwrite.hpp"
int T1()
{
   ...さっきまでの main.cpp のなかみ
}
int T2()
{
   std::filesystem::path myHome = getenv("HOME");　　//このへんは，さっきと同じ
   auto myFolder = (myHome / "Data") / "test06";　　 //このへんは，さっきと同じ
   PVDwrite Indexer(myFolder/"test.pvd");
   vtk_myIO Writer;
   Writer.ox << -3, -3, 0;
   Writer.dx << 0.01, 0.01, 0.;
   Writer.nx << 600, 600, 1;
   double time=0.0, dt=0.01;
   auto FUNC = [&time,&Writer](int i, int j) -> double
   {
       double x = Writer.ox[0] + Writer.dx[0] * i + 0.5*std::cos(2*M_PI*time);
       double y = Writer.ox[1] + Writer.dx[1] * j + 0.5*std::sin(2*M_PI*time);
       double rate= 1.0 + 0.5 * std::sin(8 * M_PI * time);
       return std::exp(-x * x - y * y)*rate;
   };
   for(int nt=0; nt<100; nt++)
   {
       time = nt * dt;
       std::ostringstream oss;
       oss << "XML_2D_" << std::setw(5) << std::setfill('0') << nt << ".vti";
       Writer.filename = (myFolder / "Contents") / oss.str();
       Writer.write_xml_vti2d(FUNC);                 //さっきと同じようにVTKファイルを出力
       Indexer.Append(Writer.filename, time);
   }
   Indexer.Write();
   std::cout << "PVD file written: " << Writer.filename << std::endl;
   return 0;
}

PVDファイルのパス名を指定して初期化しましょう. 名前をIndexerにしています. 描く関数は, まあテストなんで, 時刻timeに依存する関数を定義してみた．で，各時刻について, VTKファイルを作成します．この辺は新顔ですね！これは, 「ふんぎゃあ0023.vti」とか，番号のついたファイルをたくさん作る時に使うものです. 

• std::cout << "ほれほれ" では，画面に「ほれほれ」が送られて出力された．std::ofstream OFS であると,  OFS << "このやろう" で, ファイルOFSに「このやろう」が送られてファイルに出力された．
• std::ostringstream MOJIだと，MOJI << "もじもじ"  によって, 「もじもじ」が，文字列MOJI に送られるのだ！
• oss << "XML_2D_" によって文字列は XML_2D_ になった．std::setw(5)は，文字幅が5文字分になり，std::setfill('0') は，あきスペースは全部0で埋め尽くせ
• そこで oss << nt << ".vti" によって，文字列はXML_2D_00232.vti になった！
• こうやって作成した文字列は, oss.str() として利用できる．

そうやって作成した「ふんぎゃあ0023.vti」ファイル名と, 時刻の値を指定してIndexerにAppendで, 全時刻書き終わったら, IndexerでWriteします．すると，PVDファイルができます．

後学のために，PVDファイルの中身を見てみましょう（ただの目次のテキストファイルですので，普通にエディターで見える）

<?xml version="1.0"?><VTKFile type="Collection" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor" ><Collection>
<DataSet timestep="0" group="" part="0" file="Contents/XML_2D_00000.vti" />
<DataSet timestep="0.01" group="" part="0" file="Contents/XML_2D_00001.vti" />
<DataSet timestep="0.02" group="" part="0" file="Contents/XML_2D_00002.vti" />
<DataSet timestep="0.03" group="" part="0" file="Contents/XML_2D_00003.vti" />
<DataSet timestep="0.04" group="" part="0" file="Contents/XML_2D_00004.vti" />
...
</Collection></VTKFile>

なんかじゅげむじゅげむと書いた後，データ（ファイル名と時刻の値）をただ並べているのがわかります．groupとpartは，現在意味不明です．二つのPVDwriteを足し合わせると, たしか新しいPVDファイルになって，片っ方がpart=0, 片っ方が part=1 になった気がします．