[Unity|UE5] 論文を読んで紐物理を実装してみよう(PBD: Position Based Dynamics)

はじめに

この記事は、物理シミュレーションのPosition-Based Dynamicsの論文を読み、実装してみるというものです。一緒に勉強会をしてる知人に教えてもらった論文を読んで実装してみました。

Unityでの実装を行いましたが、UEでの実装も近々サンプルを用意してみたいと思っています。 ただ、記事のほとんどが論文に関してなので、ゲームエンジン依存の実装の話などはありません。
[追記]UE5でのサンプルを追加しました

大学生や高校生を含め、論文を読んで実装をしたことがない人にも分かるように丁寧に読み進めていくことを意識して書きました。しかしながら、大学で学ぶような内容もありますので、少々読みづらいかもしれません。

• 物理学/解析力学
• 微分・偏微分などの大学数学

数式の途中の式変形なども書きましたが、ちゃんと読もうと思うとカロリー高めかもしれません。

また、間違いなどありましたら一報いただけると嬉しいです。

PBD

以下がその論文になります。これを読みながら実装してみようと思います。
https://matthias-research.github.io/pages/publications/posBasedDyn.pdf

CEDEC2017の中川さんのスライドがとても分かりやすいです。事前にこちらに目を通しておくと理解がしやすいかもしれません。
https://github.com/nobuo-nakagawa/cedec2017/blob/master/cedec2017.pdf

1. Introduction

物理シミュレーションの方法について

• 従来の力ベース
• 従来とは異なる位置ベース

を説明しています。

力ベース

The traditional approach to simulating dynamic objects has been to work with forces.
Newton’s second law of motion relates forces to accelerations via the mass. So using the density or lumped masses of vertices, the forces are transformed into accelerations.

従来の動的オブジェクトのシミュレーション手法は、力や力積を用いたものだそうです。

その力を用いた計算の基本となるのは、ニュートンの第二法則になります。この法則は、力は質量と加速度に関係しているという式です。物理シミュレーションでも、力を与えて、加速度が生まれれば物体の位置を動かせれるということですね。

いきなり、ニュートンの第二法則と言っても結構忘れてると思います。軽く物理のおさらいをしておきましょう。

ニュートンの第二法則は以下になります。以下の式をニュートンの運動方程式と呼びます。

$$ma = F \\ \begin{aligned} \begin{cases} m &: 質量\\ a &: 加速度\\ F &: 力 \end{cases} \end{aligned}$$

At the beginning of each time step, internal and external forces are accumulated.

タイムステップごとに、内力(internal force)と外力(external force)が累積され、累積された$F$によって速度が更新されます。

これはUnityでもそうですね。 AddForce を呼ぶことで同一フレームで複数の力が加えれます。

_rigidbody.AddForce(internalForce1);
_rigidbody.AddForce(internalForce2);

AddForce する以外にも、さらにコライダの衝突によって外力も発生するでしょう。
それらが累積され、次の FixedUpdate で速度が更新され、物体が移動します。

例えば、車の実装であれば「アクセルを踏んだ状態」であれば常に前方向に AddForce で力を加えていると思います。さらに、壁にぶつかればめり込まないように制約するような外力が加わります。そして最終的にRigidbodyに累積された力 $F$ を用いて、位置が更新されます。

$F$ から位置を更新するには、運動方程式から加速度 $a$ を求めることで計算できます。加速度 $a$ が分かれば、次のタイムステップの速度を求めることができ、その速度からさらに位置を求めることができます。

$$\begin{aligned} ma &= F \\ a &= \frac{F}{m} \end{aligned}$$

この辺りについては、別の技術ブログに書いたので詳しい説明はやめておきます。
https://zenn.dev/nrdev/articles/44c5a563846e06

( $F$ によっては、そこから速度を求めて位置を更新して、となるとチラチラとめり込みそうですね。)

位置ベース

従来の方法について見ました。次は本手法で扱われている位置ベースについてです。

computer games it is often desirable to have direct control over positions of objects or vertices of a mesh. The user might want to attach a vertex to a kinematic object or make sure the vertex always stays outside a colliding object.

コンピュータゲームの場合は、衝突物体にめり込まず外側に押し出した状態にしたいので、頂点をそのまま外に直接移動させれたほうが便利そうです。この手法では、この直接位置を移動させるというアプローチによるシミュレーションだそうです。

2. Related Work

関連研究自体は割愛しますが、本手法について説明があったのでおさえておきます。

We define general constraints via a constraint function as [BW98] and [THMG04]. Instead of computing forces as the derivative of a constraint function energy, we directly solve for the equilibrium configuration and project positions.

論文 [BW98] [THMG04] のように、「制約関数」により制約(constraints)を一般化し表します。

そして、制約関数のエネルギー(a constraint function energy)の微分として力を計算する代わりに、本手法では直接的に制約を満たすような配置を解いて、位置を更新するようです。

制約関数のエネルギーが何なのかはまだ分かりませんが、本論文でやりたいこととしては 「物理学におけるエネルギーというものを考えて、そのエネルギーによって生じる力を計算して位置を更新して物理をシミュレーションしていたが、力ではなく直接的にそのエネルギーが保存されるように位置を変えてしまおう」 ということです。

「制約関数エネルギーの微分として力を計算」…。

チョットナニヲイッテイルカワカラナイ…。

本手法では、それをしないそうなのでスキップしてもいいのですが、せっかくです。
エネルギーが何なのか分からないともどかしいので、ちょっと小難しそうですが、もう少し丁寧に見てみましょう。

1. 制約関数エネルギーとは何なのか
2. それを微分すると力を計算できるのか

を見ていきたいと思います。

◆ 制約関数エネルギーとは何なのか

そもそもこれが何なのか見ていきたいと思います。まずは挙げられた論文[THMG04]で制約関数に関する記述を調べてみます。

3.1 Potential Energies
In order to represent deformations of objects, we consider constraints of the form $C(p_0,\dots, p_{n−1})$ . These scalar functions depend on mass point positions $p_i$ . They are zero if the object is undeformed. In order to compute forces based on these constraints we consider the potential energy

$$E(p_0,\dots,p_{n-1}) = \frac{1}{2}kC^2 \tag{1}$$

with k denoting a stiffness coefficient.

早速「ポテンシャルエネルギー」というのが出てきました。制約関数エネルギーのことでしょうか。

制約関数に関する記述もありますね。物体の変形を表すために制約の式 $C(p_0,\dots, p_{n−1})$ を定義しています。$p_i$ は質点の座標ですね。制約の式は、座標 $p_i$ を引数にとってスカラー値を返す関数のようです。

ポテンシャルエネルギーというものを、この制約関数 $C$ によって定義しているようです。
そしてそのポテンシャルエネルギーによって力を計算するみたいですね。

