GAMES103_mdbook

基于物理的计算机动画入门

原始课程链接

如有侵权，请联系删除

由于加载速度的限制，一个页面的内容太多会导致图片加载不出来，因此把内容分得非常细，每一个页面涉及一个很小的话题。

主要内容

说明

这个是GAMES-Webinar提供的一个课程系列。
作者王华民老师的课程是抓大放小风格。不会在公式推导的细节上深究，甚至有时有一句“你们相信我没有把公式抄错吧”就算是完成公式的推导了。但老师在关键问题上一点都含糊。喜欢把算法从大的层面进行分类和归纳。算法讲完后一定少不了关于“用不同套路解决同一问题”的比较。在传授了算法的同时，也言传身教地让听众学会了辩证地思考。
特别值得一提的是，老师在最后一课的最后的一段总结讲得很好。是整篇课程的升华。大概是关于以下内容：（1）关于图形学的学术研究。（2）关于自学的困境。（3）关于AI模拟。

本文出自CaterpillarStudyGroup，转载请注明出处。
https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P3

Prerequisites

• The course is designed for
  • Undergraduates in the 3rd or 4th year, or fresh graduates.
• Linear Algebra
  • You Should know basic linear algebra concepts, such as vectors, matrices, linear systems, SVD...
• Calculus
  • You should know how to calculate basic derivatives and integrals; you should understand chain rules, gradients, etc.
• Programming Skills
  • C, or C++, or C#, or Javascript
• Ready to learn by yourself
• The life will be much easier if you took
  • Numerical methods (numerical linear algebra, numerical PDEs), finite element analysis, fluid dynamics...

✅ 建议：读 paper 而不是教材，只读重点不读全文学知识而不是学用 Unity. 多读多写多想

P14

Graphics Pipeline

Real-Time Graphics Pipeline

P15
The number of frames sent to display in a second is called the frame rate.
For example, 24 FPS, 30 FPS, 60 FPS, …

✅ 帧率要求主要取决于交互性，因此游戏要求比电­影高。

P17

Animation Playback

✅ 由于实时比较难，可以把不需要交互的动画，例如过场动画做成离线
✅ 同理，不需要交互的场景。

P18

Movie

✅ Geometry: 离线：构造离线的3D也界
✅ 动画：渲染，实时，需要与3D世界或玩家互动
✅ 电影：离线，不需要交互，提前录下来，例如游戏中的过场动画

P19

Geometry: Three Representations

Mesh

• A mesh contains:
  • Vertices (nodes)
  • Elements (triangles, polygons, tetrahedra…)

✅ polygon 常用于 Maya,四面体常用于软体模拟。

• Triangle mesh is the foundation of graphics.

✅ 三角形不只是 Mesh 的基础，也是渲染的基础

• Problems:
  • Meshing (Delaunay triangulation)
  • Simplification/subdivision
  • Mesh optimization (smoothing, flows…)
  • Volume mesh

✅ 关于Mesh部分可以参考GAMES101
❓ 什么是 flows?
✅ Volume Mesh 的处理比普通 Mesh 要难很多

P20

✅ Structured：有规律的。Unstructured：无规律的
✅ 有些模拟算法或几何算法，可以利用 Structured 做简化或优化。但Unstructured 算法通用性更好。

P21

Point Cloud

• A point cloud is simple.
• It can be raw data from surface scan.
• Problems:
  • Mesh reconstruction from cloud
  • (Re)-Sampling
  • Neighborhood search
  • …

✅ 原始点云可能有疏有密，因些需要重采样。

P22

Volumetric Grid

• A grid partitions the space; a cell stores the physical quantities at that spot.
• Don’t confuse it with structured mesh.
• It’s often acquired from volumetric scan, e.g., CT.
• Problems:
  • Memory cost (octree?)
  • Volumetric rendering?
  • …

P23

Rendering: Non-Photorealistic vs. Photorealistic

✅ Non-Photorealistic：非真实感渲染，Photorealistic：真实感渲染
✅ 后者更主流，又分为基于光线追踪和基于传统渲染管线。见 Games 101

P27

Animation

Character and Physics-Based Animation

✅ 本课程主要是后者，但通常二者是结合的。

P38

Physics-Based Animation

Animation Paradigm

• The goal of animation is to update the state in every time step.
• The state can be:
  • Position/orientation
  • Velocity
  • Appearance
  • Density
  • …
• The time step doesn’t have to match the frame rate.
  • It’s common to animate multiple time steps then render one frame.

P39

Physics-Based Animation Topics

✅ Cloth and Hair：细的窄的一类物体
✅ Soft Bodys：包括软体，弹性体
✅ Fluids：流体，包括液体和气体

P56

✅ 刚体：常用Mesh，因为Mesh 适用于形态固定、不会剧烈拉扯或断裂的物体。用粒子模拟会破碎的刚体
衣服头发：也可以用 Grid 模拟衣服和头发可以减少碰撞处理。但计算量大，且难以处理细节
流体：烟通常使用粒子法或网格法。水波可以看作是整体，因此能用 mesh，用 mesh的好处是可以做到实时，Grid 的好处是更真实。Splashes(水花)的问题是多变，因此不能实时，通常使用粒子和网格。
Hybrid 方法：MPM = Particle + Grid,兼容二者优点，常用于模拟雪或粘滞物体
Coupling: 场景中同时有不同类别的物体，怎样模拟它们的交互。
✅ SPA 与弹性体模拟结合，可用于模拟物体破碎， 粒子法与网格法相结合，称为 MPM. 用于模拟雪、沙子。

P59

Topics in This Class

✅
Fracture 有大量的 remesh。游戏引擎中的 Fracture 通常通过预计算而不是模拟得到。
Rigid 还是Soft,看有没有形变。
Mesh 定义在物体上， Grid 定义在场景上
水波 Mesh 也会讲，图上漏掉了。

P60

My Own Expertise

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P3

Vector: Basics

定义

An (Euclidean) vector: A geometric entity endowed with magnitude and direction.

$$ \mathbf{P} =\begin{bmatrix} p_x\\ p_y\\ p_z\\ \end{bmatrix}\in \mathbf{R} ^3 $$

$$ \mathbf{o} =\begin{bmatrix} 0\\ 0\\ 0\\ \end{bmatrix} $$

The vector p is defined with respect to the origin o.

坐标系

✅ 用黑来区分，矢量：黑体小写；标量：斜体；矩阵：黑体大写；

P4

The choice of a right-hand or left-hand system is largely due to:
the convention of the screen space.

✅ 左手坐标系，E轴正方向朝屏幕内，好处是物体坐标 x、y、z 都是正值。右手系统的物体都在E轴负方向。

P5

Stacked Vector

Vectors can be stacked up to form a high-dimensional vector, commonly used for describing the state of an object.

Not a geometric vector,but a stacked vector.

P6

Vector Arithematic: Addition and Subtraction

$$ \mathbf{p\pm q=} \begin{bmatrix} p_x\pm q_x\\ p_y\pm q_y\\ p_z\pm q_z\\ \end{bmatrix} $$

$$ \mathbf{p+q=q+p} $$

P7

Example 1: Linear Representation

A (geometric) vector can represent a position, a velocity, a force, or a line/ray/segment.

✅ 图2。同一个公式，对\(t\)做不同的约束，可以定义不同的东西。
\(\mathbf{P}(t)\) 是 \(\mathbf{P}\) 和 \(\mathbf{q}\) 的 blend

P8

Vector Norm

A vector norm measures the magnitude of a vector: its length.

✅ L1-Norm 又称为曼哈顿的距离。没写下标一般默认L2-Norm

P9

Vector Norm: Usage

Distance between q and p： $$ \mathbf{||q-p||} $$

A unit vector：

$$ \mathbf{||p||} =1 $$

Normalization： $$ \mathbf{\bar{p} =p/||p||} $$

P10

Vector Arithematic: Dot Product

A dot product, also called inner product, is:

$$ \begin{array}{c} \mathbf{p\cdot q}=p_xq_x+p_yq_y+p_zq_z=\mathbf{p^Tq} \\ =||\mathbf{p} ||||\mathbf{q} ||\cos \theta \end{array} $$

• \(\mathbf{p\cdot q=q\cdot p} \)
• \(\mathbf{p\cdot (q+r)=p\cdot q+p\cdot r} \)
• \(\mathbf{p \cdot p = ||p||^2_2} \), a different way to write norm.
• If \(\mathbf{p·q} = 0\) and \(\mathbf{p,q}\ne 0\) then \(\cos \theta = 0\),then \(\mathbf{p}\) and \(\mathbf{q}\) are orthogonal.

P11

Example 2: Particle-Line Projection

✅\(X\)为物体中心点的位置，为物体上所有点的整体位移，是前面说的\(T\).
速度是加速度的积分，表示为\(V\)或\(\dot{X} \)
加速度为\(F／M\)，但\(F\)比较复杂，与时间、位置、速度都可能有关系。
位置是速度的积分。

P12

Example 3: Plane Representation

S: The signed distance to the plane

Quiz: How to test if a point is within a box?

P13

Example 4: Particle-Sphere Collision

If collision does happen, then:

$$ ||\mathbf p(t) - \mathbf{c}||^2= r^2 $$

$$ (\mathbf p-\mathbf c+t\mathbf v)·(\mathbf p-\mathbf c +t\mathbf v) =r^2 $$

$$ (\mathbf v·\mathbf v)t^2+2(\mathbf p-\mathbf c)·\mathbf vt+ (\mathbf p-\mathbf c)·(\mathbf p-\mathbf c)-r^2=0 $$

• Three possiblities:
  • No root、无碰撞
  • One root、擦边 if \(t > 0\)
  • Two roots:自碰撞 if \(t > 0 \)

P14

Vector Arithematic: Cross Product

The result of a cross product is a vector:

$$ \mathbf{r=p\times q} =\begin{bmatrix} p_yq_z-p_zq_y \\ p_zq_x-p_xq_z\\ p_xq_y-p_yq_x\\ \end{bmatrix} $$

• \(\mathbf r·\mathbf p = 0; \mathbf r·\mathbf q = 0; ||\mathbf r|| = ||\mathbf p||||\mathbf q|| \sin \theta\)
• \(\mathbf p\times \mathbf q =-\mathbf q\times \mathbf p\)
• \(\mathbf p\times (\mathbf q +\mathbf r) = \mathbf p\times \mathbf q +\mathbf p\times \mathbf r\)
• If \( \mathbf p \times \mathbf q =\mathbf 0\) and \(\mathbf p,\mathbf q\ne 0 \) then \(\sin \theta= 0\), then \(\mathbf p\) and \(\mathbf q \) are parallel (in the same or opposite direction).

P15

Example 5: Triangle Normal and Area

• Cross product gives both the normal and the area.
• The normal depends on the triangle index order, also known as topological order.

P16
Quiz: How to test if three points are on the same line (co-linear)?

P17

Example 6: Triangle Inside/Outside Test

P18

Otherwise, outside.

✅ 假设P点在三角形所在平面上
三个点的顺序很重要，不能搞反。

P19

Example 7: Barycentric Coordinates

Note that:

$$ \frac{1}{2} \mathbf{(x_0−p)×(x_1−p)\cdot n} =\begin{cases} \frac{1}{2}||\mathbf{(x_0−p)×(x_1−p)} ||& \mathrm{inside} \\
\frac{1}{2}||\mathbf{(x_0−p)×(x_1−p)} || & \mathrm{outside} \end{cases} $$

Signed areas:

$$ \mathbf{A_2=\frac{1}{2} (x_0−p)×(x_1−p)\cdot n} $$

$$ \mathbf{A_0=\frac{1}{2} (x_1−p)×(x_2−p)\cdot n} $$

$$ \mathbf{A_1=\frac{1}{2} (x_2−p)×(x_0−p)\cdot n} $$

$$ \mathbf{A_0+A_1+A_2=A} $$

Barycentric weights of p :

$$ b_0=A_0/A \quad b_1=A_1/A \quad b_2=A_2/A $$

$$ b_0+b_1+b_2=1 $$

Barycentric Interpolation

$$ \mathbf{p} =b_0\mathbf{x} _0+b_1\mathbf{x} _1+b_2\mathbf{x} _2 $$

✅ 当 \(\mathbf{p}\) 在三角形外面时，面积为负，但面积总和不变 \(b_0,b_1,b_2\) 为 \(\mathbf{p}\) 在三角形重心坐标系下的坐标

✅ \(\mathbf{p}\) 在三角形外部、重心坐标同样适用，不过权重有负数。

P20

Gouraud Shading

• Barycentric weights allows the interior points of a triangle to be interpolated.

• In a traditional graphics pipeline, pixel colors are calculated at triangle vertices first, and then interpolated within. This is known as Gouraud shading.

• It is hardware accelerated.

• It is no longer popular.

✅ 由于硬件能力提升，已经可以做到逐像素。
shading,不再需要此方法
通常也不是逐像素计算重心坐标，而是扫描线算法
例如要计算某一行，可以 ：
(1) 插值出行起点像素；
(2) 插值出行终点像素；
(3) 起点与终点间批量插值；

P21

Example 9: Tetrahedral Volume

Edge vectors:

$$ \mathbf{X_{10}=X_1-X_0 \quad \quad X_{20}=X_2-X_0 \quad \quad X_{30}=X_3-X_0} $$

Base triangle area:

$$ A=\frac{1}{2} ||\mathbf{X} _{10}\times \mathbf{X} _{20}|| $$

Height:
$$ h=\mathbf{x} _{30}\cdot\mathbf{n} =\mathbf{x} _{30}\cdot \frac{\mathbf{x} _{10}\times \mathbf{x} _{20}}{||\mathbf{x} _{10}\times \mathbf{x} _{20}||} $$

Volume:

$$ \begin{align*} V&=\frac{1}{3} ℎA=\frac{1}{6} \mathbf{x} _{30}\cdot \mathbf{x} _{10}\times \mathbf{x} _{20}\\ &=\frac{1}{6}\begin{vmatrix} \mathbf{x} _1 & \mathbf{x} _2 & \mathbf{x} _3 &\mathbf{x} _0 \\ 1& 1 & 1 &1 \end{vmatrix} \end{align*} $$

✅ 四面体
\(h\)是\(\mathbf{x}_{30}\)在 normal 上的投影
行列式是上面叉乘的另一种马法。

P22

Note that the volume \(V =\frac{1}{3}h\mathit{A} =\frac{1}{6} \mathbf{x} _ {30}\cdot (\mathbf{x} _ {10}\times \mathbf{x}_{20})\) is signed.

✅ \(\mathbf{x}_{3}\)的后面法线的同方向上，也正四面体，反之为负四面体，四点共面为零体积。

P23

Example 10: Barycentric Weights (cont.)

• p splits the tetrahedron into four sub-tetrahedra:

$$ \begin{matrix} V_0=\mathrm{Vol} (\mathbf{x}_3,\mathbf{x}_2, \mathbf{x}_1, \mathbf{p} )\\ V_1=\mathrm{Vol} (\mathbf{x}_2,\mathbf{x}_3, \mathbf{x}_0, \mathbf{p} )\\ V_2=\mathrm{Vol} (\mathbf{x}_1,\mathbf{x}_0, \mathbf{x}_3, \mathbf{p} )\\ V_3=\mathrm{Vol} (\mathbf{x}_0,\mathbf{x}_1, \mathbf{x}_2, \mathbf{p} ) \end{matrix} $$

• p is inside if and only if: \(V_0,V_1,V_2, V_3 > 0\).

• Barycentric weights:
  $$ b_0=V_0/V \quad b_1=V_1/V \quad b_2=V_2/V \quad b_3=V_3/V $$

$$ b_0+b_1+b_2+b_3=1 $$

$$ \mathbf{p} =b_0\mathbf{x} _0+b_1\mathbf{x} _1+b_2\mathbf{x} _2+b_3\mathbf{x} _3 $$

P24

Example 11: Particle-triangle Intersection

• First, we find t when the particle hits the plane:

$$ (\mathbf{p} (t)−\mathbf{x} _0)\cdot \mathbf{x} _{10}\times \mathbf{x} _{20}=0 $$

$$ (\mathbf{p}-\mathbf{x} _0+t\mathbf{v})\cdot \mathbf{x} _{10}\times \mathbf{x} _{20}=0 $$

$$ t=\frac{(\mathbf{p}−\mathbf{x}_0)\cdot  \mathbf{x} _{10}\times \mathbf{x} _{20}}{\mathbf{v}\cdot \mathbf{x} _{10}\times \mathbf{x} _{20}} $$

• We then check if \(\mathbf{p}(t)\) is inside or not.
  • See Example 6.

✅ 代入体积公式，体积为0时发生碰撞

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P26

Matrix: Basics

Matrix: Definition

A real matrix is a set of real elements arranged in rows and columns.

$$ A=\begin{bmatrix} a_{00} & a_{01} & a_{02} \\ a_{10}& a_{11} & a_{12} \\ a_{20}& a_{21} & a_{22} \end{bmatrix}=[a_{0} \quad a_{1} \quad a_{2}]\in \mathbf{R} ^{3\times 3} $$

$$ \mathbf{A^T=A} \quad \mathrm{Symmetric} $$

P27

Matrix: Multiplication

How to do matrix-vector and matrix-matrix multiplication? (Omitted)

• \(\mathbf{AB≠BA} \quad \quad \quad \quad \quad \quad \quad \quad \mathbf{(AB)x=A(Bx)} \)
• \(\mathbf{(AB)^T=B^TA^T} \quad \quad \quad \quad \quad \quad \mathbf{(A^TA)^T=A^TA}\)
• \(\mathbf{Ix=x} \quad \quad \quad \quad \quad \quad \quad \quad \quad \mathbf{AI=IA=A}\)
  \(\quad\)
• \(\mathbf{A^{−1}:  AA^{−1}=A^{−1}A=I} \quad \quad \mathrm{inverse}\)
• \(\mathbf{(AB)^{−1}=B^{−1}A^{−1}}\)
• Not every matrix is invertible, e.g., \(\mathbf{A} =\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}\)

P28

Matrix: Orthogonality

An orthogonal matrix is a matrix made of orthogonal unit vectors.

$$ \mathbf{A} =[\mathbf{a} _0\quad \mathbf{a} _1\quad \mathbf{a} _2]\quad\mathrm{such \quad that } \quad \mathbf{a}_i^\mathbf{T}\mathbf{a}_j =\begin{cases} 1，& \text{ if } i= j \text{(unit)}\\ 0.& \text{ if } i\ne j \text{(orthogonal)} \end{cases} $$

$$ \mathbf{A^TA}=\begin{bmatrix} \mathbf{a}_0^\mathbf{T} \\ \mathbf{a}_1^\mathbf{T} \\ \mathbf{a}_2^\mathbf{T} \end{bmatrix}\begin{bmatrix} \mathbf{a}_0 & \mathbf{a}_1 &\mathbf{a}_2 \end{bmatrix}=\begin{bmatrix} \mathbf{a}_0^\mathbf{T} \mathbf{a}_0 & \mathbf{a}_0^\mathbf{T} \mathbf{a}_1 & \mathbf{a}_0^\mathbf{T} \mathbf{a}_2\\ \mathbf{a}_1^\mathbf{T} \mathbf{a}_0 & \mathbf{a}_1^\mathbf{T} \mathbf{a}_1 & \mathbf{a}_1^\mathbf{T} \mathbf{a}_2\\ \mathbf{a}_2^\mathbf{T} \mathbf{a}_0 & \mathbf{a}_2^\mathbf{T} \mathbf{a}_1 & \mathbf{a}_2^\mathbf{T} \mathbf{a}_2 \end{bmatrix}=I $$

$$ \mathbf{A^T=A^{-1}} $$

P29

Matrix Transformation

A rotation can be represented by an orthogonal matrix.

✅ \(\mathbf{x、y、z}\) 是世界坐标系、 \(\mathbf{u、v、w}\) 是局部坐标系，旋转矩阵是局部坐标系在世界坐标系中的状态的描述。

P30
A scaling can be represented by a diagonal matrix.

P31

矩阵分解

Singular Value Decomposition

A matrix can be decomposed into:
\(\mathbf{A=UDV^T} \quad\)such that \(\mathbf {D}\) is diagonal,and \(\mathbf {U}\) and \(\mathbf {V}\) are orthogonal.
\(\quad \quad \quad \quad\quad\) D 的对角线元素是singular values（奇异值）

Any linear deformation can be decomposed into three steps: rotation, scaling and rotation:

✅ rotation \(\longrightarrow\) scaling \(\longrightarrow\) rotation 分别对应 \(\mathbf{V}_2^\mathbf{T},\mathbf{D}, \mathbf{U}\). 注意顺序！！！
所有 \(\mathbf{A}\) 都能做 \(\mathbf{SVD} \)

P32

Eigenvalue Decomposition

A symmetric matrix can be decomposed into:
\(\mathbf{A=UDU^{-1}}\quad\)such that \(\mathbf {D}\) is diagonal,and \(\mathbf {U}\) is orthogonal.
\(\quad \quad \quad \quad\quad\) D 的对角线元素是eigenvalues（特征值）

✅ \(\mathbf{ED}\) 看作是\(\mathbf{SVD}\)的特例，仅应用于对称矩阵，此时 \(\mathbf{U=V}\)
\(\mathbf{U}\) 是正交矩阵，因此也可写成 \(\mathbf{A = UVU^T}\)

As in the textbook
Let \(\mathbf{U} =\begin{bmatrix} \cdots & \mathbf{u} _i &\cdots \end{bmatrix}\), we have:
$$ \mathbf{Au} _i= \mathbf{UDU^T} \mathbf{u} _i=\mathbf{UD} \begin{bmatrix} \vdots \\ 0\\ 1\\ 0\\ \vdots \end{bmatrix}=\mathbf{U} \begin{bmatrix} \vdots \\ 0\\ d_i\\ 0\\ \vdots \end{bmatrix}=d_i\mathbf{u} _i $$ \(\mathbf{U}\): 是 the eigenvector of \(d_i\)
\(d_i\): 是 eigenualue

We can apply eigenvalue decomposition to asymmetric matrices too, if we allow eigenvalues and eigenvectors to be complex. Not considered here.

✅ complex：复数
图形学不考虑虚数，因此也不考虑非对称矩阵的 \(\mathbf{ED}\)

P33

Symmetric Positive Definiteness (s.p.d.)

定义

\(\mathbf{A}\) is s.p.d. if only if: \(\quad\quad\quad\quad\quad\quad\quad\quad \) \(\mathbf{v^TAv}>0\), for any \(\mathbf{v} ≠ 0. \)

\(\mathbf{A}\) is symmetric semi-definite if only if: \(\quad\quad \) \(\mathbf{v^TAv}≥0\), for any \(\mathbf{v}≠ 0\).

✅ 计算矩阵的有限元或 Hession 时会用到正定性

怎么理解SPD

\(d>0 \quad\quad\quad\quad\Leftrightarrow \quad \mathbf{v^T} d\mathbf{v} >0\), for any \(\mathbf{v} ≠ 0. \)

\(d_0, d_1,…>0 \quad\Leftrightarrow \quad \mathbf{v^TDv=v^T} \begin{bmatrix} \ddots & \Box & \Box\\ \Box & d_i & \Box\\ \Box &\Box &\ddots \end{bmatrix}\mathbf{v} >0\), for any \(\mathbf{v} ≠0.\)

✅ 一堆大于零的实数组成一个对角矩阵, 公式1的扩展

\(d_0, d_1,…>0  \quad\Leftrightarrow \quad \mathbf{v^T(UDU^T)v=v^TUU^T(UDU^T)UU^Tv}\)

\(\mathbf{U}\) orthogonal \(\quad\quad\quad\quad\quad\quad\quad\quad=\mathbf{(U^Tv)^T(D)(U^Tv)>0 } \), for any \(\mathbf{v} ≠0 \)

✅ 公式3是公式2的扩展

P34

怎么判断SPD

• A is s.p.d. if only if all of its eigenvalues are positive:
  \(\mathbf{A=UDU^T}\) and \(d_o,d_1,\cdots > 0.\)

• But eigenvalue decomposition is a stupid idea most of the time, since it takes lotsof time to compute.

✅ 实际上不会通过 \(\mathbf{ED}\) 来判断矩阵的正定性。因为ED的计算量很大。

• In practice, people often choose other ways to check if A is sp.d. For example,

\(a_{ii}>∑_{i≠j}|a_{ij}|\) for all \(i\)
A diagonally dominant matrix is p.d.

$$ \begin{bmatrix} 4&3 & 0\\ -1& 5 &3 \\ -8& 0 &9 \end{bmatrix}\begin{matrix}\quad\quad \quad4>3+0\\ \quad\quad\quad 5>1+3 \\ \quad\quad9>8 \end{matrix} $$

✅ 对角占优矩阵必定正定，正定不一定对角占优

• Finally, a s.p.d.matrix must be invertible:
  $$ \mathbf{A^{-1} =(U^T)^{-1}D^{-1}U^{-1} = UD^{-1}U^T}. $$

P35

例子

Prove that if A is s.p.d., then \(\mathbf{B} =\begin{bmatrix} \mathbf{A} &\mathbf{-A} \\ \mathbf{-A} &\mathbf{A} \end{bmatrix}\)is symmetric semi-definite.

For any \( \mathbf{x}\) and \(\mathbf{y}\), we know:

$$ \begin{bmatrix} \mathbf{ x^T}&\mathbf{ y^T} \end{bmatrix}\mathbf{B}\begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix}=\begin{bmatrix} \mathbf{ x^T}&\mathbf{ y^T} \end{bmatrix}\begin{bmatrix} \mathbf{A} &\mathbf{-A} \\ \mathbf{-A} &\mathbf{A} \end{bmatrix}\begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix} $$

$$ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\mathbf{=x^TA(x-y)-y^TA(x-y)=(x-y)^TA(x-y)} $$

Since A is sp.d., we must have:

