PyVistaで任意サーフェスメッシュを作成してSTL保存する

工学問題は複雑な形状をメッシュで分割して計算に利用する事が多々あります。この記事では形状やメッシュデータを扱う事が得意なPyVistaを使って任意サーフェスを作成し、デローニー三角形分割によりメッシュを作成します。そしてその他ソフトでの活用を目標にSTLファイルとして保存する所までを紹介します。

こんにちは。wat(@watlablog)です。ここでは任意のサーフェスからメッシュを生成してSTL保存する方法を紹介します！

目次（項目クリックでジャンプできます）

本記事の対象読者

Pythonで形状モデルを作成したい人

当ブログでも色々なプログラムを紹介していますが、Pythonは様々な工学的、科学技術的なライブラリが豊富にあります。

ここで言う形状モデルとは点や線、面や立体のデータを指します。通常これらの形状モデルはCADソフトを使って作成します。

CADソフトも大抵自動化のためのマクロ機能を有していますが、Pythonユーザからすると簡単なモデルであればPythonを使って生成した方が効率的と考えてしまいます。

CADソフトを外から自動化しようとするとGUIで操作していた機能を順番に指定して実行していくと思いますが、中には豊富な機能を持っているがゆえに動作が遅かったりします（もちろん既にCADソフトで標準化されている業務であれば、既存の汎用ソフトを自動化するメリットは大きいですが）。

Pythonで最低限の操作だけ実行するようにコーディングする事ができれば、省エネなコードになるというメリットがあります。

最近は個人で3Dプリンタを購入してDIYをしている人が増えてきました。形状モデルをSTLで作成すれば3Dプリンタ用のデータとなります。是非ご活用下さい。

PyVistaによるメッシングを覚えたい人

Pythonによる形状モデルの操作はPyVistaが大変便利です。

当ブログでは過去に「PyVistaをインストールしてPythonでSTLを扱う備忘録」や「Python/PyVistaでSTLモデルを座標変換してみた」という記事でPyVistaを使ってみましたので、こちらもご参考に。

これらの記事では主にモデルをBlenderで作成してからPyVistaで処理を行なっていました。今回はBlender等のソフトを使わずに形状を作成し、メッシング（メッシュを使って離散化する事）する事を学びます。

とは言え僕は1次情報である以下の開発者ドキュメントを参考にしてコードを書いています。まずはこちらを参照した方が良いでしょう。

外部リンク：PyVista公式ドキュメント

本記事の目標（図で説明）

任意のサーフェスメッシュを作成してSTL保存する…という記事ですが、具体的に何をするかを事前に書いて理解して頂く方が読者に優しい（~~Uターンされない~~）と思うので、ここで図を用いて説明してみます。

点列データを作成する

まずは点列データを3次元座標(x, y, z)で作成します。ここでは「PythonでNACA4桁翼の断面形状を作図する方法」という記事で描いた航空機の翼断面を使って作業をします。

今回は翼形状を奥行き方向に押し出してサーフェス（面）を作成しますが、閉じたサーフェスとするために上面(top)と下面(bottom)の2面を使う仕様とします。

点列データからスプラインを作成する

先ほど作成した点列データを使って滑らかなスプラインを作成します。…実はスプラインを使わなくても今回の目的は達成できそうだったのですが、せっかく覚えたのでPyVistaの練習としてコード内に残しました。

スプラインは点列の間隔を細かくする事によって詳細にカーブを表現する事が出来ます。最終的な面の仕上がりにも影響するので、今回のPythonコード内では細かさを調整できるようにしています。

外部リンク：PyVista：スプラインの作成

スプラインを押し出してサーフェス用の点列データを作成する

点列（ここでは閉じたサーフェスとしたいので点列で作成したプロファイル）をある方向に押し出して(Extrude)サーフェスを作る基礎とします。

後で三角形メッシュを使って面を分割するのですが、あまりにも押し出し長が長すぎると三角形1辺の比率が大きくなって品質の悪いメッシュになってしまいます。

そのため今回はDivisionsという変数を使って最大長を任意の分割数に分けるように作りました。

点列データに対してデローニー三角形分割を行う（メッシュの作成）

サーフェスの基礎となる点列データが作成できたら、次は面を貼ります。この後サーフェスデータはSTLファイル(三角形の面データで表現した3Dモデル)とするため、面を三角形で分割します。

三角形の分割アルゴリズムにはデローニー三角形分割（Delauney triangulation）が有名であり、今回はこの手法を使います（ドロネー分割とも呼ばれるけどどっちが一般的だろう？？）。