「なぜ制約関数を用いて、 $(1)$ のようにポテンシャルエネルギーを定義しているのか」はまだ謎ですが、定義の関係は把握できました。

もう少し読み進めてみましょう。

The first potential energy $E_D$ considers all pairs of mass points that are connected by tetrahedral edges. $E_D$ represents energy based on the difference of the current distance of two points and the initial or rest distance $D_0$ with $D0 \neq 0$:

$$E_D(p_i,p_j) = \frac{1}{2}k_D(\frac{|p_j-p_i| - D_0}{D_0})^2 \tag{4}$$

二点 $p_i$ と $p_j$ の頂点ペアを結んだ辺について、エネルギー $E_D$ を考えるようです。その「二点 $p_i$ と $p_j$ の頂点ペアを結んだ辺」を制約する関数 $C$ は、$(1)$ と $(4)$ から以下のようになっていることが分かります。

$$C(p_i, p_j) = \frac{|p_j-p_i| - D_0}{D_0} \\ \begin{aligned} \begin{cases} D_0 &: p_iとp_jの2点間の初期の距離\\ \end{cases} \end{aligned}$$

これは、現在の2点間の距離が、「初期の距離」からどれだけズレてしまっているか、を表しています。つまり、 $C(p_i, p_j) = 0$ であれば、ズレはなく 安定した状態 であり、これを保てば距離の制約をしていることになります。

なお、「初期の位置」ではありません。2点間の「初期の距離」です。並進移動した場合は2点間の「距離」は変わりませんので、ズレは発生しません。つまり、並進移動では $C(p_i, p_j) = 0$ が保たれます。

では、並進移動ではなく、片方だけが動いてしまった場合はどうでしょうか。

$C(p_i, p_j) = -1$ であれば縮みすぎているので元の長さになるように伸ばす必要がありますし、 $C(p_i, p_j) = 1$ であれば伸びすぎているので元の長さに縮める必要があります。これを制約として見ています。

ポテンシャルエネルギーに話を戻します。 再び式 $(1)$ を見てみます。
$k$ は弾性係数(stiffness coefficient)と書いてあります。

$$E(p_i,p_j) = \frac{1}{2}kC(p_i,p_j)^2 \tag{1} \\ \begin{aligned} \begin{cases} C(p_i,p_j) &: p_iとp_jの2点間の初期距離からのズレ\\ k &: 弾性係数(stiffness \; coefficient) \end{cases} \end{aligned}$$

この式は、弾性力の位置エネルギー (弾性エネルギー) の公式と同じですね。

https://w3e.kanazawa-it.ac.jp/math/physics/high-school_index/mechanics/energy/henkan-tex.cgi?target=/math/physics/high-school_index/mechanics/energy/elastic_potential_energy.html

https://ja.wikipedia.org/wiki/弾性エネルギー

ちなみに位置エネルギーは英語で potential energy: ポテンシャルエネルギー です。以降は、論文中でもよく出てきますので可能な限りはポテンシャルエネルギーと呼びます。

弾性力のポテンシャルエネルギー $U=\frac{1}{2}kx^2$ は「バネを自然長から $x$ だけ伸ばした」時のエネルギー量です。高校で物理を学んだ方はこの式に覚えがあるかもしれませんね。制約関数 $C(p_i, p_j)$ はこの $x$ に相当します。「辺を自然長(初期距離) $D_0$ から $C(p_i,p_j)$ だけ伸ばした(ズレた)」時のエネルギー量が $E(p_i,p_j)$ ということになります。

実際に、硬いバネをイメージしてみると、先ほどの制約関数 $C(p_i, p_j)$ との関係が想像しやすそうです。

$C(p_i, p_j) = -1$ であれば縮みすぎているので元の長さになるように伸ばす必要がありますし、 $C(p_i, p_j) = 1$ であれば伸びすぎているので元の長さに縮める必要があります。

とは言え、エネルギーは「量」の話なので、今は「元に戻す力」は忘れた方がいいですね。先ほどの記述を見るに、それはきっと微分したら出てくるのでしょう。

制約関数 $C(p_i, p_j)$ の値が大きければポテンシャルエネルギーの「量」は大きくなります。
バネを想像すると伸ばした状態ではエネルギーが大きい状態になっています。エネルギーが溜まっているなら、きっと手を離したら素早く元の位置に戻ることでしょう。

制約関数 $C(p_i, p_j)$ の値が $0$ であればポテンシャルエネルギーの「量」は $0$ になります。バネを想像すると初期位置であれば何も力は働かずにエネルギーの「量」は0です。

$$U = \frac{1}{2}k(0)^2 = 0$$

ということで、制約関数とそれからなるポテンシャルエネルギーの関係が「弾性力のポテンシャルエネルギー」としているのではないか、というところまで分かりました。

謎であった「なぜ制約関数を用いて、 $(1)$ のようにポテンシャルエネルギーを定義しているのか」は掴めた気がします。

そして、元々の疑問だった「制約関数エネルギー(a constraint function energy)」がどういうものなのかイメージできたと思います。

◆ それを微分すると力を計算できるのか

制約関数エネルギーが物理学における弾性力のポテンシャルエネルギーであることは見てきました。(実際、弾性係数を式に含んでいます)

とはいえ！！！

「弾性力のポテンシャルエネルギー」も「公式がこうだから」という、ざっくり雰囲気でしか分かってません。そもそもエネルギーって漠然としていてよく分かりません。

なので、ポテンシャルエネルギーとは何なのか調べていきたいと思います。その過程で力との関係性も見えてきて、なぜ 微分すると力になるのか がきっと分かってくるのではないかと思います。

調べてると東京大学の講義資料が非常に分かりやすかったです。
読んでもらった方が確かだと思いますが、こちらでも説明しておきます。
http://hep1.c.u-tokyo.ac.jp/~kazama/emgakushuin/em-note-2019-2.pdf

「物体が移動する」ことに着目します。
ある位置 $P_0(\vec{x_0})$ から 別の位置 $P(\vec{x})$ に物体が移動します。
( $\vec{x}$ は3次元なら $\vec{x} = (x_0,x_1,x_2)$ のような位置を表すベクトルです )

このとき、この「 $P_0$ から $P$ までの経路で加わった力 $F$ の合計」を仕事 $W[J]$ と言います。仕事をしなければ動かないのはまだ想像しやすいですね。

仕事について軽く調べておきましょう。

仕事とは

仕事(物理学)
物理学における仕事（しごと、英語: work）とは、物体に加わる力ベクトルと、物体の変位ベクトルの内積によって定義される物理量である。
仕事は、エネルギーを定義する物理量であり、物理学における種々の原理・法則に関わっている。
https://ja.wikipedia.org/wiki/仕事_(物理学)

