LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@natsuriver

gmsh を pygmsh インターフェースから利用する

Posted at

pygmsh と gmsh

gmsh https://gmsh.info/ - A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities を pygmsh インターフェースから呼び出すこともできます。gmsh 本体の Python インターフェイスと比べて、こなれたインターフェイスで、記述量が少なめな Pythonic なラッパーです。使用例を紹介します。
gmsh を使うモチベーションは、Open Cascade のエンジンを使って、ボリュームのプリミティブどうしの 和や差の演算を行いたいからです。

テスト環境

（バージョンに対してシビアではなさそうです。ざっくりした感触ですが Python 3.9 ぐらいからは問題なさそうです）

  • Windows11
  • Python 3.12.7
  • venv を作り、pip で必要なライブラリーだけを入れています。（検証はしていませんが、他の手法でも大丈夫だと思います。ただし、ライブラリー内のコード修正が必要です。）

pygmsh には、パッチ 1が必要。

  • pygmsh を手元に clone します。 https://github.com/nschloe/pygmsh
  • src/pygmsh/common/surface.py を修正 1 します。
  • VENV 環境を作る前に、この修正が必要となります。
src/pygmsh/common/surface.py
        try:
            self._id = env.addSurfaceFilling([self.curve_loop._id])
        except TypeError:       # occ backend expects single tag
            self._id = env.addSurfaceFilling(self.curve_loop._id)

(参考: 修正前の内容)

src/pygmsh/common/surface(old).py
        self._id = env.addSurfaceFilling([self.curve_loop._id])

このプルリクが反映されない、という意味ではメンテナンスされていないライブラリーという一面もあります。（でも使いやすくて良いライブラリーです）

仮想環境の作成

powershell で作業しました。

> cd ~

> ~\appData\Local\Programs\Python\Python312\python.exe -m venv .\py312_pygmsh --upgrade-deps

> ~\py312_pygmsh\Scripts\Activate.ps1

> cd ~\pygmsh

> pip install .  
# ピリオドを使って、カレントディレクトリのパッケージをインストールします。
# 他の依存ライブラリーも合わせてインストールされます

> pip install pyvista ipython

開発/動作確認のための補助ツールとして、pyvista, ipython を導入しました。必須ではありませんが、pyvista で 3D ビューを表示して、結果を素早く確認できます。

平面ポリゴンを z 方向に押し出したボリューム

ここでは、xy 面上のポリゴンを z 方向に押し出した形状を例に pygmsh を使って、ボリュームを組み立てます。三角格子のメッシュを切って、結果を pyvista で可視化して確認します。
単に押し出すだけなら、geom.extrude を使うだけでも済みます。後で、別の形状も作りたいので、
自分で全部組み立ててみましょう。

メッシュは要らない場合があるかもしれません。可視化のための便宜として、mesh_size が大きめの（軽い）メッシュは作ります。

import numpy as np
import pygmsh

from pygmsh.common.point import Point
from pygmsh.common.line import Line
from pygmsh.common.curve_loop import CurveLoop
from pygmsh.common.surface import Surface
from pygmsh.common.surface_loop import SurfaceLoop
from pygmsh.common.volume import Volume

import pyvista as pv # 結果可視化のため

geometry = pygmsh.occ.Geometry()
geom = geometry.__enter__()  #  これをしないと、初期化されていないというエラーが帰る

a, b = 0.2, 1.5
pts = np.array([(-a, a), (b, a), (b, b), (-a, b)])

z = np.ones_like(pts)[:, 0:1]

base_quadrangle = np.hstack([pts, z * 0.])
top_quadrangle = np.hstack([pts, z * 1.2])

lc = 2.0  # mesh size

pbot = [geom.add_point(botp, mesh_size=lc) for botp in base_quadrangle]
ptop = [geom.add_point(topp, mesh_size=lc) for topp in top_quadrangle]

# line を定義するときに、最初と最後をつなぐため、ひとつずらした配列を作っておく
pbot_next = np.roll(pbot, -1)  # pbot[1], pbot[2], ..., pbot[-1], pbot[0]
ptop_next = np.roll(ptop, -1)

ebot = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot, pbot_next)]
emid = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot, ptop)]
etop = [geom.add_line(phead, ptail) for phead, ptail in zip(ptop, ptop_next)]
emid_next = np.roll(emid, -1)

lbot = [geom.add_curve_loop([-e for e in ebot[::-1]])]  # 逆回り、符号反転
# 底面ループは、わざわざ逆向きにしなくても、処理できるかもしれない。
# 逆回りにするなら、マイナスも与えて、エッジの向きも反転させる

lmid = [geom.add_curve_loop([eb, emn, -et, -em]) for eb, emn, et, em in zip(ebot, emid_next, etop, emid)]
ltop = [geom.add_curve_loop(etop)]

# Create a surface for each line loop.
s = [geom.add_surface(lp) for lp in lbot + lmid + ltop]
# 表面メッシュだけで良いなら、ここですぐ `mesh = geom.generate_mesh()` でも通ります

# Create the surface loop.
surface_loop = geom.add_surface_loop(s)


vol = geom.add_volume(surface_loop)

mesh = geom.generate_mesh()

geo_filename = "./col01.vtu"


mesh.write(geo_filename)

polydata = pv.read(geo_filename)
pv.plot(polydata, show_edges=True)

