VTK メモ

2019年8月27日

はじめに

データ可視化ライブラリ VTK のメモ。

環境

• Ubuntu 16.04 LTS
• VTK 8.2.0
• CMake

VTK については、ソースからビルドする。apt-get でもよいのだが、Anaconda3 の Qt と一緒に利用しようとすると、うまく動かなかった (Anaconda3 の Qt が使いたかったというより、Qt と Anaconda3 の Python を使おうとすると、Anaconda3 の Qt が使われてしまうため)。

CMake は apt-get でインストールすればよい。

VTK を用いたプログラムのコンパイル

VTK を用いたプログラムのコンパイルには、CMake を用いる。プロジェクトディレクトリを作り、CMakeLists.txt を用意する。

cmake_minimum_required(VERSION 2.8)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

file(GLOB CXX_FILES *.cxx)

add_executable(cylinder ${CXX_FILES})
target_link_libraries(cylinder ${VTK_LIBRARIES})

コンパイルはビルドディレクトリを作り、その中で行う。

$ mkdir build
$ cd build
$ cmake ..
$ make

円筒の表示

cylinder.cxx

#include <vtkSmartPointer.h>
#include <vtkCylinderSource.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkCamera.h>

using namespace std;


int main(int argc, char **argv)
{
    // source
    auto cylinder = vtkSmartPointer<vtkCylinderSource>::New();
    cylinder->SetResolution(20);

    // mapper
    auto mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(cylinder->GetOutputPort());

    // actor
    auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);
    actor->GetProperty()->SetColor(1., .3882, .2784 );
    actor->RotateX(30.);
    actor->RotateY(-45.);

    // renderer
    auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->SetBackground(0.1, 0.2, 0.4);

    auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);

    ren->ResetCamera();
    ren->GetActiveCamera()->Zoom(1.5);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

マウスイベント

mouse.cxx

#include <vtkSmartPointer.h>
#include <vtkCylinderSource.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkCamera.h>

using namespace std;


class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void OnLeftButtonDown() noexcept override
    {
        cout << "LEFT" << endl;
        vtkInteractorStyleTrackballCamera::OnLeftButtonDown();
    }

    void OnRightButtonDown() noexcept override
    {
        cout << "RIGHT" << endl;
        vtkInteractorStyleTrackballCamera::OnRightButtonDown();
    }
};


int main(int argc, char **argv)
{
    // source
    auto cylinder = vtkSmartPointer<vtkCylinderSource>::New();
    cylinder->SetResolution(20);

    // mapper
    auto mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(cylinder->GetOutputPort());

    // actor
    auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);
    actor->GetProperty()->SetColor(1., .3882, .2784 );
    actor->RotateX(30.);
    actor->RotateY(-45.);

    // renderer
    auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->SetBackground(0.1, 0.2, 0.4);

    auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    ren->ResetCamera();
    ren->GetActiveCamera()->Zoom(1.5);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

キー入力

key.cxx

#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleTrackballCamera.h>

using namespace std;


class KeyInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void OnKeyPress() noexcept override
    {
        cout << GetInteractor()->GetKeySym() << endl;
        vtkInteractorStyleTrackballCamera::OnKeyPress();
    }
};


int main(int argc, char **argv)
{
    // renderer
    auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->SetBackground(0.1, 0.2, 0.4);

    auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new KeyInteractorStyle());

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

STL の読み書き

stl_writer.cxx

#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkSTLWriter.h>
#include <vtkSTLReader.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

using namespace std;

const char *filename = "sphere.stl";