エネルギー
物理学において、エネルギー（独: Energie）またはエナジー（英: energy）は、仕事をすることのできる能力のことを指す[1][2][3]。物体や系が持っている仕事をする能力の総称[4]。エネルギーのSI単位は、ジュール（記号：J）である。
https://ja.wikipedia.org/wiki/エネルギー

この「仕事」は「エネルギーを定義する物理量」のようです。単位がジュールで、家電とかで馴染みあがありますね。こちらも調べてみましょう。

ジュール
ジュール（英: joule、記号：J）は、エネルギー、仕事、熱量、電力量の単位である。
1 ジュールは標準重力加速度の下で 102.0 グラム（ほぼキウイフルーツまたはマンガン単1乾電池の質量）の物体を 1 メートル持ち上げる時の仕事に相当する。

定義

国際単位系 (SI) 及び日本の計量法におけるジュール J の定義は以下の通り。
・ 1 ニュートンの力がその力の方向に物体を 1メートル動かすときの仕事

仕事の単位に使われているので、その単位あたりの量がどういうものなのかが書かれてますね。 $1N$ で1メートル動かすと $1J$ のようです。たしかに、「経路上で物体にかけた力」について考えていたのでなるほどです。

これがエネルギーということでしたが、正直まだちょっとイメージしにくいです。
以下のように、いろいろ形態があるようで、「これだ！」という実態があるわけではなく、あくまで「仕事をする能力」について定量的に定義したものなのでしょうか...。

エネルギー

概説

このようにエネルギーという語・概念は「物体が仕事をなし得る能力」を意味したが、その後、自然科学の説明体系が変化し、熱・光・電磁気もエネルギーを持つことが知られるようになり、さらに、質量までがエネルギーの一形態である、と理解されるようになった[2]。

話が本題とズレて行きそうなので、「仕事がなにか」はざっくりにして、そろそろ話を戻します。
(Wikipediaを見てると無限にたどっていく、あるある現象に陥りました)

まとめると、「経路上を移動させるときに物体に加えた力の合計は仕事であり、これはエネルギーとして見れる」という感じでしょうか。

仕事は経路上で加えた力の合計ですので、この仕事 $W$ を数式で表すと以下のようになります。合計は積分ですね。

$$W = \int_{P_0}^{P} F d\vec{x}$$

次に、位置が変わっているので「位置 $P_0$ の位置エネルギー( ポテンシャルエネルギー )」と「位置 $P$ の位置エネルギー」は変化しています。これも位置が高い物体の方がエネルギーが高いというのは想像しやすいですね。なお、位置エネルギーも同じようにエネルギーですから単位も $J$ (ジュール)です。

つまり、「物体を動かす」ことに着目すると、以下のことが言えます。

• ある経路上で物体を動かすために与えた力の合計は仕事(エネルギー)
• その経路の始点と終点では、位置エネルギーが変化している

例を使って想像しておきましょう。
ジェットコースターで考えると最初は「低い所 $P_0$ 」から「高い所 $P$ 」へ行く際に、機械を使って高い所まで運びます。このときに使う力の合計が仕事です。(ロスが多いので使った電気がすべて位置エネルギーにはなりませんが)

また、低い所の位置エネルギーよりも、高い所の位置エネルギーの方が大きくなっており、位置エネルギーも変化しています。ジェットコースターはこの位置エネルギーを使って速度を上げ、下へと降りていきます。

というように、「物体が移動する」ことに着目することで、そこに働く力(合計が仕事 $W[J]$ )と、位置変化によるポテンシャルエネルギーの変化について、見ることができます。

さらに、いろんな経路を考えてみます。
これらのどの経路をたどっても仕事 $W$ が同じ場合は、 $F$ は保存力と呼ばれます。

保存力であるとき、 $P_0$ から $P$ へ行く際に失った(もしくは得た)位置エネルギーは、仕事と一致しています。数式でこの関係を表すと次のようになります。

$$U(\vec{x}) - U(\vec{x_0}) = -W$$

符号について

一応、符号についても説明しておきます。

またジェットコースターを例にします。
ジェットコースターを上まで運ぶには、その運搬機械を動かす電気がいりますね。この仕事 $W$ だけエネルギーを消費して失っています。一方で、低い位置 $P_0$ から高い位置 $P$ に動かした場合は、位置エネルギーが増えていますね。

その増えた位置エネルギー $U(\vec{x}) - U(\vec{x_0})$ だけ、代わりに機械を動かす電気を消費しています。なので、等号で表すときにはマイナスになっています。

$$貯めた位置エネルギー = 仕事で失ったエネルギー = -仕事のエネルギー$$

この移動によって生じる差分のエネルギーをポテンシャルエネルギー $U$ とします。そうすると以下のように書けます。

$$\begin{aligned} U &= U(\vec{x}) - U(\vec{x_0}) \\ &= -W \\ &= -(\int_{P_0}^{P} \vec{F}\cdot d\vec{x}) \\ &= \int_{P_0}^{P} (-\vec{F})\cdot d\vec{x} \end{aligned}$$

ポテンシャルエネルギーと力の関係式が見えてきました！！

エネルギーを微分すると力を計算できる、という話でしたね...。
び、微分をしてみましょう…！！

両辺を順番に時間微分していきます。

左辺の $U$ から微分します。なお、$U$ は位置ベクトル $\vec{x} = (x_0,x_1,x_2)$ の関数です。$t$ で全微分するとベクトルの内積の形になっています。

$$\begin{aligned} \frac{dU}{dt} &= \frac{\partial U}{\partial x_0}\frac{d x_0}{dt} + \frac{\partial U}{\partial x_1}\frac{d x_1}{dt} + \frac{\partial U}{\partial x_2}\frac{d x_2}{dt} \\ &= (\frac{\partial U}{\partial x_0}, \frac{\partial U}{\partial x_1}, \frac{\partial U}{\partial x_2}) \cdot (\frac{d x_0}{dt}, \frac{d x_1}{dt}, \frac{d x_2}{dt}) \\ &= \nabla_{\vec{x}}U \cdot \frac{d \vec{x}}{dt} \\ \end{aligned}$$

次に、右辺の積分の式を微分します。
ただ、申し訳ありませんが、途中式を書いてみたのですが自信がないです。
$P_0$ は開始地点なので、 $t=0$ での座標として $x(t_0)$ としました。$P_0$から$P$までの距離について考えていましたが、そこは任意の地点なので、$P$は任意の時間 $t$だけ移動したときの座標として $x(t)$としました。

