OpenFOAM9くらいを想定してます。
基本的なやりかたは上の記事で書いてます。
以前，
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++ -*----------------------------------*\
========= |
\\ / 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に保存します。
/*--------------------------------*- 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も入ってます。
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()