$$ \begin{bmatrix} \mathbf{ x^T} & \mathbf{y^T} \end{bmatrix}\mathbf{B} \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix}\ge 0 $$

P36

Linear Solver

Many numerical problems are ended up with solving a linear system:

It's expensive to compute \(\mathbf{A^{-1}} \), especially if \(\mathbf{A} \) is large and sparse. So we cannot simply do:\(\mathbf{x = A^{-1}b}\).

There are two popular linear solver approaches: direct and iterative.

✅ 当 \(\mathbf{A}\) 是稀疏时. \(\mathbf{A}^{-1}\)通常不是稀疏。 如果 \(\mathbf{A}\) 很大， \(\mathbf{A}^{-1}\)会占用大量空间

P37

Direct Linear Solver

方法

A direct solver is typically based LU factorization, or its variant: Cholesky, \(\mathrm{LDL^\top } \), etc…

✅ \(\mathbf{LU}\) 可用于非对称矩阵。
Cholesky 和 \( \mathbf{LDL^\top}\) 仅用于对称矩阵，但内存消耗更少。
这里不介绍如何做\(\mathbf{LU}\)分解

$$ \mathbf{A=LU=} \begin{bmatrix} l_{00} & \Box & \Box \\ l_{10} & l_{11} & \Box \\ \vdots & \cdots &\ddots \end{bmatrix}\begin{bmatrix} \ddots & \cdots &\vdots \\ \Box&u_{n−1,n−1} &u_{n−1,n} \\ \Box & \Box &u_{n,n} \end{bmatrix} $$ \(\quad\quad\quad\quad\quad\quad\quad\)lower triangular \(\quad\quad\) upper triangular

P38

分析

• When \(\mathbf{A}\) is sparse, \(\mathbf{L}\) and \(\mathbf{U}\) are not so sparse. Their sparsity depends on the permutation.(See matlab)

✅ \(\mathbf{L}、\mathbf{U}\) 和稀疏性与行列顺序有关，因此通常在\(\mathbf{LU}\) 分解之前做 permutation,使得到比较好的顺序。

• lt contains two steps: factorization and solving. lf we must solve many linear systems with the same \(\mathbf{A}\) , we can factorize it only once.

✅ \(\mathbf{LU}\) 分解是计算量的大头，只做一次 \(\mathbf{LU}\) 分解，能省去大量计算。

• Cannot be easily parallelized:Intel MKL PARDISO

P39

Iterative Linear Solver

An iterative solver has the form:

Why does it work?

$$ \begin{matrix} \mathbf{b−Ax} ^{[k+1]} =\mathbf{b−Ax} ^{[k]}−\mathbf{αAM} ^{−1}(\mathbf{b−Ax} ^{[k]}) \\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad=(\mathbf{I−αAM} ^{−1})(\mathbf{b−Ax} ^{[k]}) =(\mathbf{I−αAM} ^{−1})^{k+1}(\mathbf{b−Ax} ^ {[0]}) \end{matrix} $$

So,

\(\mathbf{b−Ax} ^{[k+1]}→0\), if \(ρ(\mathbf{I−αAM} ^{−1})<1.\)

✅\(\mathbf{b-Ax}^{[k＋1]}\) 代表下一时的残差，迭代要想收敛，\(\mathbf{b-Ax}^{[k+1]}\) 应趋于0

\(\rho\):矩阵的spectral radius (the largest absolute value of the eigenvalues)

✅ 不会真的去算 \(\rho\),而是调\(α\),试错。 因为求特征值的代价比较大

P40
\(\mathbf{M}\) must be easier to solve:

\(\quad\)

The convergence can be accelerated: Chebyshev, Conjugate Gradient, … (Omitted here.)

优点：

• simple
• fast for inexact solution
• paralleable

缺点：

• convergence condition

✅ 例如要求M是正定的或对角占优的

• slow for exact solution

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P42

Basic Concepts

1st-Order Derivatives

值是实数，变量是矢量

If \(f(\mathbf{x} )\in \mathbf{R} \), then \(df=\frac{∂f}{∂x}dx+\frac{∂f}{∂y}dy+\frac{∂f}{∂z}dz=\begin{bmatrix} \frac{∂f}{∂x} & \frac{∂f}{∂y} &\frac{∂f}{∂z} \end{bmatrix}\begin{bmatrix} dx \\ dy\\ dz \end{bmatrix}\).

$$ \frac{∂f}{∂x}=\begin{bmatrix} \frac{∂f}{∂x} & \frac{∂f}{∂y} &\frac{∂f}{∂z} \end{bmatrix}
$$

$$ \mathrm{ or } $$

✅ \(\nabla f(x)=(\frac{\partial f}{\partial x} )^T\), 重要！！！

Gradient is the steepest direction for increasing \(f\). It’s perpendicular to the isosurface.

✅ isosurface：等高面

P43

值是矢量，变量是是矢量

If \(f(\mathbf{x} )=\begin{bmatrix} f(\mathbf{x} ) \\ g(\mathbf{x} )\\ h(\mathbf{x} ) \end{bmatrix}\in \mathbf{R} ^3\),then:

✅ Divergence:散度，也是\(\mathbf{J}(\mathbf{x})\)的 trace
✅ Curl：旋度。
怎么理解 curl?把微分算子\(\nabla \)看作是个向量，让它与 \(\mathbf{f}\) 做叉乘、在流体模拟中常用。

P44

2nd-Order Derivatives

If \(f\mathbf{(x)\in R} \),then:

✅ 求导顺序不影响求导结果，因此 \(\mathbf{H}\) 是对称的
\(\mathbf{H}\)的trace称为Laplace

泰勒展开

①\(x\in R,f(x)\in R\)
$$ f(x)=f(x_0)+{f}' (x_0)(x-x_0)+\frac{1}{2} {f}'' (x_0)(x-x_0)^2+\cdots $$

②\(x\in R^n,f(x)\in R\)

$$ f(x)=f(x_0)+\rhd {f}' (x_0)\cdot (x-x_0)+\frac{1}{2}(x-x_0)^TH(x-x_0)+\cdots $$

✅ 当\(\mathbf{H}\)正定时, \(f(\mathbf{x})\)满足一些特殊的性质

P45

Quiz:

\(\frac{∂||\mathbf{x}||}{∂\mathbf{x}} = ?\)

$$ \frac{∂||\mathbf{x}||}{∂\mathbf{x} } = \frac{∂(\mathbf{\mathbf{x^Tx} } )^{1/2}}{∂\mathbf{x} }=\frac{1}{2}(\mathbf{x^{T}x} )^{−1/2} \frac{∂(\mathbf{x^Tx} )}{∂\mathbf{x} }=\frac{1}{2||\mathbf{x} ||}2\mathbf{x^T} =\frac{\mathbf{x^T} }{||\mathbf{x} ||} $$

✅ 向量梯度的物理意义：向量沿什么方向变化能最快地变短/长？答：沿它自己的当前方向。

P46

Spring Example

A Spring

🔍 Choi and Ko. 2002. Stable But Responive Cloth. TOG (SIGGRAPH)

✅ Energy：物理上的弹性势能
Force：物理上的力，是 Energy 的 gradient 的反方向;
公式后面有个 T,来源于前面的\(\nabla \)，
直观解释，前面是力的大小，后面是力的方向，推荐论文为以\(\bot\)公式推导的详细过程

P47

A Spring with Two Ends

✅ \(\nabla_0\) 代表对\(\mathbf{x}_0\)的求导

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

Rigid Bodies

Our living environment is stuffed with rigid objects.

✅ rigid：物体很硬，因此不考虑形变。

P6

Rigid Body Simulation

The goal of simulation is to update the state variable \(\mathbf{s} ^{[k]}\) over time.

P7

Rigid Body Motion

If a rigid body cannot deform, its motion consists of two parts: translation and rotation.

✅reference：参考状态,无平移,无旋转此时物体中心点在原点上，\(\mathbf{x}\) 轴向左，\(\mathbf{y}\) 轴用上\(\mathbf{z}\) 轴向前。
当前状态：旋转为\(\mathbf{R}\)，平移为\(\mathbf{T}\). 那么物体上任意点的位置为：
\(\mathbf{{x}}' = \mathbf{Rx} + \mathbf{T}\)
刚体模拟：已知当前时刻的状态和受到的力，求下一时刻的状态。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

积分

For translational motion, the state variablecontains the position \(\mathbf{x}\) and the velocity \(\mathbf{v}\).

$$ \begin{cases} \mathbf{v} (t^{[1]})=\mathbf{v} (t^{[0]})+\mathbf{M} ^{−1}\int_{t^{[0]}}^{t^{[1]}} \mathbf{f} (\mathbf{x} (t), \mathbf{v} (t), t)dt\\ \mathbf{x} (t^{[1]})=\mathbf{x} (t^{[0]})+\int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t)dt \end{cases} $$

✅ 也可以用\(\mathbf{\dot{x}} \)表示速度\(\mathbf{v} \)
速度是加速度的积分，因此\( \Delta t=\int a=\int \frac{F}{M} =M^{-1}\int F\).
位置是速度的积分
本质上是解积分
💡 积分的过程比较独立，单独放在最后，避免破坏整体的结构性。最后结论是混合式的积分方法。

P17

Types of Forces

✅ 在做模拟时，如果不要求能量守衡，出于问题简化的目的，直接对速度做衰减，代替引入阻力

P18

Rigid Body Simulation Pipeline (Translation Only)

The mass \(M\) and the time step \(\Delta t\) are user-specified variables.

✅ 实际应用中，\(\Delta t\) 要跟帧率匹配
质量 \(M\) 可以是个对角矩阵或实数

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

💡 旋转的表示比较独立，单独放在最后，避免破坏整体的结构性。最后结论是混合式的积分方法。

P27

符号定义

角度

Now we choose quaternion \(\mathbf{q}\) to represent theorientation, i.e., the rotation from the reference to the current.

角速度

We use a 3D vector \(\mathbf{\omega}\) to denote angular velocity.

$$ \begin{cases} \text{The direction of } \mathbf{\omega} \text{ is the axis.} \\
\text{The magnitude of } \mathbf{\omega} \text{ is the speed.}
\end{cases} $$

P3

力矩 torque

A torque is the rotational equivalent of a force. It describes the rotational tendency caused by a force.

✅ Torque：力矩，造成物体旋转的趋势。类比于Force：力，造成物体运动的趋势。

✅ \(\mathbf{Rr} _i\)：当前状态下质心到作用点的向量

\(\mathbf{τ} _i\) is perpendicular to both vectors: \(\mathbf{Rr} _i\) and \(\mathbf{f} _i\).

✅ 因此力矩的方向决定了旋转轴的方向，因此由叉差乘得到

\(\mathbf{τ} _i\) is porportional to ||\(\mathbf{Rr} _i\)|| and ||\(\mathbf{f} _i\)||.

\(\mathbf{τ} _i\) is porportional to \(\sin \theta\).
(\(\theta\) is the angle between two vectors.)

✅ 力矩的大小决定旋转的快慢。

P6

inertia tensor

Similar to mass, an inertia tensor describes the resistance to rotational tendency caused by torque. But different from mass, it’s not a constant.

✅ inertia 也与自身的状态相关

Which side receives greater resistance?

✅ 两图的力矩大小相同，但产生的旋转不同
inertia 看作是对运动的抵抗，其效果与力矩的方向有关，因此不是常数

P7

It’s a matrix! The mass inverse is the resistance (just like mass).

✅ 用于旋转的质量不再是实数，而是矩阵，称为 Inertia 矩阵，用 \(\mathbf{I}\) 来标记 Inertia 矩阵，其中 \(\mathbf{I}_{ref}\)为参考状态，\(\mathbf{I}\) 为当前状态，\(\mathbf{I}\) 是 \(3\times 3\) 矩阵。

reference state	current state
\(\mathbf{I} _{\mathbf{ref} }=\sum m_i(\mathbf{r} _i^\mathbf{T} \mathbf{r} _i\mathbf{1} −\mathbf{r} _i\mathbf{r} _i^\mathbf{T} )\) \(\mathbf{1}\) is the 3-by-3 identity.	\(\mathbf{I} =\sum m_i(\mathbf{r} _i^\mathbf{T}\mathbf{R} ^\mathbf{T}\mathbf{Rr} _i\mathbf{1} −\mathbf{Rr} _i\mathbf{r} _i^\mathbf{T} \mathbf{R^T} )\) \(\quad=\sum m_i(\mathbf{Rr} _i^\mathbf{T}\mathbf{r} _i\mathbf{1R} ^\mathbf{T} −\mathbf{Rr} _i\mathbf{r} _i^\mathbf{T} \mathbf{R^T} )\) \(\quad=\sum m_i\mathbf{R}(\mathbf{r}_i^\mathbf{T}\mathbf{r}_i\mathbf{1}−\mathbf{r}_i\mathbf{r}_i^\mathbf{T} ) \mathbf{R^T}\) \(\quad=\mathbf{RI _{ref}R^T}\)

✅ 不需要每次都根据当前状态计算，而是基于一个已经算好的ref状态的 inertia快速得出。

P29

更新法则

	Translational (linear)	Rotational (Angular)
Updafe
states	Velocity \(\mathbf{v}\) Position \(\mathbf{x}\)	Angular velocity \(\mathbf{ω} \) Quaternion \(\mathbf{q}\)
Physical Quantities	Mass \(\mathbf{M}\) Force \(\mathbf{f}\)	Inertia \(\mathbf{I} \) Torque \(\mathbf{τ} \)

✅ 平移： \(加速度 = \frac{力}{质量}\) ，旋转： \(加速度 =\frac{力矩}{\text{Inertia}}\)
✅ \(q\)是四元数，代表物体的旋转状态
✅ \(q_1\times q_2\)不是叉乘，而是四元数普通乘法
✅ \(\begin{bmatrix} 0 & \frac{\bigtriangleup t}{2} & w^{(1)} \end{bmatrix}\)是一个四元数，0为实部，后面为虚部
❗ 算完\(q^{[1]}\)的之后要对它 Normalize
🔎 由\(q^{[0]}\)到\(q^{[1]}\)的更新公式的推导过程见Affer Class Reading（Appendix B）
❓ 四元数相加有什么含义？

P30

Rigid Body Simulation Piplene

In practice, we update the same state variable \(\mathbf{s} =\){\(\mathbf{v,x,\omega ,q}\)} over time.

❗ Gravity doesn't cause any torque! lf your simulator does not contain any other force, there is no need to update \(\mathbf{\omega}\).

P33

After-Class Reading (Before Collision)

P35

https://graphics.pixar.com/pbm2001

❓ 建议读其中的Rigid Body Dynamics部分

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P11

补充1：Integration Methods Explained

Explicit Euler

By definition, the integral \(\mathbf{x} (t) = \int \mathbf{v} (t) dt\) is the area. Many methods estimate the area as a box.

✅ 假设\(\mathbf{x} \)和\(\mathbf{v} \)都是一维的。速度的积分就是阴影区域的面积。

✅ 近似到一阶项，因此称为一阶方法。漏掉的高阶项就是误差。

P12

Implicit Euler

✅ 使用 \(t_0\) 时刻的速度：显式积分
使用 \(t_1\) 时刻的速度：隐式积分
两种方法都只能一阶近似

P13

Mid-Point

P14

比较与混合

By definition, the integral \(\mathbf{x} (t)=∫\mathbf{v} (t) dt\) is the area. Many methods estimate the area as a box.

\(\quad\)

\(\quad\)

P15

$$ \begin{cases} \mathbf{v} (t^{[1]})=\mathbf{v} (t^{[0]})+\mathbf{M} ^{−1}\int_{t^{[0]}}^{t^{[1]}} \mathbf{f} (\mathbf{x} (t), \mathbf{v} (t), t)dt\\ \mathbf{x} (t^{[1]})=\mathbf{x} (t^{[0]})+\int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t)dt \end{cases} $$

✅ 在当前应用场景中，使用前面方法的混合

P16

Leapfrog Integration

✅ 速度和位置是错开的。上下两种写法，在计算上是一样的。

In some literature, such a approach is called semi-implicit.

It has a funnier name: the leapfrog method.

P20

补充2：Rotation Representation

Rotation Represented by Matrix

• The matrix representation is widely used for rotational motion.

• It’s friendly for applying rotation to each vertex (by matrix-vector multiplication).

• But it is not suitable for dynamics:

  • It has too much redundancy: 9 elements but only 3 DoFs.
  • It is non-intuitive.
  • Defining its time derivative (rotational velocity) is also difficult.

P21

Rotation Represented by Euler Angles

• The Euler Angles representation is also popular, often in design and control.

• It is intuitive. It uses three axial rotations to represent one general rotation. Each axial rotation uses an angle.

• In Unity, the order is rotation-by-Z, rotation-by-X, then rotation-by-Y.

• But it is not suitable for dynamics either:

  • It can lose DoFs in certain statuses: gimbal lock.
  • Defining its time derivative (rotational velocity) is difficult.

P22

Gimbal Lock

The alignment of two or more axes results in a loss of rotational DoFs.

✅ 在某些特定的情况下，自由度降低了

P23

Rotation Represented by Quaternion

Introduction

In the complex system, two numbers represent a 2D point.

What about a “complex” system for 3D point? Quaternion! Four numbers represent a 3D point (with multiplication and division).

P24

Quaternion Arithematic

Let \(\mathbf{q} = \begin{bmatrix} \mathbf{s} &\mathbf{v} \end{bmatrix} \) be a quaternion made of two parts: a scalar part \(s\) and a 3D vector part \(\mathbf{v}\), accounting for \(\mathbf{ijk}\).

✅ 在有些库里面写作： \(q = \begin{bmatrix} w & x & y &z \end{bmatrix}\)，w为实数部分

\(a\mathbf{q} =\begin{bmatrix} as &a\mathbf{v} \end{bmatrix}\quad\) Scalar-quaternion Multiplication

\(\mathbf{q} _1±\mathbf{q} _2 =\begin{bmatrix} \mathbf{s}_1±\mathbf{s}_2 & \mathbf{v} _1 ± \mathbf{v} _2 \end{bmatrix}\quad\quad\) Addition/Subtraction

\(\mathbf{q} _1×\mathbf{q} _2= \begin{bmatrix} \mathbf{s} _1\mathbf{s} _2−\mathbf{v} _1\cdot \mathbf{v} _2 & \mathbf{s} _1\mathbf{v} _2+\mathbf{s} _2\mathbf{v} _1+\mathbf{v} _1×\mathbf{v} _2 \end{bmatrix}\quad\quad\) Multiplication

\(||\mathbf{q} ||=\sqrt{\mathbf{s^2+v\cdot v} } \quad\quad\)Magnitude

\(\quad\)

P25

Rotation Represented by Quaternion

• To represent a rotation around \(\mathbf{v}\) by angle \(0\), we set the quaternion as:

  

• lt's very intuitive. lt's the built-in representation in Unity.

• Convertible to the matrix:

$$ \mathbf{R}=\begin{bmatrix} s^2+x^2-y^2-z^2 & 2(xy-sz) & 2(xz+sy)\\ 2(xy+sz) & s^2-x^2+y^2-z^2 & 2(yz-sx) \\ 2(xz-sy) & 2(yz+sx) & s^2-x^2-y^2+z^2
\end{bmatrix} $$

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P9

Topics for the Day

• Particle Collision Detection and Response

  • Penalty methods
  • Impulse methods

• Rigid Collision Detection and Response by Impulse

• Shape Matching

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P11

Particle Collision Detection --- SDF

Signed Distance Function

A signed distance function \(\phi (\mathbf{x} )\) defines the distance from \(\mathbf{x}\) to a surface with a sign. The sign indicates on which side \(\mathbf{x}\) is located.

P12

Signed Distance Function Examples

✅ 圆柱SDF基于勾股定理，\(\sqrt{\cdot }\) 内第一项为斜边长， 第二项为底边长，得出点到中轴的距离。

P13

Intersection of Signed Distance Functions

If \(\phi _0(\mathbf{x} )<0\) and \(\phi_1(\mathbf{x} )<0\) and \(\phi_2(\mathbf{x} )<0\)
then inside
\(\quad \phi (\mathbf{x} )\)=max \(⁡(\phi_0(\mathbf{x}),\phi_1(\mathbf{x}),\phi_2(\mathbf{x}))\)
Else outside
\(\quad \phi (\mathbf{x})=?\)

P14

Union of Signed Distance Functions

✅ 有时候此公式不成立，例如图中\(\mathbf{x}\) 点

Intuitively, we can consider collision detection with the union of two objects as collision detection with two separate objects.

P15

Particle Collision Response ——Penalty Method

Quadratic Penalty Method

A penalty method applies a penalty force in the next update. When the penalty potential is quadratic, the force is linear.

✅ 力的大小与距离有关，方向为normal
✅ 存在的问题：只有\(\mathbf{x}\) 进入 mesh 内部了，才会有力，但此时穿透的 artifacts 已经产生了。解决方法：使用buffor

P16

Quadratic Penalty Method with a Buffer

A buffer helps lessen the penetration issue. But it cannot strictly prevent penetration, no matter how large \(k\) is.

✅ 存在的问题：
如果 \(k\) 太小，快速的碰撞仍会产生 artifacts
如果 \(k\) 太大，碰撞的反弹过于强烈(overshooting)
解决方法：不用常数 \(k\) ，而是 \(k\) 与距离相关

P17

Log-Barrier Penalty Method

A log-barrier penalty potential ensures that the force can be large enough. But it assumes \(\phi (\mathbf{x} ) < 0\) will never happen!!! To achieve that, it needs to adjust \(\Delta t\).

✅ 用倒数关系代替线性关系。
✅ 存在的问题：
1.当\(\mathbf{x}\) 靠近物体表面时，仍然会 overshooting
2.\(\mathbf{x}\) 穿透表面后，会越陷越深。
3.本算法要求保证穿透永远不会发生，因此要仔细调节 \(\Delta t\).

P18

A Short Summary of Penalty Methods

• The use of step size adjustment is a must.

  • To avoid overshooting.
  • To avoid penetration in log-barrier methods.

• Log-barrier method can be limited within a buffer as well.

  • Li et al. 2020. Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics. TOG.
  • Wu et al. 2020. A Safe and Fast Repulsion Method for GPU-based Cloth Self Collisions. TOG.

• Frictional contacts are difficult to handle.

✅ 缺点：（1）难以模拟摩擦。（2）碰撞->施加力->调整，因此效果是滞后的。 优点：易实现
✅ 隐式积分比显式积分好，因为显式不稳定。

P19

Particle Collision Response —— Impulse Method

An impulse method assumes that collision changes the position and the velocity all of sudden.

✅ Penalty 方法是碰撞 → 力 → 下一时刻的速度和位置，效果滞后。
Impulse方法碰撞时立即更新速度和位置

✅ lmpulse 省去了力这一步，直接更新刚体状态。方法要求已经有一个比较好的\(\phi (x)\)
✅ 更新方法：N方向。更新距离：穿入的距离。

P20
Changing the position is not enough, we must change the velocity as well.

✅ \(\mathbf{v}\cdot \mathbf{N}\ge 0\)：当前速度想要让物体越陷越深, 这种情况下才需要更新速度
把\(\mathbf{v}\)分解为\(\mathbf{v_T}\)（切线方向的速度）和\(\mathbf{v_N})\)（法线方向的速度）.
✅ \(\mathbf{v_T}\)方向速度反弹， \(\mu _\mathbf{N}\) 为反弹系数。\(\mathbf{v_N}\)方向不变或由于摩擦再衰减
✅ a的约束：（1）越小越好，尽量把速度衰减掉（2）满足库仑定律（切方向的速度改变不应大于法线方向的速度改变）（3）切方向速度不能反转，即a不能为负
✅ 优点：可以精确控制摩擦力和反弹位置。缺点：计算比 Penalty 复杂
✅ 刚体常见于 Impulse； 弹性体常见于Pealty.

P21

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P23

Rigid Body Collision Detection

When the body is made of many vertices, we can detect its collision by testing each vertex

✅ 遍历 mesh 上的每个点，依次做碰撞检测。

No a perfect solution, but acceptable (will come back to this weeks later…)

P24

Rigid Body Collision Response

刚体碰撞响应与粒子碰撞响应的区别

以Vertex i为例，先分析i当前的位置和速度:

Problem: we cannot directly modify \(\mathbf{x}_i\) or \(\mathbf{v}_i\) isince they not state variables. They areindirectly determined.

✅\(x\)和\(v\)分别是刚体质心点的位置和速度,第二项为刚体上的特定点相对于质心点的位置和速度
对于粒子，可以直接用Impulse修改\(x\)和\(v\)
对于刚体，impulse只能修改\(x\)和\(v\)，不能修改\(x_i\)和\(v_i\)；其中\(x\)可以通过直接修改更新，也可以通过修改\(v\)来更新，这里选择后者。

Solution: we will find a way to modify \(\mathbf{v}\) and \(\mathbf{\omega}\).

✅ 解决方法：通过修改\(\mathbf{v}\)和\(\mathbf{\omega}\)实现修改\(x_i\)和\(v_i\)

P25

反向思考

What happens to \(\mathbf{v}_i\) when an impulse \(\mathbf{j}\) is appliedat vertex \(i\)?

✅ \(\mathbf{j}\) 是一个未知的冲量。\(\mathbf{v}_i\) 是点速度、\(\mathbf{v}\)是线速度
✅假设：此时对\(x_i\)点施加冲量\(j\)，会发生什么？

✅ 冲量 = \(Ft\) = \(m\Delta v \Rightarrow \Delta v\) = 冲量/\(m\)，由此得到\(v^{new}\)
✅ 冲量=时间 \(\cdot\)力 = 质量矩阵 * 时间 = 力矩 * t，省略公式中的时间，可得： \(Rr_i \times j\) = 冲量造成的力矩 ＝ 质量矩阵 · \(\Delta \omega \Rightarrow \Delta \omega\) ＝ 质量矩阵\(^{-1}\) · 力矩 ，由此得到\(\omega^{new}\)
❓ 为什么质量矩阵是单位阵？
✅ 由线速度\(v^{new}\)得到点速度\(\mathbf{v}_i^{new}\)