$$\begin{aligned} \frac{d}{dt} \int_{x(t_0)}^{x(t)} (-F(\overrightarrow{x(t)})) d \overrightarrow{x(t)} &= -\frac{d}{dt} \int_{x(t_0)}^{x(t)} (F(\overrightarrow{x(t)})) d \overrightarrow{x(t)} \\ &= -\frac{d}{dt} \int_{x(t_0)}^{x(t)} (F(\overrightarrow{x(t)})) d \overrightarrow{x(t)} \frac{dt}{dt} \\ &= -\frac{d}{dt} \int_{x(t_0)}^{x(t)} (F(\overrightarrow{x(t)})) \frac{d \overrightarrow{x(t)}}{dt}dt \\ &= -\frac{d}{dt} (\int_{x(t_0)}^{x(t)} (F(\overrightarrow{x(t)})) dt) \cdot \frac{d \overrightarrow{x(t)}}{dt} \\ &= -\frac{d}{dt} (\mathbb{F}(\overrightarrow{x(t)}) - \mathbb{F}(\overrightarrow{x(t_0)})) \cdot \frac{d \overrightarrow{x(t)}}{dt} \\ &= -(\frac{d \mathbb{F}(\overrightarrow{x(t)})}{dt} - \frac{d\mathbb{F}(\overrightarrow{x(t_0)})}{dt}) \cdot \frac{d \overrightarrow{x(t)}}{dt} \\ &= -(\frac{d \mathbb{F}(\overrightarrow{x(t)})}{dt} - 0) \cdot \frac{d \overrightarrow{x(t)}}{dt} \\ &= -F(\overrightarrow{x(t)}) \cdot \frac{d \overrightarrow{x(t)}}{dt} \\ \end{aligned} \\ \begin{aligned} \begin{cases} \mathbb{F}(\vec{x}) : F(\vec{x}) = \mathbb{F}(\vec{x})' となるような、F(\vec{x})の原始関数 \mathbb{F}(\vec{x}) \end{cases} \end{aligned}$$

よって

$$\begin{aligned} \nabla_{\vec{x}}U \cdot \frac{d \vec{x}}{dt} &= -F \cdot \frac{d \vec{x}}{dt} \\ \nabla_{\vec{x}}U &= -F \end{aligned}$$

ポテンシャルエネルギー $U$ の微分(勾配)が力 $F$ となっています！

運動を考えた時にエネルギーの変化を仕事を用いて関係式が得られること、
その仕事が力の積分で表されていること、これらからエネルギーを微分することで、その仕事の経路の任意の場所での力が計算できることが分かりました。

長くなりましたが、疑問だった

• 制約関数エネルギーを「弾性力のポテンシャルエネルギー」として表していること
• そのポテンシャルエネルギーを微分すると力が計算できること

が分かりました。

ついでに、弾性力のポテンシャルエネルギーがなぜ $\frac{1}{2}kx^2$ なのかも見ておきましょう。

弾性力はフックの法則よりバネ定数が $k$ であるとき、 $F(x)=kx$ です。

$$\begin{aligned} U &= -W \\ &= -\int_{0}^{x} F(x)\cdot dx \\ &= -\int_{0}^{x} (-kx)\cdot dx \\ &= k\int_{0}^{x} x\cdot dx \\ &= k \cdot (\frac{1}{2}x^2 - \frac{1}{2}0^2) \\ &= \frac{1}{2}kx^2 \end{aligned}$$

バネ定数によって強弱を変えられるような (弾性)力 を考え、その仕事を考えることでエネルギーが計算できました。

制約関数エネルギーもこのような (弾性)力 によって、距離を制約するようにモデル化していると考えられます。制約関数 $C(x_0, x_1)$ は2点間の初期距離からのズレであるスカラー値 $x$ でした。安定した初期位置の状態 $C(x_0, x_1) = 0$ (自然長)となるように弾性力が加わる、とモデル化して考えれば、上述の弾性力によるポテンシャルエネルギーで話を進めることができます。

これでやっと、制約関数エネルギーの式についても理解できました。どういう力を考えて制約をするのか、というところから来ているようです。

3. Position Based Simulation

やっと本題です。

3.1 Algorithm Overview

◆ 表記

その前に、後々読んでいく上で、どういう記号を使っているかは、混乱を回避する上で重要なのでしっかりと押さえておきたいと思います。
(私は最初に軽く流し読みで済ましてしまい、混乱しました…)

• 頂点は $N$ 個、要素は $i \in [1, \cdots, N]$
  • 後述の $n_j$ は、小文字だからといってこの要素ではないので注意
  • 他にも関係のない $n$ が稀に出てきます。文脈から別物だと分かると思いますが気を付けましょう。
• 制約は $M$ 個、要素は $j \in [1, \cdots, M]$

制約 $j$ について。

• $j$ で制約される頂点の数が $n_j$
  • 集合の基数/濃度(cardinality)
  • 多くは個々の $j$ について明示する必要はないので省略して $n$ と表記されてたりします
• $j$ の制約関数 $C_j$ は $3n_j$次元から $1$次元への写像
  • $C_j: \mathbb{R}^{3n_j} \to \mathbb{R}$
• $\{i_1, \cdots, i_{n_j}\}$ としたとき、この頂点 $i$ は、制約 $j$ に制約されている頂点の集合
  • でも、その集合の要素は $i_k \in [1, \cdots, N]$ であり、全体の頂点のどれか
  • ややこしい…。添え字、表記被りしない…？
• 剛性係数は $k_j$ で $0 \sim 1$ の間の値をとる

◆ アルゴリズム

アルゴリズムは以下の通りです。

(1)~(3)

ただの初期化です。

(4)

ここからはステップ時間ごとに進むシミュレーションのループ処理です。

(5)

外力(external forces)を加え、速度を更新します。一応、式も丁寧に見ましょう。

力 $F$ にステップ時間 $\Delta t \cdot \frac{1}{m_i}$ を掛けると、ステップ時間における速度変化 $\Delta v$ が得られます。
重力 $F=m_i \cdot g$ だけ考えるのであれば、 $\Delta v = g \cdot \Delta t$ となります。

$$\begin{aligned} F &= m_i a \\ F &= m_i \frac{\Delta v}{\Delta t} \\ F \cdot \Delta t \cdot \frac{1}{m_i} &= \Delta v \\ \end{aligned}$$

(6)

簡易化のためダンピングは飛ばします。減衰に関する処理です。

(7)

推定位置(estimates) $p_i$ を計算します。テイラー展開を用いた一次までの近似です。

$$p_i = x_i(t+1) = x_i(t) + \Delta t v_i$$

その時点での速度 $v(t)$ でざっくり $\Delta t$ だけ先の時間は $v(t) \Delta t$ だけ移動していると近似しています。

(8)

