メインコンテンツに移動

理解したくない

理解したくない.とりあえず使いたいからサンプルよこせというむけに, サンプルクラスを提供しましょう.

  • データを書き出すフォルダーを指定してインスタンスを作成
  • 等分割直交格子領域を好きなだけ追加
  • そこに格子点データを好きなだけ追加
  • ファイル名を指定して保存すると,XML-VTKで保存する

代物です.

使い方としては

#include "vtk_image.hpp"
...
VTK_IMAGE Writer("fucking_vtk");  //ホームフォルダ$HOME/fucking_vtkフォルダーに出力
Writer.Parameter=&あんたのパラメータ; //あとで値設定関数にパラメータを与える場合,こうする
...
//領域を追加していく.いくつ足しても良い
double x_r[3]={-5.,+5.,0.01};     //  X座標の範囲は[-5,5] 格子幅0.01
double y_r[3]={0.0,+5.,0.01};     //  Y座標の範囲は[0,5] 格子幅0.01
// Z座標を定義しても良い. ここでは2D問題ということでZ格子点数が1
int num[3]={int(round((x_r[1]-x_r[0])/x_r[2])),int(round((y_r[1]-y_r[0])/y_r[2])),1};
Eigen::Vector3d ORG={x_r[0],y_r[0],0.};
Eigen::Vector3d DX={x_r[2],y_r[2],0.};
Eigen::Vector3i NX={num[0],num[1],1};
int myAssHole=Writer.AddImage(ORG,DX,NX);
//  領域に整数の番号が与えられる
...
//格子点データを追加していく
{//座標を追加
  std::vector<std::string> name = { "x_1", "x_2", "x_3"};  //成分数だけ名前を定義
  Writer.AddPointData("x",name);                           //総称名,成分名で追加
  Writer.SetPointData(myAssHole,"x",[](Eigen::Vector3d X,void* p){
    return X.data(); //Eigenはdata()が配列へのポインターを与える
  }); //  データ追加関数を与えて格子データをある領域に追加
  //データ追加関数は, function(Eigen::Vector3d X,void* p)に固定である. X座標に応じて
  //値へのポインタを返すラムダ式を与える. ベクトルの場合, 配列へのポインタを返す.
  //pは,ここで与えたあんたのパラメータであるので必要なら使う. パラメータは単一の値でも良いし
  //複雑怪奇なクラスへのポインターでも良い.
}
{//別の関数も追加. ここではスカラー関数を追加してみる
  Writer.AddPointData("n");
  Writer.SetPointData(myAssHole,"n",[](Eigen::Vector3d X,void* p){
    auto& Motherfucker=*(あんたのパラメータクラス *) p;
    double th1=atan2(X[1],X[0]+0.5);
    double th2=atan2(X[1],X[0]-0.5);
    double n_val=Motherfucker.int_000(th1,th2);//あんたのクラスの計算は禁止されていない
    return &n_val;
  });
}
...
Writer.Save("baka.vtm");  //フォルダーにVTMファイルが保存される

これを使うには, 次のコードをあなたのプロジェクトに(たとえば vtk_image.hpp とかの名前で)追加します.なお

  • Eigen, Boost/FileSystemを使います.
  • libvtkCommonCore, libvtkCommonDataModel, libvtkIOXMLを使います.

ので,ライブラリーの設定は必要です.