P27

$$ \mathbf{v_i^{new}} = \mathbf{v} _i+\frac{1}{M}\mathbf{j} −(\mathbf{Rr} _i)×(\mathbf{I} ^{−1}(\mathbf{Rr} _\mathbf{i}\times \mathbf{j} )) $$

$$ \mathbf{v_i^{new}} = \mathbf{v} _i+\frac{1}{M} \mathbf{j} −(\mathbf{Rr} _i)^∗\mathbf{I} ^{−1} (\mathbf{Rr} _i)^∗\mathbf{j} $$

✅ 向量之间的点乘可以转化为矩阵与向量的乘法，方便化简。具体内容见页面最后的补充1

化简得：

$$ \mathbf{v_i^{new}}-\mathbf{v}_i=\mathbf{Kj} $$ $$ \mathbf{K} \longleftarrow \frac{1}{M} \mathbf{1} −(\mathbf{Rr} _i)^{∗}\mathbf{I} ^{−1}(\mathbf{Rr} _i)^{∗} $$

✅ 结论，当碰撞点\(i\)确定时，冲量\(j\)和其造成的速度改变量\(Δv\)是确定的，这样，可以通过施加\(j\)，精确修改\(v_i\)
✅ 已知 \(\mathbf{v}_i^{new},\mathbf{v}_i,\mathbf{K}\),可求得 \(\mathbf{j}\)

P28

Rigid Body Collision Response by Impulse

✅ i点发生碰撞 -> 算出i点碰撞后的速度 -> 算出给i点什么样的冲量能让i出现碰撞后的效果 -> 真的施加这样一个冲量 -> 更新刚体状态

P29

Some Implementation Details

• If there are many vertices in collision, we use their position average.

✅ 如果有多个顶点发生碰撞呢？ 答：方法1，问题简化，用平均值。方法2，解线性系统，见下一页

• We can decrease the restitution \(\mathbf{\mu_N} \) to reduce oscillation（抖动）.

✅ 抖动原因：重力让它往下，冲量让它往上，导致在地面上反复振荡
解决方法：接近静止时衰减 \(\mathbf{\mu_N} \)

• We don't update the position here. Why?
  • Because the problem is nonlinear.
  • We will come back to this later when we talk about constraints.

P30

多碰撞点场景

Relative velocity at joints

$$ \begin{cases} \mathbf{v} _0 ^{\mathbf{new} }− \mathbf{v} _0=\mathbf{K} _{a00 }\mathbf{j} _0+\mathbf{K} _{a01 }\mathbf{j} _1 −(−\mathbf{K} _{b00 }\mathbf{j} _0 +\mathbf{K} _{b02}\mathbf{j} _2 )\\ \mathbf{v} _1 ^{\mathbf{new} }− \mathbf{v} _1=\mathbf{K} _{a10 }\mathbf{j} _0+\mathbf{K} _{a11 }\mathbf{j} _1 −(−\mathbf{K} _{c11 }\mathbf{j} _0 +\mathbf{K} _{c13 }\mathbf{j} _3 )\\ \mathbf{v} _2 ^{\mathbf{new} }− \mathbf{v} _2=\mathbf{K} _{b20 }\mathbf{j} _0+\mathbf{K} _{b22 }\mathbf{j} _2\\ \mathbf{v} _3 ^{\mathbf{new} }− \mathbf{v} _3=\mathbf{K} _{c31 }\mathbf{j} _1+\mathbf{K} _{c33 }\mathbf{j} _3 \end{cases} $$

$$ \Downarrow $$

$$ \begin{bmatrix} \mathbf{K} _{a00 }+\mathbf{K} _{b00 } & \mathbf{K} _{a01 } & -\mathbf{K} _{b02 } & \Box \\ \mathbf{K} _{a10 } & \mathbf{K} _{a11 }+\mathbf{K} _{c11 } & \Box & -\mathbf{K} _{c13 }\\ -\mathbf{K} _{b20 } & \Box & \mathbf{K} _{b22} & \Box \\ \Box & -\mathbf{K} _{c31 } & \Box & \mathbf{K} _{c33 } \end{bmatrix}\begin{bmatrix} \mathbf{j} _{0 }\\ \mathbf{j} _{1}\\ \mathbf{j} _{2}\\ \mathbf{j} _{3} \end{bmatrix}=\begin{bmatrix} \bigtriangleup \mathbf{v} _{0}\\ \bigtriangleup \mathbf{v} _{1}\\ \bigtriangleup \mathbf{v} _{2}\\ \bigtriangleup \mathbf{v} _{3} \end{bmatrix} $$

\(\mathbf{K} _{a01}\mathbf{j} _1\) stands for the velocity change of bunny \(a\) at joint 0, caused by impulse \(\mathbf{j}_1\).

P31

After-Class Reading (Before Collision)

https://graphics.pixar.com/pbm2001
Rigid Body Dynamics

P32

Shape Matching

✅ 用粒子的方法来解决刚体的问题

P33

Basic Idea

We allow each vertex to have its own velocity, so it can move by itself.

First, move vertices independently by its velocity, with collision and friction being handled.

Second, enforce the rigidity constraint to become a rigid body again.

✅ 第二步是 Shape Matching 的关键

Rigidity

P34

更新质心位置

Now \(\mathbf{c}\) and \(\mathbf{R}\) are unknowns we want to find out from:

✅ \(\mathbf{c}\) 代表质心，即前面的 \(\mathbf{x}\)
✅ 约束：新的顶点位置与原顶点位置的距离尽量接近。
✅问题简化：用任意矩阵A代替需要满足旋转矩阵约束的\(R\)。因此\(\sum Ar_i = A \sum r_i = 0\)
✅结论：约束前后质心位置不变
❓ 优化之后的刚体可能还是与地面穿透的。

P35

更新质量速度

✅ 先假设 \(\mathbf{R}\) 是任意矩阵 \(\mathbf{A}\),再从中提取旋转成分
✅ Polar Decomposition：极性分解，把任意矩阵分解旋转部分和形变部分。

P36

极性分解

Singular value decomposition says any matrix can be decomposed into: rotation,scaling and rotation: \(\mathbf{A = UDV} ^T\).

We can rotate the object back before the final rotation: \(\mathbf{A} = (\mathbf{UV} ^T)(\mathbf{VDV} ^T)\).

✅ \(\mathbf{A} = （\mathbf{UV}^T）(\mathbf{VDV}^T) =\mathbf{RS}\)
\(\mathbf{R}\) 代表全局旋转，\(\mathbf{S}\)代表本地形变，扔掉S保留R。

$$ \mathbf{A=RS} $$

$$ \mathbf{A} ^T\mathbf{A} = \mathbf{S} ^T\mathbf{S} = \mathbf{S} ^2 $$

分解结果:unique

P39

Shape Matching Pipeline

Physical quantities are attached to each vertex, not to the entire body.

P40

算法分析

• 优点：Easy to implement and compatible with other nodal systems, i.e., cloth, soft bodies and even particle fluids.

• 局限性：Difficult to strictly enforce friction and other goals. The rigidification process will destroy them.

• 适用场景：More suitable when the friction accuracy is unimportant, i.e., buttons on clothes.

P41

After-Class Reading

Muller et al. 2005.
Meshless Deformations Based on Shape Matching. TOG (SIGGRAPH).

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

补充： Cross Product as a Matrix Product

We can convert the cross product \(\mathbf{r}\times\) into a matrix product \(\mathbf{r}^*\).

✅ \(\mathbf{r}^*\) 是 \(\mathbf{r}\) 的 cross matrix.
目的：用矩阵形式代替叉乘形式，方便公式化简

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• A Mass-Spring System
  • Explicit Integration
  • Implicit Integration
• Bending and Locking Issues
• Shape matching

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P4

弹簧结构

An Ideal Spring

An ideal spring satisfies Hooke’s law: the spring force tries to restore the rest length.

\(E(\mathbf{x})=\frac{1}{2}k (||\mathbf{x}_i −\mathbf{x}_j||−\mathbf{L} )^2\)

\(\mathbf{f} _i(\mathbf{x} )=−∇_i\mathbf{E} =−k(||\mathbf{x}_i −\mathbf{x}_j||−L)\frac{\mathbf{x}_i −\mathbf{x}_j}{||\mathbf{x}_i −\mathbf{x}_j ||}\)

\(\mathbf{f} _j(\mathbf{x})=−∇_jE=−k (||\mathbf{x}_j −\mathbf{x}_i ||−L)\frac {\mathbf{x}_j −\mathbf{x}_i}{||\mathbf{x}_j −\mathbf{x}_i||}\)

P5

Multiple Springs

When there are many springs, the energies and the forces can be simply summed up.

$$ E= {\textstyle \sum_{e=0}^{3}}E_e= {\textstyle \sum_{e=0}^{3}} (\frac{1}{2} k(||\mathbf{x} _i −\mathbf{x}_e ||−L_e)^2) $$

$$ f_i=−\nabla_iE = \textstyle \sum_{e=0}^{3}(−k(||\mathbf{x}_i−\mathbf{x}_e||−L_e)\frac{\mathbf{x}_i−\mathbf{x}_e}{||\mathbf{x}_i−\mathbf{x}_e||}) $$

✅ 能量和力都是可以叠加的

P6

Structured Spring Networks

✅ 绿线：防止斜方向的拉伸。蓝线：防止翻折。

P7

Unstructured Spring Networks

We can also turn an unstructured triangle mesh into a spring network for simulation.

✅ 蓝线：抵抗弯曲。对每条内部边，加这样一根弹簧。

P8

怎样基于三角形Mesh增加蓝线弹簧

The basic representation of a triangle mesh uses vertex and triangle lists.

✅ 已知边的信息，需找出内部边，例如\(\mathbf{x}_0\mathbf{x}_3\)，因此要基于此构造边：\(\mathbf{x}_1x 4\)
✅ Each triangle has three edges. But there are repeated ones. Repeated edges就是内部边。
✅ 1. 找出内部边。2. 找出内部边所属于的两个三角形。3. 找出两个三角形上不在这条内部边上的点。4. 连续一根弹簧。

Vertex list: {\(\mathbf{x} _0, \mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3, \mathbf{x}_4\)} (3D vectors)

Triangle list: {1, 2, 3, 0, 1, 3, 0, 3, 4} (index triples)

The key to topological construction is to sort triangle edge triples.

Each triple contains: edge vertex index 0, edge vertex index 1 and triangle index (index 0<index).

✅ 排序：基于边排序，排序后相同边会靠在一起

P11

积分系统

Explicit Integration of A Mass-Spring System

✅ 整体流程跟刚体运动很像，只是力变得复杂，每个弹簧端点上受到的力都要考虑，但没有了旋转。
✅ \( E [e] [0] ：e\)代表弹簧 ID:0或1代表弹簧两个端点
❗ 图画得不对，先提前把所有的力都算出来，再遍历所有顶点

P12

Explicit integration suffers from numerical instability caused by overshooting, when the stiffness \(k\) and/or the time step \(∆t\) is too large.

✅ Explicit：当前力 → 当前速度 → 当前位置
根据公式\(FΔt≈mv，vΔt≈Δx\)，如果\(Δt\)太大，会导致\(Δx\)太大，而导致overshooting。

A naive solution is to use a small \(∆t\) . But that slows down the simulation.

✅ 解决方法：减小\(\Delta t\)。但这个方法不解决本质问题，且会降低整个模拟系统的效率
✅ 本质上是\(Δt\)太大导致积分近似的结果与实际积分的结果有很大误差，\(k\)太大或\(Δt\)只是让这个问题更明显，减小\(k\)或\(Δt\)问题仍然存在。

P13

Implicit Integration

Implicit integration is a better solution to numerical instability. The idea is to integrate both x and v implicitly.

✅Explicit和Implicit都是用某个时刻的力代表整个\(Δt\)时间的力，就都会出现上述误差。
区别在于，Explicit用当前力，往往使结果变大，产生爆炸，Implicit用未来力，往往使结果变小，产生消失。
✅ 消失只是结果不对。但爆炸会让结果崩溃，这对是最不可接受的问题。因此用隐式代替显式。

消元得：

✅ 下面公式1通过把上面公式1代入公式2得到。下面公式2通过把上面公式写反推得到。

Assuming that \(\mathbf{f}\) is holonomic, i.e., depending on \(\mathbf{x}\) only, our question is how to solve:
$$ \mathbf{x} ^{[1]}=\mathbf{x}^{[0]}+∆t\mathbf{v} ^{[0]}+∆t^2\mathbf{M} ^{−1}\mathbf{f} (\mathbf{x}^{[1]}) $$

✅ holonomic：力的大小和方向只跟位置有关，跟速度无关。例如重力，弹力。那么 \(f\)可以写成关于位置的函数\(f(x)\)。
✅ 但\(f(x)\)不一定是线性的。因此最后转化为解非线性方程的问题。未知量为\({x} ^{[1]}\)

P14

\(\mathbf{x} ^{[1]} =\) argmin \(F(\mathbf{x})\quad\) for \(\quad F(\mathbf{x}) = \frac{1}{2∆t^2}||\mathbf{x} −\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]}||_M^2+E(\mathbf{x} )\)

✅ 前面方程解\({x} ^{[1]}\)等价于F(x)函数极小点。等价转换的推导在补充1。非线性方程问题为转化为优化问题。
其中：\(\mathbf{M}\)对角矩阵，描述质量，\(3N \times 3N\)。\(\mathbf{x}\)为 \(3N\times 1\)矢量,描述顶点信息。\(E\) 为所有的力的能量。\(\mathbf{||x||_M^2=x^TMx} \)。
✅ 只有保守力能用能量描述、非保守力（例如摩擦力）则不行。

P15

优化方法

P18

Simulation by Newton’s Method

🔎 Newton-Raphson Method见补充2. 这里直接开始Newton方向本当前场景的应用。

Specifically to simulation, we have:

$$ F (\mathbf{x} )=\frac{1}{2∆t^2}||\mathbf{x} −\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]}||_\mathbf{M} ^2+\mathbf{E} (\mathbf{x} ) $$

$$ ∇F(\mathbf{x}^{(k)})=\frac{1}{∆t^2}\mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})−\mathbf{f}(\mathbf{x}^{(k)}) $$

$$ \frac{∂^2F (\mathbf{x} ^{(k)})}{∂\mathbf{x} ^2} =\frac{1}{∆t^2} \mathbf{M} +\mathbf{H} (x^{(k)}) $$

❓ 每一个 step 都包含一次迭代优化、能做到实时？
✅ 本节课所讲的套路：分析力 → 隐式积分 → 优化问题 → 更新，对弹簧系统、有限元、弹性体等各种物理模拟同样适用
✅ 早期的方式不是用优化来做的，而是近似成线性问题后直接解方程组。这种方法相当于每一个Step做了一次牛顿法。

P19

Solve Spring Hessian

According to Lecture 2, Page 48,

✅ 弹簧系统的H是由所有弹簧的H构成的。

✅ \(H(x)\)的维度是\(3N \times 3N\)，N 是弹簧数。每个\(H_e\)的维度是\(3 \times 3\)。
✅ 课后答疑：质点的质量可以不同吗？
答：可以不同。先根据三角形的面积计算三角的质量，再把质量分配到各个顶点上。

Positive Definiteness of Hessian

✅ \(H(x)\)的正定性由\(H_e\)的正定性决定。

下面分析\(H_e\)的正定性：
For any \(\mathbf{x} _{ij}, \mathbf{v} ≠0\),

$$ \mathbf{V}^\mathbf{T}\frac{{\mathbf{x} _{ij}\mathbf{x} _{ij}}^\mathbf{T} }{||\mathbf{x} _{ij}||^2}\mathbf{V}=||\frac{{\mathbf{x} _{ij}}^\mathbf{T} \mathbf{v} }{||\mathbf{x} _{ij}||}||^2> 0 $$

$$ \mathbf{V} ^\mathbf{T} (\mathbf{I} -\frac{{\mathbf{x} _{ij}\mathbf{x} _{ij}}^\mathbf{T} }{||\mathbf{x} _{ij}||^2}) \mathbf{V} =\frac{||\mathbf{x} _{ij}||^2||\mathbf{v} ||^2-||{\mathbf{x} _{ij}}^\mathbf{T} \mathbf{v} ||^2}{||\mathbf{x} _{ij}||^2}\ge 0 $$

✅ \( \mathbf{x}_ {ij}\) 代表顶点\( \mathbf{x}_ {i}\)和顶点\( \mathbf{x}_ {j}\)的位置的差。
✅ 最后一个公式分子满足柯西不等式
✅ 结论：\(||x_{ij}||< Le\). 代表弹簧处于压缩状态。此时 He 有可能非正定，但拉伸时一定正定。
He 正定则\(H(x)\)半正定，此时弹簧系统有唯一解。

P20
When a spring is stretched, \(\mathbf{H} _e\) is s.p.d.; but when it’s compressed, \(\mathbf{H} _e\) may not be s.p.d.

As a result, \(\mathbf{H}(\mathbf{x})\) may not be s.p.d. (Lecture 2, Page 36).

\(\mathbf{A}\) may not be s.p.d. either.

✅ \(\Delta t\)越小，A越容易正定、弹簧系统越稳定。
✅为什么要讨论\(H\)矩阵是否正定？答：\(H\)矩阵相当于二阶导，正定代表开口向上，有唯一最小值。
✅ 但是A不正定，不代表没有唯一解。

P22

When a spring is compressed, the spring Hessian may not be positive definite. This means there can be multiple local minima (outcomes).

P23

Enforcement of Positive Definiteness

• Nevertheless, some linear solvers can fail to work if the matrix \(\mathbf{A}\) in \(\mathbf{A}\Delta \mathbf{x}=\mathbf{b}\) is not positive definite.

✅ 不正定最大的问题不是解不唯一，因为解出任意一个解都能让模拟系统进行下去。
非正定的主要问题，是数学计算上的不稳定，可能导致解不出来；

• One solution is to simply drop the ending term, when \({\color{Orange}{ ||\mathbf{x} _{ij}||<\mathbf{L} _e}}:\)

✅ 简单粗爆的解决方法就是把后面这项删掉。

• Other solutions exist. For example,

🔎 Choi and Ko. 2002. Stable But Responive Cloth. TOG (SIGGRAPH)

P24

Linear Solvers

The Jacobi Method

We can use the Jacobi method to solve \(\mathbf{A}∆\mathbf{x} = \mathbf{b} \).

The vanilla Jacobi method (\(α\) = 1) has a tight convergence requirement on \(\mathbf{A}\), i.e., being diagonal dominant.

The use of \(α\) allows the method to converget even when \(\mathbf{A}\) is positive definite only.

P25

An Incomplete Summary

• Direct Solvers (LU, LDLT, Cholesky, …)
  • One shot, expensive but worthy if you need exact solutions.
  • Little restriction on \(\mathbf{A}\)
  • Mostly suitable on CPUs

🔎 Intel MKL PARDISO

• Iterative Solvers
  • Expensive to solve exactly, but controllable
  • Convergence restriction on \(\mathbf{A}\), typically positive definiteness
  • Suitable on both CPUs and GPUs
  • Easy to implement
  • Accelerable: Chebyshev, Nesterov, Conjugate Gradient…

✅ 课后答疑
问题二：怎么加速？
答：用 Jacobian 可以在 GPU 上加速、直接法比迭代法慢。
问题三：共轭梯度
共轭梯度的效率很大程度上取决于 precondition,但在GPU上能使用的precondition 比较受限、 CPU 上一般选择 Incomplete LU 分解。
问题四：支持的维度
直接法比较占内存，因此支持的维度不如迭代法大。

P26

The Jacobi Method with Chebyshev Acceleration

We can use the accelerated Jacobi method to solve \(\mathbf{A}∆\mathbf{x} =\mathbf{b} \).

The Accelerated Jacobi Method
\(∆\mathbf{x} \longleftarrow \mathbf{0} \)
last_\(∆\mathbf{x} \longleftarrow \mathbf{0}\)
For \(k=0\dots \mathbf{K}\)
\(\mathbf{r} \longleftarrow \mathbf{b} −\mathbf{A} ∆\mathbf{x}\)
If \(||\mathbf{r} ||<\omega \quad\) break
If \(k=0 \quad\quad\quad \omega =1\)
Else If \( k=1 \quad \quad\quad\omega =2/(2-\rho^2)\)
Else \(\quad\quad\quad\omega =4/(4-\rho ^2\omega )\)
old_\(∆ \mathbf{x} \longleftarrow ∆ \mathbf{x}\)
\(∆\mathbf{x} ⟵∆\mathbf{x} +\mathbf{αD} ^{−1}\mathbf{r}\)
\(∆\mathbf{x} \longleftarrow \omega ∆ \mathbf{x} +(1−\omega)\)last_∆\(\mathbf{x}
\)
last_\(∆\mathbf{x} \longleftarrow \) old_\(∆\mathbf{x}\)

✅ 这一页老师没讲

P27

After-Class Reading

Baraff and Witkin. 1998. Large Step in Cloth Simulation. SIGGRAPH.

One of the first papers using implicit integration.

The paper proposes to use only one Newton iteration, i.e., solving only one linear system. This practice is fast, but can fail to converge.

✅这篇论文是衣服模拟的经典论文，第一个用隐式积分做衣服模型的论文。
论文没有用弹簧系统，而是另一套模型。
没有做非线性优化或解非线性方程，而是把非线性方程线性化，等价于做一次牛顿迭代。

补充1：非线性方程求解

求解的非线性方程如下，其中\({x} ^{[1]}\)是未知量。
$$ \mathbf{x} ^{[1]}=\mathbf{x}^{[0]}+∆t\mathbf{v} ^{[0]}+∆t^2\mathbf{M} ^{−1}\mathbf{f} (\mathbf{x}^{[1]}) $$

P14

$$ \mathbf{||x||_M^2=x^TMx} $$

✅ Note that this is applicable to every system, not just a mass-spring system.

