最近のOpenFOAMでsurfaceデータを三角形要素のVTKで出力する

# OpenFOAM9くらいを想定してます。

https://qiita.com/ishigaki/items/ce97e3a4a3088a6ce923

https://qiita.com/ishigaki/items/fbdb7697af45b4411f6b

基本的なやりかたは上の記事で書いてます。

以前，

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

と書きました。

OpenFOAM9だと，cutPlaneSurfaceってのがOpenFOAM/OpenFOAM-9/etc/caseDicts/postProcessing/surface/cutPlaneSurfaceにあります。
これつかって，

```
postProcess -func cutPlaneSurface
```

を実行すると，point dataでデータを出力したVTK（要素は四角形）ができます。

cutPlaneSurfaceの内容が下記。

``` c++

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/

#includeEtc "caseDicts/postProcessing/surface/surface.cfg"

fields      (T He U);
interpolate true;

surfaces
(
    yz
    {
        type        cuttingPlane;
        interpolate $interpolate;
        planeType   pointAndNormal;
        point       (0 0 0);
        normal      (1 0 0);
    }
);

// ************************************************************************* //
```

これ，計算領域の形状が単純なときはいいんですが，複雑な形状だと，実際にはセルのないところにコンターが書かれて紛らわしいです。なので，以前のように三角形要素のVTKを出力したいとなりました。

# 結論：typeをplaneにする。
typeにplaneがあるので，これを使います。（古いOpenFOAMにもありました。）

planeSurfaceという名前のファイルをsystemに保存します。

```cpp
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/

#includeEtc "caseDicts/postProcessing/surface/surface.cfg"

fields      (T He U);
interpolate true;

surfaces
(
    yz
    {
        type        plane;
        planeType   pointAndNormal;
        point       (0 0 0);
        normal      (1 0 0);
    }
);

// ************************************************************************* //
```

```
postProcess -func planeSurface
```
で実行します。
これで所望の三角要素で区切ったVTKが出てきます。（古いOpenFOAMで出てたVTKと同じ形式のVTK）


# matplotlibでコンターかく
三角形要素のときのVTKつかって，matplotlibでコンターかくコードの例。データ処理の都合上，図を書くのに関係ないライブラリのimportも入ってます。

```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

#plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.family'] = 'Arial'
plt.rcParams["font.size"] = 12
#plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["mathtext.fontset"] = "stix"
#plt.rcParams["mathtext.fontset"] = "dejavusans"
#from matplotlib import rc
#rc('text', usetex=False)

markers1 = [
    "d", "o", "v", "^", "<", ">", "1", "2", "3", "4", "8", "s", "p", "*", "h",
    "H", "x", "D", "|", ","
]
cmap = plt.get_cmap("tab10")
""" functions """


def he(x):
    return x / 4.0 / (x / 4.0 + (1 - x) / 29) * 1e2


def celT(x):
    return x - 273.15

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

    data = reader.GetOutput()
    cells = data.GetPolys()
    triangles = cells.GetData()
    points = data.GetPoints()
    """ mapping cell data to point data """
    mapper = vtk.vtkCellDataToPointData()
    mapper.AddInputData(data)
    mapper.Update()
    mapped_data = mapper.GetOutput()
    Valuedata = mapped_data.GetPointData().GetArray(0)

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

    ntri = triangles.GetNumberOfTuples() // 4
    npts = points.GetNumberOfPoints()
    nvls = Valuedata.GetNumberOfTuples()

    print("ntri, npts, nvls= ", ntri, npts, nvls)
    """ 
    tri: comination of triangle points 
    x, y, z: coordinate
    scalar: value
    """
    tri = np.zeros((ntri, 3))
    x = np.zeros(npts)
    y = np.zeros(npts)
    z = np.zeros(npts)
    scalar = np.zeros(nvls)

    for i in range(0, ntri):
        tri[i, 0] = triangles.GetTuple(4 * i + 1)[0]
        tri[i, 1] = triangles.GetTuple(4 * i + 2)[0]
        tri[i, 2] = triangles.GetTuple(4 * i + 3)[0]

    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, tri, scalar


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

    data = reader.GetOutput()
    cells = data.GetPolys()
    triangles = cells.GetData()
    points = data.GetPoints()
    """ mapping cell data to point data """
    mapper = vtk.vtkCellDataToPointData()
    mapper.AddInputData(data)
    mapper.Update()
    mapped_data = mapper.GetOutput()
    Valuedata = mapped_data.GetPointData().GetArray(0)

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

    ntri = triangles.GetNumberOfTuples() // 4
    npts = points.GetNumberOfPoints()
    nvls = Valuedata.GetNumberOfTuples()

    print("ntri, npts, nvls= ", ntri, npts, nvls)
    """ 
    tri: comination of triangle points 
    x, y, z: coordinate
    v1, v2, v3: vector value
    """
    tri = np.zeros((ntri, 3))
    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(0, ntri):
        tri[i, 0] = triangles.GetTuple(4 * i + 1)[0]
        tri[i, 1] = triangles.GetTuple(4 * i + 2)[0]
        tri[i, 2] = triangles.GetTuple(4 * i + 3)[0]

    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, tri, v1, v2, v3

x, y, z, tri, T = read_scalarAllVTK("T_yz.vtk")

fig_width = 6+0.3
fig_height = fig_width/(x_high-x_low)*(y_high-y_low)*aspect+0.01
figHe = plt.figure(figsize=(fig_width, fig_height))
ax1 = figHe.add_subplot(111)

image1 = ax1.tricontourf(y, z, tri, T, cmap="magma")
plt.savefig("T-yz.jpg", bbox_inches="tight", dpi=100)
plt.close()