#ifndef vtk_image_hpp
#define vtk_image_hpp
#include <Eigen/Dense>
#include <boost/filesystem.hpp>
#include <vtkSmartPointer.h>
#include <vtkMultiBlockDataSet.h> //1ケースに多数の領域が含まれる場合,これが必要
#include <vtkXMLMultiBlockDataWriter.h> //多数の領域を一発で保存する
#include <vtkImageData.h> //等分割直方格子の領域を使います
#include <vtkDoubleArray.h> //セルや格子点に取り付ける物理量を定義します. 
#include <vtkPointData.h> //格子点のデータ
class VTK_IMAGE
{
  boost::filesystem::path folder;
  vtkSmartPointer<vtkMultiBlockDataSet> multi;
  vtkSmartPointer<vtkImageData> current_image;
  vtkSmartPointer<vtkDataArray> current_array;
  vtkSmartPointer<vtkDoubleArray> current_double;
  int num_domain=0;
  int current_components=1;
  bool SelectPointData(int id,std::string name)
  {
  //指定領域を選択
    try
    {
      current_image = vtkImageData::SafeDownCast(multi->GetBlock(id));
    }
    catch (std::exception& e)
    {
      std::cout << "No such imageだぴょーん[id= " << id << "]" << std::endl;
      return false;
    }
  //指定関数を選択
    try
    {
      for(int ifunc=0;ifunc<current_image->GetPointData()->GetNumberOfArrays();ifunc++)
      {
        current_array=current_image->GetPointData()->GetArray(ifunc);
        current_components=current_array->GetNumberOfComponents();
        if (current_array->GetName()==name)
        {
          current_double=vtkDoubleArray::SafeDownCast(current_array);
          return true;
        }
        if (current_components==1) continue;
        //うっかり成分名を指定した場合に助ける
        for(int ic=0;ic<current_components;ic++)
          if (current_array->GetComponentName(ic)==name) return true;
      }
    }
    catch (std::exception& e)
    {
      std::cout << "No such componentだぴょーん[id= " << id << "]" << std::endl;
      return false;
    }
    return false;
  };
public:
  void* Parameter;
  VTK_IMAGE(std::string path):folder(path)
  {//書き出しフォルダーを指定して初期化
    boost::filesystem::path homedir=getenv("HOME");
    if (!boost::filesystem::exists(homedir/folder))
      boost::filesystem::create_directory(homedir/folder);
    boost::filesystem::current_path(homedir/folder);
    multi=vtkSmartPointer<vtkMultiBlockDataSet>::New();
  }
  int AddImage(Eigen::Vector3d origin,Eigen::Vector3d space,Eigen::Vector3i num)
  {//領域を追加
    current_image=vtkSmartPointer<vtkImageData>::New();
    current_image->SetOrigin(origin.data()); //領域の原点
    current_image->SetSpacing(space.data()); //X-Y-Z方向の格子の間隔
    int imageExtent[6]={0,num[0],0,num[1],0,num[2]}; //領域の格子点数
    current_image->SetExtent(imageExtent); //領域の原点
    int id=num_domain++;
    multi->SetBlock(id,current_image);
    return id;
  }
  bool AddPointData(std::string name)
  {//スカラー関数を追加
    for(int id=0;id<multi->GetNumberOfBlocks();id++)
    {
      current_image = vtkImageData::SafeDownCast(multi->GetBlock(id));
      auto function = vtkDoubleArray::New();
      function->SetNumberOfComponents(1);
      function->SetName(name.c_str());
      function->SetNumberOfTuples(current_image->GetNumberOfPoints());
      int ifunc=current_image->GetPointData()->AddArray(function);
    }
    return true;
  }
  bool AddPointData(std::string name,std::vector<std::string> names)
  {//ベクトル関数を追加
    for(int id=0;id<multi->GetNumberOfBlocks();id++)
    {
      current_image = vtkImageData::SafeDownCast(multi->GetBlock(id));
      auto function = vtkDoubleArray::New();
      function->SetNumberOfComponents(names.size());
      function->SetName(name.c_str());
      for(int ic=0;ic<names.size();ic++)
        function->SetComponentName(ic, names[ic].c_str());
      function->SetNumberOfTuples(current_image->GetNumberOfPoints());
      int ifunc=current_image->GetPointData()->AddArray(function);
    }
    return true;
  }
  bool Save(std::string filename)
  {//ファイル名を指定して書き出し
    boost::filesystem::path file(filename);
    auto writer=vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
    writer->SetFileName(file.c_str()); //掃除屋にファイル名を設定
    writer->SetInputData(multi); //掃除屋にボディーを引き渡す
    writer->Write(); //掃除を実行
  }
  bool SetPointData(int id,std::string name,double* (*f)(Eigen::Vector3d x,void* p))
  {
    if (!SelectPointData(id,name)) return false;
    for(vtkIdType pId=0;pId<current_image->GetNumberOfPoints();pId++)
    {
      Eigen::Vector3d XCoord;
      current_image->GetPoint(pId,XCoord.data());
      auto* fval=(*f)(XCoord,Parameter);
      if (current_components==1)
        current_double->SetValue(pId,fval[0]);
      else
        for(int ic=0;ic<current_components;ic++)
      current_double->SetComponent(pId,ic,fval[ic]);
    }
  };
};
#endif /* vtk_image_hpp */

以上