DOP Solverをつくる

# はじめに
Houdiniでは様々なシミュレーションのSolverが用意されています。
Rigid Body Solver、POP Solver、Vellum Solver、Pyro Solver...などなど。
これらのSolverの中に入っていくとわかりますが、そのほとんどもノードの組み合わせで成り立っています。
つまり、自分でもシミュレーションのSolverを作れるということです。

ということで、今回は簡単なシミュレーションのSolverを作成していきます。

~~SolverSOPを使ってしまえば済む話ですね~~

# 今回の目的
- 簡単なシミュレーションのSolverをDOPで実装する
- DOPの理解を深める

# 大まかな流れ
DOP上でシミュレーションを実装する大まかな流れとしては、以下のようになります。

1. シミュレーションオブジェクト作成
2. Solver作成
3. その他設定

1ではDOP内でオブジェクトを作成し、シミュレーションを行う対象のジオメトリやその物性値を設定します。
2では毎ステップ計算を行うSolver部分を作成し、物体が移動・変形できるようにします。
3では重力や風などシミュレーション全体に影響を与える力や、衝突判定を行う物体を設定します。

# SOPでの準備
はじめにシミュレーションを行うシーンを適当に作ります。
![SOP.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/29ef778e-fb8b-f2e5-d5cb-dae6bdcde6f3.png)
箱とそれに衝突する球を用意し、球には質量と直径、初速を設定しています。

そしてDOP Networkを作成し、以降その内部で操作していきます。

# シミュレーションオブジェクト作成
すでに実装されているシミュレーションのObject DOP (FlipObject、SmokeObject、POPObject...など)の中に入ってみるとわかりますが、やっていることはどれも変わらず、

- EmptyObjectDOPを作成
- SOPからジオメトリデータを読み込み
- 必要な物性値を設定

となります。
これらを踏まえると、ノードは以下のようになります。
![Object.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/ce94489b-dee2-5d34-9a14-0f969a1742dc.png)

EmptyObjectにジオメトリと物性値のデータをApplyDataDOPを使って設定しています。
SOPGeometryDOPのSOP Pathにはジオメトリのパスを入れ、PhysicalParametersDOPでは適当な物性値を入力します。
今回は、ジオメトリのパスには先ほどSOPで用意した球を、物性値ではFrictionを0.2に設定しました。

以上でシミュレーションオブジェクトは完成です。

ちなみに、DOPにおいて灰色の線は**オブジェクトとデータ**、緑色の線は**データ**を受け渡ししていることを示しているそうです。

# Solver作成
オイラー法を使った簡単なシミュレーションのSolverを組んでいきます。
オイラー法というと難しいように聞こえますが、簡単に言えば、力で速度を更新し、速度で位置を更新する方法です。


ノードは以下のようになります。
![Solver.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/9ed42f16-d509-903d-71e8-a5a8ec94f006.png)
先ほど作ったオブジェクトをMultiSolverDOPの左に入力し、
GasExternalForcesDOPとGasIntegratorDOPをMergeしてMultiSolverDOPの右に入力しています。

ここで、GasExternalForcesDOPではオブジェクトに影響している外力で速度を更新しています。
外力が影響するのはSOPGeometryDOPで設定したジオメトリなので、Geometryパラメータにはそのジオメトリの名前を入力します。
![gasexternalforces.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/9b7929e2-a2e4-22d9-5b3e-eb1156b89d22.png)

また、GasIntegratorDOPでは速度をもとに位置を更新しています。
このとき、球の直径@pscaleを考慮した衝突処理を行いたいので、Use Pscale for SDFにチェックを入れておきます。
![gasintegrator.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/22b4259f-bcc7-61d3-f2bf-1a2853b1276b.png)

以上でSolverも完成です。

# その他設定
最後にGroundPlaneやStaticObject、Gravityを追加し、完成です。
![wholenodes.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/6f9e083e-ea35-964a-3d74-ecc3ac532946.png)

# 結果
![Euler.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/1b7de9c2-9900-078e-c42c-b157c384df20.gif)

球が壁に衝突して、跳ねる様子をシミュレーションすることができました。
球と壁の衝突判定やその処理はGasIntegratorDOPが行っているため、簡単にできました。

# 応用(シンプレクティック法)
応用として、シンプレクティック法のソルバーを実装します。
シンプレクティック法については書ききれないので、ネットや本で調べてみてください。
ざっくり言えば、オイラー法よりも精度が良い方法といった感じです。
(原理は全く異なるものですが...)

参考
- https://www.kspub.co.jp/book/detail/1572058.html
- http://cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/symplectic.html

今回は1次のバイラテラル・シンプレクティック数値積分でシミュレーションを行います。
計算式で表すと、以下のようになります。