把公式处理一下得，
$$ x^{[0]}+Δtv^{[0]}+Δt^2M^{-1}f(x^{[1]})-x^{[1]}=0 $$ 左右两边同时乘以\(\frac{M}{Δt^2}\)得
$$ \frac{1}{Δt^2} M(x^{[1]}-x^{[0]}-Δtv^{[0]})-f(x^{[1]})=0 $$
这里面唯一的未知量是\(x^{[1]}\)，定义函数 $$ y=\frac{1}{Δt^2} M(x-x^{[0]}-Δtv^{[0]})-f(x) $$
当\(x = x^{[1]}\) 时，\(y = 0\), 即 \(y(x^{[1]}) = 0\)
从另一个角度讲， $$ \begin{eqnarray} x^{[1]} & = \mathrm{argmin}& F(x)\Rightarrow {F}' (x^{[1]}) & = & 0 \end{eqnarray} $$
因此, \({F}' (x) = y. \quad F(x) = \int ydx \)

补充2：Newton-Raphson Method

x是值的F(x)函数

The Newton-Raphson method, commonly known as Newton’s method, solves the optimization problem: \(x^{[1]}\) = argmin \(F(x)\).

Given a current \(x^{(k)}\), we approximate our goal by:

$$ 0={F}' (x)≈{F}'(x^{(k)})+{F}'' (x^{(k)})(x−x^{(k)}) $$

✅ \(a = \min F(x)⇒ F'(a)= 0\)，\({F}' (x)\) 是非线性函数，直接解\({F}' (x)=0\) 很难解
✅ 对\({F}'(x)\) 做一阶泰勒展开，保留到二阶项。
假设\(x^{[k]}\)为任意已知值，就变成了解线性方程，很容易解出\(x\).
因为\({F}'(x)\) 是一个近似的，\(x\) 也是一个近似解。但\(x^{[k]}\) 越接近真实解，\(x\) 也会越接近真实解。因此，选代是\(x^{[k]}\)和\(x\) 都不断逼近真实解的过程。
✅ 普通的梯度下降是把\({F}' (x)\) 近似到一阶，牛顿法是近似到二阶，因此下降更快。

✅ Overshooting 的本质：误差会积累和放大

P16
Newton’s method finds an extremum, but it can be a minimum or maximum.

• At a minimum \(x^∗, {F}'' (x^∗)>0\).
• At a maximum \(x^∗, {F}''(x^∗)<0\).
• If \({F}''(x)>0\) is everywhere, \(F(x)\) has no maximum. \(=> F(x)\) has only one minimum.

✅ \(F'(a)= 0,a\) 有可能是最大值或最小值，因此要判定解是否合理。判定方法： \({F}''(x)\)

P17

x是向量的F(x)函数

Now we can apply Newton’s method to: \(x^{[1]} \)= argmin \(F(x)\). Given a current \(x^{(k)}\), we approximate our goal by:

$$ 0=\nabla F( \mathbf{x}) ≈\nabla F (\mathbf{x} ^{(k)})+\frac{∂F ^2(\mathbf{x} ^{(k)})}{∂\mathbf{x} ^2} (\mathbf{x−x} ^{(k)}) $$

✅ 按照 \(\Delta x\) 的更新公式，只需要用到\(F'(x)\) 和 \({F}''(x)\)， 不需要知道 \(F(x)\).
✅ 此处\(x\)是向量，因此\(F'(x)\)是向量，\({F}''(x)\)是 Hession 矩阵

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P29

The Bending Spring Issue

A bending spring offers little resistance when cloth is nearly planar, since its length barely changes.

✅黑线为三角形面片，每条边一根弹簧，并增加一根蓝线弹簧，构成弯曲弹簧，阻止两个面片弯折。
存在的问题：小的弯折，弹簧长度几乎不变，抵抗弯曲的力量非常弱。（不适用于类似于纸的弯折效果）。

P30

A Dihedral Angle Model

A dihedral angle model defines bending forces as a function of \(\theta : \mathbf{f} _i=f (\theta )\mathbf{u} _i\).

✅ Dihedarl Angel:二面角
✅ 把弯曲的力写成关于二面角的函数

✅ \(x_1, x_2, x_3, x_4\) 都会受到 bending force. 力的大小相同但方向不同，但都是关于二面角的函数。
✅\(u_i\)：描述力的方向，与\(\theta\)大小无关。\(f(\theta)\)：描述力的大小，是关于\(\theta\)的函数。

• First, \(\mathbf{u}_1\) and \(\mathbf{u}_2\) should be in the normal directions \(\mathbf{n}_1\) and \(\mathbf{n}_2\).

• Second, bending doesn’t stretch the edge, so \(\mathbf{u}_4\)−\(\mathbf{u}_3\) should be orthogonal to the edge, i.e., in the span of \(\mathbf{n}_1\) and \(\mathbf{n}_2\).

• Finally, \(\mathbf{u}_1+\mathbf{u}_2+\mathbf{u}_3+\mathbf{u}_4=\mathbf{0}\), which means \(\mathbf{u}_3\) and \(\mathbf{u}_4\) are in the span of \(\mathbf{n}_1\) and \(\mathbf{n}_2\).

✅ 合力为0。

P31
Conclusion:

✅ N是未归一化的 normal. N 的方向与 normal 相同。大小为三角形的面积。
✅ 重要的不是结果，而是根据观察进行合理假设的思考过程。

P32

Planar case:

$$ \mathbf{f} _i=k\frac{||\mathbf{E}||^2}{||\mathbf{N}_1||+||\mathbf{N}_2||} \sin(\frac{π−\theta}{2})\mathbf{u} _i $$

Non-planar case:

$$ \mathbf{f} _i=k\frac{||\mathbf{E} ||^2}{||\mathbf{N} _1||+||\mathbf{N} _2||}(\sin(\frac{π−\theta}{2})-\sin(\frac{π−\theta_0}{2}))\mathbf{u}_i $$

✅ Non-planar case：不是指弯曲时的力，而是指静止状态(reference state)为非平面的场景下，弯曲为\(\theta\)时的力。\(\theta_0\)表示 reference state. ✅ 老师没解释公式怎么来的
🔎 Bridson et al. 2003. Simulation of Clothing with Folds and Wrinkles. SCA.
✅ 此论文适合读完。除了弯曲模型，还有一些有意思的设计。

Explicit integration.
Derivative is difficult to compute.

✅ 由于完全基于力而不考虑能量，适合用显式积分。

P34

A Quadratic Bending Model

✅二面角方法是纯分析力的方法，比较复杂。此处是Bending issue的另一个方法。

A quadratic bending model has two assumptions: 1) planar case; 2) little stretching.

$$ E(\mathbf{x} )=\frac{1}{2} \begin{bmatrix} \mathbf{x}_0 & \mathbf{x}_1 & \mathbf{x}_2 & \mathbf{x}_3 \end{bmatrix}\mathbf{Q} \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \\ \mathbf{x}_2\\ \mathbf{x}_3 \end{bmatrix} $$

$$ \mathbf{Q} =\frac{3}{\mathbf{A} _0+\mathbf{A} _1}\mathbf{qq^T} $$

✅ \({\mathbf{A} _0}\)和\({\mathbf{A} _1}\)是两个三角形在reference状态下的面积。

$$ \mathbf{q} = \begin{bmatrix} (\cot\theta _1+ \cot\theta _3)\mathbf{I} \\ (\cot\theta _0+ \cot\theta _2)\mathbf{I} \\ (-\cot\theta _0- \cot\theta _1)\mathbf{I} \\ (-\cot\theta _2- \cot\theta _3)\mathbf{I} \end{bmatrix} $$

\(\mathbf{I}\) is 3-by-3 identity.

✅ \(\mathbf{Q}\)只与\(\mathbf{\theta}\)有关，因此是一个定值。

It’s not hard to see that: \(E (\mathbf{x} )=\frac{3||\mathbf{q} ^\mathbf{T}\mathbf{x} ||^2}{2(A_0+A_1)}\). Also, \(E (\mathbf{x} )=0\) when the triangles are flat.

✅ \(\mathbf{q^T}\mathbf{x}\)在估算两个三角形的拉普拉斯，即两个三角的曲率、当两个三角形共面时， \(E(\mathbf{x})=0\)
🔎 离散曲面的拉普拉斯，见GAMES102
✅ \(E(\mathbf{x})\) 来自数学上曲率的推导，而不是来自物理意义的推导。
✅ 问题：能量的思想能用在刚体上吗？
答：这里的能量是弹性能量、刚体无弹性，因此也无所谓能量。

Pros of The Quadratic Bending Model

• Easy to implement:

✅ \(E(\mathbf{x})\)是关于\(\mathbf{x}\)的二次函数，很容易计算\(E(\mathbf{x})\)的一阶导（力）和二阶导\(\mathbf{H} \)

$$ \mathbf{f} (\mathbf{x} )=−\nabla \mathbf{E} (x)= −\mathbf{Q} \begin{bmatrix} \mathbf{x} _0\\ \mathbf{x} _1\\ \mathbf{x} _2 \\ \mathbf{x} _3 \end{bmatrix} $$

$$ \mathbf{H} (\mathbf{x} )=\frac{∂^2E(\mathbf{x} )}{∂\mathbf{x} ^2}=\mathbf{Q} $$

• Compatible with implicit integration.

Cons of The Quadratic Bending Model

• No longer valid if cloth stretches much.

✅方法假设面料拉伸比较小，当面料拉伸太大，\(\mathbf{\theta}\)就会改变，\(\mathbf{Q}\)就不准了。

• Not suitable if the rest configuration is not planar.
  • Cubic shell model.
  • Projective dynamics model.

After Class Reading

🔎 Bergou et al. 2006. A Quadratic Bending Model for Inextensible Surfaces. SCA.
✅ 这篇论文是在本算法上的进一步工作。

P37

The Locking Issue

So far we talked about the mass-spring model and other bending models, assuming cloth planar deformation and cloth bending deformation are independent.

Is it true? Think about a zero bending case. Can a simulator fold cloth freely?

✅ 正常来讲拉伸和弯曲是两件独立的事情。但在弹簧模型系统中，把它们耦合了。
例如纸这种无弹性的面料，会把它的弹性系数调得很大，达到无弹性的效果。但导致了它无法弯折的artifacts。
✅ 在K很大或网格分辨率低时， locking issue 会特别明显。

P38

The fundamental reason is due to a short of degrees of freedoms (DoFs).
For a manifold mesh, Euler’s formula says:#edges=3#vertices-3-#boundary_edges.
So if edges are all hard constraints, the DoFs are only: 3+ #boundary_edges.

✅ 自由度 = 变量数 - 约束数。
每个顶点有3个自由度、每条边是一个约束，因此单纯加点不会改善，但让点变密可以改善
✅ 实操套路：1. 弹簧压缩时让k比较小；2. 假设弹簧在一定长度范围内可自由活动，不受力，以上方法都不解决根本问题；3. 把自由度定义在边上不­是顶点上，但把问题搞得更复杂了。

P43

A Summary For the Day

• A mass-spring system

  • Planar springs against stretching/compression \(\quad\)- replaceable by co-rotational model
  • Bending springs \(\quad\)- replaceable by dihedral or quadratic bending
  • Regardless of the models, as long as we have \(E (\mathbf{x})\), we can calculate force \(\mathbf{f} (\mathbf{x} )=−∇ \mathbf{E} (\mathbf{x})\) and Hessian \(\mathbf{H} (\mathbf{x} )=∂E^2(\mathbf{x} )/∂\mathbf{x} ^2\). Forces and Hessians are stackable.

• Two integration approaches

  • Explicit integration, just need force. Instability
  • Implicit integration, as a nonlinear optimization problem
  • One way is to use Newton’s method, which solves a linear system in every iteration:

$$ (\frac{1}{∆t^2}\mathbf{M} +\mathbf{H} (\mathbf{x} ^{(k)}))∆\mathbf{x} =− \frac{1}{∆t^2} \mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})+\mathbf{f} (\mathbf{x} ^{(k)}) $$

• There are a variety of linear solvers (beyond the scope of this class).
• Some simulators choose to solve only one Newton iteration, i.e., one linear system per time step.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P39

Shape Matching

✅ Shape Matching 跳过了。

P40

Shape Matching

The basic idea is to define a quadratic energy based on the rotated reference element. To do so, we split transformation into deformation + rotation.

P41

Shape Matching

The basic idea is to define a quadratic energy based on the rotated reference element. To do so, we split transformation into deformation + rotation.

P42

Shape Matching

We can then define the quadratic energy as:

$$ E (\mathbf{x} )=\frac{1}{2}||\mathbf{F−R} ||^2 $$

(\(\mathbf{R}\) is the rotation inside of \(\mathbf{F}\). This energy tries to penalize the existence of \(\mathbf{S}\)).

Assuming that \(\mathbf{R}\) is constant, this \(E(\mathbf{x})\) becomes a quadratic function. We can then derive the force and the Hessian.

$$ E(\mathbf{x} ) =\frac{1}{2} ||\begin{bmatrix} \mathbf{x} _1-\mathbf{x} _0 &\mathbf{x} _2-\mathbf{x} _0 \end{bmatrix}\begin{bmatrix} \mathbf{r} _1-\mathbf{r} _0 &\mathbf{r} _2-\mathbf{r} _0 \end{bmatrix}^{−1}−\mathbf{R}||^2 $$

P43

A Summary For the Day

• A mass-spring system

  • Planar springs against stretching/compression \(\quad\)- replaceable by co-rotational model
  • Bending springs \(\quad\)- replaceable by dihedral or quadratic bending
  • Regardless of the models, as long as we have \(E (\mathbf{x})\), we can calculate force \(\mathbf{f} (\mathbf{x} )=−∇ \mathbf{E} (\mathbf{x})\) and Hessian \(\mathbf{H} (\mathbf{x} )=∂E^2(\mathbf{x} )/∂\mathbf{x} ^2\). Forces and Hessians are stackable.

• Two integration approaches

  • Explicit integration, just need force. Instability
  • Implicit integration, as a nonlinear optimization problem
  • One way is to use Newton’s method, which solves a linear system in every iteration:

$$ (\frac{1}{∆t^2}\mathbf{M} +\mathbf{H} (\mathbf{x} ^{(k)}))∆\mathbf{x} =− \frac{1}{∆t^2} \mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})+\mathbf{f} (\mathbf{x} ^{(k)}) $$

• There are a variety of linear solvers (beyond the scope of this class).
• Some simulators choose to solve only one Newton iteration, i.e., one linear system per time step.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• Strain Limiting and Position Based Dynamics
• Projective Dynamics
• Constrained Dynamics

P4

The Stiffness Issue

• Real-world fabrics resist strongly to stretching, once they stretch beyond certain limits.

• But, increasing the stiffness can cause problems.

  • Explicit integrators will be unstable
    • Solution: smaller time steps and more computational time.
  • The linear systems involved in Implicit integrators will be ill-conditioned.
    • Solution: more iterations and computational time.

• Can we achieve high stiffness, with a low computational cost?

✅当弹簧系数大的情况下，仍能保证系统稳定。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

投影函数 Projection Function

P5

A Single Spring

If a spring is infinitely stiff, we can treat the length as a constraint and define a projection function.

\(\mathbf{ϕ} (\mathbf{x} )=||\mathbf{x} _i− \mathbf{x} _j||−L=0\)

✅ projection function.投影函数、会移动端点使弹簧满足约束。

P6

✅ 把\(\mathbf{x}_ i\)和\(\mathbf{x}_ j\)拼成6维空间中的点\(\mathbf{x}\)，满足约束的\(\mathbf{x}\)构成6D空间中的一块区域；

{\(\mathbf{x} _i^{\mathbf{new}},\mathbf{x} _j^{\mathbf{new} }\)}= argmin \( \frac{1}{2}\){\(m_i||\mathbf{x} _i^{\mathbf{new} }−\mathbf{x} _i||^2+m_j||\mathbf{x} _j^{\mathbf{new}} −\mathbf{x} _j||^2\)}

such that \(\mathbf{ϕ} (\mathbf{x} )=0\)

✅ 投影函数的目标：(1)把\(\mathbf{x}\)移到区域内。 (2)移动距离最短，因此构成优化问题。
✅优化问题，但不是通过选代解决，而是数值求解，直接算出最优的\(\mathbf{x}_i\)和\(\mathbf{x}_j\).

P7

$$ \mathbf{x} ^{\mathbf{new} } \longleftarrow \mathrm{Projection} (\mathbf{x}) $$

$$ \mathbf{x} _i^{\mathbf{new} }\longleftarrow \mathbf{x} _i−\frac{m_j}{m_i+m_j} (||\mathbf{x} _i−\mathbf{x} _j||−L)\frac{\mathbf{x} _i−\mathbf{x}_j}{||\mathbf{x} _i−\mathbf{x} _j||} $$

$$ \mathbf{x} _j^{\mathbf{new} }\longleftarrow \mathbf{x} _j+\frac{m_i}{m_i+m_j} (||\mathbf{x} _i−\mathbf{x} _j||−L)\frac{\mathbf{x} _i−\mathbf{x}_j}{||\mathbf{x} _i−\mathbf{x} _j||} $$

$$ \quad $$

$$ \mathbf{ϕ} (\mathbf{x} ^{\mathbf{new} })=||\mathbf{x} _i^{\mathbf{new} }− \mathbf{x} _j^{\mathrm{new} }||−L=||\mathbf{x} _i−\mathbf{x} _j−\mathbf{x} _i+\mathbf{x} _j+L||−L=0 $$

✅ 对推导结果的合理性解释：(1) 移到前后质心不变。(2) 移到方向为沿着或远离质心。(3) 移到距离与自身重量有关。

By default, \(m_i=m_j\), but we can also set \(m_i=\infty\) for stationary nodes.

✅ 对于固定点，不更新就好了。

P8

Multiple Springs – A Gauss-Seidel Approach

What about multiple springs? The Gauss-Seidel approach projects each spring sequentially in a certain order. Imagine two springs with unit rest lengths…

P9

• We cannot ensure the satisfaction of every constraint. But the more iterations we use, the better those constraints are satisfied.

• Although the name is related to Gauss-Seidel, it differs from Gauss-Seidel. It is more relevant to stochastic gradient descent (in machine learning).

• The order matters. The order can cause bias and affect convergence behavior.

✅顺序影响结果的偏向性和收敛速度。

P10

Multiple Springs – A Jacobi Approach

• To avoid bias, the Jacobi approach projects all of the edges simultaneously and then linearly blend the results.

• The problem is an even lower convergence rate.

• Again, the more iterations it uses, the better the constraints are enforced.

  

✅ 基于每条边计算\(\mathbf{x}\)的更新但不真的更新、每个点会得到多种更新方案，最后取更新方案的均值。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

Position Based Dynamics (PBD)

算法过程

Position based dynamics (PBD) is based on the projection function.

• The stiffness behavior, i.e., how tightly constraints are enforced, is subject to non-physical factors.

  • The number of iterations
  • The mesh resolution

• The velocity update following projection is important to dynamic effects.

• This method is applicable to other constraints as well, including triangle constraints, volume constraints, and collision constraints.

  • To implement these constraints, simply define their projection functions.

A PBD Simulator
//Do Simulation, update \(\mathbf{x}\) and \(\mathbf{v}\)
$$ \mathbf{v}\longleftarrow\dots $$

$$ \mathbf{x}\longleftarrow\dots
$$

✅ 第一步：不考虑约束，基于粒子运动方法更新 \(\mathbf{v}\) 和 \(\mathbf{x}\)；

//Now PBD starts.

$$ \mathbf{x} ^{\mathbf{new} } \longleftarrow \mathrm{Projection} (\mathbf{x} ) $$

$$ \mathbf{v}\longleftarrow \mathbf{v} +(\mathbf{x} ^{\mathbf{new} }−\mathbf{x})/∆t $$

$$ \mathbf{x}\longleftarrow \mathbf{x} ^{\mathbf{new}} $$

✅ 第二步：基于约束和投影函数更新 \(\mathbf{v}\) 和 \(t\).
✅第二步中速度更新很重要，会影响dynamic模拟的效果。
✅\(\mathbf{v}\)的更新不是直接覆盖，而是叠加。

P12

Pros and Cons of PBD

• Pros
  • Parallelable on GPUs (PhysX)
  • Easy to implement
  • Fast in low resolutions
  • Generic, can handle other coupling and constraints, including fluids

✅ 一般来说，少于1000个点时能实时，多于1000个点时效率明显下降
✅ PBD 适用于低分辨率场景、常见的低精度实时模拟的套路。
❗ 模拟真正的时间开销不在计算 (虽然有很多计算公式) 而是在内存的访问上。
PBD 的优点是内存访问少、因为它没有太多物理变量。
因此，对追求效率的场景，主要优化内存访问而不是计算。

• Cons
  • Not physically correct
  • Low performance in high resolutions
    • Hierarchical approaches (can cause oscillation and other issues…)
    • Acceleration approaches, like Chebyshev

✅ 问：PBD 非物理方法，怎么体现出弹性效果？
答：迭代数多则弹性差、网格顶点少则弹性差。
✅ 弹性表现受网格数量影响(难以控制)、没有所谓的精确解，哪怕迭代数足够多、迭代数过多会导致locking issue.

P13

After-Class Reading

Muller. 2008. Hierarchical Position Based Dynamics. VRIPHYS.

✅ NVIDIA的很多物理引擎都是基于PBD的

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P14

Strain Limiting

原理

Strain limiting aims at using the projection function for correction only.

✅ 投影函数作为模拟过程的后处理，防止模拟后产生大的形变，使模拟结果更稳定。

P15

例一： Spring Strain Limit

We can set the spring strain, i.e., the stretching ratio σ, to be within a limit.

✅ 这一页以弹簧为例子说明 Strain Limit
✅ Strain：物理上描述形变的量，即本页的\(\sigma \)

$$ \sigma ^\mathrm{{min}}≤\frac{1}{L}||\mathbf{x} _i− \mathbf{x} _j||≤\sigma^\mathrm{{max}} $$

✅ 仅要求弹簧长度满足某不比例，不要求一定到达某个位置。

Constraint

P16

$$ \quad $$

$$ \mathbf{x}^{\mathbf{new}}\longleftarrow \mathrm{Projection} (\mathbf{x} ) $$

$$ \sigma\longleftarrow \frac{1}{L}||\mathbf{x}_i− \mathbf{x}_j|| $$

✅ 计算当前拉伸比

$$ \sigma _0\longleftarrow \mathrm{min} (\mathrm{max} (\sigma,\sigma^{\mathrm{min} }),\sigma^{\mathrm{max} }) $$

✅ 计算期望的拉伸比

$$ \mathbf{x} _i^{\mathrm {new} }⟵\mathbf{x} _i−\frac{m_j}{m_i+m_j}(||\mathbf{x} _i− \mathbf{x} _j||−σ_0L)\frac{\mathbf{x} _i− \mathbf{x} _j}{||\mathbf{x} _i− \mathbf{x} _j||} $$

$$ \mathbf{x} _j^{\mathrm {new} }⟵\mathbf{x} _j+\frac{m_j}{m_i+m_j}(||\mathbf{x} _i− \mathbf{x} _j||−σ_0L)\frac{\mathbf{x} _i− \mathbf{x} _j}{||\mathbf{x} _i− \mathbf{x} _j||} $$

✅ 用\(\sigma _0L\)代替原长\(L\)．

\(\mathrm{PBD}: \sigma _0≡1;\quad\quad\)

✅ PBD可以看作是Strain Limit的特例。

No limit: \(\sigma ^{\mathrm{min} } = 0, σ^{\mathrm{max} } = \infty\)

✅ Strain Limit 的应用场景：(1) 模拟布料：“拉伸到一定范围后变得非常 stiff” 的效果 (2) 防止“形变大发生数值不稳定”。

P17

例二：Triangle Area Limit

We can limit the triangle area as well. To do so, we define a scaling factor.

✅这是另一个例子。希望顶点移动尽量少，因此定义约束：三角形面积变化在一定范围内。

{\(\mathbf{x}_i^{\mathrm{new} },\mathbf{x}_i^{\mathrm{new} },\mathbf{x}_k^{\mathrm{new} }\)}  = \(\mathrm{argmin} \frac{1}{2} \){\(m_i||\mathbf{x}_i^{\mathrm{new} }−\mathbf{x}_i||^2+m_j||\mathbf{x}_j^{\mathrm{new} }−\mathbf{x}_j||^2+m_j||\mathbf{x}_k^{\mathrm{new} }−\mathbf{x}_k||^2\)} 

such that the constraint is satisfied.

✅ strain s为面积的缩放量

P18

$$ \quad $$

$$ \mathbf{x} ^{\mathbf{new}} ⟵\mathrm{Projection} (\mathbf{x} ) $$

$$ \mathbf{A}\longleftarrow \frac{1}{2} ||(\mathbf{x} _j− \mathbf{x} _i)\mathbf{×} (\mathbf{x}_k− \mathbf{x} _i)|| $$

✅ 计算当前三角形的面积

$$ \mathbf{s} \longleftarrow \sqrt{\mathrm{\mathrm{min}} (\mathrm{\mathrm{max}} (A,A^{\mathrm{min}}),A^{\mathrm{max}})/A} $$

✅ 计算期望的面积缩放比

$$ \mathbf{c} \longleftarrow \frac{1}{m_i+m_j+m_k} (m_i\mathbf{x} _i+m_j\mathbf{x} _j+m_k\mathbf{x} _k) $$

✅ C为质心，要求缩放前后质心不变。数学上、质心不变，点的移动最少；物理上，质心变了代表物体运动了，scale 不应该导致物理运动。

$$ \mathbf{x} _i^{\mathrm{new}}\longleftarrow \mathbf{c} +s(\mathbf{x} _i−\mathbf{c} ) $$

$$ \mathbf{x} _j^{\mathrm{new}}\longleftarrow \mathbf{c} +s(\mathbf{x} _j−\mathbf{c} ) $$

$$ \mathbf{x} _k^{\mathrm{new}}\longleftarrow \mathbf{c} +s(\mathbf{x} _k−\mathbf{c} ) $$

✅ 通过对顶点到质心的距离的缩放，得到顶点的新的位置

P19

Strain Limiting在Simulation中的作用

• Strain limiting is widely used in physics- based simulation, typically for avoiding instability and artifacts due to large deformation.

• Strain limiting is useful for nonlinear effects, in a biphasic way.

✅ 例如布料一般一开始抵抗比较小，拉到一定程度后抵抗迅速变大。对这种非线性的表现，可以把模拟分布两个阶段，前面用普通模拟，后面用strain limiting。

• Strain limiting also helps address the locking issue.

✅两个阶段有不同的算法，针对两个阶段的不同特点，可以分别解决两个阶段的问题。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P22

Projective Dynamics

原理

Instead of blending projections in a Jacobi or Gauss-Seidel fashion as in PBD, projective dynamics uses projection to define a quadratic energy.

✅ PBD方法直接拿约束来修复顶点位置，没有物理含义。而本页方法把projection方法跟物拟模拟结合起来，主要体现在用约束来做什么。

能量和力

\(E(\mathbf{x} ) = {\textstyle \sum_{e = (i,j)}}\frac{1}{2} ||(\mathbf{x} _i−\mathbf{x} _j)−(\mathbf{x} _{e,i}^{\mathrm{new} }−\mathbf{x} _{e,j}^{\mathrm{new} })||^2\)

{\(\mathbf{x} _{e,i}^{\mathrm{new} },\mathbf{x} _{e,j}^{\mathrm{new} }\)} = \(\mathrm{Projection} _e(\mathbf{x}_i,\mathbf{x}_j)\) for every edge \(e\)

✅ 本文基于约束定义能量。{\(\mathbf{x} _{e,i}^{\mathrm{new} },\mathbf{x} _{e,j}^{\mathrm{new} }\)}为期望的顶点位置。不直接把顶点从当前位置移到期望位置。而是把当前位置和期望位置的距离转化为能量，通过能量推动顶点从当前位置移到目标位置。

$$ E(\mathbf{x})=\sum _{e=(i,j)}\frac{k}{2}(||\mathbf{x}_i-\mathbf{x}_j||-L_e)^2 $$

$$ \mathbf{f} _i=−\nabla_iE(\mathbf{x} )=−{\textstyle \sum _{e:i\in e}}(\mathbf{x} _i−\mathbf{x} _j)−(\mathbf{x} _{e,i}^{\mathrm{new}} −\mathbf{x} _{e,j}^{\mathrm{new} }) $$

✅ 基于 \(E(\mathbf{x})、\mathbf{x}_i^{\mathrm{new} } 、\mathbf{x}_j^{\mathrm{new} }\) 计算力，此时假设\(\mathbf{x}_i^{\mathrm{new} }\)和 \(\mathbf{x}_j^{\mathrm{new} }\)都是定值，\(\mathbf{x}_i\)和 \(\mathbf{x}_j\)是变量。
✅ 本文基于约束定义能量和力，得到的结果与基于弹簧能量计算的能量和力相同。
✅ 既然 \(E\) 和 \(F\) 是一样的，何必多此一举? 答：H不同。

P23

Hessian 矩阵

Instead of blending projections in a Jacobi or Gauss-Seidel fashion as in PBD, projective dynamics uses projection to define a quadratic energy.

✅ 同一个顶点在三个不同边上的投影是不同的。
✅ 可以直接根据Mesh的拓扑关系构造H矩阵。
✅ 为什么能简化\(\mathbf{H}\)的计算？答：在计算某一个端点时，假设另一个端点不动（常量），那么能量就是只关于这个端点的二次函数

P24

Shape Matching

Shape matching is also projective dynamics, if we view rotation as projection:

|

The 2D Space	The 3D Space

Assuming that \(\mathbf{{\color{Orange} R} }\) is constant,
$$ \begin{matrix} \mathbf{f} _0=−\nabla_0E(\mathbf{x} )\\ \mathbf{f} _1=−\nabla_1E(\mathbf{x} ) \\ \mathbf{f} _2=−\nabla_2E(\mathbf{x} )\\ \mathbf{H} =\frac{∂E^2(\mathbf{x} )}{∂x^2} \quad \text{is a constant !}
\end{matrix} $$

P25

Simulation by Projective Dynamics

• According to implicit integration and Newton’s method, a projective dynamics simulator looks as follows, with matrix \(\mathbf{A} =\frac{1}{∆t^2}\mathbf{M+}\mathbf{H} \) being constant.

• We can use a direct solver with only one factorization of A.

✅ 解线性系统的主要耗时在LU分解，而这个算法中\(\mathrm{H}\)是常数矩阵，只需要做一次LU分解，简化了对\(\mathrm{H}\)分解的计算量。

Initialize \(\mathbf{x} ^{(0)}\), often as\( \mathbf{x} ^{[0]} \)or \(\mathbf{x} ^{[0]} +∆t\mathbf{v} ^{[0]} \)

For \(k=0\dots K\)
\(\quad\) Recalculate projection

✅ 对于弹簧系统，Recaculate projection 这一步实际上不需要，因为直接用弹簧系统的公式算力，得到的\(f\)是一样的。
✅ 如果是做 shape matching, 还是需要这一步，用于算 \(f\)

\(\quad\) Solve \((\frac{1}{∆t^2}\mathbf{M} +\mathbf{H} )∆\mathbf{x} =−\frac{1}{∆t^2}\mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})+\mathbf{f} (\mathbf{x} ^{(k)})\)