衝突制約条件(collision constraints)による位置の変更です。ポジションベースの解決なので、推定位置 $p_i$ を制約にもとづき変更します。

これも簡易化のため、今回は省略します。

(9)~(11)

ソルバーの処理です。ソルバーとはシミュレーションなどにおいて、数値計算を実行して方程式を解くシステムです。この論文では反復法によって、方程式の解を近似計算しているようです。

ここがメインの処理です。制約 $C$ を満たすように推定位置 $p_i$ を修正します。

例として、単純な距離制約関数は $C(\mathbf{p_1} , \mathbf{p_2}) = |\mathbf{p_1} - \mathbf{p_2}| - d$ です。

(12)~(15)

現在位置 $x_i$ から推定位置 $p_i$ へ頂点を直接移動させます。

そして、頂点の移動から速度を決定します。外力を加えた速度から、位置を推定し、制約を満たすように位置を加工し、最終的な推定位置に移動しました。その移動を行うような速度に設定します。

(16)

衝突のあった頂点では、摩擦(friction)や反発係数(restitution coefficients)によって速度を更新します。

(8)と同様に簡易化のため衝突に関しては省略します。

(17)

ループ終了。以上がアルゴリズムの概要です。

ソルバーの制約 $C$ による推定位置 $p_i$ の修正を抑えれば実装できそうです。次にそれを見ていきましょう。

3.3 Constraint Projection

制約射影(constraint projection)を行います。

The method we propose for constraint projection conserves both momenta for internal constraints. Again, the point based approach is more direct in that we can directly use the constraint function while force based methods derive forces via an energy term (see [BW98, THMG04]).

力ベースの手法では [BW98, THMG04] にあるようなエネルギー項によって得られる力を使うのに対して、位置ベース手法では制約関数から直接位置を更新します。

参照論文にて、エネルギーと力についてしっかりおさえたのでイメージしやすいですね。
ですが、上述の通り、使いません…！！直接位置を更新します！ただ、徒労ではありません。このPBDを改善したXPBDではポテンシャルエネルギーが出てくるので無駄にはなりません…！！そのうちXPBDの論文読みも書きたいと思います。

では、論文に沿って読んでいきたいと思います。

◆ 制約を満たせば運動量は保存される

省略

シミュレーションループ内で直接位置を移動させる際に最も重要なことは、線形運動量(linear momentum)と角運動量(angular momentum)の保存です。

それぞれ保存されるには以下であることです。

$$\sum_{i} m_i \Delta \mathbf{p_i} = 0 \tag{1}$$

$$\sum_{i} \mathbf{r_i} \times m_i \Delta \mathbf{p_i} = 0 \tag{2}$$

運動量は「質量と速度の積」、角運動量は「位置座標と運動量の外積」になります。

$\Delta \mathbf{p_i}$ はステップ時間 $\Delta t$ における推定位置までの移動距離です。ゆえに速度と言えます。なので、上式のように質量に掛けて運動量としています。

この運動量の保存の関係を崩すような移動を考えるとき、ゴースト力(ghost forces)と呼ぶ外力を導入します。

例えば、オブジェクトを引っ張ったりした場合などです。

ただ、運動量を保存するには、個々の頂点間に与えられている内部制約(internal constraints)だけでOKです。

そしてその内部制約として与えられる接続制約(attachment constraints)などによって、オブジェクト全体の運動量は保存されるので、全体が帳尻合わせで動いたりします。

要は、頂点の連結をした紐を考えた場合、紐をつかんで移動させるにはゴースト力を与えて動かせばよく、運動量の保存は与えた紐の頂点間の接続制約によって保存される、ということだと思います。

◆制約関数 $C$ の値変化

制約関数 $C$ の値は初期距離からのズレという話はしました。なので、おさらいのような感じになります。

まず表記について。

頂点 $\mathbf{p_1}, \cdots, \mathbf{p_n}$ であるとき $\mathbf{p} = [ \mathbf{p_1}^T, \cdots, \mathbf{p_n}^T]^T$ と表記しています。

$n$ はある制約 $j$ において、それが制約する頂点の数でした。

論文中では文中で表記するために転置記号 $T$ を用いていますが、以下のように単に座標をまとめてベクトル表記しただけです。

$$\mathbf{p} = \begin{bmatrix} [x_1, y_1, z_1] \\ \vdots \\ [x_n, y_n, z_n] \\ \end{bmatrix}$$

For internal constraints, $C$ is independent of rigid body modes, i.e. translation and rotation. This means that rotating or translating the points does not change the value of the constraint function.

(rigid body modesってどういう訳なんでしょう？)

制約関数 $C$ の値は並進運動や回転運動によって変化しない、独立な値です。2点間の距離による弾性エネルギーに関する制約関数であれば、2点間が並進運動しようが回転しようが距離は変わりません。

ですが並進運動ではなく、片方の頂点だけが移動した場合は変わってしまいます。

図にすると以下のようになります。実線が元の位置、そこから片方の頂点だけ移動した場合を赤色で図示しています。制約関数 $C$ の値は初期距離からのズレでしたね。図で見て分かるように、 制約関数 $C$ の値が変化するのは、2点間を結んだ辺と同じ方向に移動したときです。 移動してもそれが回転運動であれば値は変化しません。

◆初期距離からズレがない状態に位置を修正する

今やりたいこととしては、制約によって、修正される推定位置 $\mathbf{p_i}$ を求めることでした。

頂点が外力によって移動したとき、この制約を満たすように図の赤点をもとの位置、つまり $C=0$ となる位置まで動かしたいわけです。(図を見ると辺に沿って動かせば良さそうではあります)

$C=0$ になるまで、 $C$ が変化するような方向へ頂点を動かせば いいわけです。

それはどの方向でしょうか？

「頂点を動かしたとき $C$ が変化するような方向」 は、つまり $C$の傾き であり、微分を用いて得ることができます。

二次元で考えてみます。
$y=3x$ であれば、 $x$ を $1$ だけ動かしたら $y$ が正の方向へ $3$ 増えます。
現在が $x = -3$ で $y=3x = -9$ のとき、 $y=0$ にしたいなら正の方向に変化するように $x$ を加えていくようなイメージです。

制約関数 $C$ は、3次元座標 $\mathbf{p}$ の関数 $C(\mathbf{p})$ であり、スカラー値を返します。そのため、微分を考えるには、3次元空間に関して微分する必要があります。

このように三次元空間に対してスカラー値が割り当てられるような空間をスカラー場と言いますが、このスカラー場のスカラーをベクトルで微分するような演算を「勾配(gradient)」を取ると言います。

定義は以下のようになっています。$\nabla$ はナブラと読みます。