ポップアップウィンドウ内に、以下のような構造が表示されます。

quadangular_extruded.png

オブジェクトの構造

>>> lmid[0]._id
 2

>>> lmid[0].dim_tag
 (1, 2)

>>> lmid[0].dim_tags
 [(1, 2)]

>>> s[0]._id
 1

>>> s[0].dim_tags
 [(2, 1)]

ユーザは、変数名でオブジェクトを管理しますが、内部的には、_id や、これをオブジェクトの次元と組み合わせた dim_tag で管理しているようです。これは Gmsh の仕様で、pygmsh はこの内部情報を透過させているのでしょう。dim_tags の方では配列にしてありますが、これが必要になる状況には、遭遇していません。

任意の多角形

上に示したコードで、配列 pts を書き換えて、所望の多角形の形状にすれば、 押し出した 3D が生成できます。

a, b = 0.2, 1.5
pts = np.array([(a, -a), (a, -b), (b, -b), (b, b), (a, b), (a, a), (-a, a), (-a, b), (-b, b), (-b, -b), (-a, -b), (-a, -a)])

image.png

立体図形どうしの演算

図形同士の差などが計算できます。 
import numpy as np
import pygmsh

import pyvista as pv # 結果可視化のため

geometry = pygmsh.occ.Geometry()
geom = geometry.__enter__()  #  これをしないと、初期化されていないというエラーが返る

a, b = 0.2, 1.1
pts_0 = np.array([(-a, a), (b, a), (b, b), (-a, b)])

z = np.ones_like(pts_0)[:, 0:1]

base_quadrangle = np.hstack([pts_0, z * 0.])
top_quadrangle = np.hstack([pts_0, z * 1.2])

lc = 2.0  # mesh size

pbot = [geom.add_point(botp, mesh_size=lc) for botp in base_quadrangle]
ptop = [geom.add_point(topp, mesh_size=lc) for topp in top_quadrangle]

pbot_next = np.roll(pbot, -1)  # pbot[1], pbot[2], ..., pbot[-1], pbot[0]
ptop_next = np.roll(ptop, -1)

ebot = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot, pbot_next)]
emid = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot, ptop)]
etop = [geom.add_line(phead, ptail) for phead, ptail in zip(ptop, ptop_next)]
emid_next = np.roll(emid, -1)
lbot = [geom.add_curve_loop([e for e in ebot])]  # NG 例
# lbot = [geom.add_curve_loop([-e for e in ebot[::-1]])]  # 逆回り、符号反転
lmid = [geom.add_curve_loop([eb, emn, -et, -em]) for eb, emn, et, em in zip(ebot, emid_next, etop, emid)]
ltop = [geom.add_curve_loop(etop)]

# Create a surface for each line loop.
s = [geom.add_surface(lp) for lp in lbot + lmid + ltop]

# Create the surface loop.
surface_loop = geom.add_surface_loop(s)

vol = geom.add_volume(surface_loop)

a, b = 0.05, 0.6
pts_1 = np.array([(a, -a), (a, -b), (b, -b), (b, b), (a, b), (a, a), (-a, a), (-a, b), (-b, b), (-b, -b), (-a, -b), (-a, -a)])
z_1 = np.ones_like(pts_1)[:, 0:1]

base_quadrangle_1 = np.hstack([pts_1, z_1 * 0.])
top_quadrangle_1 = np.hstack([pts_1, z_1 * 1.5])

pbot_1 = [geom.add_point(botp, mesh_size=lc) for botp in base_quadrangle_1]
ptop_1 = [geom.add_point(topp, mesh_size=lc) for topp in top_quadrangle_1]

pbot_next_1 = np.roll(pbot_1, -1)  # pbot[1], pbot[2], ..., pbot[-1], pbot[0]
ptop_next_1 = np.roll(ptop_1, -1)

ebot_1 = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot_1, pbot_next_1)]
emid_1 = [geom.add_line(phead, ptail) for phead, ptail in zip(pbot_1, ptop_1)]
etop_1 = [geom.add_line(phead, ptail) for phead, ptail in zip(ptop_1, ptop_next_1)]
emid_next_1 = np.roll(emid_1, -1)
lbot_1 = [geom.add_curve_loop([e for e in ebot_1])]  # NG 例
# lbot = [geom.add_curve_loop([-e for e in ebot[::-1]])]  # 逆回り、符号反転
lmid_1 = [geom.add_curve_loop([eb, emn, -et, -em]) for eb, emn, et, em in zip(ebot_1, emid_next_1, etop_1, emid_1)]
ltop_1 = [geom.add_curve_loop(etop_1)]

# Create a surface for each line loop.
s_1 = [geom.add_surface(lp) for lp in lbot_1 + lmid_1 + ltop_1]

# Create the surface loop.
surface_loop_1 = geom.add_surface_loop(s_1)

vol_1 = geom.add_volume(surface_loop_1)

# vd = geom.boolean_union([vol, vol_1])
vd = geom.boolean_difference(vol_1, vol)

mesh = geom.generate_mesh()

geo_filename = "./col01.vtu"

mesh.write(geo_filename)

polydata = pv.read(geo_filename)
pv.plot(polydata, show_edges=True)

image.png

  1. https://github.com/nschloe/pygmsh/pull/579/files ↩ ↩2

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?