matplotlibで非構造格子（2次元, 三角形分割）のデータを可視化する

0. 概要

　この記事では、Pythonの可視化ライブラリであるmatplotlibを使用して2次元で三角形分割された非構造格子上のデータを図示する方法を紹介します。例として単純支持梁の曲げ問題を取り上げると、下の図のように変形後の形状と応力分布を可視化することができます。

1. はじめに

　流体や構造、電磁界などを対象とした科学技術計算には、有限要素法や有限体積法などの数値解析手法があり、これらの方法では空間の離散化に非構造格子が使用されています。計算において要素ごとに求まる物理量を可視化するには、メッシュを物理量の値に応じた色で塗り分ける必要があります。とくに、構造体では変形後の形状に合わせて応力分布を表示したいのではないでしょうか。
　pythonの可視化ライブラリであるmatplotlibのtripcolorを用いると、三角形分割の非構造格子上のデータを図示できます。筆者の調べたところ、tripcolorに関する日本語の記事はほとんど存在しませんでした。そこで、この記事では使い方と実行例を解説します。

前提条件

以下の知識を前提とします。

• Pythonの基本的な使い方
• matplotlibの基本的な使い方

動作環境

• OS: Ubuntu 20.04.4 LTS on Windows 10 X86_64
• Python 3.8.8
• Numpy 1.20.2
• matplotlib 3.3.4

2. tripcolor

2.1 関数の概要

　tripcolorは三角形要素で分割された領域上のデータを色で可視化します。matplotlibの中には以下の二つのインターフェースが用意されています。

• matplotlib.axes.Axes.tripcolor
• matplotlib.pyplot.tripcolor
  上はAxesクラスのメソッド、下はpyplotモジュール内の関数です。筆者は普段Axesのメソッドで使用しているので、デモプログラムではこちらを使用します。tripcolorの引数は以下のように渡します（これ以外の引数の渡し方も存在するので詳細は公式ドキュメントを参照してください）。

tripcolor(x, y, triangles, C, ...)

xとyは節点のx, y座標を表すarray-likeなオブジェクト（numpy.ndarrayなど）であり、(節点数,)の形状を持ちます。trianglesは各三角形を構成する節点の番号で、こちらは(三角形要素の数, 3)の形状を持つint型のarray-likeなオブジェクトです。Cは可視化するデータのarray-likeなオブジェクトであり、その形状は(三角形要素の数,)もしくは(節点数,)である必要があります。また、...の箇所にはその他のキーワード引数を渡すことができます。

2.2 実行例

2.2.1 簡単な例

とりあえず簡単な図を作成するプログラム例を示します。

simple_example.py

# import libraries
import numpy as np
import matplotlib.pyplot as plt

# 5つの節点のx,y座標. 配列要素の順番がそのまま節点番号に対応する.
x = np.array([0.0, 0.0, 0.5, 1.0, 1.0])
y = np.array([0.0, 1.0, 0.5, 0.0, 1.0])
# 4つの三角形要素を構成する節点の番号, e.g., 要素番号: 0 -> 構成する節点: 0, 1, 2
triangles = np.array([[0, 1, 2],
                      [1, 4, 2],
                      [4, 3, 2],
                      [3, 0, 2]])
# 各三角形要素のデータの値
C = np.array([0.0, 1.0, 3.0, 5.0])

fig:plt.Figure = plt.figure(figsize=(4.0,3.0))
ax:plt.Axes = fig.add_subplot()

tpc = ax.tripcolor(x, y, triangles, C)

# カラーバーを表示
fig.colorbar(tpc)
# アスペクト比を1対1に, レイアウトを調整
ax.set_aspect('equal')
fig.tight_layout()

fig.savefig('simple_example.png')

