gmsh を使って .stl ファイルのリメッシュを行う手順

## gmsh を使って .stl ファイルのリメッシュを行う手順

### 背景

gmsh: https://gmsh.info/ を使えば .stl ファイルをリメッシュできる，という記載 [^1], [^2] はあるものの，具体的な手順の説明がなくて，あれこれ試行錯誤しました．分かってみれば簡単なことだったのでまとめます．
ひとことで言えば，チュートリアルの `t13.geo` ファイルを読んだ後で、キーボードの `2` を押すとリメッシュされた結果が得られます．

### 環境構築

gmsh: https://gmsh.info/ 4.8.4 をインストールしました

* OS: Windows10 Home 20H2

* gmsh は、Python の venv 仮想環境の中に pip でインストールしたものを利用
  仮想環境は venv を用いて `C:\Users\water\py39\gmsh` 以下に作成

  ```bash
  C:\Users\mizut\py39>python -m venv gmsh
  C:\Users\mizut\py39>gmsh\Scripts\activate.bat
  (gmsh) C:\Users\mizut\py39>python -m pip install --upgrade gmsh
    ...
  (gmsh) C:\Users\mizut\py39>gmsh\Scripts\gmsh.bat 
  ```

  

  * （趣味の問題） `.venv` を使わないで `C:\Users\water\sample1> python -m venv gmsh` としている
  * gmsh のコマンドライン起動は，この仮想環境内にインストールされているバッチファイル`C:\Users\water\py39\meshtools\Scripts\gmsh.bat` で起動します．仮想環境内に入っている必要があります．（バッチファイルから python を起動して gmsh を起動させるからです）

### 検討の方針

チュートリアルフォルダ `C:\Users\water\py39\gmsh\Lib\site-packages\gmsh-4.8.4-Windows64-sdk\share\doc\gmsh\tutorial` から例題 `t13.geo` と `t13_data.stl` ファイルを手元にコピーして，これを修正しながら，`gmsh` の操作を理解します． 


### ステップ1: チュートリアルの起動と操作

```bash
C:\Users\mizut\py39>gmsh\Scripts\activate.bat
(gmsh) C:\Users\mizut\py39>gmsh\Scripts\gmsh.bat 
```

#### 1.1 チュートリアルの `t13.geo` 読み込み

起動した gmsh の GUI から チュートリアルの `t13.geo` を読み込みます．
[File] - [Open...] から .geo を読み込むことができます．


![image-001.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/4af782fd-d93a-d65c-fe95-976bcc387ed0.png)


モデルはマウス操作でぐるぐる回せます．

* ***左ボタン*** : 回転
* **中央ボタン** : 拡大縮小
* **右ボタン** : 視点中心の移動

  ちなみに `t13_data_stl` は，こんな感じで箱に丸穴が 2 つと四角の穴が開いた物体です。

  ![stlモデル](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/013f9bfe-4da0-3ccd-d913-b33e3f5284d3.png)

   次の項目に書いたように 2D element faces を表示させた上でモデルを回転させれば、gmsh の GUI 上で構造を確認できます．

#### 1.2 gmsh で扱うデータの構造を眺める

geo ファイルの記述内容によって， stl ファイルの表面メッシュは複数の surface に分割されています．この様子を見るために，2 つの補助ウインドウを表示します．

[Tools] - [Options]

![Options](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/5e191ca1-92dd-d223-c7eb-270cd3d09db9.png)

[Tools] - [Visibility] 

![Visibility](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/545da182-5a2b-733e-ce42-5ce075a2761f.png)


Visibilty ウィンドウで List タブを選んだ状態で、どれか 1 つ Surface を選択して右下の Apply ボタンを押します．また矢印キー上下で選択内容を切り替えると、順次表示内容が切り替わるので、Surface 2 から Surface 12 までの内容が確認できます．Visibility のリスト内で複数の要素を選んでから このウィンドウ右下の Apply ボタンを押すことで，選択した Element だけを表示させて確認できます．

![element 表示](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/720611cf-b6d1-48b4-3fa1-21a6d7ee7a27.png)


このとき表示されているのは、Mesh の中の 2D element だということが、以下の操作で確かめられます．

Option ウィンドウ － Mesh を選んで、2D element faces にチェック

![2D element faces](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/8a4901e8-158d-dc34-fac2-0348439de0e5.png)

#### 1.3 リメッシュの実行

メインウィンドウの Mesh の下にある 2D をクリックすると，あらかじめセットされている条件にしたがって、リメッシュが実行されます．（キーボードショートカット は `2` ）

![remesh](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/2fd95efa-ca26-7981-0dbb-146e7dd0f012.png)


#### 1.4 結果の保存

[File] [Save mesh] で  `t13.msh` ファイルにこのメッシュが保存できます．また [Export] で名前をつけて，stl ファイルとして保存できます．

#### 1.5 3D メッシュ

メインウィンドウの Mesh の下にある 3D をクリックすると，表面メッシュをもとに 3D メッシュが作成できます．（キーボードショートカット は `3` ） 

[Export] で .vtk ファイルに書き出すと paraview を使って 3D メッシュの状態を確認できます．

* paraview の Filter Clip で切るときに Crinkle Clip のチェックボックスの状態に応じて、メッシュを維持したり，メッシュをさらに slice したりする．ON/OFF どちらかの設定で、メッシュ状態が確認できるはずです． 

### ステップ2: パラメータをあれこれ調整する

#### 2.1 エッジ認識の角度パラメータ

stl 形状からエッジを認識する条件の角度設定は、GUI の Parameters の項目に表示されている `40` という数字で調整できます．この効果を確認するために，まず geometry ファイルをリロードします． Geometry の下の Reload script をクリックします．（キーボードショートカット は `0` ） 

