OpenFOAMのpostProcessで出力したVTKをつかってコンターをmatplotlibで描く

# 概要
https://qiita.com/ishigaki/items/ce97e3a4a3088a6ce923

以前，OpenFOAMのpostProcessでVTKを出力する方法について書きました。
そのデータからmatplotlibでコンターを描く方法です。OpenFOAM 8で確かめました。

# VTKのデータ
## OpenFOAM 2.3くらいまでの話
OpenFOAM 2.3くらいまでは，sampleコマンドで出力したVTKでは，物理量はcell上の値が保存されてました。下のCELL_DATA 14324の部分。
また，polygonsは3点使った，三角形になってました。

なので，matplotlibでコンターにするとき，cell dataをpoint dataにマッピングしなおす作業が必要で，tricontourとか使ってました。
こことか参考にして，プログラム書いてました。
https://gist.github.com/BHCLL/0791889ec09686541e3aa38e1942ca21

### 補足
OpenFOAM 6は2.3と同じVTKのフォーマットで出力されてました。

``` sample.vtk
# vtk DataFile Version 2.0
sampleSurface
ASCII
DATASET POLYDATA
POINTS 7401 float
0 0.065 0
0 0.056875 0
0 0.04875 0
中略

POLYGONS 14324 57296
3 17 18 1
3 1 0 17
3 18 19 2
中略
CELL_DATA 14324
FIELD attributes 1
T 1 14324 float
398.985 398.985 398.958 398.958 398.959 398.959 398.973 398.973 398.931 398.931
略
```

## 新し目のOpenFOAM 
OpenFOAM 8で出力したVTKははじめからpoint dataで物理量も出力されるようになってます。（今日知ってびっくりしました）
polygonsも4点使っていて，四角形になってます。

なので，tricontourを使わず，contourでかけます。

``` sample2.vtk
# vtk DataFile Version 2.0
sampleSurface
ASCII
DATASET POLYDATA
POINTS 16900 float
0.00232558 0.00232558 0 0 0.00232558 0 0 0 0 0.00232558 0
0 0.00465116 0.00232558 0 0.00465116 0 0 0.00697674 0.00232558 0
中略
POLYGONS 16641 83205
4 0 1 2 3 4 4 0 3 5 4
6 4 5 7 4 8 6 7 9 4
10 8 9 11 4 12 10 11 13 4
14 12 13 15 4 16 14 15 17 4
中略
POINT_DATA 16900
FIELD attributes 1
T 1 16900 float
309.072 309.073 309.153 309.18 309.066 309.18 309.055 309.18 309.038 309.18 309.016
309.18 308.988 309.18 308.957 309.18 308.926 309.18 308.895 309.18 308.867
略
```

# matplotlibでコンターかく
手順は

1. VTKから点と物理量を読み込む。
2. VTKから読み込んだ座標は順番がへんなので，一旦座標を作り直す
3. contourを描画

以下，例のプログラム。

``` python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import pandas as pd
import glob
import os
from scipy import interpolate
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from matplotlib.colors import Normalize
import matplotlib.colors as colors
import matplotlib

def read_scalarAllVTK(filename):
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.ReadAllScalarsOn()
    reader.ReadAllVectorsOn()
    reader.Update()

    data = reader.GetOutput()
    points = data.GetPoints()
    Valuedata = data.GetPointData().GetArray(0)

    npts = points.GetNumberOfPoints()
    nvls = Valuedata.GetNumberOfTuples()

    print("npts, nvls= ", npts, nvls)
    """ 
    x, y, z: coordinate
    scalar: value
    """
    x = np.zeros(npts)
    y = np.zeros(npts)
    z = np.zeros(npts)
    scalar = np.zeros(nvls)

    for i in range(npts):
        pt = points.GetPoint(i)
        x[i] = pt[0]
        y[i] = pt[1]
        z[i] = pt[2]

    for i in range(0, nvls):
        Value = Valuedata.GetTuple(i)
        scalar[i] = Value[0]

    return x, y, z, scalar

def read_vectorAllVTK(filename):
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.ReadAllScalarsOn()
    reader.ReadAllVectorsOn()
    reader.Update()

    data = reader.GetOutput()
    points = data.GetPoints()

    Valuedata = data.GetPointData().GetArray(0)

    #print("puints=\n",points)
    #print("mapped_data \n", mapped_data)

    npts = points.GetNumberOfPoints()
    nvls = Valuedata.GetNumberOfTuples()

    print("npts, nvls= ",  npts, nvls)
    """ 
    x, y, z: coordinate
    v1, v2, v3: vector value
    """
    x = np.zeros(npts)
    y = np.zeros(npts)
    z = np.zeros(npts)
    v1 = np.zeros(nvls)
    v2 = np.zeros(nvls)
    v3 = np.zeros(nvls)

    for i in range(npts):
        pt = points.GetPoint(i)
        x[i] = pt[0]
        y[i] = pt[1]
        z[i] = pt[2]

    for i in range(0, nvls):
        Value = Valuedata.GetTuple(i)
        v1[i] = Value[0]
        v2[i] = Value[1]
        v3[i] = Value[2]

    return x, y, z, v1, v2, v3

# x,y,z座標と温度TをVTKから読み込み
x, y, z, T = read_scalarAllVTK("T_xy.vtk")

# 座標を作り直すときの分割数。OpenFOAMのpointと同じ個数でOK。129x129セルなので，点は130x130個。
divx=130
divy=130

# xig, yigが作り直した座標。
xi = np.linspace(x.min(), x.max(), divx)
yi = np.linspace(y.min(), y.max(), divy)
xig, yig= np.meshgrid(xi, yi)

# 物理量も座標に合わせて作り直す
Ti = interpolate.griddata((x, y), T, (xi[None,:], yi[:,None]), method='linear')

aspect = 1
fig_width = 6.5 + 0.3
fig_height = fig_width / (x_high - x_low) * (y_high - y_low) * aspect + 0.01

fig = plt.figure(figsize=(fig_width, fig_height))
ax1 = figHe.add_subplot(111)

image1 = ax1.contourf(xig, yig, Ti)                            
                                
```

適宜カラーバーなどを調整してください。