OpenFOAMのCavityチュートリアルのメッシュを台形にする

PyGEM

• Python Geometrical Morphing.
• 1400万格子でも大丈夫！らしい
• MITライセンス
• ちょっとしたGUIもある

pip show PyGEM

Name: pygem
Version: 0.2
Summary: Tools to apply FFD.
Home-page: https://github.com/mathLab/PyGeM
Author: Filippo Salmoiraghi, Marco Tezzele
Author-email: [e-mail], [e-mail]
License: MIT
Location: /home/ubuntu/.pyenv/versions/anaconda3-4.3.1/envs/py35/lib/python3.5/site-packages
Requires: numpy, numpy-stl, scipy, matplotlib, enum34, Sphinx, sphinx-rtd-theme

Python3系利用する場合は2to3を使う

---

メッシュモーフィングの準備

OpenFOAMのCavityケースでメッシュモーフィングを行う
blockMeshコマンドを実行後、boundingBoxを見る

blockMesh|  grep boundingBox

boundingBox: (0 0 0) (0.1 0.1 0.01)

格子点データをコピーしておく。

import shutil
from PyFoam.IPythonHelpers.Case import Case

case = Case("../cavity")
org_points_path = case.sol.polyMeshDir()+"/points_org"
points_path = case.sol.polyMeshDir()+"/points"
shutil.move(points_path, org_points_path)

---

boundingBoxを定義する

格子点が完全にコントロールボックス内に入るように少し大きめのコントロールボックスにする。

import numpy as np
import pygem as pg
import pygem.affine as at

params = pg.params.FFDParameters()

e = 1.0e-6

params.rot_angle = np.array([0,0,0])
params.lenght_box =np.array([0.1,0.1,0.01]) + e
params.origin_box = np.array([0.,0.,0.]) -e

params.rotation_matrix = at.angles2matrix(
    params.rot_angle[2] * np.pi / 180,
    params.rot_angle[1] * np.pi / 180,
    params.rot_angle[0] * np.pi / 180)

params.position_vertex_0 = params.origin_box
params.position_vertex_1 = params.position_vertex_0 + np.dot(params.rotation_matrix, [params.lenght_box[0], 0, 0])
params.position_vertex_2 = params.position_vertex_0 + np.dot(params.rotation_matrix, [0, params.lenght_box[1], 0])
params.position_vertex_3 = params.position_vertex_0 + np.dot(params.rotation_matrix, [0, 0, params.lenght_box[2]])

params.psi_mapping  = np.diag(1. / params.lenght_box)
params.inv_psi_mapping = np.diag(params.lenght_box)

---

台形に変形するように設定する

cavityメッシュの正方形(xy平面)を台形にする
コントロールポイントと変形量をparamsに追加する

params.n_control_points=np.array([2,2,2])
params.array_mu_x = np.zeros(params.n_control_points.cumprod()[-1]).reshape(params.n_control_points)
params.array_mu_y = np.zeros(params.n_control_points.cumprod()[-1]).reshape(params.n_control_points)
params.array_mu_z = np.zeros(params.n_control_points.cumprod()[-1]).reshape(params.n_control_points)

params.array_mu_x[0,1,:]  = 0.25
params.array_mu_x[1,1,:]  = -0.25
params.array_mu_x

array([[[ 0.  ,  0.  ],
        [ 0.25,  0.25]],

       [[ 0.  ,  0.  ],
        [-0.25, -0.25]]])

---

台形にする

格子点を書きこむ（元のファイルを上書きするような破壊的な使い方はできない）

of_handler = pg.openfhandler.OpenFoamHandler()
points = of_handler.parse(org_points_path)

free_form = pg.freeform.FFD(params, points)
free_form.perform()
new_points = free_form.modified_mesh_points

of_handler.write(new_points, points_path)

---

台形になっていることをParaViewで確認する

台形になった