vtkUnstructuredGridの内部構造の探索例

はじめに

OpenFOAMの結果をVTKに渡す際にセルの情報を知る必要がありましたので，VTKの非構造ボリュームメッシュのセルの情報をC++を使って取り出したプログラムを作成しました．

#例
以下の例は，非構造メッシュを読み込んで，各セルを少し縮退させたポリゴンをobjで書き出した例です．

• (参考)vtkUnstructuredGridReaderクラスの定義
• (参考)vtkCellクラスの定義

C++

main.cpp

#include <vtkUnstructuredGridReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCell.h>
#include <fstream>

int main( int argc, char** argv )
{

        vtkUnstructuredGridReader* reader = vtkUnstructuredGridReader::New();
        reader->SetFileName( argv[1] );
        reader->Update();

        std::ofstream fout( "out.obj" );
        if ( !fout ) return -1;

        vtkUnstructuredGrid* grid = reader->GetOutput();
        const int nblock = grid->GetNumberOfCells();

        for ( int i = 0; i < nblock; i++ ) {
                vtkCell* cell = grid->GetCell( i );

                vtkPoints* points = cell->GetPoints();
                const int npoints = points->GetNumberOfPoints();
                double cx = 0, cy = 0, cz = 0;
                for ( int j = 0; j < npoints; ++j ) {
                        double *p = points->GetPoint( j );
                        cx += p[0];
                        cy += p[1];
                        cz += p[2];
                }
                cx *= 1.0 / npoints;
                cy *= 1.0 / npoints;
                cz *= 1.0 / npoints;

                const int nfaces = cell->GetNumberOfFaces();
                for ( int j = 0; j < nfaces; ++j ) {
                        vtkCell* face = cell->GetFace( j ); // face is represented as vtkCell
                        const int nfpoints = face->GetNumberOfPoints();
                        for ( int k = 0; k < nfpoints; ++k ) {
                                double *p0 = face->GetPoints()->GetPoint( k );
                                p0[0] = 0.9 * p0[0] + 0.1 * cx;
                                p0[1] = 0.9 * p0[1] + 0.1 * cy;
                                p0[2] = 0.9 * p0[2] + 0.1 * cz;
                                fout << "v " << p0[0] << " " << p0[1] << " " << p0[2] << std::endl;
                        }
                }

        }

        int count = 1;
        for ( int i = 0; i < nblock; i++ ) {
                vtkCell* cell = grid->GetCell( i );
                const int nfaces = cell->GetNumberOfFaces();
                for ( int j = 0; j < nfaces; ++j ) {
                        fout << "f";
                        vtkCell* face = cell->GetFace( j );
                        const int nfpoints = face->GetNumberOfPoints();
                        for( int k = 0 ; k < nfpoints; ++k ) {
                                fout << " " << count;
                                count++;
                        }
                        fout << std::endl;
                }
        }
        return 0;
}

CMakeLists.txt

cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
project(vtkusg)
set(CMAKE_MODULE_PATH /usr/local/Cellar/vtk/6.3.0/lib/cmake/vtk-6.3)
find_package(VTK 6.3 REQUIRED)

include_directories(${VTK_INCLUDE_DIRS})
link_directories(${VTK_LIBRARY_DIRS})
#add_definitions(${VTK_DEFINITIONS})
add_executable (vtkusg main.cpp)
target_link_libraries (vtkusg ${VTK_LIBRARIES})

###実行

$ ./vtkusg input.vtk

入力ファイルにかかわらず，out.objが出力されます．
実験では，https://github.com/naucoin/VTKData/blob/master/Data/uGridEx.vtk を使いました．

結果は以下の通りです．ファイル内には四面体と６面体がそれぞれ２個ありました．エッジだけのセルもありましたが，それらは出力の際に出てきませんでした．

##C#版(抜粋)
もともとはC#で作っていました．関数はほぼ１対１で対応しています．違いは，C++ではポインタだったのが参照渡しになっているくらいだと思います．残念ながらC#版のリファレンスの整備はまだまだらしいです（学生談）．そのため利用は，C++とにらめっこする必要がありそうです．

main.cs

vtkUnstructuredGridReader reader = vtkUnstructuredGridReader.New();
reader.SetFileName( filePath );
reader.Update();

System.IO.StreamWriter sw = new System.IO.StreamWriter( "test.obj" );
vtkUnstructuredGrid grid = reader.GetOutput();
int nblock = grid.GetNumberOfCells();

for ( int i = 0; i < nblock; i++ )
{
        int cellId = i;
        vtkCell cell = grid.GetCell( cellId );

        vtkPoints points = cell.GetPoints();
        int npoints = points.GetNumberOfPoints();
        double x = 0, y = 0, z = 0;
        for ( int j = 0; j < npoints; ++j ) {
                double[] p = points.GetPoint( j );
                x += p[0];
                y += p[1];
                z += p[2];
        }
        x *= 1.0 / npoints;
        y *= 1.0 / npoints;
        z *= 1.0 / npoints;
        //x,y, z center of points.

        int nfaces = cell.GetNumberOfFaces();
        for ( int j = 0; j < nfaces; ++j ) {
                vtkCell face = cell.GetFace( j );
                for ( int k = 0; k < face.GetNumberOfPoints(); ++k ) {
                        double[] p0 = face.GetPoints().GetPoint( k );
                        p0[0] = 0.9 * p0[0] + 0.1 * x;
                        p0[1] = 0.9 * p0[1] + 0.1 * y;
                        p0[2] = 0.9 * p0[2] + 0.1 * z;

                        sw.WriteLine( "v " + p0[0] + " " + p0[1] + " " + p0[2] );

                }
        }

}

int count = 1;
for ( int i = 0; i < nblock; i++ )
{
        int cellId = i;
        vtkCell cell = grid.GetCell( cellId );

        int nfaces = cell.GetNumberOfFaces();
        for ( int j = 0; j < nfaces; ++j ) {
                vtkCell face = cell.GetFace( j );
                sw.Write( "f " );
                for ( int k = 0; k < face.GetNumberOfPoints(); ++k ) {
                        sw.Write( " " + count );
                        count++;
                }
                sw.WriteLine( " " );

        }
}
sw.Close();