外部リンク：PyVista：三角形サーフェスを作成

それでは早速Pythonでコーディングしてみましょう！

Pythonで任意サーフェスメッシュを作成してSTL保存するコード

動作環境

本記事のプログラムは以下の環境で動作を確認しています。

Windows	OS	Windows10 64bit
	CPU	2.4[GHz]
	メモリ	4[GB]

Mac	OS	macOS Catalina 10.15.7
	CPU	1.4[GHz]
	メモリ	8[GB]

Python	Python 3.9.6
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.21.1
pyvista	0.31.3

全コード（コピペ用）

以下が全コードです。「def NACA()」関数は「PythonでNACA4桁翼の断面形状を作図する方法」で断面を作図した関数ですが、一部変更しています。

デローニー三角形分割はpyvista.delaunay_2dで行いますが、この関数はどうやらXY方向に拡がっている面に対して分割が実行されるようで、全体的に座標変換をしています。

サーフェスはNACA関数の「x_resolution」を大きくする事で細かくなります。

「def generate_surf_from_spline()」関数のlengthに押し出し長、divisionsに分割数を渡します。round_numはデータの四捨五入を調整します。

import numpy as np
import pyvista as pv

# 最大反り量m, 最大反り位置p, 最大厚み比tからNACA4桁翼形状を作成する関数
def NACA4(m, p, t, x_resolution):
    # 無次元横軸(0<=x<=1)を作成
    x = np.linspace(0, 1, x_resolution)

    # 公称パラメータを実際の桁に合わせる
    t = t / 100
    m = m / 100
    p = p / 10
    print('t=', t, 'm=', m, 'p=', p)

    # 翼厚分布
    yt = (t / 0.20) * \
         (0.29690 * np.sqrt(x)
          - 0.12600 * x
          - 0.35160 * (x ** 2)
          + 0.28430 * (x ** 3)
          - 0.10150 * (x ** 4))

    # キャンバー線(反り)
    yc_front = (m / (p ** 2)) * ((2 * p * x) - (x ** 2))
    yc_back = (m / (1 - p) ** 2) * ((1 - 2 * p) + (2 * p * x) - (x ** 2))

    # 翼厚分布とキャンバー線から翼形状を作成
    yc = []
    for i in range(len(x)):
        if x[i] < p:
            yc.append(yc_front[i])
        else:
            yc.append(yc_back[i])

    # 翼厚分布と反りを形状に反映
    y_top = yt + yc
    y_bottom = - yt + yc

    # 後でサーフェスを閉じるために末尾を0にする
    y_top[-1] = 0
    y_bottom[-1] = 0

    # Z座標を設定
    z = np.zeros_like(x)

    # 点列をまとめる
    points_top = np.column_stack((x, z, y_top))
    points_bottom = np.column_stack((x, z, y_bottom))

    # 点列データからスプラインを作成
    spline_top = pv.Spline(points_top, x_resolution)
    spline_bottom = pv.Spline(points_bottom, x_resolution)

    ''' 表示用
    cloud1 = pv.PolyData(spline_top.points)
    cloud2 = pv.PolyData(spline_bottom.points)
    add = spline_top + spline_bottom + cloud1 + cloud2
    add.plot(point_size=15, color='tan')
    '''

    return spline_top, spline_bottom

def generate_surf_from_spline(spline_top, spline_bottom, length, divisions, round_num):
    # splineの点データをコピーする（コピーしないと元データが書き換わる）
    base_top = np.copy(spline_top.points)
    base_bottom = np.copy(spline_bottom.points)

    # 押し出し方向の初期値を準備
    dy = 0

    # divisions数の間隔で点列データを押し出しコピー
    for i in range(int(divisions)):
        dy += length / round(divisions, round_num)

        # y方向の座標に増分値を代入
        base_top[:, -2] = np.round(dy, round_num)
        base_bottom[:, -2] = np.round(dy, round_num)

        # splineの点列データに行方向に追加
        spline_top.points = np.vstack([spline_top.points, base_top])
        spline_bottom.points = np.vstack([spline_bottom.points, base_bottom])

    # 点列データからXY方向にデローニー三角形分割を行う
    points_top = pv.PolyData(spline_top.points)
    points_bottom = pv.PolyData(spline_bottom.points)
    surf_top = points_top.delaunay_2d()
    surf_bottom = points_bottom.delaunay_2d()

    '''表示用
    add_points = points_top + points_bottom
    add_points.plot()
    '''

    # ブーリアン加算演算をして面を一つにする
    surf_add = surf_top + surf_bottom
    surf_add = surf_add.clean(tolerance=1e-2)   # 面同士が結合している点を重複削除する
    surf_add.plot(show_edges=True)

    # 点の数を確認（重複がなくなっていればtop≠bottomとなる）
    print('top=', len(surf_top.points))
    print('bottom=', len(surf_bottom.points))
    print('add=', len(surf_add.points))

    return surf_add