このプログラムでは下の図のように、5つの節点で構成される合計4つの三角形要素のデータを色で表示します。上のコードではx、yには節点の座標を左下(0,0)、左上(0,1)、中央(0.5,0.5)、右下(1,0)、右上(1,1)の順に与え、trianglesには三角形を構成する3つの節点の番号を左、上、右、下の要素の順番に格納しています。最後に、要素毎のデータをCに与え、以上の変数をtripcolor()メソッドに渡しています。tripcolor()はmatplotlib.collections.Collectionオブジェクトを返すので、これをtpcに代入しています。このオブジェクトをFigure.colorbar()に渡すとカラーバーを表示できます。また、x, y軸のスケール比が異なると図が引き伸ばされた形になるので、Axes.set_aspect('equal')でスケール比を揃えています。

2.2.2 節点データの場合

　データが節点毎に与えられている場合もtripcolorで図示できます。2.2.1のコードで、Cに節点と同じ数のデータを与えると自動的に節点毎のデータだと解釈します。

node_data.py

C = np.array([0.0, 1.0, 3.0, 5.0, 6.0])

上記のように変更すると、三角形を構成する3節点の平均値で要素が色付けされます。（もしも節点と要素の数が同一だった場合は節点毎のデータだとみなされます。要素毎のデータであることを指定したいときは、キーワード引数でfacecolors=Cと与えます。）

また、節点データの場合はtripcolorのキーワード引数にshading='gouraud'を渡すとGouraudシェーディングを使用します。この機能を用いると、グラデーションのように描画されます。

gouraud_shading.py

tpc = ax.tripcolor(x, y, triangles, data, shading='gouraud')

2.2.3 その他の機能

(a) カラーマップの範囲

キーワード引数vmin、vmaxでカラーマップの最小値・最大値を指定します。

change_vmax.py

tpc = ax.tripcolor(x, y, triangles, data, vmin=0.0, vmax=10.0)

(b) 枠線の描画

キーワード引数edgecolorsに渡した色で、三角形の枠線を描画します。

edgecolors.py

tpc = ax.tripcolor(x, y, triangles, data, edgecolors='black')

(c) 透明度

キーワード引数alphaで透明度（0.0~1.0）を指定します。また、edgecolors='face'で三角形内部と同色（alphaの影響は受けない）の枠線を描画します。

alpha.py

tpc = ax.tripcolor(x, y, triangles, data, alpha=0.5, edgecolors='face')

(d) カラーマップの変更

キーワード引数cmapにmatplotlib.colors.Colormapオブジェクトを渡して、カラーマップを変更します。デフォルトで用意されているカラーマップは、matplotlib.pyplot.get_cmapで取得できます。カラーマップの名前はこちらを参照してください。https://matplotlib.org/stable/tutorials/colors/colormaps.html

cmap.py

tpc = ax.tripcolor(x, y, triangles, data, cmap=plt.get_cmap('jet'))

3. 単純支持梁の応力分布の可視化

　実用的な例として、弾性論の平面問題で解が得られている単純支持梁の応力分布を可視化を行います。ここでは解析解を用いていますが、もちろんFEMやFVMなどの数値計算結果も可視化できます。こちらは、ソースコードと出力画像だけ置かせていただきます。

ソースコード

visualize_beam_stress.py

#!/usr/bin/env python3
""" Demo script to visualize the unstructured triangular grid data
    using matplotlib.
"""
# import libraries
from collections import namedtuple
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation as Tri

param_names = ['length', 'height', 'NX', 'NY', 'force',
                   'youngs_modulus', 'poissons_ratio']
Params = namedtuple('Params', param_names)

# create mesh
def create_mesh(params:Params) -> tuple:
    """ Create triangular mesh.
    """
    # 節点座標を作成
    l = 0.5*params.length
    c = 0.5*params.height
    x = np.linspace(-l, l, params.NX+1)
    y = np.linspace(-c, c, params.NY+1)
    x, y = np.meshgrid(x, y)

    # 配列形状を二次元から一次元に変更
    x = x.ravel()
    y = y.ravel()

    # triangulationオブジェクトを作成
    # 第三引数 triangles : (ntri, 3) array-like of int, optional
    #  で三角形の節点番号を指定可能
    # NoneならDelaunay三角形分割で自動的に作成
    triangles = Tri(x, y)

    return x, y, triangles