\(\quad\) \(\mathbf{x} ^{(k+1)}\longleftarrow \mathbf{x} ^{(k)}+∆\mathbf{x} \)

\(\quad\) If \(||∆\mathbf{x}||\) is small \(\quad\) then break

\(\mathbf{x} ^{[1]}\longleftarrow \mathbf{x} ^{(k+1)}\)

\(\mathbf{v} ^{[1]}\longleftarrow (\mathbf{x} ^{[1]}-\mathbf{x} ^{[0]})/∆t\)

“Newton’s Method”

$$ \quad $$

P26

Preconditioned Steepest Descent

• Mathematically, this approach is preconditioned steepest descent, in which:

$$ F(\mathbf{x} )=\frac{1}{2∆t^2} ||\mathbf{x} −\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]}||_\mathbf{M} ^2+E(\mathbf{x} ) $$

The performance depends on how well \(\mathbf{{\color{Orange} H} }\) approximates the real Hessian.

✅\(\mathrm{H}\)不需要很精确，一个近似的正定的矩阵，就能让结果收敛。

P27

Pros and Cons of Projective Dynamics

• By building constraints into energy, the simulation now has a theoretical solution with physical meaning.
• Fast on CPUs with a direct solver. No more factorization!

✅ Fast on CPU,因为它只作一次\(\mathbf{LU}\)分解。

• Fast convergence in the first few iterations.
• Slow on GPUs. (GPUs don’t support direct solver wells.)

✅ Slow on GPU，因为\(\mathbf{LU}\)分解不适用于 \(\mathbf{GPU}\)

• Slow convergence over time, as it fails to consider Hessian caused by projection.
  • Still suffering from high stiffness
• Cannot easily handle constraint changes.

✅ constraint changes: 网格关系改变导至弹簧结构改变，原来的\(\mathbf{H}\)将不再适用。
- Contacts
- Remeshing due to fracture, etc.

✅ 课后答疑：
能量优化的方法很少用于刚体，主要是有限元、弹性体、衣服模拟。

P28

After-Class Reading

Bouaziz et al. 2014. Projective Dynamics: Fusing Constraint Projections for Fast Simulation. TOG (SIGGRAPH).

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P30

Constrained Dynamics

要解决的问题

A critical problem exists: what if constraints/forces are very very stiff? Or infinitely stiff?

✅ 此算法用于处理 very very stiff 的场景即约束必须严格满足，而前面算法需要做很多次迭代才能产生这种效果（计算量大）。
✅ 此算法常用于衣服、刚体、人体。比如人体的关节联结，是一种非常stiff的约束。

根据约束建立模型

Compliant constraint

$$ \mathbf{ϕ} _e(\mathbf{x} )=||\mathbf{x} _{ei}− \mathbf{x} _{ej}||−L_e $$

$$ E(\mathbf{x} )={\textstyle \sum_e}\frac{1}{2} k(||\mathbf{x} _{ei} −\mathbf{x} _{ej}||−L_e)^2=\frac{1}{2} \mathbf{\mathbf{ϕ}^T(x)C} ^{−1}\mathbf{ϕ} (\mathbf{x} ) $$

✅ \(E(X)=\sum _ {e} \frac{1}{2}k( \phi _ e(x))^2\)

$$ \mathbf{f} (\mathbf{x} )=−∇E=-\begin{pmatrix} \frac{∂E}{∂\mathbf{ϕ}} & \frac{∂\mathbf{ϕ}}{∂x} \end{pmatrix}^\mathbf{T} =−\mathbf{J^TC} ^{−1}\mathbf{ϕ} =\mathbf{J^Tλ} $$

Let \(N\) be the number of vertices and E be the number of constraints,

✅ \(E\) 和 \(f\) 变成了关于两个变量\((x、\lambda )\)的函数。
✅把能量写成约束的形式
$$ E(x)=\frac{1}{2}\phi ^T (x)C^{-1}\lambda \\ f(x)=J^T \lambda
$$ 根据\(E(x)\)和\(f(x)\)做隐式积分

转化为隐式积分问题

P31

By implicit integration, we get:

$$ \mathbf{Mv} ^{\mathrm{new} }−∆t\mathbf{J^Tλ} ^{\mathrm{new} }=\mathbf{Mv} $$

✅ 动量守衡公式：\(Mv'- Mv = Ft = \)冲量
✅ 此处新 \(\lambda^{\mathrm{new}} \)来计算 F. 说明是 Implicit

Meanwhile, $$ \mathbf{Cλ} ^{\mathrm{new} }=−\mathbf{ϕ} ^{\mathrm{new} }≈−\mathbf{ϕ} −\mathbf{J} (\mathbf{x} ^{\mathrm{new} }−\mathbf{x} )≈−\mathbf{ϕ} −∆t\mathbf{Jv} ^{\mathrm{new} } $$

✅ 对\(-\phi ^{\mathrm{new}}\) 的泰勒展开
✅\(J\)是上页中的Jacobian.

$$ \begin{bmatrix} \mathbf{M} & −∆t\mathbf{J^T} \\ ∆t\mathbf{J} & \mathbf{C} \end{bmatrix}\begin{bmatrix} \mathbf{v} ^{\mathrm{new} } \\ \mathbf{λ} ^{\mathrm{new} } \end{bmatrix}\begin{bmatrix} \mathbf{Mv} \\ -\mathbf{ϕ} \end{bmatrix} $$

✅ 最后的矩阵公式由上面两个公式整理合并得到。
\(x ^{\mathrm{new}}-x =\bigtriangleup t\cdot v \)

解隐式积分

P32
Now we have a system with two sets of variables: the primal variable \(\mathbf {x}\) (or \(\mathbf {v=x}\) ̇) and the dual variable \(\mathbf {λ}\).

✅ 之前的隐式积分只有一个变量，此处的隐式积分有两个变量。需要同时解出两个变量。

• Method 1: We can solve the two variables by a direct solver together, in a primal-dual fashion:

$$ \begin{bmatrix} \mathbf{M} & −∆t\mathbf{J^T} \\ ∆t\mathbf{J} & \mathbf{C} \end{bmatrix}\begin{bmatrix} \mathbf{v} ^{\mathrm{new} } \\ \mathbf{λ} ^{\mathrm{new} } \end{bmatrix} = \begin{bmatrix} \mathbf{Mv} \\ -\mathbf{ϕ} \end{bmatrix} $$

✅ 注意：Method 1 中的矩阵有可能不正定、因此很多数学方法用不了。

• Method 2: We can reduce the system by Schur complement and solve \(\mathbf {λ}^{\mathrm{new} }\) first.

✅ Method 2 的消元过程不容易构造、尤其是当矩阵比较复杂时。

$$ (∆t^2\mathbf{JM} ^{−1}\mathbf{J} ^\mathbf{T} +\mathbf{C} )\mathbf{λ} ^{\mathrm{new} } =−\mathbf{ϕ} −∆t\mathbf{Jv} $$

$$ \mathbf{v} ^{\mathrm{new}}\longleftarrow \mathbf{v} +−∆t\mathbf{M} ^{−1}\mathbf{J^Tλ} ^{\mathrm{new}} $$

✅ 用哪种方法取决于应用场景

优点

Infinite stiffness? \(\mathbf{C \longrightarrow 0}\).

✅ 此方法将软硬度量解耦出来，并用矩阵C来表示，使得可以方便控制软硬度，例如让\(c=0\) 来表示 infinite stiffness.

应用

Articulated Rigid Bodies (ragdoll animation)

P34

Stable Constrained Dynamics

✅ 没有展开讲，见after reading中的论文

From a mass-spring system, we know spring Hessian (tangent stiffness) is:

$$ \mathbf{H} (\mathbf{x} )=∑_{e=(i,j)}\begin{bmatrix} \Box & \Box & \Box & \Box \\ \Box & \mathbf{H} _e & -\mathbf{H} _e & \Box \\ \Box & -\mathbf{H} _e &\mathbf{H} _e & \Box \\ \Box & \Box & \Box & \Box
\end{bmatrix} $$

According to constrained dynamics:\(\mathbf{f} (\mathbf{x} )=\mathbf{J^Tλ}\) and \(\mathbf{λ} =−\mathbf{C} ^{−1}\mathbf{ϕ} \), so:

$$ \mathbf{J}_e=\frac{∂\phi _e}{∂\mathbf{x} }=\begin{bmatrix} \frac{\mathbf{x} _{ij}^\mathbf{T} }{||\mathbf{x} _{ij}||} & -\frac{\mathbf{x} _{ij}^\mathbf{T} }{||\mathbf{x} _{ij}||} \end{bmatrix} $$

P35

Stable Constrained Dynamics

According Lecture 5, Page 16, implicit integration is:

$$ (\frac{1}{∆t^2}\mathbf{M+H} (\mathbf{x} ^{[0]}))∆\mathbf{x} = \frac{1}{∆t^2}\mathbf{M} (∆t\mathbf{v} ^{[0]})+\mathbf{f} (\mathbf{x} ^{[0]}) $$

$$ \Downarrow $$

$$ (\mathbf{M} +∆t^2\mathbf{H} (\mathbf{x} ^{[0]}))\mathbf{v} ^{\mathrm{new} }= \mathbf{Mv} ^{[0]}+∆t\mathbf{f} (\mathbf{x} ^{[0]}) $$

Missing geometric stiffness matrix here…

P36

After-Class Reading (optional)

Tournier et al. 2015. Stable Constrained Dynamics. TOG (SIGGRAPH).

✅ Method 2 Gauss 消元，如果把\(\lambda\)消掉，会得到一个基本上与隐式积分相似的公式，唯一的区别是\(\mathbf{H}\)上略有不同 —— 隐式积分多了一项。如果把用这一项加回去，会使constrain dynamic 变得稳定。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P37

A Summary For the Day

• Position-based dynamics and strain limiting
  • The key is to build a projection function for every constraint.
  • Two approaches for integration: Jacobi and Gauss-Seidel.
  • Fast in low resolutions, but problematic in high resolutions.
  • Not physically correct.
• Projective Dynamics
  • Also uses projection functions, but they are now built into energies.
  • In every iteration, projections are first updated, and then treated as constants in implicit formulation.
  • The matrix in the system becomes constant, can be pre-factorized for fast simulation.
  • Converges fast only in the first few iterations, slow afterwards. CPU friendly.
• Constrained Dynamics
  • Focused on very stiff constraints. Introduces dual variables.
  • Also built upon implicit integration. Two methods: primal-dual, pure dual.
  • Restrictions on the solvers.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• The linear finite element method (FEM)
• The finite volume method (FVM)
• Hyperelastic models

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P3

Linear Finite Element Method

P4

The Linear FEM Assumption

✅ 假设：三角形的形变是均匀的

In a nutshell, linear FEM assumes that for any point \(\mathbf{X}\) in the reference triangle, its deformed correspondence is: \(\mathbf{x=FX+c}\).

✅ reference triangle：三角形处于没有发生形变的静止的状态。
✅ \(\mathbf{X}\)和\(\mathbf{x}\)可以分别是 reference 和 deformed 三角形的顶点或内部点，公式都同样适用。

For any vector between two points, we can use F to convert it from reference to deformed:
$$ \mathbf{x} _{ba}=\mathbf{x} _b−\mathbf{x} _a=\mathbf{FX} _b+\mathbf{c} −\mathbf{FX} _a−\mathbf{c} =\mathbf{FX} _{ba}. $$

P5

计算Deformation Gradient

Therefore, we can calculate the deformation gradient by edge vectors.

Problem: \(\mathbf{F}\) is related to deformation, but it contains rotation.

✅ 期望\(\mathbf{F}\)只包含形变量、不包含平移和旋转、因为刚体运动不应该有形变，所以要把形变提取出来。
✅平移已经在\(\mathbf{c}\)里面了，所以只需考虑旋转。

P6

从F中去除旋转

回顾SVD分解

Ideally, we need a tensor to describe shape deformation only. Recall that SVD gives \(\mathbf{F=UDV^T}\), where only \(\mathbf{V^T}\) and \(\mathbf{D}\) are relevant to deformation.

✅ \(\mathbf{V^T}\) 看上去是旋转、实际上是为了确定形变的方向、 \(\mathbf{U}\) 才是真正的旋转
✅ 目的：把\(\mathbf{F}\)中的\(\mathbf{U}\)去掉、可以先做 \(\mathbf{SVD}\) 分解再把\(\mathbf{U}\)去掉。但本文使用了更简单的方法

Green Strain

So we get rid of \(\mathbf{U}\) as: \(\mathbf{G} =\frac{1}{2} (\mathbf{F^TF−I} )=\frac{1}{2} (\mathbf{VD} ^2\mathbf{V} ^\mathbf{T} −\mathbf{I} )=\begin{bmatrix} \varepsilon _{uu} & \varepsilon _{uv}\\ \varepsilon _{uv} & \varepsilon _{vv} \end{bmatrix}\), Green strain.

✅ \(\mathbf{G}\) 是一个描述物体形变的有无和大小矩阵，且与关旋转

• If no deformation, \(\mathbf{G=0}\); if deformation increases, ||\(\mathbf{G}\)|| increases.
• Three deformation modes: \(\varepsilon _{uu}\), \(\varepsilon _{vv}\) and \(\varepsilon _{uv}\).
• \(\mathbf{G}\) is rotation invariant: if additional rotation \(\mathbf{R}\), then deformation gradient is \(\mathbf{RF}\) but green strain is the same: \(\mathbf{G} =\frac{1}{2} (\mathbf{F^TR^TRF−I} )=\frac{1}{2} (\mathbf{VD} ^2\mathbf{V} ^\mathbf{T} −\mathbf{I} )\).

P7

计算能量

Let \(\mathbf{G}\) be the the green strain describing deformation. We consider the energy density per reference area as: \(W (\mathbf{G})\).

✅用形变程度来定义能量。\(W\)代表单位面积上的能量，因此称为能量密度。总能量为单位能量\(\mathbf{X}\)面积.
✅ \(A^{ref}\) 为 reference status 下三角形的面积
✅ StVK是一种经典的能量密度函数（Strain Energy Density Function）， 在力学中不常用，但在图形学中很常用、原因是简单

✅ S 是一个与力有关的物理量。会在FVM内容中进一步解释。
✅ 能量对位移求导是力。形变是一种位移。能量密度对位移求导是一种类似于力的密度的量。

P8

计算Forces

Given everything we have, we can now calculate the forces.

✅ 力是形变施加到顶点上的力

✅ 绿色部分是上一页S中的内容、灰色部分将在下一页推导。

P9

方法一

Recall that,

❗ \(\mathbf{F}\)不是力，是deformation gradient．
✅ 假设a,b,c,d是形变后的顶点。

By definition,
$$ \mathbf{G} =\frac{1}{2} (\mathbf{F^TF−I} )=\begin{bmatrix} \frac{1}{2}(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} (a\mathbf{x} _{10}+c\mathbf{x} _{20})−\frac{1}{2} & \frac{1}{2}(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} (b\mathbf{x} _{10}+d\mathbf{x} _{20})\\ \frac{1}{2}(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} (b\mathbf{x} _{10}+d\mathbf{x} _{20}) & \frac{1}{2}(b\mathbf{x} _{10}+d\mathbf{x} _{20})^\mathbf{T} (b\mathbf{x} _{10}+d\mathbf{x} _{20})−\frac{1}{2} \end{bmatrix} $$

So:

$$ \frac{∂\varepsilon _{uu}}{∂\mathbf{x} _1}=a(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} \quad\quad \frac{∂\varepsilon _{vv}}{∂\mathbf{x} _1}=b(b\mathbf{x} _{10}+d\mathbf{x} _{20})^\mathbf{T} \quad\quad \frac{∂\varepsilon _{uv}}{∂\mathbf{x} _1}=\frac{1}{2} a(b\mathbf{x} _{10}+d\mathbf{x} _{20})^\mathbf{T} +\frac{1}{2} b(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} $$

$$ \frac{∂\varepsilon _{uu}}{∂\mathbf{x} _2}=c(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} \quad\quad \frac{∂\varepsilon _{vv}}{∂\mathbf{x} _2}=d(b\mathbf{x} _{10}+d\mathbf{x} _{20})^\mathbf{T} \quad\quad \frac{∂\varepsilon _{uv}}{∂\mathbf{x} _2}=\frac{1}{2} c(b\mathbf{x} _{10}+d\mathbf{x} _{20})^\mathbf{T} +\frac{1}{2} d(a\mathbf{x} _{10}+c\mathbf{x} _{20})^\mathbf{T} $$

✅ \(\mathbf{x}\)为current边的矩阵，\(\mathbf{r}\)为reference边的矩阵。

P10

方法二

✅ 把 P9 代入 P8 得到 P10
✅ 上一页推导方法从定义出来，过程简单，但很容易出错。这里用矩阵来简化计算，得到同样的结果。

P11

结论

In conclusion, we have:

$$ \mathbf{f} _1=−A^{\mathrm{ref} }\mathbf{FS} \begin{bmatrix} a\\ b \end{bmatrix} \quad\quad \mathbf{f} _2=−A^{\mathrm{ref} }\mathbf{FS} \begin{bmatrix} c\\ d \end{bmatrix} $$

$$ \begin{bmatrix} \mathbf{f} _1 &\mathbf{f} _2 \end{bmatrix}= − A ^{\mathrm{ref} }\mathbf{FS} \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} \end{bmatrix}^\mathbf{−T} $$

✅ \(f_0=-f_1-f_2\)

P12

Implementations

🔎 Volino et al. 2009. A simple approach to nonlinear tensile stiffness for accurate cloth simulation. TOG

Only talks about cloth (2D reference -> 3D deformation)

• What about tetrahedron (3D reference -> 3D deformation)?
  • Same idea, but everything is now in 3D.
  • Deformation gradient \(\mathbf{F} \in \mathbf{R} ^{3×3}\)
  • Green strain \(\mathbf{G} \in \mathbf{R} ^{3×3}\)
  • Stress tensor \(\mathbf{S} \in \mathbf{R} ^{3×3}\)
  • Forces \(\mathbf{F}_i \in \mathbf{R} ^3\)

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P13

Finite Volume Method

✅ FEM求导，FVM积分。在线性场景下，FEM和FVM本质上等价的。
✅ FVM基于力从何而来的思想。

符号定义

P14

Traction

First, let’s consider traction t: the internal force per unit length (area).

✅ 两个弹性体被界面 \(L\) 分开、求 \(L\) 上的力、

Total interface force:

$$ f\mathbf{} =\oint _L \mathbf{t} dl $$

✅ \(\mathbf{t}\)是\(L\) 上的单位面积/长度上的力。那么总的力是t的积分。

Stress

Stress tensor \(\mathbf{σ} \) (interface normal -> traction):

$$ \mathbf{t=σn} $$

So,

$$ \mathbf{f} =\oint_{L} \mathbf{σn} dl $$

✅ \(\mathbf{σ}\) 一个矩阵

P15

计算FVM中的力

2D

FVM considers force calculation in an integration perspective, not a differentiation perspective.

✅ \(\mathbf{X}_0\)是顶点\(\mathbf{X}_0\)附近邻域面积上的力。

Force contributed by an element:

$$ \mathbf{f}_0 =\oint _L \mathbf{σn} dl $$

✅ \(\mathbf{X}_0\)上的力是邻域面边界\(L\)上的力的积分、不考虑边界内部的力，因为认为内力为0。
✅ 仅看其中一个三角形、假设曲线经过 \(\mathbf{X}_0\mathbf{X}_1\)和 \(\mathbf{X}_0\mathbf{X}_2\)的中点。因为三角形的力对三个顶点是平均的。

Since \(\mathbf{σ}\) is constant within the element,

$$ \oint_L \mathbf{σn} dl + \oint_{L_{20}} \mathbf{σn} dl+\oint_{L_{10}}\mathbf{σn} dl=0 $$

✅ 对于封闭曲线L+L10+L20做积分，\(\int _n=0\)，因此 \(\sigma \int _n=0\)

(Divergence Theorem)

We know the force is:

$$ \mathbf{f}_0 = - \oint _ {L _{20}} \mathbf{σn} _ {20} dl - \oint _ {L _{10}} \mathbf{σn} _ {10} dl = - \mathbf{σ}(\frac{||\mathbf{X} _ {20}||}{2}\mathbf{n} _ {20}+ \frac{||\mathbf{X} _ {20}||}{2}\mathbf{n} _ {10}) $$

P16

3D

✅ 三维场景是对四面体的四个面积分。
✅ 每个三角形的 stress 都不同、同一个三角形内部 stress 是常数。

In 3D, FVM works in the same way.

Force:

$$ \begin{array}{l} \mathbf{f} _ 0 = −\oint _ Ω \mathbf{σn} dA=−\mathbf{σ} (\frac{A _ {012}}{3}\mathbf{n} _ {012} + \frac{A _ {023}}{3}\mathbf{n} _ {023} + \frac{A _ {031}}{3}\mathbf{n} _ {031})\\ =−\frac{σ}{3}(\frac{||\mathbf{x} _ {10} \times \mathbf{x} _ {20}||}{2} \frac{\mathbf{x} _ {10} × \mathbf{x} _ {20}}{||\mathbf{x} _ {10} × \mathbf{x} _ {20}||} + \frac{||\mathbf{x} _ {20} × \mathbf{x} _ {30}||}{2}\frac{\mathbf{x} _ {20} × \mathbf{x} _ {30}}{||\mathbf{x} _ {20} × \mathbf{x} _ {30}||} +\frac{||\mathbf{x} _ {30} × \mathbf{x} _ {10}||}{2}\frac{\mathbf{x} _ {30} × \mathbf{x} _ {10}}{||\mathbf{x} _ {30} × \mathbf{x} _ {10}||})\\ =−\frac{\mathbf{σ}}{6} (\mathbf{x} _ {10} × \mathbf{x} _ {20} + \mathbf{x} _ {20} × \mathbf{x} _ {30} + \mathbf{x} _ {30} × \mathbf{x} _ {10}) \end{array} $$

❓ 遗留问题， stress 如何计算？

✅\(f_0\)是\(\sigma n\)在绿色体截面上的积分。
✅类似于上一页合力为零的原理，\( \oint 截面＋\oint 表面=0\)
✅“面积／3”是因为面上的贡献均匀地分布到三个点上。

P17

计算FVM中的stress

Although the use of stress tensor is the same: mapping from the interface normal to the traction, it can be defined by different configurations.

In FEM, we define the energy density \(W\) in the reference state. Therefore, this stress \(\mathbf{S}\) is a mapping from the normal \(\mathbf{N}\) to the traction \(\mathbf{T}\), both in the reference state.	In FVM, we need \(\mathbf{σ}\) to convert the normal into \(\mathbf{t}\) for force calculation. Therefore, this stress assumes the normal \(\mathbf{n}\) and the traction \(\mathbf{t}\) are in the deformed state.

✅ 在 reference 状态下有 normal. traction 和 stress.在形变状态下也有 normal traction 和 stress.FEM 使用的是 reference 空间下的量。

P18

Different Stresses

We can now have different stresses, serving the same purpose but in different forms.

✅ FVM 需要的是 Cauchy Stress.（\(\sigma \)）、上节课讲了(S)的 计算方法，需要根据(S)求(\(\sigma \)).
\(P → \sigma \) 的过程没有展开讲，结论在P21

P19

P与\(\sigma \)的关系：Area Weighted Normals

Now let’s figure out the relationship between \(A^{\mathrm{ref} }\mathbf{N}\) and \(A\mathbf{n}\), the two area weighted normals.

$$ 2A^{\mathrm{ref} }\mathbf{N} =\mathbf{X} _ {a0}×\mathbf{X} _ {b0} $$