$$\begin{aligned} \nabla f &= (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z} ) \\ &= (\Delta xだけ動かしたときの傾き、\Delta yだけ動かしたときの傾き, \Delta zだけ動かしたときの傾き) \end{aligned}$$

それぞれの軸に対して偏微分をし、それぞれの軸方向への変位を考えたベクトルになっています。勾配というだけあって、各座標地点におけるスカラー値 $C$ の傾きを表しています。

スカラー場とその勾配の例としては部屋の温度です。部屋の空間座標ごとに温度というスカラー値があります。この温度の変化の傾きが、座標ごとに決まるはずです。きっと結露するような季節は窓際の「勾配」は大きいことでしょう。

話を戻します。

制約関数 $C$ の勾配を考えます。$C$ は二点間の辺に沿って移動した場合にだけ値が変わりました。つまり、その方向に勾配ができているはずです。その勾配ベクトル $\nabla_{\mathbf{p}} C$ こそ、 $C=0$ を満たすように頂点を移動させるべき方向です。

今考えるべき修正量(correction) $\Delta \mathbf{p}$ は、 $\mathbf{p}$ に関する勾配 $\nabla_{\mathbf{p}} C$ に沿った方向のベクトルだと言えます。よって、この修正量 $\Delta \mathbf{p}$ は以下のように表せると論文に書かれています。

$$\Delta \mathbf{p} = \lambda \nabla_{\mathbf{p}}C(\mathbf{p}) \tag{4}$$

外力によって移動される推定位置が $\mathbf{p}$。それを制約を満たす位置に戻す修正量 $\Delta \mathbf{p}$。そして修正後の位置 $\mathbf{p} + \Delta \mathbf{p}$ となります。よって、 $C(\mathbf{p} + \Delta \mathbf{p})=0$ となるような修正量 $\Delta \mathbf{p}$ を探すことが目的となります。その式を、論文では以下のように近似しています。

$$\begin{aligned} C(\mathbf{p} + \Delta \mathbf{p}) &= 0 \\ &\approx C(\mathbf{p}) + \nabla_{\mathbf{p}}C(\mathbf{p}) \cdot \Delta \mathbf{p} \tag{3} \end{aligned}$$

これはどういった近似なのでしょうか？
現在の位置での $C(\mathbf{p})$ に何かを足しています。$\Delta \mathbf{p}$ だけ移動したら $C(\mathbf{p})$ を打ち消すような何かを足しています。

$$\begin{aligned} C(\mathbf{p}) + A &= 0\\ C(\mathbf{p}) &= -A \\ \end{aligned}$$

その位置$\mathbf{p}$での勾配 $\nabla_{\mathbf{p}}C(\mathbf{p})$ と 未知の修正量 $\Delta \mathbf{p}$ の内積を足しているようです。この二つは同じ向きでしたね。 $\Delta \mathbf{p} = \lambda \nabla_{\mathbf{p}}C(\mathbf{p})$ でした。これを式 $(3)$ に代入してみます。未知なのはスカラー $\lambda$ だけです。

$$\begin{aligned} C(\mathbf{p}) + \nabla_{\mathbf{p}}C(\mathbf{p}) \cdot \Delta \mathbf{p} &= C(\mathbf{p}) + \nabla_{\mathbf{p}}C(\mathbf{p}) \cdot (\lambda \nabla_{\mathbf{p}}C(\mathbf{p})) \\ &= C(\mathbf{p}) + \lambda (\nabla_{\mathbf{p}}C(\mathbf{p}) \cdot \nabla_{\mathbf{p}}C(\mathbf{p})) \\ &= C(\mathbf{p}) + \lambda |\nabla_{\mathbf{p}}C(\mathbf{p})|^2 \\ \end{aligned}$$

ちょっと図示してイメージしてみます。

勾配 $\nabla_{\mathbf{p}}C(\mathbf{p})$ の大きさが何かはよく分かりませんが、 未知数スカラー$\lambda$ 倍でどうとでも調整できますね。そして、初期距離からのズレである $C(\mathbf{p})$ を打ち消すような値になるような $\lambda$ であればいいわけです。未知数 $\lambda$ が1つだけの方程式だったので、式変形だけで $\lambda$ は求めれます。

$$\begin{aligned} C(\mathbf{p} + \Delta \mathbf{p}) &= 0 \\ C(\mathbf{p}) + \lambda |\nabla_{\mathbf{p}}C(\mathbf{p})|^2 &= 0 \\ \lambda |\nabla_{\mathbf{p}}C(\mathbf{p})|^2 &= - C(\mathbf{p})\\ \lambda &= - \frac{C(\mathbf{p})}{|\nabla_{\mathbf{p}}C(\mathbf{p})|^2} \end{aligned}$$

よって、式$(4)$にこれを代入し、本来求めたかった未知の修正量 $\Delta \mathbf{p}$ は次式で表されます。

$$\begin{aligned} \Delta \mathbf{p} &= - \frac{C(\mathbf{p})}{|\nabla_{\mathbf{p}}C(\mathbf{p})|^2} \nabla_{\mathbf{p}}C(\mathbf{p}) \tag{5} \end{aligned}$$

なんかいつのまにか要素の少ない式で表現されましたね…。
式$(3)$にて近似式で内積を使って式を組み立てた部分が、数学センスがないのかピンと来ていません…。ただ内積を使うというのがうまそうですね。任意のベクトル $\Delta \mathbf{p}$ に対して、修正方向ベクトル(勾配)との内積を取ることで、修正方向の成分に関する値だけに着目した$\lambda$ で修正量を表現しつつ、勾配の具体的な値については無視することができるようになっています。

修正量 $\Delta \mathbf{p}$ を式で表すことができました。 $\mathbf{p}$は制約する頂点をまとめたベクトルでした。個々の頂点の修正量について見ていきます。個々の頂点 $\mathbf{p_i}$ の修正ベクトル $\Delta \mathbf{p_i}$ について記述すると以下のようになります。

$$\Delta \mathbf{p_i} = -s \nabla_{\mathbf{p_i}}C(\mathbf{p_i}, \cdots, \mathbf{p_n}) \tag{6}$$

$$s = \frac{C(\mathbf{p_i}, \cdots, \mathbf{p_n})}{\sum_{j} |\nabla_{\mathbf{p_j}} C(\mathbf{p_i}, \cdots, \mathbf{p_n})|^2} \tag{7}$$

なお、$s$ は個々の頂点 $\mathbf{p_i}$ に依存しない式になっているので、別変数に置き換えて式が見やすくなるようにされています。$s$ の式に急に $\sum$ が出てきましたね。もう少し丁寧に見ていきましょう。