Parameters の箱の中の数字をキーボード入力かマウス操作で変更しましょう．数字を小さくすると， Surface の分割数が増えます．どういう面ができたかを，1.2 で説明したように Visibility を切り替えて確認することができます．

30 度にすると，小さな三角形が，別の Surface として分離されてしまう箇所がありました．青い三角形の部分は、外側の面（surface 2）に穴が開いていて、内側の表面（surface 9）が見えています．

![image-008](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/86b66abf-8f41-c0b9-4573-a7b105b17a7a.png)


このときも， Mesh 2D を実行すれば，リメッシュは実行されますがこの穴 （surface 3）近傍の三角形はいびつになっています．

![image-009](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/f428da46-ae83-1ebd-eecc-e68622c22a5b.png)


また角度パラメータを 60 にすると、上側の四角い開口のエッジ面の「ロ」型が独立な面として認識されなくなります．リメッシュ後の形状においても、開口エッジ面の一部は斜めに傾いて側面とつながります．

![image-20210807214145808](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/3be0c332-4189-1bfd-d2b6-cfe0f7586852.png)


#### 2.2 フィールドメッシュサイズ

##### 2.2.a Parameters の調整例

Parameters の □ Apply funny mesh size をチェックすると、メッシュサイズを 座標 x と y の関数として調整したメッシュが作成されます．

![image-010.png]
(https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/283cd272-5a87-01e3-2ae6-2d517fbf2c2c.png)

これはチュートリアルファイル `t13.geo` の末尾のコードと対応し，変数 `funny` を `1` と設定することに相当します．

```sh
funny = DefineNumber[0, Choices{0,1},
  Name "Parameters/Apply funny mesh size field?" ];

Field[1] = MathEval;
If(funny)
  Field[1].F = "2*Sin((x+y)/5) + 3";
Else
  Field[1].F = "4";
EndIf
```



##### 2.2.b GUI からのメッシュサイズの調整

左側ツリーの中で Mesh / Define / Size fields を選ぶと Size fields というパネルが出現します．
MathEval を選んで （これも 上記 Field[1] に対応）Options を書き換えることができます．

例えば `6` に書き換えて，Set as background field のチェックを確認したうえで Apply すると修正が反映されて，このとき入力ファイル末尾に、新たな行が記録されます．

```bash
//+
Field[1].F = "6";

```

![image-012.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/967528a1-d320-cbc6-d8a4-3eac544d6e48.png)

こういう修正をしたときに、生成されるメッシュの例です．

![image-013.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/9f11bf21-58a4-f654-e9ae-642262467aac.png)


element 数は 4129 で，`Field[1].F = "4";` の場合の 8949 element に対して半減しています．一方で当然ながらフィレット部分や r のついた曲がり角付近の形状表現は低下しています．

さらに F を大きくするとどうなるでしょうか．

![image-014.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/dec3a890-90e3-3594-314f-958b0b174a81.png)


右下で緑系の角錐が見えているのは，内壁が外壁を貫通して外側に出てしまっている部分です．
もともとこの stl 構造では，壁の厚さが 3 mm 程度なので，`Field[1].F = "12";` の設定では、バックグランドフィールドの `12` という値が大きすぎることを示唆しています．

他のチュートリアルファイルも参考に、必要なパラメータの調整方法を探っていくことで、所望のリメッシュを追求できそうです． t13.geo ファイルのコメント部分には，以下のように CAD モデルが手に入るなら，そちらを使ったほうが正確なジオメトリーが得られるので好ましいと記載されています．

```
// Note that if a CAD model (e.g. as a STEP file, see `t20.geo') is available
// instead of an STL mesh, it is usually better to use that CAD model instead of
// the geometry created by reparametrizing the mesh. Indeed, CAD geometries will
// in general be more accurate, with smoother parametrizations, and will lead to
// more efficient and higher quality meshing. Discrete surface remeshing in Gmsh
// is optimized to handle dense STL meshes coming from e.g. imaging systems
// where no CAD is available; it is less well suited for the poor quality STL
// triangulations (optimized for size, with e.g. very elongated triangles) that
// are usually generated by CAD tools for e.g. 3D printing.

```

#### 2.3 形状分割ロジックの変更

Parameters の 3 つめのチェックボックス Create surfaces guaranteed to be parameterizable をチェックすると、開発チームの論文 [^2] に記載された parameterizable な分割が適用されます．

実質的に特徴線抽出が働かない条件として，角度パラメータを 160 として適用してみます．

![image-015.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/477a5ef8-0950-e07a-c776-531957a8d036.png)


リメッシュしてみるとこんな形状になりました．四角い開口の周辺の構造が落ちてしまっているので，特徴線も活用したほうが良さそうです．

![image-016.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/c4dc1c45-ddd0-571f-fb43-2dd4b94f4cf7.png)


角度パラメータ 80 ぐらいを選ぶと次の図のようになり，顕著な欠落はなくなりました．

![image-017.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/60b14017-d215-d5c3-8ff9-4c6c2d5d4db7.png)


### まとめ
GUI で設定したパラメータがジオメトリを規定する .geo ファイルに反映されること，キーストロークの `2` や `0` でメッシュ作成，リロードができるようになっているなど，gmsh には素早い試行錯誤をしやすい仕組みが用意されていることがわかりました．
その他のチュートリアルも試していくことで， gmsh を必要に応じて活用できるようになると思います．


[^1]: https://gmsh.info/doc/course/general_overview.pdf 42ページ以降 
[^2]: https://doi.org/10.1016/j.jcp.2020.109575  https://arxiv.org/pdf/2001.02542.pdf (arXiv)