```math
\left\{
\begin{array}{ll}
q(t+\frac{\Delta t}{2}) = q(t) + \Delta t \frac{\partial H}{\partial p}_{p=p(t)} \\
p(t+\frac{\Delta t}{2}) = p(t) - \Delta t \frac{\partial H}{\partial q}_{q=q(t+\frac{\Delta t}{2})} \\
p(t+\Delta t) = p(t + \frac{\Delta t}{2}) - \Delta t \frac{\partial H}{\partial q}_{q=q(t+\frac{\Delta t}{2})} \\
q(t+\Delta t) = q(t + \frac{\Delta t}{2}) + \Delta t \frac{\partial H}{\partial p}_{p=p(t+\Delta t)} \\
\end{array}
\right.
```
q:変位、p:運動量、H:ハミルトニアン

衝突処理についてはGasIntegratorDOPから衝突時のデータを取得し、反映させるようにします。

### ネットワーク
ネットワーク構成は以下のようにしました。
![symplectic.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/7f417d6e-3dd3-13c5-0a55-a629f27394b2.png)

先ほどのネットワークでのSolver部分だけ変更しています。
Solver部分を拡大すると以下のようにしています。
![symp_solver.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/8641e8cb-425a-bdbd-265d-3733e6cc1ee4.png)
各ノードは画像内のメモに記した通りの計算を行っています。

### 外力の読み込み
![symp_gasex.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/e745294f-259c-9c5c-13f3-02c8384639fd.png)
GasExternalForceDOPのパラメータでは先ほどと同様に外力を読み込みますが、今回は速度の更新を自分で行いたいため、ComputeNewParticleVelocitiesのチェックは外しておきます。

### 衝突判定前の位置と速度を保存
![symp_store_prev.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/e07c063d-7ad2-f89e-293a-a748dda62127.png)
シンプレクティック積分で位置と速度を更新したいので、念のためGasIntegratorDOPの計算の前に位置と速度を別アトリビュートとして保存します。
GeometryWrangleDOPを使っています。

### 衝突判定
![symp_collision.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/ace5cd35-324c-0649-16e2-fb4b6f360afc.png)

先ほどのGasIntegratorDOPでは位置の更新を行っていましたが、今回は衝突判定だけを行います。
MoveOutofSDFColliderは物体が衝突時に貫入しないためにチェックをつけます。
そして、StoreFinalPointPositionsは必要ないのでチェックを外しています。
また、AddImpactDataのチェックをつけ、衝突時のデータを使えるようにします。
最後に、IntegrateVelocityのチェックを外すことで、位置の更新を行わないようにします。

ここで、AddImpactDataについて追記します。
この項目にチェックをつけると、下図のように衝突時にImpactsというデータがシミュレーションオブジェクトに追加されます。
![impacts3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/019b89ff-278e-a679-0adc-818a6c04e0b4.png)
このImpactsには衝突衝突時の力やその方向、衝突したオブジェクトの情報が含まれています。
GeometrySpreadSheetで見るどんなデータが入っているかがわ分かりますが、このままだと扱いが難しいです。
そこで、SOPSolverDOPを使うことで簡単に扱えるようにします。



### 変位・運動量の設定
![symp_setpq.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/b3166d40-141c-da7a-4dc1-a945953c2989.png)
シンプレクティック積分では変位と運動量を使うので、それぞれ計算し、@q,@pのアトリビュートに納めます。

### シンプレクティック積分
積分を行う部分はSOPSolverDOPで行っています。
SOPSolverDOPを使うのは先ほど記述したように、Impactsのデータを簡単に扱うためです。
SOPSolverDOPの内部は以下のようにしています。
![solverinside.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/2c303bb8-c40d-6ce1-de18-0f7a676014dd.png)
赤枠で囲ったノードから衝突時の情報を得られるので、＠idを使ってアトリビュートをコピーしています。
衝突時の力については、＠impulseから衝突時の力積の大きさ、＠normalから力積が働く方向を得られます。
これらの積をとり、力積のベクトルを求め、運動量＠pに加えることで衝突処理を行います。

### force,impulseのリセット
GasLinearCombinationDOPをつかうことで、forceとimpulseのアトリビュートをリセットしています。
リセットしないと、外力が毎フレーム加算されてしまいます。
これは、GeometryWrangleDopを使ってコードを書いても問題ないかと思います。
![symp_wipe_f.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/ae3e994e-acce-8545-6488-e0270f1b1387.png)
![symp_wipe_imp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/0039e121-86c7-4c71-e2c4-3447a9581cc5.png)

# 結果
![Symplectic.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186607/39d7681c-74c9-d5a6-3cf1-d1991efbafb4.gif)

当然ながら、オイラー法とさほど変わりません。
強いて言えば、静止直前まで跳ねているように見えるくらいでしょうか。

正直、映像として使う分にはオイラー法で十分だと思います。
しかし、Houdiniのその柔軟性から、高精度のシミュレーションも実装できることが分かりました。
つまり、Houdiniを力学解析ソフトとして活用することもできるのではないかと思いました。(暴論)

# まとめ
シミュレーションオブジェクトとSolverを自作し、簡単なシミュレーションを行うことができました。
また、Solverが自作できるようなったことで、より精度の良いアルゴリズムや最新の研究を実装し、既存のシミュレーションと組み合わせるといったことが簡単にできるようになったと思います。

普段気にすることのないDOP Solverたちの中身を見ると新しい発見もあるかと思います。
是非、Solverの自作を試してみてください。