$$ \begin{array}{l} 2A\mathbf{n} =\mathbf{x} _ {a0}×\mathbf{x} _ {b0}=\mathbf{FX} _{a0} × \mathbf{FX} _ {b0} = (\mathbf{UDV^TX} _ {a0}) × (\mathbf{UDV^TX} _ {b0}) \\ \quad\quad=\mathbf{U} ((\mathbf{DV^TX} _ {a0}) × (\mathbf{DV^TX} _ {b0})) \\ \quad\quad=\mathbf{U} \begin{bmatrix} d_1d_2& \Box & \Box \\ \Box & d_0d_2 & \Box \\ \Box & \Box &d_0d_1 \end{bmatrix} ((\mathbf{V^TX} _ {a0})×(\mathbf{V^TX} _ {b0}))\\ \quad\quad=\mathbf{U} \begin{bmatrix} d_1d_2& \Box & \Box \\ \Box & d_0d_2 & \Box \\ \Box & \Box &d_0d_1 \end{bmatrix} \mathbf{V^T} (\mathbf{X} _ {a0} × \mathbf{X} _ {b0})\quad=d_0d_1d_2\mathbf{U} \begin{bmatrix} 1/d_0& \Box & \Box \\ \Box & 1/d_1 & \Box \\ \Box & \Box &1/d_2 \end{bmatrix} \mathbf{V^T} (\mathbf{X} _ {a0}×\mathbf{X} _ {b0}) \\ \quad\quad=\mathrm{det} (\mathbf{F} )\mathbf{F} ^{−\mathbf{T}}(\mathbf{X} _ {a0}×\mathbf{X} _ {b0})=\mathrm{det} (\mathbf{F} )\mathbf{F} ^{−\mathbf{T}}(2A^\mathrm{ref}\mathbf{N}) \end{array} $$

P20

Now we know: \(A\mathbf{n} =\mathrm{det} (\mathbf{F})\mathbf{F^{−T}} (A^{\mathrm{ref}}\mathbf{N} )\).

We also know the force can be calculated using two different stresses:

$$ \mathbf{f} =−\frac{1}{3} {\textstyle \sum {A^{\mathrm{ref} }}}\mathbf{PN} =−\frac{1}{3}\sum A\mathbf{σn} $$

Therefore, we get:

$$ \mathbf{P} (A^{\mathrm{ref} }\mathbf{N} )=\mathbf{σ} \mathrm{det} (\mathbf{F} )\mathbf{F^{−T}} (A^{\mathrm{ref} }\mathbf{N} ) $$

$$ \mathrm{det} ^{−1}(\mathbf{F} )\mathbf{PF^T=σ} $$

P21

结论

We can now have different stresses, serving the same purpose but in different forms.

P22

根据stress算出力

The previous analysis suggests we can use reference normals instead.

Second Piola–Kirchhoff stress:
\(\mathbf{S} =\frac{∂W}{∂\mathbf{G}}\), as in previous FEM formulation

✅ 第一行公式：用 deformed position 计算 deformed position. 第二行公式：用 ref position 计算 deformed position, 因此直接把\(\sigma \)换成 \(\mathbf{P} \) 就可以。
✅ 好处：ref position 是常数，可以做预计算、并存储为\(b_1\).
✅ F：deformation gradient.见P5
✅ 用三种不同定义的 stress 来算力、目的是得到计­算最友好的公式
✅ 此处内容涉及材料力学、

P23

关于b1

✅ 问： \(\mathbf{X}_ {20}^\mathbf{T} b_1\)的计算公式中、绿色的\(\mathbf{X}_ {01}×\mathbf{X} _ {21}\)怎么变成了\(\mathbf{X}_ {20}\times \mathbf{X}_ {10}\)？下面的\(\mathbf{X}_ {30}^\mathbf{T} b_1\),也一样。
答：因为\(\mathbf{X}_o，\mathbf{X}_1，\mathbf{X}_2\)是同一个三角形上的顶点，任意两条边做cross都是一样的，处理好正负就好了。
也可以用cross的乘法分配律得出相同的结论

P24
Therefore, We get:

$$ \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^\mathbf{T} \mathbf{b} _1=\begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^\mathbf{T}(\mathbf{X} _{01}×\mathbf{X} _{21}+\mathbf{X} _{21}×\mathbf{X} _{31}+\mathbf{X} _{31}×\mathbf{X} _{01})=\begin{bmatrix} 6Vol\\ 0\\ 0 \end{bmatrix} $$

$$ \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^\mathbf{T} \mathbf{b} _2=\begin{bmatrix} 0\\ 6Vol\\ 0 \end{bmatrix} $$

$$ \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^\mathbf{T} \mathbf{b} _3=\begin{bmatrix} 0\\ 0\\ 6Vol \end{bmatrix} $$

$$ \begin{bmatrix} \mathbf{b} _{1} & \mathbf{b} _{2} &\mathbf{b} _{3} \end{bmatrix}^\mathbf{T} =6Vol\begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^{-\mathbf{T}}\\ \quad\quad=\frac{1}{\mathrm{det}( \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^{-1})} \begin{bmatrix} \mathbf{X} _{10} & \mathbf{X} _{20} &\mathbf{X} _{30} \end{bmatrix}^{-\mathbf{T}} $$

P25

A Quick Summary

✅reference 状态下，\(F＝I,G=0,P=0,f=0\)

P26

After-Class Reading

重点推荐：

Teran et al. 2003. Finite Volume Methods for the Simulation of Skeleton Muscles. SCA.

✅ 这篇论文重点推荐，但论文中的公式与 PPT 中的不完全一样. PPT 上的又进一步简化。

Optional

Volino et al. 2009. A Simple Approach to Nonlinear Tensile Stiffness for Accurate Cloth Simulation. TOG.

✅2D有限元

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P28

Hyperelastic Models

✅ 前面的内容，都假设使用 StVK 材料、优点是简单；缺点是无法处理反转。因此在材料力学中不常用。
Hyperplasia 利用能量密度(W)、提供一个从 Strain (G) 到 Stress (S)的映射

跳过

P29

First Piola–Kirchhoff stress

We treat the first Piola–Kirchhoff stress tensor \(\mathbf{P}\) as a function of deformation gradient \(\mathbf{F}\):

$$ \mathbf{f} _0= −\frac{\mathbf{P} (\mathbf{F} )}{6}(\mathbf{X} _{10}×\mathbf{X} _{20}+\mathbf{X} _{20}×\mathbf{X} _{30}+\mathbf{X} _{30}×\mathbf{X} _{10}) $$

It converts an interface normal \(\mathbf{N}\) in the reference state to a traction \(\mathbf{t}\) in the deformed state.

$$ \mathbf{t}=\mathbf{P} (\mathbf{UDV^T} )\mathbf{N} $$

✅ 没讲，

P30

Rotation-Invariance

The stress tensor \(\mathbf{P}\) is rotation-invariant to \(\mathbf{U}\):

$$ \mathbf{P} (\mathbf{UDV^T} )=\mathbf{UP} (\mathbf{DV^T} ) $$

✅ 没讲，

P31

Isotropic Materials

The stress tensor \(\mathbf{P}\) is rotation-invariant to \(\mathbf{U}\):

$$ \mathbf{P} (\mathbf{DV^T} )=\mathbf{P} \mathbf{(D)V^T} $$

✅ 没讲，

P32

Isotropic Materials

Isotropic Materials：各向同性材料

✅ 符号解释：\(\mathbf{P}\)：First Piola Stress、 \(\mathbf{F}\)：Deformation Gradient
✅ 各向同性公式认为：\(\mathbf{P}\) 是关于 \(\mathbf{F}\) 的函数
✅ 对F做 \(\mathbf{SVD}\) 分解可得到 \(\mathbf{UDV^T}\)，其中\(D\)是对角矩阵、其对角元素描述了三个方向的拉伸的量、把公式中的旋转分量剔除掉、 \(\mathbf{P}\) 只与 Principal stretches 有关。

In many literatures, people parameterize \(\mathbf{P} (I_\mathbf{C},II_\mathbf{C},III_\mathbf{C} )\) by principal invariants, for:

$$ I_\mathbf{C} =\mathrm{trace} (\mathbf{C} )=λ_0^2+λ_1^2+λ_2^2 $$

$$ III_\mathbf{C} =\mathrm{det} (\mathbf{C} ^2)=λ_0^4+λ_1^4+λ_2^4 $$

$$ II_\mathbf{C} =\frac{1}{2} (\mathrm{trace} ^2(\mathbf{C} )−\mathrm{trace} (\mathbf{C} ^2))=λ_0^2λ_1^2+λ_0^2λ_2^2+λ_1^2λ_2^2 $$

\(\mathbf{C=U^TU}\) is the right Cauchy-Green deformation tensor.

✅ \(I_c、 II_c、 III_c\) 的定义是基于材料学、数学的先验知识

P33

Isotropic Models

✅ 材料力学中更常用 neo-Hookean
✅Fung常用来模拟人体组织。

P34

计算P

Anyway, we still use the principal stretches for computation:

$$ \mathbf{P} (λ_0,λ_1,λ_2)=\begin{bmatrix} \frac{∂W}{∂λ_0} & \Box &\Box \\ \Box & \frac{∂W}{∂λ_1} & \Box \\ \Box & \Box &\frac{∂W}{∂λ_2} \end{bmatrix} $$

And we compute the first Piola-Kirchhoff stress as:

$$ \mathbf{P} = \mathbf{UP} (λ_0,λ_1,λ_2)\mathbf{V} ^\mathbf{T} $$

P35

A Quick Summary

✅ 主要是算P的方法不同

P36

The Limitation of StVK

Irving et al. 2004. Invertible Finite Elements For Robust Simulation of Large Deformation. SCA

✅ 纵轴是力、横轴长度为弹簧长度、参考长度是1， 因此横轴为1时纵轴为0. 横轴 > 1 代表拉伸、拉伸越大代表力越大。但压缩时， \(StVK\) 表现出的力不对，且当弹簧（四面体）反转以后，力也会反转，这种表现也不对，因为最后会停在横轴-1的状态上。

P37

✅ Poison Effect： 弹性体往上拉时两边会凹进去，本质原因是保体积。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P39

A Summary For the Day

• FEM uses the derivates of the strain energy function to obtain the force.

• FVM uses the integral of the interface traction to obtain the force.

• The two approaches lead to the identical outcome, in different formulations

• Hyperelastic models define the strain energy function by principal stretches, i.e., the singular values of the deformation gradient.

• For isotropic materials, we can calculate the stress through diagonalization.

✅ Level：1. 了解，会用；2. 理解、举一反三；3. 跳出图形学；
✅ 图形学关注的不是数学模型，而是快。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P3

Topics for the Day

• Hessian of Elastic Energy

• Implicit Integration

• Nonlinear optimization.

✅ 隐式积分不讲了。推荐P12页的论文
直接跳到 P18、非线性优化。

P4

Hessian of Elastic Energy

P5

Recall that…

P6

Energy Hessian:

P7

SVD Derivative

Since \(\mathbf{F=UΛV^T}\)is singular value decomposition \((\mathbf{Λ} =\mathrm{diag} (λ_0,λ_1,λ_2))\), we can do:

$$ \mathbf{U^T} \frac{∂\mathbf{F} }{∂\mathbf{F} _{kl}}\mathbf{V} =\mathbf{U^T} (\frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}}\mathbf{ΛV^T} +\mathbf{U}\frac{∂\mathbf{Λ} }{∂\mathbf{F} _{kl}}\mathbf{V^T}+\mathbf{UΛ}\frac{∂\mathbf{V^T} }{∂\mathbf{F} _{kl}})\mathbf{V} $$

$$ \mathbf{U^T} \frac{∂\mathbf{F} }{∂\mathbf{F} _{kl}}\mathbf{V} =(\mathbf{U^T} \frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}})\mathbf{Λ} +\frac{∂\mathbf{Λ} }{∂\mathbf{F} _{kl}}+ \mathbf{Λ}(\frac{∂\mathbf{V^T} }{∂\mathbf{F} _{kl}}\mathbf{V}) $$

$$ \mathbf{U^T} \frac{∂\mathbf{F} }{∂\mathbf{F} _{kl}}\mathbf{V} = \mathbf{AΛ} +\frac{∂\mathbf{Λ} }{∂\mathbf{F} _{kl}}+ \mathbf{ΛB} $$

P8

Skew-Symmetric Matrix

Matrix \(\mathbf{A}\) is skew-symmetric (or anti-symmetric), if \(\mathbf{A=−A^T}\):

$$ \mathbf{A} =\begin{bmatrix} 0 & a & b \\ -a & 0 & c \\ -b & -c & 0 \end{bmatrix} $$

If \(\mathbf{D}\) is diagonal, then:

$$ \mathbf{AD} =\begin{bmatrix} 0 & a & b \\ -a & 0 & c \\ -b & -c & 0 \end{bmatrix}\begin{bmatrix} d & 0 &0 \\ 0 & e & 0\\ 0 & 0 &f \end{bmatrix}=\begin{bmatrix} 0 & ? & ?\\ ? & 0 &? \\ ? & ? &0 \end{bmatrix} $$

$$ \mathbf{DA} = \begin{bmatrix} d & 0 & 0\\ 0 & e & 0\\ 0 & 0 &f \end{bmatrix}\begin{bmatrix} 0 & a & b\\ -a & 0 &c \\ -b & -c &0 \end{bmatrix}=\begin{bmatrix} 0 & ? &? \\ ? & 0& ?\\ ? & ? &0 \end{bmatrix} $$

When \(\mathbf{U}\) is orthogonal, we have:

$$ \mathbf{0} =\frac{∂(\mathbf{U^TU)}}{∂\mathbf{F} _{kl}} =\mathbf{U^T} \frac{ ∂\mathbf{U} }{∂\mathbf{F} _{kl}}+\frac{∂\mathbf{U^T}}{∂\mathbf{F} _{kl}}\mathbf{U} =\mathbf{U^T}\frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}}+(\mathbf{U^T} \frac{∂\mathbf{U}}{∂\mathbf{F} _{kl}})^\mathbf{T} $$

Therefore, \(\mathbf{A=U^T}\frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}}\) is skew-symmetric. So is \(\mathbf{B} =\frac{∂\mathbf{V^T} }{∂\mathbf{F} _{kl}} \mathbf{V}\).

P9

SVD Derivative (cont.)

Since \(\mathbf{F=UΛV^T}\) is singular value decomposition \(\mathbf{Λ} =\mathrm{diag} (λ_0,λ_1,λ_2)\), we can do:

$$ \mathbf{U^T} \frac{∂\mathbf{F} }{∂\mathbf{F} _{kl}}\mathbf{V} = \mathbf{AΛ} +\frac{∂\mathbf{Λ} }{∂\mathbf{F} _{kl}}+\mathbf{ΛB} $$

After expansion, we get:

$$ \mathbf{U^T} \frac{∂\mathbf{F}}{∂\mathbf{F} _{kl}} \mathbf{V} =
\begin{bmatrix} 0 & a_0 & a_1\\ -a_0 & 0 & a_2\\ -a_1& -a_2 &0 \end{bmatrix} \begin{bmatrix} λ_0 & \Box & \Box\\ \Box & λ_1 &\Box \\ \Box & \Box &λ_2 \end{bmatrix} \begin{bmatrix} \frac{∂λ_0}{∂\mathbf{F} _{kl}} & \Box & \Box\\ \Box & \frac{∂λ_1}{∂\mathbf{F} _{kl}} & \Box\\ \Box & \Box &\frac{∂λ_2}{∂\mathbf{F} _{kl}} \end{bmatrix} \begin{bmatrix} λ_0 & \Box & \Box \\ \Box & λ_1 &\Box \\ \Box & \Box &λ_2 \end{bmatrix} \begin{bmatrix} 0 & b_0& b_1\\ -b_0 & 0 &b_2 \\ -b_1 & -b_2 &0 \end{bmatrix} $$

P10

SVD Derivative (cont.)

Since \(\mathbf{F=UΛV^T}\) is singular value decomposition \(\mathbf{Λ}=\mathrm{diag} (λ_0,λ_1,λ_2)\), we can do:

$$ \mathbf{U^T} \frac{∂\mathbf{F} }{∂\mathbf{F} _{kl}}\mathbf{V} =\mathbf{AΛ} +\frac{∂\mathbf{Λ} }{∂\mathbf{F} _{kl}}+\mathbf{ΛB} $$

After expansion, we get:

$$ \mathbf{U^T} \frac{∂\mathbf{F}}{∂\mathbf{F} _ {kl}}\mathbf{V}
= \begin{bmatrix} ∂λ _ 0/∂\mathbf{F} _ {kl} & λ _ 1 a _ 0 + λ _ 0 b _ 0 & λ _ 2 a_ 1 + λ _ 0 b _ 1 \\ −λ _ 0 a_ 0 − λ _ 1 b _ 0 & ∂λ _ 1/∂ \mathbf{F} _ {kl} & λ _ 2 a _ 2 + λ _ 1 b _ 2\\ −λ _ 0 a_ 1−λ_ 2b_ 1 & −λ_ 1a_ 0a_ 2−λ_ 2b_ 2 &∂λ_ 2/∂ \mathbf{F} _ {kl} \end{bmatrix} $$

$$ \begin{bmatrix} m_ {00} & m_ {01} & m_ {02} \\ m_ {10} & m_ {11} & m_ {12} \\ m_ {20} & m_ {21} & m_ {22} \end{bmatrix} = \begin{bmatrix} ∂λ_0/∂\mathbf{F} _{kl} & λ_1a_0+λ_0b_0 & λ_2a_1+λ_0b_1\\ −λ_0a_0−λ_1b_0 & ∂λ_1/∂ \mathbf{F} _{kl} & λ_2a_2+λ_1b_2\\ −λ_0a_1−λ_2b_1 & −λ_1a_0a_2−λ_2b_2 &∂λ_2/∂ \mathbf{F} _ {kl} \end{bmatrix} $$

Eventually, we get: \(\mathbf{A=U^T} \frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}}, \mathbf{B} = \frac{∂\mathbf{V^T} }{∂\mathbf{F} _{kl}}\mathbf{V} \) and \(∂λ_0/∂\mathbf{F} _{kl}, ∂λ_1/∂\mathbf{F} _{kl}, ∂λ_2/∂\mathbf{F} _{kl}\).

P11

A Quick Summary

• Step 1: By SVD derivatives, we get: \(\frac{∂\mathbf{U} }{∂\mathbf{F} _{kl}},\frac{∂λ_d}{∂\mathbf{F} _{kl}},\frac{∂\mathbf{V^T} }{∂\mathbf{F} _{kl}}\).

• Step 3: Finally, we reach our goal in Hessian matrix:

$$ \frac{∂\mathbf{f} _i}{∂\mathbf{x} _j}=\sum _ {k,l}\frac{∂\mathbf{f} _i}{∂\mathbf{F} _{kl}} \frac{∂\mathbf{F} _{kl}}{∂\mathbf{x} _j}=\sum _ {k,l} \frac{∂(\mathrm{Udiag} (\frac{∂W} {∂λ_0},\frac{∂W} {∂λ_1},\frac{∂W} {∂λ_2})\mathbf{V^T}) }{∂\mathbf{F} _{kl}}−V\mathbf{d} _i\frac{∂\mathbf{F} _{kl}}{∂\mathbf{x} _j} $$

P12

After-Class Reading

Xu et al. 2015. Nonlinear Material Design Using Principal Stretches. TOG (SIGGRAPH).

Definitely read this paper if you decide to implement it.

P13

Implicit Integration

P14

Implicit Integration

Recall that we need implicit integration to avoid numerical instability (Class 5, page 13):

$$ \begin{cases} \mathbf{v} ^{[1]}=\mathbf{v} ^{[0]}+∆t\mathbf{M} ^{−1}\mathbf{f} ^{[1]} \\
\mathbf{x} ^{[1]}=\mathbf{x} ^{[0]}+∆t\mathbf{v} ^{[1]} \end{cases} $$

$$\mathrm{or} $$

$$ \begin{cases} \mathbf{x} ^{[1]}=\mathbf{x} ^{[0]}+∆t\mathbf{v} ^{[0]}+∆t^2\mathbf{M} ^{−1}\mathbf{f} ^{[1]} \\
\mathbf{v} ^{[1]}=(\mathbf{x} ^{[1]}−\mathbf{x} ^{[0]})/∆t
\end{cases} $$

We also said that:

$$ \mathbf{x} ^{[1]}=\mathbf{x} ^{[0]}+∆t\mathbf{v} ^{[0]}+∆t^2\mathbf{M} ^{−1}\mathbf{f}(\mathbf{x}^{[1]}) $$

$$ \Updownarrow $$

$$ \mathbf{x} ^{[1]}=\mathrm{argmin } F (\mathbf{x} ) \quad \mathrm{for} \quad\mathbf{F} (\mathbf{x} )=\frac{1}{2∆t^2}||\mathbf{x} −\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]}||_M^2+E(\mathbf{x} ) $$

P15

Newton-Raphson Method

The Newton-Raphson method, commonly known as Newton’s method, solves the optimization problem: \(x^{[1]}\)=argmin \(F(x)\). \(\quad\quad\) \((F(x)\) is Lipschitz continuous.)

Given a current \(x^{(k)}\), we approximate our goal by:

$$ 0={F}'(x)≈{F}'(x^{(k)})+{F}'' (x^{(k)})(x−x^{(k)}) $$

P16

Newton-Raphson Method

Now we can apply Newton’s method to: \(\mathbf{x} ^{[1]}\)= argmin \( F (\mathbf{x} )\).
Given a current \(\mathbf{x}^{(k)}\), we approximate our goal by:

$$ \mathbf{0} =∇F(\mathbf{x} )≈∇F(\mathbf{x} ^{(k)})+\frac{∂F^2(\mathbf{x} ^{(k)})}{∂\mathbf{x} ^2} (\mathbf{x} −\mathbf{x} ^{(k)}) $$

Newton’s Method
Initialize \(x^{(0)}\)
For \(k=0…K\)
$$ ∆\mathbf{x} \longleftarrow −(\frac{∂F^2(\mathbf{x} ^{(k)})}{∂\mathbf{x} ^2})^{−1}∇F(\mathbf{x} ^{(k)}) $$ $$ \mathbf{x} ^{(k+1)}\longleftarrow \mathbf{x} ^{(k)}+∆\mathbf{x} $$ If \(||∆\mathbf{x}||\) is small \(\quad\quad\)then break
$$ \mathbf{x} ^{[1]}\longleftarrow \mathbf{x} ^{(k+1)} $$

P17

Simulation by Newton’s Method

Specifically to simulation, we have:

$$ F (\mathbf{x} )=\frac{1}{2∆t^2} ||\mathbf{x} −\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]}||_M^2+E(\mathbf{x} ) $$

$$ ∇F(\mathbf{x} ^{(k)})=\frac{1}{∆t^2}\mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})−\mathbf{f} (\mathbf{x} ^{(k)}) $$

$$ \frac{∂^2F(\mathbf{x} ^{(k)})}{∂\mathbf{x} ^2} =\frac{1}{∆t^2}\mathbf{M} +\mathbf{H} (\mathbf{x} ^{(k)}) $$

Initialize \(\mathbf{x}^{(0)}\), often as \(\mathbf{x} ^{[0]}\) or \(\mathbf{x} ^{[0]} +∆t\mathbf{v} ^{[0]}\)
For \(k=0…K\)
$$ \mathrm{Solve}\quad (\frac{1}{∆t^2} \mathbf{M+H} (\mathbf{x} ^{(k)}))∆\mathbf{x} =− \frac{1}{∆t^2}\mathbf{M} (\mathbf{x} ^{(k)}−\mathbf{x} ^{[0]}−∆t\mathbf{v} ^{[0]})+\mathbf{f} (\mathbf{x} ^{(k)}) $$ $$ \mathbf{x} ^{(k+1)}\longleftarrow \mathbf{x} ^{(k)}+∆\mathbf{x} $$ If ||\(∆\mathbf{x}\)|| is small \(\quad\quad\) then break
$$ \mathbf{x} ^{[1]}\longleftarrow \mathbf{x} ^{(k+1)} $$ $$ \mathbf{v} ^{[1]}←(\mathbf{x} ^{[1]}−\mathbf{x} ^{[0]})/ ∆t $$

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P18

Nonlinear Optimization

P19

Gradient Descent

Another way to solve \(\mathbf{x}^∗\)=argmin \(F(\mathbf{x})\) is the gradient descent method.

How to find the optimal step size becomes a critical question.

P20

step size adjustment

优点：simple, Low overhead

P21

Descent Directions

The direction \(\mathbf{d(x)}\) is descending, if a sufficiently small step size \(α\) exists for:

$$ F(\mathbf{x} )>F(\mathbf{x} +α\mathbf{d} (\mathbf{x} )) $$

✅沿负梯度方向可以下降，但不一定是最好的方向。怎样判断一个方向是否可以下降？答：看与负梯度方向是否在同侧。

P22
With line search, we can use any search direction as long as it’s descending:

$$ F(\mathbf{x} ^{(0)})>F(\mathbf{x} ^{(1)})>F(\mathbf{x} ^{(2)})>F(\mathbf{x} ^{(3)})>… $$

P23

Descent Methods

• Gradient descent is a descent method, since:

$$ \mathbf{d} (\mathbf{x} )=−∇F(\mathbf{x} )\quad \Rightarrow \quad −∇F(\mathbf{x} )\cdot (−∇F(\mathbf{x} ))>0 $$

• Newton’s method is also a descent method, if the Hessian is always positive definite:

$$ \mathbf{d} (\mathbf{x} )=−(\frac{∂^2F(\mathbf{x} )}{∂\mathbf{x} ^2})^{−1}∇F(\mathbf{x} ) \quad \Rightarrow \quad −∇F(\mathbf{x} )\cdot (−(\frac{∂^2F(\mathbf{x} )}{∂\mathbf{x} ^2})^{−1}∇F(\mathbf{x} ))>0 $$

✅牛顿法不一定收敛，\(\mathbf{H}\)正定场景牛顿法一定收敛。

• Any method using a positive definite matrix P to modify the gradient yields a descent method:

$$\mathbf{d} (\mathbf{x} )=−\mathbf{P} ^{−1}∇F(\mathbf{x} ) \quad \Rightarrow \quad −∇F(\mathbf{x} )\cdot (−\mathbf{P} ^{−1}∇F(\mathbf{x} ))>0 $$