def calculate_exact_solution(x:np.ndarray, y:np.ndarray,
        triangles:Tri, params:Params) -> tuple:
    """ Calculate displacement and stress of 2-dimensional
        linear elasticity problem. """
    l = 0.5 * params.length
    c = 0.5 * params.height
    E = params.youngs_modulus
    nu = params.poissons_ratio
    q = params.force
    I = c**3 / 96.0   # 断面二次モーメント
    
    # 文献ではy座標の向きが逆
    y = -y

    # 中央部分のたわみ
    delta = -5.0*q*l**4/(24.0*E*I)\
             * (1.0 + 12.0/5.0 *(c/l)**2 *(0.8+0.5*nu))
    # x方向変位
    u = q/(2.0*E*I) *( (l**2*x - x**2/3.0)*y\
                       + x*(2.0*y**3/3.0 - 0.4*c**2*y)\
                       + nu*x*(y**3/3.0 - c**2*y + 2.0*c**3/3.0) )
    # y方向変位
    v = -q/(2.0*E*I) *( y**4/12.0 - 0.5*(c*y)**2 + 2.0*c**3*y/3.0\
                        + nu*( 0.5*(l**2-x**2)*y**2 + y**4/6.0\
                               -c**2*y**2/5.0)\
                       - 0.5*(l*x)**2 - x**4/12.0 -(c*x)**2/5.0\
                       +(1.0+0.5*nu)*(c*x)**2 ) + delta

    # 要素中心の座標を計算 (FEMでは三角形要素の中心で応力が求まることに合わせている)
    tris = triangles.triangles
    x_c = x[tris].mean(axis=1)
    y_c = y[tris].mean(axis=1)
    # x方向垂直応力
    stress_x = -0.75*q/c**3 * (x_c**2 * y_c - 2.0/3.0 * y_c**3) \
                + 0.75*q/c*((l/c)**2 - 0.4)*y_c

    return u, v, stress_x

# main function
def main():
    """ Main function.
    """
    # パラメータの設定
    params = Params(
        length = 10.0,
        height = 1.0,
        NX = 300,
        NY = 30,
        force = 1.0,
        youngs_modulus = 1e6,
        poissons_ratio = 0.3
    )

    x, y, tris = create_mesh(params)
    u, v, stress = calculate_exact_solution(x, y, tris, params)

    fig:plt.Figure = plt.figure(figsize=(6.0,1.0))
    ax:plt.Axes = fig.add_subplot()

    # tripcolor(x, y, triangles, C)
    # x, y : (nnode,) array-like
    #   x, y 座標
    # triangles : (ntri, 3) array-like of int
    #   各要素を構成する節点の番号
    # C : (ntri,) or (nnode,) array-like
    #   可視化を行う要素or節点上のデータ
    tpc = ax.tripcolor(x+u, y+v, tris.triangles, stress,
                       cmap=plt.get_cmap('jet'))

    # カラーバーを表示
    fig.colorbar(tpc)

    ax.set_aspect('equal')
    fig.tight_layout()
    # 以下のコメントを外すと枠と目盛りが消える
    # ax.set_frame_on(False)
    # ax.set_xticks([])
    # ax.set_yticks([])
    fig.savefig('beam_stress.png', dpi=600)

# main
if __name__ == "__main__":
    main()

出力画像

4. おわりに

この記事は筆者のQiita初投稿です。誰かの役に立てば嬉しいです。

5. 参考文献

• matplotlib 公式ドキュメント https://matplotlib.org/stable/index.html
• matplotlib.axes.Axes.tripcolor https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tripcolor.html
• matplotlib.axes.Axes.pcolor https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.pcolor.html#matplotlib.axes.Axes.pcolor
• matplotlib.colors.Colormap https://matplotlib.org/stable/api/_as_gen/matplotlib.colors.Colormap.html#matplotlib.colors.Colormap
• S. Timoshenko, and J. N. Goodier, "Theory of Elasticity," McGraw-Hill, (1951).