# NACA4桁翼の断面3D点列座標を計算する関数を実行
spline_top, spline_bottom = NACA4(m=3, p=3, t=15, x_resolution=30)

# スプラインから断面一定のNACA翼サーフェスメッシュを作成
surf = generate_surf_from_spline(spline_top, spline_bottom, 1, 10, 1)

# STLファイルを作成
pv.save_meshio('surf.stl', surf)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

import numpy as np

import pyvista as pv

# 最大反り量m, 最大反り位置p, 最大厚み比tからNACA4桁翼形状を作成する関数

def NACA4(m, p, t, x_resolution):

# 無次元横軸(0<=x<=1)を作成

x = np.linspace(0, 1, x_resolution)

# 公称パラメータを実際の桁に合わせる

t = t / 100

m = m / 100

p = p / 10

print('t=', t, 'm=', m, 'p=', p)

# 翼厚分布

yt = (t / 0.20) * \

(0.29690 * np.sqrt(x)

- 0.12600 * x

- 0.35160 * (x ** 2)

+ 0.28430 * (x ** 3)

- 0.10150 * (x ** 4))

# キャンバー線(反り)

yc_front = (m / (p ** 2)) * ((2 * p * x) - (x ** 2))

yc_back = (m / (1 - p) ** 2) * ((1 - 2 * p) + (2 * p * x) - (x ** 2))

# 翼厚分布とキャンバー線から翼形状を作成

yc = []

for i in range(len(x)):

if x[i] < p:

yc.append(yc_front[i])

else:

yc.append(yc_back[i])

# 翼厚分布と反りを形状に反映

y_top = yt + yc

y_bottom = - yt + yc

# 後でサーフェスを閉じるために末尾を0にする

y_top[-1] = 0

y_bottom[-1] = 0

# Z座標を設定

z = np.zeros_like(x)

# 点列をまとめる

points_top = np.column_stack((x, z, y_top))

points_bottom = np.column_stack((x, z, y_bottom))

# 点列データからスプラインを作成

spline_top = pv.Spline(points_top, x_resolution)

spline_bottom = pv.Spline(points_bottom, x_resolution)

''' 表示用

cloud1 = pv.PolyData(spline_top.points)

cloud2 = pv.PolyData(spline_bottom.points)

add = spline_top + spline_bottom + cloud1 + cloud2

add.plot(point_size=15, color='tan')

'''

return spline_top, spline_bottom

def generate_surf_from_spline(spline_top, spline_bottom, length, divisions, round_num):

# splineの点データをコピーする（コピーしないと元データが書き換わる）

base_top = np.copy(spline_top.points)

base_bottom = np.copy(spline_bottom.points)

# 押し出し方向の初期値を準備

dy = 0

# divisions数の間隔で点列データを押し出しコピー

for i in range(int(divisions)):

dy += length / round(divisions, round_num)

# y方向の座標に増分値を代入

base_top[:, -2] = np.round(dy, round_num)

base_bottom[:, -2] = np.round(dy, round_num)

# splineの点列データに行方向に追加

spline_top.points = np.vstack([spline_top.points, base_top])

spline_bottom.points = np.vstack([spline_bottom.points, base_bottom])

# 点列データからXY方向にデローニー三角形分割を行う

points_top = pv.PolyData(spline_top.points)

points_bottom = pv.PolyData(spline_bottom.points)

surf_top = points_top.delaunay_2d()

surf_bottom = points_bottom.delaunay_2d()

'''表示用

add_points = points_top + points_bottom

add_points.plot()

'''

# ブーリアン加算演算をして面を一つにする

surf_add = surf_top + surf_bottom

surf_add = surf_add.clean(tolerance=1e-2) # 面同士が結合している点を重複削除する

surf_add.plot(show_edges=True)

# 点の数を確認（重複がなくなっていればtop≠bottomとなる）

print('top=', len(surf_top.points))

print('bottom=', len(surf_bottom.points))

print('add=', len(surf_add.points))

return surf_add

# NACA4桁翼の断面3D点列座標を計算する関数を実行

spline_top, spline_bottom = NACA4(m=3, p=3, t=15, x_resolution=30)

# スプラインから断面一定のNACA翼サーフェスメッシュを作成