int main(int argc, char **argv)
{
    // source
    auto sphere = vtkSmartPointer<vtkSphereSource>();
    sphere->SetThetaResolution(50);
    sphere->SetPhiResolution(50);

    // writer
    auto writer = vtkSmartPointer<vtkSTLWriter>();
    writer->SetFileName(filename);
    writer->SetInputConnection(sphere->GetOutputPort());
    writer->Write();

    // reader
    auto reader = vtkSmartPointer<vtkSTLReader>();
    reader->SetFileName(filename);

    // mapper
    auto mapper = vtkSmartPointer<vtkPolyDataMapper>();
    mapper->SetInputConnection(reader->GetOutputPort());

    // actor
    auto actor = vtkSmartPointer<vtkActor>();
    actor->SetMapper(mapper);

    // renderer
    auto ren = vtkSmartPointer<vtkRenderer>();
    ren->AddActor(actor);
    ren->SetBackground(0.1, 0.2, 0.4);

    auto renWin = vtkSmartPointer<vtkRenderWindow>();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    auto inter = vtkSmartPointer<vtkRenderWindowInteractor>();
    inter->SetRenderWindow(renWin);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

VTK の読み書き

VTK ファイル (Unstructure Grid) を読み込む。

scalar.vtk

# vtk DataFile Version 2.0
scalar
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
0 0 0
1 0 0
1 1 0
0 1 0
CELLS 1 5
4 0 1 2 3
CELL_TYPES 1
9
POINT_DATA 4
SCALARS point_scalars float
LOOKUP_TABLE default
1
2
3
4
CELL_DATA 1
SCALARS cell_scalars float
LOOKUP_TABLE default
5

vtk_reader.cxx

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkUnstructuredGridReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

using namespace std;

const char *filename = "scalar.vtk";


int main(int argc, char **argv)
{
    // reader
    auto reader = vtkSmartPointer<vtkUnstructuredGridReader>();
    reader->SetFileName(filename);
    reader->Update();

    auto output = reader->GetOutput();
    auto scalarRange = output->GetScalarRange();
    cout << scalarRange << endl;

    // mapper
    auto mapper = vtkSmartPointer<vtkDataSetMapper>();
    mapper->SetInputData(output);
    mapper->SetScalarRange(scalarRange);

    // actor
    auto actor = vtkSmartPointer<vtkActor>();
    actor->SetMapper(mapper);

    // renderer
    auto ren = vtkSmartPointer<vtkRenderer>();
    ren->AddActor(actor);
    ren->SetBackground(0.1, 0.2, 0.4);

    auto renWin = vtkSmartPointer<vtkRenderWindow>();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    auto inter = vtkSmartPointer<vtkRenderWindowInteractor>();
    inter->SetRenderWindow(renWin);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

OpenFOAM のケースの読み込み

OpenFOAM の計算結果を読み込む。ここでは cavity ケースを用いる。

• cavity.tar.gz

of_reader.cxx

#include <iostream>
#include <string>
#include <vtkSmartPointer.h>
#include <vtkDoubleArray.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOpenFOAMReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkGeometryFilter.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkLookupTable.h>
#include <vtkCompositePolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkScalarBarActor.h>
#include <vtkTextProperty.h>
#include <vtkTextActor.h>

using namespace std;

const string caseName = "cavity";
const string filename = caseName + "/system/controlDict";
const string field = "U";
const string dataType = "cell";


class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void Dolly() override
    {
        MotionFactor *= -1;
        vtkInteractorStyleTrackballCamera::Dolly();
        MotionFactor *= -1;
    }
};


int main(int argc, char **argv)
{
    // reader
    const auto reader = vtkSmartPointer<vtkOpenFOAMReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();

    //const auto n = reader->GetTimeValues()->GetNumberOfValues();
    //const auto latestTime = reader->GetTimeValues()->GetValue(n - 1);
    const auto latestTime = reader->GetTimeValues()->GetRange()[1];

    reader->UpdateTimeStep(latestTime);
    reader->Update();

    const auto block = reader->GetOutput();
    //cout << block->GetBlock(0)->GetClassName() << endl;
    const auto grid = vtkUnstructuredGrid::SafeDownCast(block->GetBlock(0));

    string arrayName;
    vtkDataArray *array;

    if(field == "p"){
        vtkDataArray *p;
        if(dataType == "cell"){
            p = grid->GetCellData()->GetArray("p");
        }else if(dataType == "point"){
            p = grid->GetPointData()->GetArray("p");
        }

        arrayName = "p";
        array = p;
    }else if(field == "U"){
        const auto magU = vtkSmartPointer<vtkDoubleArray>::New();
        magU->SetName("magU");
        magU->SetNumberOfComponents(1);

        vtkDataArray *U;
        if(dataType == "cell"){
            U = grid->GetCellData()->GetArray("U");
        }else if(dataType == "point"){
            U = grid->GetPointData()->GetArray("U");
        }else{
            assert(0);
        }

        for(int i = 0; i < U->GetNumberOfTuples(); i++){
            const auto Ux = U->GetComponent(i, 0);
            const auto Uy = U->GetComponent(i, 1);
            const auto Uz = U->GetComponent(i, 2);
            magU->InsertNextValue(sqrt(Ux*Ux + Uy*Uy + Uz*Uz));
        }

        if(dataType == "cell"){
            grid->GetCellData()->AddArray(magU);
        }else if(dataType == "point"){
            grid->GetPointData()->AddArray(magU);
        }else{
            assert(0);
        }

        arrayName = "magU";
        array = magU;
    }else{
        assert(0);
    }

    // filter
    const auto geomFilter = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter->SetInputData(block);

    // lookup table
    const auto lut = vtkSmartPointer<vtkLookupTable>::New();
    lut->SetHueRange(0.667, 0.);
    lut->Build();

    // mapper
    const auto mapper = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper->SetInputConnection(geomFilter->GetOutputPort());
    if(dataType == "cell"){
        mapper->SetScalarModeToUseCellFieldData();
    }else if(dataType == "point"){
        mapper->SetScalarModeToUsePointFieldData();
    }
    mapper->SelectColorArray(arrayName.c_str());
    mapper->SetScalarRange(array->GetRange());
    mapper->SetLookupTable(lut);

    // actor
    const auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    const auto prop = actor->GetProperty();
    prop->SetAmbient(0.5);

    prop->EdgeVisibilityOn();
    prop->SetEdgeColor(0, 0, 0);
    prop->SetLineWidth(2);

    prop->SetRepresentationToSurface();

    // renderer
    const auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->GradientBackgroundOn();
    ren->SetBackground(1., 1., 1.);
    ren->SetBackground2(0.4, 0.55, 0.75);

    const auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    const auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    // sclar bar
    const auto scalarBar = vtkSmartPointer<vtkScalarBarActor>::New();
    scalarBar->SetOrientationToVertical();
    scalarBar->SetLookupTable(lut);
    const auto textProp = vtkSmartPointer<vtkTextProperty>::New();
    textProp->SetColor(0., 0., 0.);
    textProp->SetFontSize(20);
    textProp->SetFontFamilyToArial();
    textProp->ItalicOff();
    textProp->BoldOff();
    scalarBar->UnconstrainedFontSizeOn();
    scalarBar->SetTitleTextProperty(textProp);
    scalarBar->SetLabelTextProperty(textProp);
    scalarBar->SetLabelFormat("%5.2f");
    scalarBar->SetTitle(arrayName.c_str());
    ren->AddActor(scalarBar);

    // text
    const auto text = vtkSmartPointer<vtkTextActor>::New();
    text->SetInput(caseName.c_str());
    text->GetTextProperty()->SetColor(0., 0., 0.);
    text->GetTextProperty()->SetFontSize(24);
    text->SetPosition(300, 30);
    ren->AddActor(text);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

ここでは、マウスの操作を ParaView に合わせるために dolly の方向を逆にしている。

OpenFOAM のマルチリージョンケースの読み込み

OpenFOAM のマルチリージョンケースの計算結果を読み込む。ここでは multiRegionHeater ケースを用いる。

• multiRegionHeater.tar.gz

of_reader_multi.cxx

#include <iostream>
#include <string>
#include <limits>
#include <vtkSmartPointer.h>
#include <vtkDoubleArray.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOpenFOAMReader.h>
#include <vtkCompositeDataIterator.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkGeometryFilter.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkLookupTable.h>
#include <vtkCompositePolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkScalarBarActor.h>
#include <vtkTextProperty.h>
#include <vtkTextActor.h>

using namespace std;

const string caseName = "multiRegionHeater";
const string filename = caseName + "/system/controlDict";

const auto enabledBlock = {
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
};

const string field = "T";
const string dataType = "point";


class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void Dolly() override
    {
        MotionFactor *= -1;
        vtkInteractorStyleTrackballCamera::Dolly();
        MotionFactor *= -1;
    }
};


int main(int argc, char **argv)
{
    // reader
    const auto reader = vtkSmartPointer<vtkOpenFOAMReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();

    const auto latestTime = reader->GetTimeValues()->GetRange()[1];
    cout << "Latest time: " << latestTime << endl;

    reader->UpdateTimeStep(latestTime);
    reader->Update();


    const auto n = reader->GetNumberOfPatchArrays();
    for(int i = 0; i < n; i++){
        cout << reader->GetPatchArrayName(i) << endl;
    }

    reader->DisableAllPatchArrays();
    for(const auto& b : enabledBlock){
        reader->SetPatchArrayStatus(b, 1);
    }
    reader->Update();


    const auto block = reader->GetOutput();

    const auto arrayName = field;
    auto arrayMin = numeric_limits<double>::max();
    auto arrayMax = numeric_limits<double>::min();

    const auto itr = block->NewIterator();
    itr->InitTraversal();
    while(!itr->IsDoneWithTraversal()){
        //cout << itr->GetCurrentDataObject()->GetClassName() << endl;
        const auto grid = vtkUnstructuredGrid::SafeDownCast(
                itr->GetCurrentDataObject()
        );

        vtkDataArray *T;
        if(dataType == "cell"){
            T = grid->GetCellData()->GetArray(field.c_str());
        }else if(dataType == "point"){
            T = grid->GetPointData()->GetArray(field.c_str());
        }else{
            assert(0);
        }

        auto range = T->GetRange();
        if(range[0] < arrayMin){
            arrayMin = range[0];
        }
        if(range[1] > arrayMax){
            arrayMax = range[1];
        }

        itr->GoToNextItem();
    }


    // filter
    const auto geomFilter = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter->SetInputData(reader->GetOutput());

    // lookup table
    const auto lut = vtkSmartPointer<vtkLookupTable>::New();
    lut->SetHueRange(0.667, 0.);
    lut->Build();

    // mapper
    const auto mapper = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper->SetInputConnection(geomFilter->GetOutputPort());
    if(dataType == "cell"){
        mapper->SetScalarModeToUseCellFieldData();
    }else if(dataType == "point"){
        mapper->SetScalarModeToUsePointFieldData();
    }
    mapper->SelectColorArray(arrayName.c_str());
    mapper->SetScalarRange(arrayMin, arrayMax);
    mapper->SetLookupTable(lut);

    // actor
    const auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    const auto prop = actor->GetProperty();
    prop->SetAmbient(0.5);

    prop->EdgeVisibilityOn();
    prop->SetEdgeColor(0., 0., 0.);
    prop->SetLineWidth(2);

    prop->SetRepresentationToSurface();

    // renderer
    const auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->GradientBackgroundOn();
    ren->SetBackground(1., 1., 1.);
    ren->SetBackground2(0.4, 0.55, 0.75);

    const auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    const auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    // sclar bar
    const auto scalarBar = vtkSmartPointer<vtkScalarBarActor>::New();
    scalarBar->SetOrientationToVertical();
    scalarBar->SetLookupTable(lut);
    const auto textProp = vtkSmartPointer<vtkTextProperty>::New();
    textProp->SetColor(0., 0., 0.);
    textProp->SetFontSize(20);
    textProp->SetFontFamilyToArial();
    textProp->ItalicOff();
    textProp->BoldOff();
    scalarBar->UnconstrainedFontSizeOn();
    scalarBar->SetTitleTextProperty(textProp);
    scalarBar->SetLabelTextProperty(textProp);
    scalarBar->SetLabelFormat("%5.2f");
    scalarBar->SetTitle(arrayName.c_str());
    ren->AddActor(scalarBar);

    // text
    const auto text = vtkSmartPointer<vtkTextActor>::New();
    text->SetInput(caseName.c_str());
    text->GetTextProperty()->SetColor(0., 0., 0.);
    text->GetTextProperty()->SetFontSize(24);
    text->SetPosition(200, 30);
    ren->AddActor(text);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

ブロックの抽出

上の例では、ブロックの選択を reader の SetPatchArrayStatus() で行なったが、これを vtkExtractBlock で行うこともできる。ParaView で、OpenFOAM リーダーの設定でブロックを選択するか、Extract Block フィルタで選択するかというのと同じである。

extract_block.cxx

#include <iostream>
#include <string>
#include <regex>
#include <limits>
#include <vtkSmartPointer.h>
#include <vtkDoubleArray.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOpenFOAMReader.h>
#include <vtkExtractBlock.h>
#include <vtkCompositeDataIterator.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkGeometryFilter.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkLookupTable.h>
#include <vtkCompositePolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkScalarBarActor.h>
#include <vtkTextProperty.h>
#include <vtkTextActor.h>

using namespace std;

const string caseName = "multiRegionHeater";
const string filename = caseName + "/system/controlDict";

const auto enabledBlock = {
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
};

const string field = "T";
const string dataType = "point";


class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void Dolly() override
    {
        MotionFactor *= -1;
        vtkInteractorStyleTrackballCamera::Dolly();
        MotionFactor *= -1;
    }
};


int main(int argc, char **argv)
{
    // reader
    const auto reader = vtkSmartPointer<vtkOpenFOAMReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();

    const auto latestTime = reader->GetTimeValues()->GetRange()[1];
    cout << "Latest time: " << latestTime << endl;

    reader->UpdateTimeStep(latestTime);
    reader->Update();


    //
    // create patch to index table
    //
    reader->EnableAllPatchArrays();
    reader->Update();

    map<string, int> blockIndex;
    auto index = 1;
    const auto n = reader->GetNumberOfPatchArrays();
    for(int i = 0; i < n; i++){
        const auto name = reader->GetPatchArrayName(i);

        cmatch m;
        regex_match(name, m, regex(".*internalMesh$"));
        if(!m.empty()){
            index++;
        }
        blockIndex[name] = index;
        if(!m.empty()){
            index++;
        }
        index++;
    }

    for(int i = 0; i < n; i++){
        const auto name = reader->GetPatchArrayName(i);
        cout << name << " " << blockIndex[name] << endl;
    }


    const auto extractBlock = vtkSmartPointer<vtkExtractBlock>::New();
    extractBlock->SetInputData(reader->GetOutput());
    for(const auto& b : enabledBlock){
        extractBlock->AddIndex(blockIndex[b]);
    }
    extractBlock->Update();


    const auto block = extractBlock->GetOutput();

    const auto arrayName = field;
    auto arrayMin = numeric_limits<double>::max();
    auto arrayMax = numeric_limits<double>::min();

    const auto itr = block->NewIterator();
    itr->InitTraversal();
    while(!itr->IsDoneWithTraversal()){
        //cout << itr->GetCurrentDataObject()->GetClassName() << endl;
        const auto grid = vtkUnstructuredGrid::SafeDownCast(
                itr->GetCurrentDataObject()
        );

        vtkDataArray *T;
        if(dataType == "cell"){
            T = grid->GetCellData()->GetArray(field.c_str());
        }else if(dataType == "point"){
            T = grid->GetPointData()->GetArray(field.c_str());
        }else{
            assert(0);
        }

        auto range = T->GetRange();
        if(range[0] < arrayMin){
            arrayMin = range[0];
        }
        if(range[1] > arrayMax){
            arrayMax = range[1];
        }

        itr->GoToNextItem();
    }


    // filter
    const auto geomFilter = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter->SetInputConnection(extractBlock->GetOutputPort());

    // lookup table
    const auto lut = vtkSmartPointer<vtkLookupTable>::New();
    lut->SetHueRange(0.667, 0.);
    lut->Build();

    // mapper
    const auto mapper = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper->SetInputConnection(geomFilter->GetOutputPort());
    if(dataType == "cell"){
        mapper->SetScalarModeToUseCellFieldData();
    }else if(dataType == "point"){
        mapper->SetScalarModeToUsePointFieldData();
    }
    mapper->SelectColorArray(arrayName.c_str());
    mapper->SetScalarRange(arrayMin, arrayMax);
    mapper->SetLookupTable(lut);

    // actor
    const auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    const auto prop = actor->GetProperty();
    prop->SetAmbient(0.5);

    prop->EdgeVisibilityOn();
    prop->SetEdgeColor(0., 0., 0.);
    prop->SetLineWidth(2);

    prop->SetRepresentationToSurface();

    // renderer
    const auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->GradientBackgroundOn();
    ren->SetBackground(1., 1., 1.);
    ren->SetBackground2(0.4, 0.55, 0.75);

    const auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    const auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    // sclar bar
    const auto scalarBar = vtkSmartPointer<vtkScalarBarActor>::New();
    scalarBar->SetOrientationToVertical();
    scalarBar->SetLookupTable(lut);
    const auto textProp = vtkSmartPointer<vtkTextProperty>::New();
    textProp->SetColor(0., 0., 0.);
    textProp->SetFontSize(20);
    textProp->SetFontFamilyToArial();
    textProp->ItalicOff();
    textProp->BoldOff();
    scalarBar->UnconstrainedFontSizeOn();
    scalarBar->SetTitleTextProperty(textProp);
    scalarBar->SetLabelTextProperty(textProp);
    scalarBar->SetLabelFormat("%5.2f");
    scalarBar->SetTitle(arrayName.c_str());
    ren->AddActor(scalarBar);

    // text
    const auto text = vtkSmartPointer<vtkTextActor>::New();
    text->SetInput(caseName.c_str());
    text->GetTextProperty()->SetColor(0., 0., 0.);
    text->GetTextProperty()->SetFontSize(24);
    text->SetPosition(200, 30);
    ren->AddActor(text);

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

実行結果は of_reader_multi.py と同じである。

ハイライト表示

注目しているブロックをハイライト表示したい。対象のブロックだけ Extract Block で取り出して、色をつければよい。

highlight.cxx

#include <iostream>
#include <string>
#include <regex>
#include <limits>
#include <vtkSmartPointer.h>
#include <vtkIntArray.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOpenFOAMReader.h>
#include <vtkExtractBlock.h>
#include <vtkCompositeDataIterator.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkGeometryFilter.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkLookupTable.h>
#include <vtkCompositePolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

using namespace std;

const string caseName = "multiRegionHeater";
const string filename = caseName + "/system/controlDict";

const vector<string> enabledBlock = {
    "bottomWater/internalMesh",
    "topAir/internalMesh",
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
};

const auto highlightBlock = {
    "topAir/minX",
};


class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void Dolly() override
    {
        MotionFactor *= -1;
        vtkInteractorStyleTrackballCamera::Dolly();
        MotionFactor *= -1;
    }
};


int main(int argc, char **argv)
{
    // reader
    const auto reader = vtkSmartPointer<vtkOpenFOAMReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();


    //
    // create patch to index table
    //
    reader->EnableAllPatchArrays();
    reader->Update();

    map<string, int> blockIndex;
    auto index = 1;
    const auto n = reader->GetNumberOfPatchArrays();
    for(int i = 0; i < n; i++){
        auto name = reader->GetPatchArrayName(i);

        cmatch m;
        regex_match(name, m, regex(".*internalMesh$"));
        if(!m.empty()){
            index++;
        }
        blockIndex[name] = index;
        if(!m.empty()){
            index++;
        }
        index++;
    }

    for(int i = 0; i < n; i++){
        auto name = reader->GetPatchArrayName(i);
        cout << name << " " << blockIndex[name] << endl;
    }


    const auto extractBlock = vtkSmartPointer<vtkExtractBlock>::New();
    extractBlock->SetInputData(reader->GetOutput());
    for(const auto& b : enabledBlock){
        extractBlock->AddIndex(blockIndex[b]);
    }
    extractBlock->Update();

    const auto extractBlock2 = vtkSmartPointer<vtkExtractBlock>::New();
    extractBlock2->SetInputData(reader->GetOutput());
    for(const auto& b : highlightBlock){
        extractBlock2->AddIndex(blockIndex[b]);
    }
    extractBlock2->Update();


    //
    // block index
    //
    const auto block = extractBlock->GetOutput();

    const string arrayName = "index";
    auto arrayMin = numeric_limits<int>::max();
    auto arrayMax = numeric_limits<int>::min();

    for(const auto& item : blockIndex){
        if(item.second < arrayMin){
            arrayMin = item.second;
        }else if(item.second > arrayMax){
            arrayMax = item.second;
        }
    }

    const auto itr = block->NewIterator();
    itr->InitTraversal();
    auto i = 0;
    while(!itr->IsDoneWithTraversal()){
        //cout << itr->GetCurrentDataObject()->GetClassName() << endl;
        const auto grid = vtkUnstructuredGrid::SafeDownCast(
                itr->GetCurrentDataObject()
        );

        const auto index = vtkSmartPointer<vtkIntArray>::New();
        index->SetName(arrayName.c_str());
        index->SetNumberOfComponents(1);

        for(int j = 0; j < grid->GetNumberOfCells(); j++){
            index->InsertNextValue(blockIndex[enabledBlock[i]]);
        }
        grid->GetCellData()->AddArray(index);

        itr->GoToNextItem();
        i++;
    }


    // filter
    const auto geomFilter = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter->SetInputConnection(extractBlock->GetOutputPort());

    const auto geomFilter2 = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter2->SetInputConnection(extractBlock2->GetOutputPort());

    // lookup table
    const auto lut = vtkSmartPointer<vtkLookupTable>::New();
    lut->SetHueRange(0.667, 0.);
    lut->Build();

    // mapper
    const auto mapper = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper->SetInputConnection(geomFilter->GetOutputPort());
    mapper->SetScalarModeToUseCellFieldData();
    mapper->SelectColorArray(arrayName.c_str());
    mapper->SetScalarRange(arrayMin, arrayMax);
    mapper->SetLookupTable(lut);

    const auto mapper2 = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper2->SetInputConnection(geomFilter2->GetOutputPort());
    mapper2->SetScalarModeToUseCellFieldData();

    // actor
    const auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    const auto prop = actor->GetProperty();
    prop->SetAmbient(0.5);

    prop->EdgeVisibilityOn();
    prop->SetEdgeColor(0., 0., 0.);
    prop->SetLineWidth(2);

    prop->SetRepresentationToSurface();

    const auto actor2 = vtkSmartPointer<vtkActor>::New();
    actor2->SetMapper(mapper2);

    const auto prop2 = actor2->GetProperty();
    prop2->SetColor(1., 0., 1.);
    prop2->SetOpacity(0.8);

    // renderer
    const auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->AddActor(actor2);
    ren->GradientBackgroundOn();
    ren->SetBackground(1., 1., 1.);
    ren->SetBackground2(0.4, 0.55, 0.75);

    const auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    const auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}

特徴線

特徴線は vtkFeatureEdges で描画できる。OpenFOAM のマルチリージョンケースで特徴線を描画してみる。

feature_edge.cxx

#include <iostream>
#include <string>
#include <vtkSmartPointer.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOpenFOAMReader.h>
#include <vtkGeometryFilter.h>
#include <vtkFeatureEdges.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkCompositePolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

using namespace std;

const string caseName = "multiRegionHeater";
const string filename = caseName + "/system/controlDict";

const auto enabledBlock ={
    "bottomWater/internalMesh",
    "topAir/internalMesh",
    "rightSolid/internalMesh",
    "leftSolid/internalMesh",
    "heater/internalMesh",
};

class MouseInteractorStyle
    : public vtkInteractorStyleTrackballCamera
{
public:
    void Dolly() override
    {
        MotionFactor *= -1;
        vtkInteractorStyleTrackballCamera::Dolly();
        MotionFactor *= -1;
    }
};


int main(int argc, char **argv)
{
    // reader
    const auto reader = vtkSmartPointer<vtkOpenFOAMReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();

    reader->DisableAllPatchArrays();
    for(const auto& b : enabledBlock){
        reader->SetPatchArrayStatus(b, 1);
    }
    reader->Update();

    // filter
    const auto geomFilter = vtkSmartPointer<vtkGeometryFilter>::New();
    geomFilter->SetInputData(reader->GetOutput());

    const auto featureEdge = vtkSmartPointer<vtkFeatureEdges>::New();
    featureEdge->SetInputConnection(geomFilter->GetOutputPort());

    // mapper
    const auto mapper = vtkSmartPointer<vtkCompositePolyDataMapper>::New();
    mapper->SetInputConnection(featureEdge->GetOutputPort());
    mapper->SetScalarModeToUseCellFieldData();

    // actor
    const auto actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    const auto prop = actor->GetProperty();
    prop->SetColor(0., 0., 0.);
    prop->SetLineWidth(2);
    prop->SetRepresentationToSurface();

    // renderer
    const auto ren = vtkSmartPointer<vtkRenderer>::New();
    ren->AddActor(actor);
    ren->GradientBackgroundOn();
    ren->SetBackground(1., 1., 1.);
    ren->SetBackground2(0.4, 0.55, 0.75);

    const auto renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren);
    renWin->SetSize(640, 480);

    // interactor
    const auto inter = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    inter->SetRenderWindow(renWin);
    inter->SetInteractorStyle(new MouseInteractorStyle());

    renWin->Render();
    inter->Initialize();
    inter->Start();

    return 0;
}