$\mathbf{p}$ は 拘束 $C$ に関する頂点をまとめた $\mathbf{p} = [ \mathbf{p_1}^T, \cdots, \mathbf{p_n}^T]^T$ でした。 式 $(5)$をこの表記で見てみます。

$$\Delta \begin{bmatrix} \mathbf{p_1} \\ \vdots \\ \mathbf{p_i} \\ \vdots \\ \mathbf{p_n} \\ \end{bmatrix} = -s \begin{bmatrix} \nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \cdots, \mathbf{p_n}) \\ \vdots \\ \nabla_{\mathbf{p_i}} C(\mathbf{p_1}, \cdots, \mathbf{p_n}) \\ \vdots \\ \nabla_{\mathbf{p_n}} C(\mathbf{p_1}, \cdots, \mathbf{p_n}) \\ \end{bmatrix}$$

$\nabla_{\mathbf{p}} C(\mathbf{p})$ の $\nabla_{\mathbf{p}}$ が、個々の頂点の式にしたときに、どのように表記されるかちょっと自信なかったですが、真ん中がそのまま式 $(6)$ になっており、この表記の仕方で合っていそうですね。

次に、式 $(6)$もこの表記で表して、どこで $\sum$ がでてくるのか見ていきます。分子はそのままなので、 $\sum$ のある分母だけ見ていきましょう。

$$\begin{aligned} |\nabla_{\mathbf{p}}C(\mathbf{p})|^2 &= \begin{vmatrix} \begin{bmatrix} \nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \cdots, \mathbf{p_n}) \\ \vdots \\ \nabla_{\mathbf{p_n}} C(\mathbf{p_1}, \cdots, \mathbf{p_n}) \\ \end{bmatrix} \end{vmatrix}^2 \\ &= \begin{pmatrix} \sqrt{ |\nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \cdots, \mathbf{p_n})|^2 + \cdots + |\nabla_{\mathbf{p_n}} C(\mathbf{p_1}, \cdots, \mathbf{p_n})|^2 } \end{pmatrix} ^2 \\ &= |\nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \cdots, \mathbf{p_n})|^2 + \cdots + |\nabla_{\mathbf{p_n}} C(\mathbf{p_1}, \cdots, \mathbf{p_n})|^2 \\ &= \sum_{\alpha = 1}^{n} |\nabla_{\mathbf{p_\alpha}} C(\mathbf{p_1}, \cdots, \mathbf{p_n})|^2 \end{aligned}$$

$\nabla_{\mathbf{p}}C(\mathbf{p})$ は、先ほどのように各頂点についてベクトル表記でまとめたものでした。ベクトルですので、ノルムを計算すると平方根が出てきます。出てきた平方根は2乗で打ち消されます。残るはベクトルの各要素を二乗したものの総和です。よって、 $\sum$ を用いて、 式 $(7)$ の分母の形に変形できました。

表記について

$p_j$ という表記が式$(7)$でしか出ておらずちょっとよく分かりませんでした。 $j \in [1, \cdots, M]$ で既に $j$ は制約を表す記号で使われていますし、和の数は制約 $j$ で拘束している頂点数 $n_j$ (制約頂点集合の濃度)のはずです。ちょっと $p_j$ のどういう表記の仕方が分からなかったので、上記の式変形過程で $\alpha$ という適当な変数を用いているのはご了承ください。(制約$j$の要素$[1, \cdots, n]$の集合の和という意味で $\sum_{j} \mathbf{p_j}$と表記しているのでしょうか？でも $j$ は集合ではなく制約としているのでよく分かりませんでした)

あとは、質量によって修正量の重み付けをしています。重いほうが抵抗があり、移動しにくくなるので直感にもあっています。

やっと終わりです！！各頂点 $\mathbf{p_i}$ における、制約を満たすために移動させる修正ベクトル $\Delta \mathbf{p_i}$ の式を得ることができました。

最後に具体例が書かれているので見てみましょう。

◆具体例

2点間の初期距離を維持するような、距離拘束関数 $C(\mathbf{p_1}, \mathbf{p_2}) = | \mathbf{p_1} - \mathbf{p_2} | -d$ を考えます。この距離拘束関数 $C(\mathbf{p_1}, \mathbf{p_2})$ の勾配は次のようになります。

$$\nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \mathbf{p_2})= \frac{\mathbf{p_1} - \mathbf{p_2}}{|\mathbf{p_1} - \mathbf{p_2}|}$$

これはベクトル $\overrightarrow{\mathbf{p_{2}p_{1}}}$ を正規化したベクトルです。正規化とは長さが1となるようなベクトルにすることで、向きのみを表すような方向ベクトルになります。

勾配の大きさなんて分からないのに、なぜ方向ベクトルと置くことができるのでしょうか？説明がなかったので最初はちょっと悩んだのですが、おそらく下記の理由だと思います。

勾配 $\nabla_{\mathbf{p}}C(\mathbf{p})$ の大きさが何かはよく分かりませんが、 未知数スカラー$\lambda$ 倍でどうとでも調整できますね。

という話がありました。

大事なのは 勾配の方向 。修正のために動かす向きです。勾配が方向ベクトルより多い少ないがあっても、方向が同じで単なるスカラー倍の違いであれば、それは未知数 $\lambda$ に含めてしまっても同じです。なので、 $\mathbf{p_2}$ から $\mathbf{p_1}$ への 方向ベクトル を勾配 $\nabla_{\mathbf{p_1}}C(\mathbf{p})$ として考えてしまおうということかと思われます。

加えて、そう考えると式も少し簡単になり便利になります。長さが $1$ ですから、式中に出てくる $|\nabla_{\mathbf{p_1}} C(\mathbf{p_1}, \mathbf{p_2})|^2$ は $1$ になり、消すことができます。

$$\begin{aligned} s &= \frac{C(\mathbf{p_1},\mathbf{p_2})}{w_1|\nabla_{\mathbf{p_1}} C(\mathbf{p_1},\mathbf{p_2})|^2 + w_2|\nabla_{\mathbf{p_2}} C(\mathbf{p_1},\mathbf{p_2})|^2} \\ &= \frac{| \mathbf{p_1} - \mathbf{p_2} | -d}{w_1 + w_2} \\ \end{aligned}$$

勾配が消えて式が単純になりました。

$$\Delta \mathbf{p_1} = - s w_1 \frac{\mathbf{p_1} - \mathbf{p_2}}{|\mathbf{p_1} - \mathbf{p_2}|} \tag{10}$$

$$\Delta \mathbf{p_2} = - s w_2 \frac{\mathbf{p_1} - \mathbf{p_2}}{|\mathbf{p_1} - \mathbf{p_2}|} \tag{10}$$

すごく簡単な式ですね！もう実装に入れそうです！試しに実装してみます。2点間の距離束縛の物理といえば紐が実装できるはずです。