surf = generate_surf_from_spline(spline_top, spline_bottom, 1, 10, 1)

# STLファイルを作成

pv.save_meshio('surf.stl', surf)

結果

実行するとグラフウィンドウが表示されますが、ウィンドウを閉じると「surf.stl」というSTLファイルがプログラム実行フォルダ直下に作られます。

PyVistaでSTLファイルを読み込んで確認する方法

作成したSTLファイルがちゃんとしているかどうか、確認してみましょう。

確認もPyVistaで可能です。以下のコードでSTLファイルを読み込んでメッシュの状態を確認してみましょう。

import pyvista as pv

mesh = pv.read('surf.stl')
mesh.plot(show_edges=True)

import pyvista as pv

mesh = pv.read('surf.stl')

mesh.plot(show_edges=True)

下図のように作成した形状が読み込まれればOKです。

今回はサーフェスデータをSTLに保存しましたが、次回記事は面を全て閉じたSTLにしようと思います。

おまけ：面を作成してSTLにするPythonコード(簡易版)

先ほどまではちょっと凝ったサーフェスデータを作っていましたが、ここでは簡単に平面を作ってSTLファイルに保存するサンプルコードを紹介します。

import pyvista as pv
import numpy as np

def generate_surface(filename):
    ''' サーフェスを作成してSTLファイルに保存する '''
    # サーフェスデータを作成
    x = np.arange(0, 10, 0.2)
    y = np.arange(0, 10, 0.2)
    X, Y = np.meshgrid(x, y)
    #Z = np.sin(X) * np.cos(Y)
    Z = 0.0 * X + 0.0 * Y

    # メッシュデータに変換
    points = np.column_stack((X.ravel(), Y.ravel(), Z.ravel()))
    points = pv.PolyData(points)
    mesh = points.delaunay_2d()

    # 表示させて確認
    mesh.plot(show_edges=True)

    # STLファイルに保存
    pv.save_meshio(filename, mesh)

    return

if __name__ == '__main__':
    generate_surface('mesh.stl')

import pyvista as pv

import numpy as np

def generate_surface(filename):

''' サーフェスを作成してSTLファイルに保存する '''

# サーフェスデータを作成

x = np.arange(0, 10, 0.2)

y = np.arange(0, 10, 0.2)

X, Y = np.meshgrid(x, y)

#Z = np.sin(X) * np.cos(Y)

Z = 0.0 * X + 0.0 * Y

# メッシュデータに変換

points = np.column_stack((X.ravel(), Y.ravel(), Z.ravel()))

points = pv.PolyData(points)

mesh = points.delaunay_2d()

# 表示させて確認

mesh.plot(show_edges=True)

# STLファイルに保存

pv.save_meshio(filename, mesh)

return

if __name__ == '__main__':

generate_surface('mesh.stl')

まとめ

本記事ではPyVistaというライブラリを使って点列に対してスプラインを作ったり、三角形メッシュを作ったり、STLファイルにしたりしました。

STLファイルを作る所は自前でやろうとするとテキストファイルで色々細かい指定をしなければならず、かなり骨が折れる作業と思います（理解のために一度やってみたいですが）。

普段はCA○IAとか使っているので、最初は「CUIでモデルを作成する」という事に少し抵抗を感じていましたが、このくらいの規模であればやってみるとそんなにストレスはありませんでした。

これもライブラリのおかげと思います。

数式から作成した形状を実際にSTLモデルにまでする事ができました！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

本記事の対象読者

Pythonで形状モデルを作成したい人

PyVistaによるメッシングを覚えたい人

本記事の目標（図で説明）

点列データを作成する

点列データからスプラインを作成する

スプラインを押し出してサーフェス用の点列データを作成する

点列データに対してデローニー三角形分割を行う（メッシュの作成）

Pythonで任意サーフェスメッシュを作成してSTL保存するコード

動作環境

全コード（コピペ用）

結果

PyVistaでSTLファイルを読み込んで確認する方法

おまけ：面を作成してSTLにするPythonコード(簡易版)

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

本記事の対象読者

Pythonで形状モデルを作成したい人

PyVistaによるメッシングを覚えたい人

本記事の目標（図で説明）

点列データを作成する

点列データからスプラインを作成する

スプラインを押し出してサーフェス用の点列データを作成する

点列データに対してデローニー三角形分割を行う（メッシュの作成）

Pythonで任意サーフェスメッシュを作成してSTL保存するコード

動作環境

全コード（コピペ用）

結果

PyVistaでSTLファイルを読み込んで確認する方法

おまけ：面を作成してSTLにするPythonコード(簡易版)

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル