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

# 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])
```

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

[^1]: https://github.com/nschloe/pygmsh/pull/579/files

## 仮想環境の作成

powershell で作業しました。

```pwsh
> 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` が大きめの（軽い）メッシュは作ります。

```python
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/98be7fe2-9667-caa5-1cdb-89c9c274b114.png)

## オブジェクトの構造

```python
>>> 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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/16b9a000-4feb-7d46-2e09-b94442883d18.png)

## 立体図形どうしの演算

<details><summary>図形同士の差などが計算できます。</summary>

```python
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)
```
</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/38a5733a-248d-27fa-7def-ce63e4c2530d.png)