P24

A unified descent framework

A unified descent framework

P25

✅ 图形学中更关注 Total Cost. 让 P 更加接近 H，可以减少迭代数，让 P 更容易得到，减少迭代成本。
Traction：物体表面上的力的密度，有点像压强

P27

After-Class Reading

Wang. 2016. Descent Methods for Elastic Body Simulation on the GPU. TOG (SIGGRAPH Asia).

P28

A Summary For the Day

• We can calculate the Hessian of the FEM elastic energy based on SVD derivatives.

• The goal of doing this is for implicit time integration.

• Fundamentally, the goal is to solve a nonlinear optimization.

  • Gradient Descent, Newton’s method, and others can all be considered as descent methods.
  • The key question is the matrix for calculating the search direction.
  • We need both the per-iteration cost and the number of iterations to be small.

✅ 模拟的公式通常都固定，很难有突破、瓶颈在于计算量、随着分辨率的提升，模拟的计算量几乎是无止境的。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• Collision Detection
• Interior Point Methods
• Impact Zone Optimization
• Untangling Cloth

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P4

Collision Detection Pipeline

P5

Spatial Partitioning

Spatial partitioning divides the space by a grid and stores objects into grid cells.

P6

静止场景

To find pair candidates for collision test, we just have to check the grid cells.

P7

运动场景

If we need to consider moving objects, we just expand the object region.

P8

数据存储

Instead of allocating memories to cells, we can build an object-cell list and then sort them. This avoids memories wasted in empty cells.

✅ 用 hash 链表代替数组
✅要解决的问题：3D空间需要划分出大量的小格。有的格子可能包含很多object. 大多数格子可能没有object.
✅方法：3D数组转为list表示法。
✅缺点：内存访问不连续。

P9

Morton Code

One question is how to define the cell ID. Using the grid order is not optimal, since it cannot be easily extended and it is lack of locality. Morton code uses a Z-pattern instead.

✅ 希望内存访问尽量连续。也就是下一次访问的内存地址在上次的附近
✅ Giid Order：横向访问连续、纵向访问不连续、三维情况会更严重。

✅ Morton Code：一种对格子编号的顺序。

P10

After-Class Reading

GPU Gems 3

Chapter 32. Broad-Phase Collision Detectionwith CUDA

P14

Bounding Volume Hierarchy

Bounding volume hierarchy is built on geometric/topological proximity of objects.

✅空间划分→物体划分

构造BVH

P15

外部物体与身体相交检测

To find elements potentially in collision with an object, we just traverse the tree.

P16

Bounding Box

Axis-aligned bounding box (AABB) is the most popular bounding volume. Besides that, there are also spheres and oriented bounding box (OBB).

Two AABBs intersect if and only if they intersect in every axis.

P17

自相交检测

To process self collisions by BVH, we define two procedures.

P18

Bounding Volume Hierarchy

The performance depends on the effectiveness of culling.

Zheng and James. 2012. Energy-based Self-Collision Culling for Arbitrary Mesh Deformations. TOG (SIGGRAPH)

✅ 对每个区域计算能量，根据形变能量的大小来判断有没有可能相交，此方法不适用于衣服，因为在衣服模拟中大形变很常见、不代表有相交。

Changxi Zheng, Doug L.James, Cornell University

P19

Comparison between SH and BVH

• Spatial Hashing

  • Easy to implement
  • GPU friendly
  • Needs to recompute after updating objects

• Bounding Volume Hierarchy

  • More involved
  • Not GPU friendly
  • To update BVH, just update bounding volumes

✅ CUDA 代码. INVIDI 代码通常使用SH。GPU 喜欢简单粗爆的数据结构，BVH相对于GPU过于复杂。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P20

Collision Detection Pipeline

✅ 图画得不对。 DCD 和CCD对输入的需求不同，但SH和 BVH 跟这两种输入不是 一一 对应的关系。
✅ 但DCD和CCD对输入的要求是不一样的。 P20页所谓的 D 还是 C,是以时间角度说的。

P21

Discrete Collision Detection (DCD)

DCD tests if any intersection exists in each state at discrete time instant: \(\mathbf{x}^{[0]}\), \(\mathbf{x}^{[1]}\), …

✅ 准确来说。DCD检测的不是碰撞，而是相交

edge-triangle intersection

To a triangle mesh, the basic test is edge-triangle intersection test.

\(t\) 代表相交位置对应 \(\mathbf{x}_a\) 和\(\mathbf{x}_b\)的插值量

✅ 检测在特定状态下是否相交，每一帧都不相交就认为无碰撞。
✅ 相交和碰撞的区别：相交分析的是运动前后的状态、碰撞检测的是运动的过程、未相交不一定无碰撞、

P22

Tunneling

DCD is simple and robust, but it suffers from the tunneling problem: objects penetrating through each other without being detected.

tunneling problem：当物体运动特别快时，有可能的穿透另一物体而没有被检测到，常见于细薄物体、例如衣服
这种情况无相交但是有碰撞

P23

Continuous Collision Detection (CCD)

CCD tests if any intersection exists between two states: \(\mathbf{x} ^{[0]}\) and \(\mathbf{x} ^{[1]}\).

To a triangle mesh, there two basic tests: vertex-triangle and edge-edge tests.

vertex-triangle tests

✅ 当四点共面时，构成的四面体体积为0、利用四面体的体积公式，可求出四点共面的时间 \(t\) . 这里的\(\mathbf{t}\)是时间
✅ 假设运动是匀速的，\( \mathbf{x}_ {30}(t)、 \mathbf{x}_ {10}(t)、\mathbf{x}_ {20}(t)\)都是关于\(t\)的线性函数。
✅ 一元三次方程有公式解，但用到\(\sqrt[3]{\cdot}\)，因此不建议使用，建议用牛顿法。因为\(\sqrt[3]{\cdot }\) 的误差非常大。

P24

edge-edge tests.

✅ 为什么要检测边边相交，因为有可能三角形相交但点面没有相交。

✅ 先求四点共面的 \(t\)
✅ 解一元三次方程也不建议牛顿法，而是二分法，因为\(t\)的范围是[0,1]

P25

Issues with CCD

• Floating-point errors, especially due to root finding of a cubic equation
  • Buffering epsilons, but that causes false positives.
  • Gaming GPUs often use single floating-point precision.

✅ 游戏 GPU 以单精度为主，因此要注意浮点误差问题。

• Computational costs: more expensive than DCD.

  • Some argue that broad-phase collision culling is the bottleneck.

• Difficulty in implementation.

P26

After-Class Reading

Bridson et al. 2002. Robust Treatment of Collisions, Contact and Friction for Cloth Animation. TOG (SIGGRAPH).

Relative simple explicit integration of cloth dynamics

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P27

Interior Point Methods and Impact Zone Optimization

✅ 发现碰撞的pairs后如何处理。

P28

Two Continuous Collision Response Approaches

✅ 这是两个大的套路，不是具体的方法。

Given the calculated next state \(\mathbf{x} ^{[1]}\), we want to update it into \(\bar{\mathbf{x} } ^{[1]}\), such that the path from \(\mathbf{x} ^{[0]}\) to \(\bar{\mathbf{x} } ^{[1]}\) is intersection-free.

内点法：

	内点法	Impact Zone 法
			✅ 蓝色区域为安全区域
	✅ 从\(\mathbf{x}^{[0]}\)出来，朝\(\mathbf{x}^{[1]}\)走，并永远保证只在安全区域走，直到不能走为止。	✅ 从\(\mathbf{x}^{[1]}\)出发，反复优化结果（投影），直到回到安全区域为止。
优点	Always succeed	Fast. 1. Close to solution. 2. Only vertices in collision (impact zones). 3. Can take large step sizes.	✅ Impact Zone：\(\mathbf{x}^{[1]}\)通常离安全区域不太远，且优化时只针对 Impact Zone 优化，因此快。
局限性	Slow. 1. Far from solution. 2. All of the vertices. 3. Cautiously by small step sizes.	May not succeed.	✅ 内点：为保证每一步安全，步长不能太大，因此慢、哪怕\(\mathbf{\bar{x}}^{[1]}\)最终没有到最佳位置，但能保证一定在安全区域，因此一定成功。\(\mathbf{x}^{[0]}\)和\(\mathbf{x}^{[1]}\)可能比较远，也导致慢。 ✅ 内点法慢的原因2没有解释。

P30

Log-Barrier Interior Point Methods

算法过程

For simplicity, let’s consider the Log-barrier repulsion between two vertices.

$$E(\mathbf{x} )=−\rho \text{ log} ||\mathbf{x} _{ij}||$$

$$ \mathbf{f} _i(\mathbf{x} )=−∇_iE=ρ\frac{\mathbf{x} _ {ij}}{||\mathbf{x} _ {ij}||^2} \\ \mathbf{f} _j(\mathbf{x} )=−∇_jE=−ρ\frac{\mathbf{x} _ {ij}}{||\mathbf{x} _{ij}||^2} $$

✅ 用 Log 定义能量、前面某一节课讲过。距离 → 能量 → 斥力

✅ 不需要互斥力一直存在，因此做了一个截断（IPC）

P31

算法实现

We can then formulate the problem as:

$$ \bar{\mathbf{x} }^ {[1]}\longleftarrow \mathrm{argmin} _\mathbf{x} (\frac{1}{2} ||\mathbf{x} −\mathbf{x} ^{[1]}||^2−ρ\sum \mathrm{log} ||\mathbf{x} _{ij}||) $$

✅ 优化目标：点的位置与目标位置（穿模）尽量接近，然后优化。

Gradient Descent:

\(\mathbf{x} ^{(0)}\longleftarrow \mathbf{x} ^{[0]}\)
For \(k=0…K\)
$$\mathbf{x} ^{(k+1)}\longleftarrow \mathbf{x} ^{(k)}+α(\mathbf{x} ^{[1]}−\mathbf{x} ^{(k)}+ρ\sum \frac{\mathbf{x} _{ij}}{||\mathbf{x} _{ij}||^2})$$ \(\bar{\mathbf{x} }^ {[1]}\longleftarrow \mathbf{x} ^{(k+1)}\)

The step size \({\color{Red} α}\) must be adjusted to ensure that no collision happens on the way. To find \({\color{Red} α}\), we need collision tests.

✅ 绿色是来自\(\mathbf{x}^{[1]}\)的引力，黄色是来自边界的斥力、关键是步长\(\alpha \)， 每走一小步都需要反复的碰撞检测。

P32

Impact Zone Optimization

The goal of impact zone optimization is to optimize \(\mathbf{x}^{[1]}\) until it becomes intersection-free. (This potentially suffers from the tunneling issue, but it’s uncommon.)

$$ \bar{\mathbf{x} }^{[1]}\longleftarrow \mathrm{argmin} _\mathbf{x}  \frac{1}{2} ||\bar{\mathbf{x} }-\mathbf{x}^{[1]}||^2 $$

$$ \text{such that} \begin{cases} C(\mathbf{x} )=−(\mathbf{x} _3−b_0\mathbf{x} _0−b_1\mathbf{x} _1−b_2\mathbf{x} _1)\cdot \mathbf{N} ≤0 & \text{ For each detected vertex-triangle pair } \\ C(\mathbf{x} )=−(b_2\mathbf{x} _2+b_3\mathbf{x} _3−b_0\mathbf{x} _0−b_1\mathbf{x} _1)\cdot \mathbf{N}≤0 & \text{ For each detected edge-edge pair } \end{cases} $$

✅ 利用 constrain（不是能量）转化成优化问题，具体没讲。

P33

Geometric Impulse

The goal of impact zone optimization is to optimize \(\mathbf{x}^{[1]}\) until it becomes intersection-free. (This potentially suffers from the tunneling issue, but it’s uncommon.)

Every pair gives new positions to the involved vertices. We can combine them together in a Jacobi, or Gauss-Seidel fashion, just like position-based dynamics.

P34

After-Class Reading (Cont.)

Bridson et al. 2002. Robust Treatment of Collisions, Contact and Friction for Cloth Animation. TOG (SIGGRAPH).

Relative simple explicit integration of cloth dynamics

P35

Augmented Lagrangian

$$ \bar{\mathbf{x} }^{[1]}\longleftarrow \text{argmin} _\mathbf{x}  \frac{1}{2} ||{\mathbf{x} }-\mathbf{x}^{[1]}||^2 $$

$$ \text{such that} \begin{cases} C(\mathbf{x} )=−(\mathbf{x} _3−b_0\mathbf{x} _0−b_1\mathbf{x} _1−b_2\mathbf{x} _1)\cdot \mathbf{N} ≤0 & \text{ For each detected vertex-triangle pair } \\ C(\mathbf{x} )=−(b_2\mathbf{x} _2+b_3\mathbf{x} _3−b_0\mathbf{x} _0−b_1\mathbf{x} _1)\cdot \mathbf{N}≤0 & \text{ For each detected edge-edge pair } \end{cases} $$

P36

Augmented Lagrangian

We can then convert it into an unconstrained form:

$$\bar{\mathbf{x} } {[1]}\longleftarrow \text{argmin}_{x,λ}(\frac{1}{2} ||\mathbf{x} −\mathbf{x} ^{[1]}||^2+\frac{ρ}{2} ||\text{max}(\tilde{C} (\mathbf{x} ))||^2−\frac{1}{2ρ}||\mathbf{λ} ||^2) $$

$$ \tilde{C} (\mathbf{x})= \text{max}(\mathbf{C} (\mathbf{x} )+\mathbf{λ} /ρ) $$

Augmented Lagrangian:

\(\mathbf{x} ^{(0)} \longleftarrow \mathbf{x} ^{[0]}\)
\(\mathbf{λ \longleftarrow 0} \)
For \(k=0…K\)
$$\mathbf{x} ^{(k+1)} \longleftarrow \mathbf{x} ^{(k)}−α∇(\frac{1}{2} ||\mathbf{x} −\mathbf{x} ^{[1]}||^2+\frac{ρ}{2} ||\text{max}(\tilde{C} (\mathbf{x} ))||^2−\frac{1}{2ρ}||\mathbf{λ} ||^2)$$
\(λ\longleftarrow ρ\tilde{C} (\mathbf{x} )\)
\(\bar{\mathbf{x} } ^{[1]}\longleftarrow \mathbf{x} ^{(k+1)}\)

Tang et al. 2018. I-Cloth: Incremental Collision Handling for GPU-Based Interactive Cloth Simulation. TOG. (SIGGRAPH Asia)

P37

About Impact Zone Optimization

• Fast per iteration

  • Only have to deal with vertices in collision.

• Convergence sensitive to \(||\mathbf{x} ^{[0]}−\mathbf{x} ^{[1]}||^2\), or the time step \(∆t\)

  • Can take many iterations to, or never achieve intersection-free.
  • Easy solution is to reduce \(∆t\), but that increases total costs.

P38

Rigid Impact Zones

The rigid impact zone method simply freezes vertices in collision from moving in their pre-collision state. It’s simple and safe, but has noticeable artifacts.

✅ 检测到碰撞，则把这个区域退回到上一帧。

P39

A Practical System

✅ 有碰撞，先做 Impact Zone. 因为这个快、不能解决再用后面方法、计算量不允许则选择 Rigid Impact.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P41

Intersection Elimination

• Let’s consider how to eliminate existing intersections, but without using any collision history.

• Such a method is useful when there are already intersections in simulation, due to:

  • Past collision handling failures
  • Intense user interaction

• In this case, we don’t require the simulation is to always intersection-free.

P42

对于有体积的物体

Eliminating cloth-volume and volume-volume intersections is straightforward: simply pushing vertices/edges in the volume out.

P43

对于没有体积的物体，Untangling Cloth问题

The situation is complicated in cloth-cloth intersection, since we don’t have a clear definition of inside and outside.

方法一

算法过程

Baraff et al. used flood-fill to segment cloth into regions and decided which region is in intersection. (Cannot handle boundary well.)

Baraff et al. 2003. Untangling Cloth. TOG (SIGGRAPH)

✅P42适用于有体积的物体，但布没有封闭体积，两根线没有里面外面之分，因此相交时不知道哪一段是正确的。
✅ 方法：对布分段，根据分段区域决定谁在上谁在下，以此为依据推动顶点。
✅ 此方法缺点：1. 无法处理边界；2. 难以在 GPU 上实 现；

P44

算法效果

Baraff et al. 2003. Untangling Cloth. TOG (SIGGRAPH)

✅缺点：1. 难以处理边界；2. 对整个面进行评估，难以用于GPU.

P45

方法二

Volino and Magnenat-Thalmann proposed to untangle cloth by reducing the intersection contour.
Their method can handle boundaries, but it doesn’t always work.

✅ 两个面相交会产生一条曲线，目标是让曲线变短。优点：可以处理边界；缺点：基于局部优化、可用于 GPU。
✅可以处理边界情况，缩短边界也能解除相交。

P46

After-Class Reading

Volino and Magnenat-Thalmann et al. 2006. Resolving Surface Collisions through Intersection Contour Minimization. TOG (SIGGRAPH).

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P47

A Summary For the Day

• Collision handling involves two steps: collision detection and collision response.

• Collision detection contains two phases: broad-phase culling and narrow-phase test.

• There are two types of collision detection tests: discrete and continuous.

• Similarly, there are discrete and continuous collision responses.

• For continuous collision responses, we must update the state to become collisionfree state. There are two approaches: interior point method and impact zone optimization. Rigid impact zone is also a method, but it’s problematic.

• For discrete collision responses, we allow intersections to stay and hope to remove them in long turn. Cloth-cloth intersections are difficult to handle.

✅ 如果考虑摩擦，通常把摩擦做为后处理，但这样结果不精确。如果同时处理摩擦和碰撞、会很复杂。
✅Impulse方法的碰撞检测通常用SDF．但很多形变体无法使用SDF.
✅Impulse响应方式是离散响应方式，无法处理穿透问题。
✅碰撞问题通常不使用物理方法，因为使用物理方法需要小步长，效率非常低。
✅碰撞开源代码：bullet. physics X

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Fluid Effects

Unlike other bodies, Fluids exhibit highly volatile behaviors. As a result, it’s difficult to come up with a single, efficient way for simulating all of fluid effects.

✅ 流体的形态很多，例如水滴、水花、水浪，对应的模拟方法也不同。难以用通用的方法高效地模拟所有场景。

P3

Two Types of Simulation Approaches

Lagrangian Approach (dynamic particles or mesh) Node movement carries physical quantities (mass, velocity, …).	Eulerian Approach (static grid or mesh) Grid/Mesh doesn’t move. Stored physical quantities change.
✅ 左：无 Grid. 物理量附加在粒子上，粒子运动时更新自身物理量。	✅ 右：固定 Grid. 物理量固定在 Grid 上。粒子运动后统一新格子的物理量。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P4

A Height Field Model

高度场（Height Field）和速度场

P5

高度场（Height Field）和速度场的定义

✅ 利用高度函数来表达波的平面，通过更新\(h(x)\)来表达水面随时间波动的效果。
✅ 由于用函数表达，无法描述大海浪的效果，因为这种情况下一个位置对应多个高度。
✅ 速度场描述水流的速度和方向， 速度< 0 则右往左， > 0 则左往右。

P6

高度场的更新

✅ 高度场更新公式第一项：高度场随时间的变化。
✅ 高度场更新公式第二项，根据微分的定义：
\(h(x)u(x)\): 单位时间内流过x线的水量。
\(d(h(x)u(x))\) ：单位时间内区域 \([x \quad x+dx]\) 的水量变化、
\(d(h(x)u(x))1/dx\) ：单位时间内区域区 \([x \quad x+dx]\) 的水位高度变化

P7

速度场的更新

The velocity is also a function of \(x:u(x)\).

✅ 第一项：当水在流动时，速度应该跟水一起流动，下节课再讲。
✅ 第三项：外力，当前也不考虑。

P8

Ignoring advection and external acceleration, we get a simple form:

$$ \begin{matrix} \frac{du(x)}{dt}=−\frac{1}{ρ} \frac{dP(x)}{dx} \quad \quad & ρ: \text{density} \quad \quad & P(x):\text{pressure} \end{matrix} $$

✅ 第二项： 在短时间内、速度变化由左右压强差决定。同样的压强下，密度大则难推，密度小则好推。

P9

Shallow Wave Equation

✅ 为什么叫 Shallow Wave, 因为该算法假设水波很小，因此 \(dh / dx\) 可忽略不计。

We now have two equations:

✅ 公式化简的目的：不需要关心速度场\(u\)、仅关注高度场就可以。
✅ 第一个公式：（1）对\(d(hu)\)展开（2）再求一次\(dt\)
✅ 第二个公式：对\(x\)求导

We can then eliminate \(u\) and formulate the shallow wave equation:

✅ 合并同类项，得到最终方程
✅ 但引擎无法直接处理微分程，因此要离散化开求解。

P10

高度场离散化

We discretize a continuous height field into a discrete set of height columns.

✅ 高度场离散又化为多个水柱，微分算子也要离散化。

P11

微分算子离散化

前向差分与后向差分

The idea of finite differencing is to use the difference to approximate the derivative.

$$ f(t_0+∆t)=f(t_0)+∆t\frac{df(t_0)}{dt} +\frac{∆t^2}{2} \frac{d^2f(t_0)}{dt^2} +… $$

Forward differencing (first-order)

$$ f(t_0−∆t)=f(t_0)−∆t\frac{df(t_0)}{dt}+\frac{∆t^2}{2}\frac{d^2f(t_0)}{dt^2} +… $$

Backward differencing (first-order)

P12

Central Differencing

The idea of finite differencing is to use the difference to approximate the derivative.

$$ f(t_0+∆t)=f(t_0)+∆t\frac{df(t_0)}{dt}+\frac{∆t^2}{2}\frac{d^2f(t_0)}{dt^2} +… $$

$$ f(t_0−∆t)=f(t_0)−∆t\frac{df(t_0)}{dt}+\frac{∆t^2}{2}\frac{d^2f(t_0)}{dt^2} +… $$

Central differencing (second-order)

P13

二阶微分算子离散化

高度

We apply central differencing twice to estimate \(d^2ℎ_i/dt^2\).

$$ \begin{matrix} \frac{dℎ_i(t_0+0.5∆t)}{dt}≈\frac{ℎ_i(t_0+∆t)−ℎ_i(t_0)}{∆t} \quad\quad& \frac{dℎ_i(t_0−0.5∆t)}{dt}≈\frac{ℎ_i(t_0)−ℎ_i(t_0−∆t)}{∆t} \end{matrix} $$

✅ 先用 central difference 求出两个中点的一阶导数，再基于此计算 \(t_0\) 处的二阶导。这种操作又称为一维Laplace 算子。

P14

压强

Similarly, we apply central differencing twice to estimate \(d^2P/dx^2\).

$$ \begin{matrix} \frac{dP_{i+0.5}}{dt} ≈\frac{P_{i+1}−P_i}{∆x} \quad\quad & \frac{dP_{i−0.5}}{dx} ≈\frac{P_i−P_{i−1}}{∆x} \end{matrix} $$

✅ 二维情况用周围4个元素，见 Games102 离散拉普拉斯算子。

P15

Discretized Shallow Wave Equation

We can now discretize the shallow wave equation \(\frac{d^2ℎ}{dt^2}=\frac{ℎ}{ρ}\frac{d^2P}{dx^2}\).

\(\quad\)

\(\quad\)

✅ 更新目标：下一个时刻的水柱的高度，即 \(h_i(t_0 + ∆t)\)
✅ 但按此公式模拟可能出现水的体积变多或变少的问题。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P16

Volume Preservation

We want the volume to stay the same. Suppose that \(\sum ℎ_i(t)=\sum ℎ_i(t−∆t)=V\). But,

$$ ℎ_i(t_0+∆t)=2ℎ_i(t_0)−ℎ_i(t_0−∆t)+\frac{∆t^2ℎ_i}{∆x^2ρ}(P_{i+1}+P_{i−1}−2P_i) $$

✅ 体积会变大还是变小，取决于桔色项，但很难保证这一项是0.

P17

Volume Preservation – Solution 1

✅ 保证 \(h_i\) 和 \(h_{i+1}\)的交换的水量相等、因此保体积
✅ 把\(h_{i-1}\)与\(h_i\)的交换和\(h_i\)与\(h_{i+1}\)的交换拆开。即：
（1）把\((P_{i+1}+P_{i−1}−2P_i)\)拆成\(P_{i−1}−P_i\)和\(P_{i+1}−P_i\)
（2）把\(h_i\)拆成\(\frac{h_{i-1}+h_i}{2}\)和\(\frac{h_{i+1}+h_i}{2}\)
✅ 直观理解：对每个水柱而言，流入的量和流出的量是等价的。

P18

🔎 Kass and Miller. 1990. Rapid, Stable Fluid Dynamics for Computer Graphics. Computer Graphics.

P19

Volume Preservation – Solution 2

An easier way to preserve volume is to simply assume \(h_i\) in the right term is constant.

P20

Pressure

P21

Viscosity

Like damping, viscosity tries to slow down the waves.

✅ Viscosity: 粘滞，相当于流体的阻尼。

P22

Algorithm

$$\text{A Shallow Wave Simulator}$$ For every cell \(i\) $$ℎ_i^{new}←ℎ_i+β(ℎ_i−ℎ_i^{old})\\ ℎ_i^{new}←ℎ_i^{new}+α(ℎ_{i−1}−ℎ_i)\\ ℎ_i^{new} ←ℎ_i^{new}+α(ℎ_{i+1}−ℎ_i)\\$$ For every cell \(i\) $$ℎ_i^{old}←ℎ_i\\
ℎ_i←ℎ_i^{new}$$

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P23

Boundary Conditions

Dirichlet boundary

A Dirichlet boundary assumes that the boundary height \(H_{i+1}\) is constant. It’s considered as an open boundary.