◆剛性(stiffness)について

その前に、剛性についても記述があるので読んでおきます。拘束を弾性力のポテンシャルエネルギーでモデル化していましたね。バネ定数 $k$ も剛性の一種です。

一番シンプルなのは剛性 $k$ を修正 $\Delta \mathbf{p}$ に掛けることです。しかし、ソルバーのループ内で何度も掛けると非線形な影響を与えてしまいます。

まずは、線形な場合について考えます。

剛性 $k = 1$ であれば固いので、伸びはなく修正はそのまま $\Delta \mathbf{p} \cdot k = \Delta \mathbf{p} \cdot 1$ だけ行われます。この剛性によって修正されなかった距離は残りの$(1-k)$ですので、 $\Delta \mathbf{p} \cdot (1 - k) = \Delta \mathbf{p} \cdot (1 - 1) = 0$ です。つまり、完全に制約を満たす位置まで修正しています。

剛性 $k = 0.3$ であれば柔らかいので、修正がしっかりとされず修正される距離は $\Delta \mathbf{p} \cdot 0.3$ だけです。この剛性によって修正されなかった距離は $\Delta \mathbf{p} \cdot (1 - 0.3) = \Delta \mathbf{p} \cdot 0.7$ です。本来戻る距離の $0.7$ 倍の距離だけまだ残っており、制約を満たす位置まではまだ戻っていません。

このようになっていれば、剛性に対して線形で扱いやすいのですが、ループ内でやるとこの修正されなかった距離が非線形な $\Delta \mathbf{p} \cdot (1 - k)^{n_{loop}}$ になってしまいます。しかし、直接$k$を使わずに、少し加工した $k’ = 1 - (1-k)^{\frac{1}{n_{loop}}}$ をループ内で掛けることで線形な $\Delta \mathbf{p} \cdot (1 - k)$ になるとのことです！

剛性$k$によって弱まり、修正されなかった距離は $\Delta \mathbf{p} \cdot (1 - k)^{n_{loop}}$ になるとのことでしたが、実際に $k’$ をこの式にいれて、そうなるのか計算してみます。

$$\begin{aligned} \Delta \mathbf{p} \cdot (1 - k')^{n_{loop}} &= \Delta \mathbf{p} \cdot (1 - (1 - (1-k)^{\frac{1}{n_{loop}}}))^{n_{loop}} \\ &= \Delta \mathbf{p} \cdot (1 - 1 + (1-k)^{\frac{1}{n_{loop}}}))^{n_{loop}} \\ &= \Delta \mathbf{p} \cdot ((1-k)^{\frac{1}{n_{loop}}})^{n_{loop}} \\ &= \Delta \mathbf{p} \cdot (1-k)^{\frac{n_{loop}}{n_{loop}}} \\ &= \Delta \mathbf{p} \cdot (1-k) \\ \end{aligned}$$

なりました。一応、愚直に $k = 0.3$ 、 $n_{loop} = 10$ を入れて計算して、先ほどの $0.7$ になるか見てみましょう。$0.7$ になりますね。期待する線形な結果が得られました。

実装

簡約のため飛ばした部分もありますが、実装に必要な部分の論文は読み終わりました。提示されたアルゴリズムを愚直に実装すれば大丈夫です。計算の効率は考えず、理解のために極力アルゴリズムや数式のままコードに落とし込んでみます。以下のような感じになりました。

本記事では、アルゴリズムや制約処理の部分のみを載せています。
全体のサンプルコードはUnityProjectをGithubに上げてあります。
https://github.com/ner-develop/PositionBasedDynamics

void Simulate(float dt)
{
	// (5)
	foreach (var (_, vertex) in _vertexes)
	{
		if (vertex.isFixed) { continue; }
		vertex.v = vertex.v + dt * _parameter.gravity;
	}

	// (6) スキップ
	// Damp Velocities

	// (7)
	foreach (var (_, vertex) in _vertexes)
	{
		vertex.p = vertex.x + dt * vertex.v;
	}

	// (8) スキップ
	// generateCollisionConstraints

	// (9) ~ (11)
	float k = 1 - Mathf.Pow(1 - Mathf.Clamp01(_parameter.stiffness), 1f / _parameter.iterations);
	for (int iterationCount = 0; iterationCount < _parameter.iterations; iterationCount++)
	{
		ProjectConstraints(k);
	}

	// (12) ~ (15)
	foreach (var (_, vertex) in _vertexes)
	{
		if (vertex.isFixed) { continue; }
		vertex.v = (vertex.p - vertex.x) / dt;
		vertex.x = vertex.p;
	}

	// (16) スキップ
	// velocityUpdate
}

void ProjectConstraints(float k)
{
	// 拘束
	for (int i = 0; i < _vertexes.Length - 1; i++)
	{
		var (_, vertex1) = _vertexes[i];
		var (_, vertex2) = _vertexes[i + 1];

		var p1_p2 = vertex1.p - vertex2.p;
		var p1_p2_length = p1_p2.magnitude;
		var w1 = vertex1.w;
		var w2 = vertex2.w;

		var gradC = p1_p2 / p1_p2_length;
		var C = Constraint(vertex1, vertex2);
		var s = C / ((w1 + w2) * Vector3.Dot(gradC, gradC));

		var Δp1 = -s * w1 * gradC;
		var Δp2 = +s * w2 * gradC;

		vertex1.p += Δp1 * k;
		vertex2.p += Δp2 * k;
	}
}

float Constraint(Vertex v1, Vertex v2)
{
	float d = (v1.initialPosition - v2.initialPosition).magnitude;
	return (v1.p - v2.p).magnitude - d;
}

Unreal Engine サンプルプロジェクト

UE5でのサンプルプロジェクトも用意しました。C++ですがほとんどコードは同じです。
https://github.com/ner-develop/PositionBasedDynamicsForUE

おわりに

論文に示されたアルゴリズム通りに実装するだけならそのままでいい場合は簡単ですが、ちゃんと示された数式を読み取り、意図を理解するとなるとだいぶ大変だなと思いました。技術ブログを書く上で、説明するために文章を書き始めると理解していない部分が多くあったことに気づきかなり勉強になりました。足りない知識の分野も見えてきたのでその点でも論文を読んで実装してみるというのは非常に学びがあります。物理は大学で専攻していないので理解するのに苦戦しました。

XPBDについても論文を読んで実装してみたのですが、同様に技術ブログにする際には苦戦するかもしれません。勉強になるので書いてみたいなとは思っております。

有識者の方がもし目を通して間違いなどにお気づきになったら、お手数ですが一報頂けると喜びます。大学であれば教授や先輩などに間違えを教えてもらえますが、社会人になってからの独学はその点で大変ですね。