$$ ℎ_{i+1}−ℎ_i+ℎ_{i−1}−ℎ_i=H_{i+1}−ℎ_i+ℎ_{i−1}−ℎ_i $$

✅ 这种方法用于模拟开放的水面，例如大海的区域、假设被模拟的区域外是静止的水面、高度为常数，(Dirichlet)
✅ \(h\)为边界内，\(H\)为边界外。

P25

Algorithm with Neumann Boundaries

Extending the simulator to 3D is also straightforward.

Neumann boundary

A Neumann boundary specifies the boundary derivatives. For example, a zero-derivative boundary means \(ℎ_{i+1}≡ℎ_i\). It’s considered as a closed boundary.

$$ ℎ_{i+1}−ℎ_i+ℎ_{i−1}−ℎ_i=ℎ_{i−1}−ℎ_i $$

✅ Neuman 用于模拟有边界水域，例如池堂、假设边界上没有水流交换。

P24

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P26

Two-Way Coupling

✅ 水和水中的物体相互作用，物体可以是刚体、弹性体能各种类型的物体。

The coupling between a solid and a liquid should be two-way, i.e., liquid->solid and solid->liquid.

✅ 水 → 物体：浮力。物体 → 水，会把这个水柱的水排出去，此处只讲 “物体 → 水” 部分

P27

关键问题

The coupling between solid and water should be two-way, i.e., water>solid and solid- >water.

The key question is how to expel water out of the gray cell regions???

P28

Virtual Height

✅ 在要排的水柱上面增加一个虚拟的高度，然后正常模拟，关键是求出要加多少虚拟高度，能正好达到排出那么多水的效果。

The idea is to set up a virtual height \(v_i\), so that \(ℎ_i^\text{real_new}=ℎ_i−e_i\).

$$ ℎ_i−e_i=ℎ_i+β(ℎ_i−ℎ_i^{old})+α(v_{i+1}+ℎ_{i+1}+ℎ_{i−1}−2v_i−2{ℎ_i})=ℎ_i^{new}+α(v_{i+1}−2v_i) $$

✅ \(ℎ_i^\text{real_new}=ℎ_i−e_i\)：下图左边格子的理想高度。

$$ ℎ_{i+1}−e_{i+1}=ℎ_{i+1}+β(ℎ_{i+1}−ℎ_{i+1}^{old})+α(ℎ_{i+2}+v_i+ℎ_i−2v_{i+1}−2ℎ_{i+1})=ℎ_{i+1}^{new}+α(v_i−2v_{i+1}) $$

✅ 公式2对应右边格子。

P29

Poisson’s Equation

The outcome is Poisson’s equation, with \(v_i\) and \(v_{i+1}\) being unknowns.

$$ 2v_i−v_{i+1}=\frac{1}{α}(ℎ_i^{new}−ℎ_i+e_i)=b_i $$

$$ −v_i+2v_{i+1}=\frac{1}{α}(ℎ_{i+1}^{new}−ℎ_{i+1}+e_{i+1})=b_{i+1} $$

✅ 通过公式化简提取出其中的线性关系。

P30

The outcome is Poisson’s equation, with \(v_i\) and \(v_{i+1}\) being unknowns.

✅ 由于木块位置会变，需要解的\(v_i\)也要改变。 为了让公式统一方便计算，把图左矩阵乘法写成右边形式。公式结果不变，只是工程实现上的简化。

P31

Algorithm with Coupling

✅ \(\gamma \) 的作用：本算法显式积分，不稳定、\(\gamma \) 会让水波小很多。

P32

Rigid Body Update

We estimate the floating force by the actual water expelled in every column.

$$ f_i=ρg∆x(ℎ_i−ℎ_i^{new}) $$

Or in 3D,

$$ f_{i,j}=ρg∆A(ℎ_{i,j}−ℎ_{i,j}^{new}) $$

✅ 阿基米得定律：物体受到的浮力 = 排出去的水的重力
✅ 同时要考虑旋转和力矩。但目前旋转的效果不太好，可以考虑改进为隐式积分。
✅ 流体对方块的效果。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P33

A Summary For the Day

• The shallow wave model simulates waves over a height field.

• It’s based on a lot of simplification. We will discuss what fluid dynamics really looks like without simplification.

• The strength of the shallow wave model is its simplicity and efficiency. It can easily simulate water-solid coupling too.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• A Grid Representation and Finite Differencing

• Incompressible, Viscous Navier Stokes’ equations

✅ Incompressible：粒子不可压缩
✅ Viscous：有粘滞

• Air and liquid

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P3

A Grid Representation and Finite Differencing

P4

A Regular Grid Representation

✅ 把场定义在标准格子上的好处：(1)把物理量定义在格子的中心（2）计算导数或利用导数进行微分计算变得容易了。
✅ 上节课grid用1D来表示2D，2D表示3D，不是真正的grid方法。
🔎 Central Differencing：L10.

P6

Finite Differencing on Grid

一阶层数

The grid is very friendly with central differencing.

P7

二阶导数

The grid is very friendly with central differencing.

P8

Discretized Laplacian

We can then obtain the discretized Laplacian operator on grid.

$$ \frac{∂^2f_{i,j}}{∂x^2}≈\frac{\frac{∂f_{i−0.5,j}}{∂x}−\frac{∂f_{i+0.5,j}}{∂x}}{ℎ}≈\frac{f_{i−1,j}+f_{i+1,j}−2f_{i,j}}{ℎ^2} $$

$$ \frac{∂^2f_{i,j}}{∂y^2}≈\frac{\frac{∂f_{i,j+0.5}}{∂y}−\frac{∂f_{i,j−0.5}}{∂y}}{ℎ} ≈\frac{f_{i,j−1}+f_{i,j+1}−2f_{i,j}}{ℎ^2} $$

✅ 网格上的Laplace算子。

P9

Boundary Conditions

The boundary condition specifies \(f_{i−1,j}\) if it’s outside.

A Dirichlet boundary: \(f_{i−1,j}=C\)

A Neumann boundary: \(f_{i−1,j}=f_{i,j}\)

✅ 至少有一个边界使用Dirithlet．否则会全部收缩为0．
✅ Neumann是约束相对关系，没有绝对数值，会有无穷多解。

P12

Diffusion

The process of applying Laplacian smoothing is called diffusion.

✅ Laplace的本质是与邻居做平均。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P13

Problem with Central Differencing

Central differencing gives the derivative in the middle.

• The cell doesn’t exist at (i+0.5, j).

• To get \( \frac{∂f_ {i,j}}{∂x} \), we need \(f_{i−1,j}\) and \(f_{i+1,j}\). But this is weird, because \(f_{i,j}\) is unused.

✅ 前面假设所有物理量定义在格子的中间。但此处算出来的一阶微分量不在格子中间。

P14

Solution: Staggered Grid

✅ 不规定所有物理量都定义在格子中间，也可以定义在墙上。

We define some physical quantities on faces, specifically velocities.

• The x-part of the velocity is defined on vertical faces.

• The y-part of the velocity is defined on horizonal faces.

✅ 把速度定义在墙上的好处量，速度是矢量、可以用不同方向的墙表达不同方向上的速度、直观。

• Intuitively, they represent the flow speed between two cells. For example, we write the volume changing speed at cell (i,j) as:

✅ 通过四面墙上的速度计算当前格子的净流出（注意正负号）

P15

Divergence-Free Condition

No volume change is equal to say the fluid is incompressible. This can be formally written as a divergence-free velocity field.

❓ 这一页没听懂、净流入流出为0，水面还怎么动呢？
✅ 由于格子不可压，每个格子的净流出（入）应该为0．
✅ \(\nabla\)为散度符号，见前面课程。
✅ 公式1为直观理解，公式2为数学推导，本质上是一致的。

P16

Bilinear Interpolation

🔎 双线性插值：见GAMES 101

P17

We use bilinear interpolation to interpolate staggered velocities as well.

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P18

Incompressible, Viscous Navier-Stokes Equations

P19

Equation Fomulation

✅ 这是一个描述了速度场的公式，它可以告许你速度如何更新。公式1限制不可压，公式2 diffusion 的目的是粘滞。

Method of Characteristics: solving a long partial differential equation (PDE) in steps

• Step 1: Update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=\mathbf{g}\)
• Step 2: Update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=−(\mathbf{u}\cdot ∇)\mathbf{u}\)
• Step 3: Update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=υ∆\mathbf{u}\)
• Step 4: Update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=−∇\mathbf{p}\)

✅ 把偏微分方程分解几个小块，依次轮流优化每一小块。
❓ 这种方法为什么可行？

P20

Step 1: External Acceleration

The Update of \(\mathbf{u}\) by \(∂\mathbf{u}∕∂t=\mathbf{g}\) is straightforward, just add acceleration to \(u\) and \(v\).

✅ \(v_{i,j}\)代表向下的速度，对所有格子更新\(v_{i,j}\).
✅ 其它外部速度同理。

P21

Step 2: Advection

✅ Advection,代表流动。即速度会跟着粒子移动，基于欧拉的方法才需要考虑这个问题。因为固定的格子无法描述水的流动。
✅ 基于拉格朗日的方法，变量定义在粒子上，天然满足这个特点。

数学模型

Next we need to update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=−(\mathbf{u}\cdot ∇)\mathbf{u}\).

Solving this in an Eulerian way can be a source of instability.

✅ Eulerian way： \(\mathbf{u}^{\mathrm{new} }=\frac{\partial u}{\partial t} ·Δt＋\mathbf{u}\) 不稳定

To solve this problem, we come to realize that advection means to carry physical quantities by velocity.

P22

Solution: Semi-Lagrangian Method

The solution is to trace a virtual particle backward over time.

✅ 例如要求\(\mathbf{x}_0\)的速度，倒推哪个粒子会运动到\(\mathbf{x}_0\)处;因此找到\(\mathbf{x}_1\)，从\(\mathbf{x}_1\)的下一刻速度来更新\(\mathbf{x}_0\)的速度。

• Define \(\mathbf{x}_0←(i−0.5, j)\)
• Compute \(\mathbf{u}(\mathbf{x}_0)\)
• \(\mathbf{x}_1←\mathbf{x}_0−∆t \mathbf{u}(\mathbf{x}_0)\)

✅ 假设短时间内速度不变，根据当前速度猜测上一帧的位置。

• Compute \(\mathbf{u}(\mathbf{x}_1)\)
• \(u_{i,j}^{new}←u(\mathbf{x}_1)\)

Note that if the velocities are staggered, we need to do staggered bilinear interpolation.

P23

✅ 对每个墙上的速度都以相同的方式更新。

P24

We could also subdivided the time step for better tracing.

✅ 反推找\(\mathbf{x}_1\)时 step 细一点，这样能找得准一点
✅ 怎么计算每个\(\mathbf{x}\)的\(\mathbf{u}\)?答：双线性插值方法、
✅ 做模拟通常更在乎稳定而不是误差，此方法更稳定，但会有模糊的 artifacts.

P25

Step 3: Diffusion

Next we need to update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=\upsilon ∆\mathbf{u}\).

✅ 分别对\(u\)和 \(v\) 做 laplacian.
✅ 注意公式中\(v\)和\(\nu \)的不同，后者为粘滞系数。

We could also use even smaller sub-steps…

P27

Step 4: Pressure Projection

Finally, we need to update \(\mathbf{u}\) by solving \(∂\mathbf{u}∕∂t=−∇\mathbf{p}\).

Staggering makes this very straightforward:

$$ u_{i,j}^{new}←u_{i,j}−\frac{∆t}{ℎ}(p_{i,j}−p_{i−1,j}) $$

$$ v_{i,j}^{new}←v_{i,j}−\frac{∆t}{ℎ}(p_{i,j}−p_{i,j−1}) $$

But what is \(\mathbf{p}\)?

✅ 公式第二项离散化后在特定方向上的压强。
✅ \(u\)和\(v\)分别为两个方向上的速度。

P28

压强的来源

The pressure is caused by incompressibility.

✅ 压强的原因：由于流体不可压缩、对于流体的压力会传导到每个点上。
✅ 每个点都有压强，虽然压强未知，但可以根据不可压条件构造方程组。
✅ 不可压的表现为有压强，产生的效果是散度为0．

In other words, after this update by pressure, we should achieve:

which means

$$ \Downarrow $$

P29

压强的数学模型

The pressure is caused by incompressibility. Eventually, we get a Poisson equation:

Eventually, we get a Poisson equation:

$$ 4p_{i,j}−p_{i−1,j}−p_{i+1,j}−p_{i,j−1}−p_{i,j+1}= \\ \\ ℎ(−u_{i+1,j}−v_{i,j+1}+u_{i,j}+v_{i,j}) $$

with boundary conditions:

$$ \text{Dirichlet boundary (open) } p_{i−1,j}=P \\ \text{Neumann boundary (close) } p_{i−1,j}=p_{i,j}$$

Once we solve \(\mathbf{p}\), we update \(\mathbf{u}\) and done.

P30

After-Class Reading

Jos Stam. 1999. Stable Fluids. TOG (SIGGRAPH).

✅ 这篇论文主要讨论了step2，但也包含了全部过程

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P31

Air and Smoke

✅ 前面讲的是怎么更新速度；后面讲怎么利用速度做出效果。

P32

Air Simulation

• Air simulation is done in two steps.
• In Step 1, we update the flow (the velocity field) \(\mathbf{u}\).
• In Step 2, we use semi-Lagrangian (page 22) advect all of the other physical quantities, i.e., density, temperature…
• Typically we use Dirichlet boundaries for an open space (or Neumann boundaries for a container.)
• We can use it to simulate underwater as well.

P33

Water Simulation

✅ 要渲染的不是水，而是水与空气的接触面。但通常只模拟水不模拟空气。

• Two representations
  • Volume-of-fluid (as the name suggests…)
  • A signed distance function defined over the grid.

✅ 表示1：例如一个格子存储水的体积的百分化。用于早期，无法描述水的界面，因此不精准。

• How to advect(更新)?
  • Semi-Lagrangian (volume loss)
  • Level set method (volume loss)
  • Needs corrections.

✅ advect 2：专用于更新 SDF 的方法。
✅ 水变少是常见问题，两种advect都存在。

P35

After-Class Reading

Osher and Fedkiw.
Level Set Methods and Dynamic Implicit Surfaces.

✅ 介绍流体模拟的很好的书。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P36

A Summary For the Day

• The Eulerian grid presentation is very friendly with finite differencing. This makes calculus a lot easier.

• For velocity fields, we can use staggered grid.

• For low-speed, incompressible, viscous flow, we need to solve the Navier-Stokes equations.

• To solve the equations, we can do this in step-by-step (method of characteristics).

• To simulate air and water, we need to advect some physical quantities.

  • Smoke (density); water (volume-of-fluid, or signed distance function)
  • Volume loss issue in water (how to fix it?)
  • If you need to create a mesh from grid for rendering, you need something like marching cube.

✅ 用有限元方法也能模拟水，但难以模拟流动性。流动需要对四面体重新构造，成本很高。

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P2

Topics for the Day

• A SPH model

✅ SPH = Smoothed Particle Hydrodynamics

• SPH-based fluids

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

A SPH Model

✅ 搞一个模型，能够用于计算偏微分，把微分向量应用到方程上进行求解。
✅ SP:smooth particle

P4

A SPH Model

Consider a (Lagrangian) particle system: each water molecule is a particle with physical quantities attached, such as position \(\mathbf{x}_i\), velocity \(\mathbf{v}_i\), and mass \(m_i\).

✅ 用粒子来表达流体，物理变量附着在粒子上，粒子转化为三角网格再渲染，或直接渲染带透明贴图的粒子(游戏)。

P5

原理

• Suppose each particle j has a physical quantity \(A_j\).
• The quantity can be: velocity, pressure, density, temperature….
• How to estimate the quantity at a new location \(\mathbf{x}_i\)?

✅ 空间中有很多带有物理量的粒子，求任意位置上的物理量。这是插值问题，关键是要插值结果平滑。

模型

A Simple Model

$$ \begin{matrix} A_i^{\mathbf{smooth}}=\frac{1}{n}\sum _jA_j & \text{ For } ||\mathbf{x}_i−\mathbf{x}_j||<R \end{matrix} $$

存在的问题

✅ 取平均的方式没有考虑粒子的分布。

P7

A Better Model

• Let us assume each one represents a volume \(V_j\).
• So a better solution is:

$$ \begin{matrix} A_i^{\mathbf{smooth} }=\frac{1}{n}\sum_jV_jA_j & \text{ For } ||\mathbf{x} _i−\mathbf{x} _j||<R \end{matrix} $$

✅ 公式假设总球的体积是1，球内的粒子瓜分这些体积。所以\(\sum _jV_j=1\)

P8

存在的问题

• One problem of this solution:

$$ \begin{matrix} A_i^{\mathbf{smooth} }=\frac{1}{n}\sum_jV_jA_j & \text{ For } ||\mathbf{x} _i−\mathbf{x} _j||<R \end{matrix} $$

• Not smooth! (7 -> 9!)

✅ 微小的移动，圆内多了两个点，导致结果突变。

P9

Final Solution

• Final solution:

$$ \begin{matrix} A_i^{\mathbf{smooth}}=\sum _ j V_jA_jW_{ij} & \text { For } ||\mathbf{x} _ i− \mathbf{x} _j||< R \end{matrix} $$

• \(W_{ij}\) is called smoothing kernel.
• When \(||\mathbf{x} _ i − \mathbf{x} _ j||\) is large, \(W_{ij}\) is small.
• When \(||\mathbf{x} _ i−\mathbf{x} _ j||\) is small, \(W_{ij}\) is large.

P10

Particle Volume Estimation

• But how do we get the volume of particle \(i\)?
  $$ V_i=\frac{m_i}{ρ_i} $$

$$ ρ_i^ \mathbf{smooth} =\sum _ j V_ j ρ_ j W _ {ij}= \sum _ jm_jW_{ij} $$

✅ 粒子在运动过程中，疏密会有变化，因此体积不是常数，要实时计算。
✅ 公式中的\(\rho \)不是指水的密度，而是粒子分布的密度。
✅ 把密度当作粒子的物理量。用同样的方法插出某个点的密度。

P11

Smoothed Interpolation – Final Solution

• So the actual solution is:

P12

Kernal函数

Kernal函数的作用

• We can easily compute its derivatives:

  • Gradient

  $$ \begin{matrix} A_i^ \mathbf{smooth} = \sum _ jV_jA_ jW_ {ij} \quad & ∇A_i ^\mathbf{smooth} = \sum_jV_jA_j∇W_ {ij} \end{matrix} $$

  • Laplacian

  $$ \begin{matrix} A_i^ \mathbf{smooth} = \sum _ j V_ j A_ jW_ {ij} \quad & ∇A_i^\mathbf{smooth} = \sum_ jV_ jA_ j∇W_ {ij} \end{matrix} $$

❓ 为什么认为体积是常数？答：假设一个点的运动不影向周围邻居的体积。
✅ 对于当前点来说，周围粒子的物理量是常数，只有\(W_{ij}\)与当前点有关。
✅ 而\(W_{ij}\)来自于已知的kernel函数，其derivative也是已知的。

P13

A Smoothing Kernel Example

$$ W_{ij}=\frac{3}{2\pi h^3} \begin{cases} \frac{2}{3}-q^2+\frac{1}{2} q^3 \quad & (0\le q<1) \\ \frac{1}{6}(2-q)^3 \quad& (1\le q<2) \\ 0 \quad & (2\le q) \end{cases} $$

$$ q=\frac{||\mathbf{x} _i-\mathbf{x} _j||}{h} $$

\(h\) is called smoothing length

✅ smooth Kernal 有很多种，这种最常见。

P14

Kernel Derivatives

• Gradient at particle i (a vector)

$$ \nabla _ i W _ {ij} = \begin{bmatrix} \frac{\partial W _ {ij}}{\partial x _ i} \\ \frac{\partial W _ {ij}}{\partial y _ i} \\ \frac{\partial W _ {ij}}{\partial z _ i} \end{bmatrix} = \frac{\partial W_ {ij}}{\partial q} \nabla _ iq= \frac{\partial W _ {ij}}{\partial q} \frac{\mathbf{x} _ i-\mathbf{x} _ j}{|| \mathbf{x} _ i - \mathbf{x} _ j||h} $$

$$ q=\frac{||\mathbf{x} _i-\mathbf{x} _j||}{h} $$

$$ W_{ij}=\frac{3}{2\pi h^3} \begin{cases} \frac{2}{3}-q^2+\frac{1}{2} q^3 \quad & (0\le q<1) \\ \frac{1}{6}(2-q)^3 \quad& (1\le q<2) \\ 0 \quad & (2\le q) \end{cases} $$

$$ \frac{\partial W_{ij}}{\partial q} =\frac{3}{2\pi h^3} \begin{cases} -2q+\frac{3}{2}q^2 \quad & (0\le q<1) \\ -\frac{1}{2}(2-q)^2 \quad& (1\le q<2) \\ 0 \quad & (2\le q) \end{cases} $$

P15

Kernal Laplacian

$$ \frac{\partial W_{ij}}{\partial q} =\frac{3}{2\pi h^3} \begin{cases} -2q+\frac{3}{2}q^2 \quad & (0\le q<1) \\ -\frac{1}{2}(2-q)^2 \quad& (1\le q<2) \\ 0 \quad & (2\le q) \end{cases} $$

$$ \frac{\partial^2 W_{ij}}{\partial q^2} =\frac{3}{2\pi h^3} \begin{cases} -2+3q \quad & (0\le q<1) \\ 2-q \quad& (1\le q<2) \\ 0 \quad & (2\le q) \end{cases} $$

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P16

SPH-Based Fluids

P17

Fluid Dynamics

• We model fluid dynamics by applying three forces on particle i.
  • Gravity
  • Fluid Pressure
  • Fluid Viscosity

P18

Gravity Force

• Gravity Force is:

$$ \mathbf{F} _ \mathbf{i}^ \mathbf{gravity} = m _i \mathbf{g} $$

P19

Pressure Force

✅ 要解决的问题：怎么计算压强？怎么把压强转化为力？

怎么计算压强

• Pressure is related to the density

  • First compute the density of Particle i:

  $$ \rho _ i = \sum _ j m _ j W _ {ij} $$

  • Convert it into pressure (some empirical function):

  $$ P_i=k((\frac{\rho _i}{\rho _\mathrm{constant } } )^7-1) $$

✅ 密度到压强的计算是一个经验公式。

压强转化为力

P20

• Pressure force depends on the difference of pressure:

P21

• Mathematically, the difference of pressure => Gradient of pressure.

$$ \mathbf{F} _i^{pressure}=-V_i\nabla _iP^{smooth} $$

✅ 体积为粒子在空间中占有的体积，体积越大受到的压力越大、\(\nabla\)代表压强的差。

• To compute this pressure gradient, we assume that the pressure is also smoothly represented:

$$ P_i^{smooth}= \sum _ j V_jP_j W_{ij} $$

✅ 假设空间是一个压强场、粒子是空间中的采样。\(P^{smooth}\)是通过周粒子\(P\)的插值得到的采样点压强。

• So:

$$ \mathbf{F} _ i^{pressure} = - V _ i \sum _ j V _ j P _ j \nabla _ i W _ {ij} $$

P22

Viscosity Force

粘滞所产生的效果

• Viscosity effect means: particles should move together in the same velocity.
• In other words, minimize the difference between the particle velocity and the velocities of its neighbors.

✅ Viscosity (粘滞)类似于 damping (阻尼)，但有些区别，后者的目标是让粒子的运动停下来，前者的目的是让所有粒子的运动整齐划一，即速度差趋于0.

P23

粘滞力 Viscosity Force

• Mathematically, it means:
  $$ \mathbf{F} _i^{vis \cos ity}=-\nu m_i\Delta _i\mathbf{V} ^{smooth} $$

✅ \(V\)：粘滞系数， \(\nabla V\)：速度的 Laplacian.注意速度是3D矢量。

• To compute this Laplacian, we assume that the velocity is also smoothly represented:

$$ \mathbf{V} _i^{smooth}= \sum_jV_j \mathbf{v} _ j W _ {ij} $$

• So:

$$ \mathbf{F} _i^{vis \cos ity}=-\nu m_i\sum _jV_j\mathbf{v} _j\Delta _iW _{ij} $$

✅ smooth会产生粘滞的效果。

P24

Algorithm

• For every particle i
  • Compute its neighborhood set
  • Using the neighborhood, compute:
    • Force = 0
    • Force + = The gravity force
    • Force + = The pressure force
    • Force + = The viscosity force
• Update \(v_i = v_i + t * \text{ Force } / m_i\);
• Update \(x_i = x_i + t * v_i\);

✅ 性能瓶颈：计算邻居，因为总粒子数为百万级。

Spatial Partition加速求最近邻

P25

Exhaustive Neighborhood Search

• Search over every particle pair? O(\(N^2\))
• 10M particles means: 100 Trillion pairs…

P26

Solution: Spatial Partition

• Separate the space into cells
• Each cell stores the particles in it
• To find the neighborhood of i, just look at the surrounding cells

P27

Spatial Partition

• What if particles are not uniformly distributed?

✅ 例如水花喷溅的效果，通常靠近水面的粒子小一点，更利于表现细节。

• Solution: Octree, Binary Spatial Partitioning tree…

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/

P28

Fluid Display

• Need to reconstruct the water surface from particles!

✅ 点云转成三角面片用于渲染也是一个比较复杂的问题。
✅（1）平滑方法：bias kemal（见GAMES 102）
✅（2）把球转为SDF，SDF转为Mesh

P29

Ongoing Research

• How to make the simulation more efficient?

• How to make fluids incompressible?

• When simulating water, only use water particles, no air particles. So particles are sparse on the water-air boundary. How to avoid artifacts there?

• Using AI, not physics, to predict particle movement?

本文出自CaterpillarStudyGroup，转载请注明出处。

https://caterpillarstudygroup.github.io/GAMES103_mdbook/