GAMES 图形学系列笔记（九）
for 每个内部顶点 v_i:
Vector3f delta(0, 0, 0);
float weight_sum = 0;
for 每个邻接顶点 v_j:
float cot_alpha = cot(角 alpha_ij);
float cot_beta = cot(角 beta_ij);
float weight = cot_alpha + cot_beta;
delta += weight * (v_j - v_i);
weight_sum += weight;
}
// 简单起见，这里使用均匀拉普拉斯平滑作为演示
// 实际应使用cot权重公式
Vector3f laplacian = delta / weight_sum;
v_i_new = v_i + lambda * laplacian;
}
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_101.png>
## Utopia框架使用指南 🛠️
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_103.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_105.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_107.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_109.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_111.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_113.png>
以下是Utopia框架的核心架构与使用方法简介。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_115.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_117.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_119.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_121.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_123.png>
### 框架架构
框架主要分为四层：
1. **基础层**：提供数学库、反射系统等基础功能。
2. **渲染层**：基于DirectX 12的渲染管线，管理着色器、材质、网格等渲染资源。
3. **逻辑层**：采用ECS（实体-组件-系统）架构组织代码和数据。
* **实体**：世界的对象，只是一个ID。
* **组件**：附加在实体上的数据（如位置、网格、材质）。
* **系统**：处理具有特定组件组合的实体的逻辑。
4. **编辑器层**：提供场景编辑、属性调试的可视化界面。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_125.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_127.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_129.png>
### 完成作业的关键步骤
作业目标是实现离散平均曲率流算法来生成极小曲面。
1. **创建数据与系统**：在ECS框架下，创建一个组件（如`DenoiseData`）来存储网格数据，创建一个系统（如`DenoiseSystem`）来执行算法逻辑。
2. **操作半边网格**：框架提供了`HalfedgeMesh`库。你需要：
* 将渲染网格转换为半边网格结构。
* 遍历内部顶点，计算其一环邻域重心（使用cot权重公式）。
* 根据更新公式调整顶点位置。
* 将修改后的半边网格数据传回渲染网格。
3. **可视化曲率（可选）**：计算网格的高斯曲率或平均曲率，并将曲率值映射为顶点颜色，在材质中显示。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_131.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_133.png>
在编辑器中，你可以将网格对象拖拽到你的数据组件中，并通过系统控制面板触发算法的执行和迭代。
## 课程总结 🎯
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_135.png>
本节课我们一起学习了以下内容：
1. 三角网格作为离散曲面的两种理解方式。
2. 用于表示网格拓扑的**半边数据结构**及其基本操作。
3. 微分几何的基础概念：切平面、法向、主曲率、高斯曲率与平均曲率。
4. 如何将光滑曲面的微分性质**离散化**到三角网格上，进行法向和曲率的估计。
5. **极小曲面**的概念，以及利用**离散平均曲率流**迭代算法来生成极小曲面的原理与步骤。
6. 用于完成作业的**Utopia框架**的基本架构、ECS设计模式和使用方法。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_137.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_139.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_141.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_143.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_145.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/d4185f64a49599b0e4f858aaabb3a4b5_146.png>
通过本课的学习，你将能够理解离散曲面的微分性质计算方法，并利用提供的框架动手实现一个简单的几何处理算法。
# GAMES102-几何建模与处理---P9-微分坐标---GAMES-Webinar---BV1NA411E7Yr_note
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_0.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_2.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_3.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_5.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_7.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_9.png>
在本节课中，我们将学习微分坐标（拉普拉斯坐标）的概念及其在几何处理中的应用。我们将从回顾作业六开始，逐步理解拉普拉斯算子的几何意义、离散形式，并探讨其在网格光滑化、参数化和编辑等任务中的核心作用。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_11.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_13.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_15.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_17.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_19.png>
---
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_21.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_23.png>
## 作业六回顾与总结 📊
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_25.png>
上一节我们介绍了网格的基本概念和数据结构。本节中，我们来看看作业六的完成情况，它为我们理解微分坐标打下了实践基础。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_27.png>
作业六要求对三角网格进行编程，通过迭代修改顶点坐标，使其朝着拉普拉斯方向移动，最终使网格变得光滑。
以下是作业完成情况的要点：
* **核心操作**：将每个顶点向其**一邻域**（直接相连的顶点）的加权平均位置移动。
* **迭代过程**：通过不断迭代上述操作，网格逐渐光滑，最终在边界固定的条件下逼近**极小曲面**。
* **观察现象**：
* 顶点密集的区域（如兔子耳朵）收敛速度较慢。
* 使用均匀权重可能导致网格自交（面片翻转）。
* 使用与几何相关的权重（如余切权重）通常能获得更好的性质。
通过这个作业，我们直观地体验了拉普拉斯坐标作为“几何细节”或“尖锐度”度量的作用：移动顶点以减少其拉普拉斯坐标的长度，相当于磨平该处的几何特征。
---
## 离散曲面与图结构 🔗
在深入微分坐标之前，我们需要巩固对离散曲面的理解。离散三角网格本质上是一个图（Graph）在三维空间中的嵌入。
* **两种观点**：
1. **参数曲面观点**：将网格视为从二维参数域到三维空间的映射。每个顶点拥有三维坐标。
2. **图嵌入观点**：将网格视为一个二维图，其顶点被赋予了三维坐标。其本质是二维流形。
* **数据结构的重要性**：处理网格的核心是掌握图的数据结构（如半边结构）。数据结构的选择需要在**存储空间**和**计算时间**之间权衡。没有最好的结构，只有最适合当前应用的结构。
* **邻域概念**：要研究一个顶点的局部性质（如曲率），我们需要考察其邻域。
* **一邻域**：与该顶点直接通过边相连的所有顶点集合。
* **k邻域**：可通过不超过k条边到达该顶点的所有顶点集合。通常，一邻域足以近似“无穷小”邻域的性质。
---
## 拉普拉斯坐标（微分坐标）的定义与几何意义 🧮
理解了局部邻域后，我们就可以正式引入本节课的核心概念——拉普拉斯坐标，也称为微分坐标。
对于一个顶点 **v_i**，其拉普拉斯坐标 **δ_i** 定义为该顶点与其一邻域加权平均位置的差向量：
**δ_i = v_i - Σ_{j∈N(i)} ω_{ij} v_j**
其中，**N(i)** 是顶点 **v_i** 的一邻域，**ω_{ij}** 是权重，满足 Σ ω_{ij} = 1。
* **几何意义**：这个向量 **δ_i** 的长度和方向刻画了该顶点偏离其局部邻域所构成平面的程度。长度越大，该点越“尖锐”；长度为零，则该点与其邻域共面，局部平坦。
* **与平均曲率的关系**：在连续曲面的离散化中，可以证明，当网格足够精细时，拉普拉斯坐标 **δ_i** 近似等于该点的**平均曲率向量**，即：**δ_i ≈ H_i * n_i**，其中 **H_i** 是平均曲率，**n_i** 是法向。
* **权重的选择**：权重 **ω_{ij}** 的选择至关重要。从离散微分几何推导出的**余切权重**（Cotangent Weight）具有优良的几何性质，通常比均匀权重效果更好，能减少网格变形和自交。
---
## 拉普拉斯光滑与全局方法 ⚙️
上一节我们通过迭代局部操作实现了光滑化。本节中，我们来看看如何从全局角度一次性求解光滑结果。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_29.png>
* **拉普拉斯光滑（局部迭代法）**：即作业六使用的方法，公式为：
**v_i‘ = v_i + λ * δ_i**
通过不断迭代，逐渐减小 **δ_i**，实现光滑。λ 是步长参数。
* **全局方法（一次性求解）**：目标是找到所有顶点的新位置 **V‘**，使得每个内部顶点的拉普拉斯坐标为零（或按比例缩小），同时固定边界顶点。这可以归结为求解一个线性方程组：
**L * V‘ = B**
其中：
* **L** 是**拉普拉斯矩阵**，一个大型稀疏矩阵，编码了网格的邻接关系和权重。
* **V‘** 是未知的顶点坐标向量（x, y, z分量需分别求解）。
* **B** 是由边界条件构成的右侧向量。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_31.png>
* **优势对比**：
* **局部迭代法**：实现简单，但可能收敛慢，且易产生自交。
* **全局求解法**：通过求解稀疏线性方程组，可直接得到结果，数值上更稳定。这是作业七的核心内容，需要学习使用数值计算库（如Eigen）来求解此类方程组。
---
## 参数化初步应用 🗺️
参数化是将三维网格映射到二维平面的过程，是纹理映射等应用的基础。利用我们刚学的全局拉普拉斯方法，可以实现一种简单的参数化。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_33.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_35.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_37.png>
* **基本思想**：
1. 将三维网格的**边界顶点**映射到二维平面上的一个凸多边形（如正方形）边界上。
2. 对于内部顶点，在二维平面上求解拉普拉斯方程 **L * U = B**，其中 **U** 是内部顶点的二维坐标，**B** 由映射后的边界顶点决定。
3. 解出的 **U** 即为每个顶点在二维参数平面上的坐标。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_39.png>
* **特点**：
* 这种方法称为**保角参数化**的线性近似，计算简单高效。
* 可以证明，如果边界映射到凸多边形，则生成的参数化不会发生三角形翻转。
* 缺点是可能产生较大的面积扭曲，特别是当三维网格形状与平面区域差异大时。
* **应用**：生成的二维坐标（UV坐标）即可用于纹理映射，将二维图像贴合到三维模型表面。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_41.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_43.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_45.png>
---
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_47.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_49.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_51.png>
## 约束与网格编辑 ✏️
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_53.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_55.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_57.png>
最后，我们探讨如何利用拉普拉斯坐标进行受约束的网格变形和编辑。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_59.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_61.png>
* **硬约束与软约束**：在全局求解框架中，我们可以方便地加入约束。
* **硬约束**：某些顶点必须固定在指定位置。这通过将这些顶点的坐标从变量变为已知量（移到方程右侧）来实现。
* **软约束**：希望某些顶点尽量靠近但不必须到达指定位置。这通过将其作为最小二乘项加入目标函数来实现。
* **网格编辑的应用**：例如，用户拖动模型上的一个点，希望模型其他部分随之自然变形。
* **核心思想**：在变形过程中，尽量保持每个顶点的**拉普拉斯坐标**（即局部细节）不变。但为了适应旋转，需要对拉普拉斯坐标的方向进行相应的旋转变换。
* **数学形式**：这通常转化为一个带约束的**最小二乘优化问题**，目标是最小化拉普拉斯坐标的变化，约束条件是用户指定的顶点位移。
---
## 总结与作业预告 📚
本节课中，我们一起学习了微分坐标（拉普拉斯坐标）这一核心概念。
* **核心概念**：拉普拉斯坐标 **δ_i = v_i - Σ ω_{ij} v_j** 是顶点与其局部邻域平均的差，是局部几何细节的度量。
* **关键应用**：
1. **网格光滑**：通过减小 δ_i 实现。
2. **参数化**：在二维平面上求解 L * U = B，将网格展开。
3. **网格编辑**：保持 δ_i 的变换，实现细节保持的变形。
* **方法演进**：从**局部迭代**的直观方法，到构建**全局拉普拉斯矩阵**并求解线性方程组的稳健方法。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_63.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/f7a007474d09cc323264b272f2f53de2_64.png>
**作业七预告**：请使用全局求解方法（构建拉普拉斯矩阵并调用Eigen等库求解方程组）重新实现网格光滑化，并进一步将其扩展，实现简单的网格参数化。重点在于掌握稀疏矩阵的构建与线性方程组的求解。
# GAMES103-基于物理的计算机动画入门---P1-Lecture-01-Intro-to-Physics-Based-Animation---GAMES-Webinar---BV12Q4y1S73g_note
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_0.png>
在本节课中，我们将要学习这门课程的整体框架、学习目标、课程机制以及基于物理的计算机动画的基本概念。我们将从图形学的基础知识开始，逐步深入到物理模拟的核心领域。
## 课程概述与机制 📋
本课程名为“基于物理的计算机动画入门”。课程主要讨论如何将物理模拟技术应用于计算机动画中，并介绍相关技术的基本原理和算法。课程包含相应的编程作业以巩固所学知识。
我是课程讲师王华明。课程相关问题可以发送至课程邮箱，会有助教或老师解答。课程资料、回放、PPT及作业将发布在GAMES论坛上。
课程时间为每周一16:00至17:30，时长约一个半小时，具体可能根据内容调整。课程后半段会预留时间进行答疑。
从第三周开始，课程将转为线下授课与线上转播结合的形式。线下地点初步定在林地科技会议室（浙江大学紫金港校区附近），方便杭州及周边地区的同学参与。
课程将持续12周，大约到春节前两周结束。作业批改和课程辅助将由助教团队负责。
## 预备知识要求 📚
为了降低学习门槛，本课程仅要求学员具备基础的数学知识和编程能力。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_2.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_4.png>
以下是所需的核心知识：
* **数学基础**：需要掌握线性代数（矢量、矩阵、线性系统、特征值分解等概念）和微积分（求导、积分、链式法则、梯度、泰勒展开等基本概念和计算）。
* **编程基础**：需要具备C、C++、C#或Java等语言的编程能力。课程作业将使用Unity引擎，其脚本语言为C#。有相关语言经验即可快速上手。
* **图形学基础**：了解简单的图形学概念（如变换、旋转）和渲染基本原理即可，无需深入复杂的图形学知识。
课程会涵盖数值算法、偏微分方程、有限元分析、流体力学等进阶内容，但即使没有这些背景知识，也可以通过课程学习掌握。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_6.png>
## 课程工具与环境 🛠️
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_8.png>
课程作业将使用Unity引擎完成。Unity对学生和个人用户免费，对硬件要求较低，大部分计算基于CPU。
以下是关于工具的重点说明：
* **Unity版本**：对版本要求不严格，较新的版本均可。支持Windows和macOS系统。
* **学习重点**：使用Unity主要是将其作为学习和实践的工具，课程核心是学习物理模拟算法本身，而非Unity引擎的使用。大部分作业将替代Unity原有的物理引擎，通过编写脚本来实现模拟。
* **资源获取**：Unity软件可从官网下载，其官方论坛可以解答通用的引擎使用问题。
## 课程大纲与作业 📅
课程总计12周，涵盖物理模拟的多个核心方向，内容相对独立，便于学习。
以下是课程的大致安排：
* **第1周**：课程简介（本周）。
* **第2周**：数学基础回顾，结合图形学实例讲解。
* **第3-4周**：刚体动力学与刚体碰撞处理。
* **第5-7周**：布料模拟，引入PBD、Projective Dynamics等算法。
* **第8-9周**：软体动力学与有限元方法。
* **第10-12周**：流体模拟，包括表面波、网格法和粒子法。
课程包含四次编程作业，分别对应刚体、布料、软体和流体四个方向。每次作业包含基础任务和可选的高级任务。完整完成所有作业的同学将获得课程纪念品。本课程没有考试。
课程没有固定教材，但每节课会提供相关的论文或文献作为选读材料，建议有精力的同学课后阅读以加深理解。
## 什么是计算机图形学？ 🖥️
计算机图形学简而言之，就是研究如何构建三维虚拟世界，并将其以二维图像形式呈现出来的学科。它与计算机视觉方向相反，后者是从二维图像理解三维世界。
图形学主要包含三个方向：
1. **几何**：研究如何构造和表达三维虚拟世界。
2. **渲染**：研究如何将三维世界转化为二维图像并显示出来。
3. **动画**：研究如何让三维世界中的物体运动起来。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_10.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_12.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_14.png>
一个理想的实时图形学管线是：几何处理（可离线） -> 动画模拟（实时） -> 渲染输出（实时） -> 显示。帧率是衡量实时性的关键指标，例如电影通常为24帧/秒，而交互性强的游戏可能需要60帧/秒或更高。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_16.png>
## 图形学的应用领域 🌐
计算机图形学技术已广泛应用于多个行业：
* **娱乐产业**：如电子游戏、电影特效、社交媒体滤镜和虚拟数字人等。
* **设计与工程**：如计算机辅助设计、建筑设计、时尚设计等。
* **电子商务与智能制造**：如虚拟试衣、数字产品展示，连接设计与生产。
* **前沿领域**：如元宇宙、VR/AR/MR，构建沉浸式虚拟世界的基石。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_18.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_20.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_22.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_24.png>
## 什么是基于物理的动画？ ⚙️
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_26.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_28.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_30.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_32.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_34.png>
动画的本质是在离散的时间点上更新物体的状态。两个时间点之间的间隔称为时间步长 `Δt`。物理模拟的核心问题就是在每个时间步长内，如何根据物理定律更新物体的状态（如位置、速度、形状等）。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_36.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_38.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_40.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_42.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_44.png>
基于物理的动画主要模拟四大类物质：
1. **刚体**：假设物体无形变。
2. **布料与头发**：属于细薄物体。
3. **软体/弹性体**：物体可发生弹性形变。
4. **流体**：包括液体和气体。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_46.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_48.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_50.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_52.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_54.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_56.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_58.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_60.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_62.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_64.png>
相应地，有三种主要的模拟表达方式：
* **网格法**：用三角形或四面体网格表示物体表面或体积。适用于形态相对固定的物体，如刚体、布料、弹性体。
* **粒子法**：用一堆离散的点云表示物体。优点是无须处理网格拓扑，适用于流体、碎裂等效果。
* **网格法**：将空间划分为规则的小格子，在每个格子存储物理量。常用于流体、烟雾的模拟，内存消耗较大。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_66.png>
此外，还有**混合方法**，如物质点法，它结合了粒子和网格的优点，常用于模拟雪、沙等物质。不同物质间的耦合交互也是一个重要的研究课题。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_68.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_70.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_72.png>
## 本课程涵盖内容 🎯
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_74.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_76.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_78.png>
结合课程安排与讲师专长，本课程将重点讲解以下内容：
* **刚体**：刚体动力学和碰撞处理，不涉及破碎模拟。
* **布料**：主要讲解布料模拟，头发模拟暂不涉及。
* **软体**：弹性体模拟、有限元方法及超弹性模型。
* **流体**：涵盖表面波、基于网格的不可压缩流体模拟以及基于粒子的流体模拟。
* **通用技术**：专门讲解碰撞检测与处理，以及基于约束的动力学方法。
## 总结与学习建议 💡
本节课我们一起学习了课程的基本信息、预备要求、工具使用以及基于物理的动画的核心概念。
成功的学习需要做到以下几点：
1. **做好准备**：掌握必要的数学和编程基础，提前熟悉Unity的基本操作。
2. **积极参与**：按时参加课程（或观看回放），认真完成编程作业。
3. **主动拓展**：利用课余时间阅读推荐的文献资料，深入理解算法原理。
4. **多读多写多想**：这是掌握任何技术的关键。通过阅读积累知识，通过编程实践深化理解，通过思考融会贯通。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_80.png>
希望本课程能帮助大家对物理模拟算法和计算机动画建立一个扎实的基础。未来如果大家兴趣浓厚，我们也可能开设相关的高级课程。
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_82.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/447852504d67eea817fbb890cb3a6fcc_84.png>
我们下周将从必要的数学基础开始，并结合图形学中的具体问题展开讲解。
# GAMES103-基于物理的计算机动画入门---P10-Lecture-10-Surface-Waves--Lab-4----GAMES-Webinar---BV12Q4y1S73g_note
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/0bbffc89557c994fb814abda28dacfef_0.png>
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/0bbffc89557c994fb814abda28dacfef_2.png>
在本节课中，我们将学习如何模拟水面上的波浪效果。我们将介绍一种名为“浅水波”的简化模型，它非常适合在游戏等实时应用中模拟水面波动。课程内容将涵盖从基本概念到具体实现的完整流程，并最终与我们的实验作业相结合。
---
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/0bbffc89557c994fb814abda28dacfef_4.png>
## 概述：两种模拟方法
在深入浅水波模型之前，我们需要了解模拟流体的两种基本方法：拉格朗日法和欧拉法。
上一节我们介绍了弹性体和刚体的模拟，它们都属于拉格朗日法。本节中我们来看看这两种方法的区别。
* **拉格朗日法**：将物理属性（如速度、密度）定义在随物质一起运动的“质点”上。例如，模拟一群水分子，每个分子都有自己的属性，并随水流移动。
* **欧拉法**：将空间划分为固定的网格（格子），物理属性定义在这些固定的空间位置上。当流体流过时，我们需要更新每个格子中的属性值。
我们今天要学习的浅水波模型，就是基于欧拉法的一种应用。
---
## 浅水波模型与高度场 🌊
我们的目标是模拟一个水面。一个直观的方法是使用**高度场**来描述它。
我们可以把水面想象成一个二维网格，每个网格点都有一个高度值 `h`，代表该处水面的海拔。通过更新整个网格的高度值，就能模拟出波浪传播的效果。这种描述被称为**2.5维高度场**。
除了高度，我们还需要知道水流的运动情况，因此引入**速度场** `u`。它定义了在每个水平位置 `x` 上，水流穿过该垂直截面的水平速度。
现在我们有了描述状态的**高度场** `h(x)` 和描述运动的**速度场** `u(x)`。接下来，我们需要找到更新它们的物理规则。
---
## 从物理方程到更新公式
根据流体力学的基本原理，我们可以推导出高度和速度随时间变化的方程。经过一系列假设和简化（特别是假设波浪很“浅”，即高度变化平缓），我们可以得到一个核心方程，它只与高度场和压强有关：
**∂²h/∂t² = (h / ρ) * (∂²p/∂x²)**
其中：
* `h` 是高度。
* `ρ` 是水的密度（常数）。
* `p` 是压强。
这个方程就是**浅水波方程**。它告诉我们，水面高度的二阶时间导数（即加速度）与压强的二阶空间导数成正比。
在只考虑重力的情况下，水下压强 `p` 可以简化为：**p = ρ g h**。其中 `g` 是重力加速度。代入上式，并合并常数，我们可以得到一个更简洁的、只关于高度 `h` 的方程：
**∂²h/∂t² = α * (∂²h/∂x²)**， 其中 **α = g h**（通常近似为常数）。
---
## 离散化与数值求解 🔢
计算机无法直接处理连续方程，我们需要将其**离散化**。我们将空间划分为许多小格子，时间也划分为小步长 `Δt`。
首先，我们需要用离散的格子值来近似方程中的导数。这里我们使用**中心差分法**。
对于时间的二阶导数，我们有近似公式：
**(hᵢⁿ⁺¹ - 2hᵢⁿ + hᵢⁿ⁻¹) / Δt² ≈ ∂²h/∂t²**
对于空间的二阶导数，我们有近似公式：
**(hᵢ₊₁ⁿ - 2hᵢⁿ + hᵢ₋₁ⁿ) / Δx² ≈ ∂²h/∂x²**
将这两个近似代入简化后的浅水波方程，并进行整理，我们就可以得到每个格子在下一时刻的高度更新公式：
**hᵢⁿ⁺¹ = 2hᵢⁿ - hᵢⁿ⁻¹ + α * (Δt²/Δx²) * (hᵢ₊₁ⁿ - 2hᵢⁿ + hᵢ₋₁ⁿ)**
这个公式非常直观：下一个时刻的高度，由当前时刻的高度、前一时刻的高度以及左右邻居的高度共同决定。
---
## 保持体积守恒与添加阻尼
直接使用上述公式模拟，可能会导致水的总体积（所有格子高度之和）发生变化，这不符合物理规律。为了保证体积守恒，我们需要对公式进行修正。
以下是两种常用方法：
1. **修改耦合系数**：将更新公式中与邻居交互的系数，从 `hᵢ` 改为 `(hᵢ + h邻居)/2`。这样可以保证从格子A流到格子B的水量，等于从格子B流到格子A的水量。
2. **常数化处理（作业采用）**：直接将公式中的 `α` 视为常数。这样在求和时，邻居项会相互抵消，从而自动保持体积守恒。
此外，真实的水有粘滞性，波浪会逐渐衰减。我们可以在更新公式中引入一个**阻尼系数 β**（0 < β ≤ 1）：
**hᵢⁿ⁺¹ = 2hᵢⁿ - hᵢⁿ⁻¹ + β * [ α * (Δt²/Δx²) * (hᵢ₊₁ⁿ - 2hᵢⁿ + hᵢ₋₁ⁿ) ]**
当 `β = 1` 时，无阻尼；`β` 越小，阻尼越大，波浪衰减越快。
---
## 处理边界条件 🧱
水面不可能无限大，我们需要定义边界处的行为。主要有两种边界条件：
* **狄利克雷边界条件**：固定边界处的高度为一个常数值（如 `H`）。这用于模拟开放水域（如海洋），边界外是静止的。
* 在代码中，将边界外虚拟格子的高度固定为 `H` 即可。
* **诺伊曼边界条件**：固定边界处高度的导数为零（即边界两侧高度相等）。这用于模拟无法穿越的墙壁（如水池边）。
* 在代码中，当计算边界格子时，忽略其越界一侧的邻居（不进行水流交换）即可实现。
以下是诺伊曼边界条件的简化代码逻辑示意：
```cpp
// 假设 h_new 已初始化为 2*h_curr - h_old
for each grid i {
if (i-1 是有效格子) {
h_new[i] += beta * alpha * (h_curr[i-1] - h_curr[i]);
}
if (i+1 是有效格子) {
h_new[i] += beta * alpha * (h_curr[i+1] - h_curr[i]);
}
}
// 更新状态：h_old = h_curr; h_curr = h_new;
流体与刚体的耦合 ⚙️
在作业中，我们需要模拟方块在水面移动并激起波浪的效果。这涉及到流体-刚体耦合。耦合是双向的：
- 刚体对流体：方块排开水体，从而扰动水面。
- 流体对刚体：被排开的水体产生浮力，作用于方块。
我们重点解决第一个问题：如何模拟方块排水？
一个巧妙的方法是引入虚拟高度 v。假设在方块占据的格子上，我们不是直接移除水，而是临时给这些格子一个额外的虚拟高度。然后，通过求解一个线性方程组，计算出需要多少虚拟高度 v，才能使得在下一个模拟步长后，这些格子达到我们预期的“被排空”后的目标高度。
求解这个方程组可以使用共轭梯度法等数值方法。在作业框架中，已经提供了相关的求解器（PCG）。你需要做的是：
- 设置好需要求解的格子（
mask）。 - 根据方块的位置和目标排水量，计算方程组的右侧项（
b）。 - 调用求解器得到虚拟高度
v。 - 将
v乘以一个衰减系数（用于稳定模拟，避免因拖动过快产生过大波浪），然后加入到高度场的更新计算中。
对于第二个问题（浮力），可以根据阿基米德原理计算：每个被方块覆盖的格子所产生的浮力，等于其排开的水的重量（F = ρ * g * 被排开体积）。将所有格子的浮力向量相加，并计算它们对刚体质心产生的力矩，即可得到作用在方块上的总浮力和扭矩，进而影响其刚体运动。
总结
本节课中，我们一起学习了基于物理的水面波浪模拟。
我们从拉格朗日与欧拉两种模拟思路出发，引入了用于描述水面的高度场概念。通过物理推导和大量简化，得到了核心的浅水波方程。为了在计算机中求解，我们利用中心差分法对方程进行离散化，得到了直观的格子更新公式。为了保证模拟的合理性，我们探讨了体积守恒的方法、添加了阻尼效果，并介绍了两种边界条件的处理方式。
最后，我们深入探讨了本次作业的核心：流体-刚体耦合。通过引入虚拟高度并求解线性方程组，来模拟刚体（方块）排水并激发波浪的过程，同时也简要说明了如何计算水流对刚体产生的浮力。
这套基于浅水波模型的模拟方法，效率高、实现相对简单，是游戏中实现实时水面效果的常用技术。希望大家通过实验，能更深入地理解和掌握这些概念。
GAMES103-基于物理的计算机动画入门—P11-Lecture-11-Incompressible-Fluid-Dynamics-and-Eulerian-Fluids—GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将要学习基于网格（欧拉方法）的流体模拟技术。我们将首先介绍网格的表达方式，然后利用有限差分法计算微分算子，接着讨论不可压缩粘性流体的纳维-斯托克斯方程及其数值解法，最后探讨如何描述流体状态（如烟和水）。
1. 网格表达与有限差分法 🕸️
上一节我们介绍了基于粒子的模拟方法，本节中我们来看看基于固定网格的欧拉方法。这种方法将空间划分为规则的网格，并将物理量（如密度、速度）定义在网格上。
1.1 规则网格与物理场
规则网格（Regular Grid）将空间均匀划分为正方形（二维）或立方体（三维）的格子。每个格子中心可以存储一个物理量，例如标量（密度、压强、温度）或矢量（速度）。整个网格构成了一个物理场（标量场或矢量场）。
这种规则结构带来了一个巨大优势：计算导数变得非常容易。
1.2 利用中心差分法计算导数
有限差分法是数值计算导数的核心工具。对于定义在网格中心的值，我们可以使用中心差分法。
- 一阶导数：对于函数
f在i位置（假设网格间距为h）的一阶导数，公式为：(f(i+1) - f(i-1)) / (2h) - 二阶导数：二阶导数可以通过连续两次应用一阶导数公式得到，最终形式为：
(f(i-1) + f(i+1) - 2*f(i)) / (h^2)
在二维网格中，一个点 (i, j) 的二阶导数计算如下：
- x方向二阶导数：
(f(i-1, j) + f(i+1, j) - 2*f(i, j)) / (h^2) - y方向二阶导数：
(f(i, j-1) + f(i, j+1) - 2*f(i, j)) / (h^2)
1.3 拉普拉斯算子
拉普拉斯算子是梯度的散度，在模拟中极为重要。在二维网格上，离散化的拉普拉斯算子 ∇²f 计算公式为：∇²f(i, j) = (f(i-1, j) + f(i+1, j) + f(i, j-1) + f(i, j+1) - 4*f(i, j)) / (h^2)
直观理解是：一个点的拉普拉斯值等于其所有邻居值的和减去四倍自身值，再除以网格间距的平方。
1.4 边界条件
与求解其他偏微分方程一样，流体模拟也需要定义边界条件，主要有两种：
- 狄利克雷边界条件：指定边界上的函数值为已知常数。
- 诺伊曼边界条件：指定边界上函数导数的值，例如规定边界值与内部相邻值相等。
需要注意的是，在求解某些方程（如拉普拉斯方程 ∇²f = 0）时，不能全部使用诺伊曼边界条件，否则会导致系统矩阵奇异，解不唯一。
2. 交错网格与散度 🌀
上一节我们介绍了将物理量定义在格子中心的常规网格。但在流体模拟中，对速度场采用一种特殊的网格——交错网格（Staggered Grid）会更加方便。
2.1 交错网格的定义
在交错网格中，速度矢量并不定义在格子中心，而是定义在格子的面上：
- x方向速度
u：定义在垂直网格面的中心。 - y方向速度
v：定义在水平网格面的中心。
这样定义非常直观：速度值直接代表了流体通过该网格面的流量。
2.2 散度与不可压缩条件
对于一个格子，单位时间内流体的净流出量可以通过其四个面上的速度计算：净流出量 = u(i+1, j) + v(i, j+1) - u(i, j) - v(i, j)
在物理学中，不可压缩流体意味着每个格子内流体的净流入/流出量为零，即体积不变。这等价于速度场的散度为零：∇ · U = ∂u/∂x + ∂v/∂y = 0
在交错网格上，这个条件的离散形式恰好就是上面净流出量为零的公式。因此，使用交错网格能非常自然且精确地实施不可压缩约束。
2.3 双线性插值
由于物理量定义在离散的网格点上，当需要获取任意位置（非网格点）的值时，就需要进行插值。对于定义在格子中心的量（如压强），使用标准的双线性插值。对于定义在面上的速度，插值前需要对坐标进行0.5个网格的偏移，以对齐到正确的存储位置。
3. 纳维-斯托克斯方程与解法 ⚙️
本节中我们来看看描述流体运动的核心方程——纳维-斯托克斯方程，以及如何数值求解它。
3.1 方程概述
对于不可压缩粘性流体，纳维-斯托克斯方程描述了速度场的演化，包含两个部分：
- 动量方程：描述速度如何受外力、粘性、对流和压强梯度影响。
- 不可压缩条件：
∇ · U = 0，确保流体体积不变。
图形学中通常求解的是不可压缩形式。
3.2 分裂法求解策略
直接求解完整的纳维-斯托克斯方程很复杂。我们采用分裂法，将一个时间步内的更新分解为几个连续的、更简单的步骤，逐步更新速度场 U：
- 外力项：添加重力等外力。
U1 = U0 + Δt * g - 对流项：处理流体自身流动导致的速度迁移。
- 粘性项：处理流体的扩散（粘性）效应。
- 投影步：调整速度场，使其满足不可压缩条件。
接下来，我们详细看看每一步如何实现。
3.2.1 对流项：半拉格朗日法
对流项是欧拉方法中的难点。由于网格固定，而流体在流动，我们需要计算流体微团从上一时刻移动到当前位置所携带的速度。半拉格朗日法是解决此问题的稳定方案。
其核心思想是反向追踪：要得到当前网格点 X 在新时刻的速度，就沿着当前速度场反向追溯 Δt 时间，找到上一时刻该流体微团的位置 X_prev，然后将 X_prev 处的速度（通过插值得到）作为 X 的新速度。U_new(X) = U_old(X_prev)，其中 X_prev = X - Δt * U_old(X)
这种方法稳定，但可能引入数值耗散（模糊）。
3.2.2 粘性项：扩散求解
粘性项的形式是 ν∇²U，其中 ν 是粘性系数。这本质上是一个扩散过程，可以用显式方法求解：U_new = U_old + Δt * ν * ∇²U_old
其中拉普拉斯算子 ∇² 的计算方法已在第1.3节介绍。若 Δt 较大，此步可能不稳定，可采用小步长多次迭代或隐式方法提高稳定性。
3.2.3 投影步：压强求解与速度修正
这是确保不可压缩条件的关键一步。我们引入压强 p，并通过压强梯度来修正速度：U_final = U_before_projection - Δt * ∇p
修正后的速度必须满足散度为零 ∇ · U_final = 0。将上式代入，可得到一个关于压强 p 的泊松方程：∇²p = (∇ · U_before_projection) / Δt
求解这个泊松方程得到压强场 p，再用 p 去修正速度，最终得到既满足动量守恒又不可压缩的速度场。这个步骤被称为投影，因为它将速度场投影到散度为零的子空间。
4. 流体状态模拟：烟与水的表达 🌫️💧
上一节我们解决了速度场的更新问题，本节中我们来看看如何更新流体的可见状态，例如烟的密度和水的形状。
4.1 烟的模拟
烟的模拟相对直接。除了速度场，我们还需要模拟一个标量场，例如密度 ρ 或温度 T。这些被动标量场的更新也遵循类似的物理过程：
- 对流：使用半拉格朗日法，让密度随速度场移动。
- 扩散：密度自身也会扩散。
- 源：在烟源位置添加密度值。
烟的渲染通常使用体渲染技术，将密度场转换为视觉效果。
4.2 水的模拟与界面捕捉
水的模拟更复杂，因为水有明确的自由表面（界面）。我们需要额外的方法来刻画这个界面随速度场的演化。两种主流方法是：
- 体积分数法：每个网格存储一个值，表示该网格被水占据的体积百分比。这种方法简单，但界面模糊、不精确。
- 水平集法：每个网格存储一个有符号距离函数，其绝对值表示到水面的最短距离，符号表示内外（正为空气，负为水）。界面就是零等值面。水平集函数本身也通过对流方程更新。
水平集法能提供更清晰的界面，但面临一个挑战：体积损失。在数值模拟中，水的总体积可能无法严格保持，导致水面下降或上升。需要额外的算法（如粒子水平集法）来进行体积修正。
将模拟得到的水平集场（或密度场）转换为可用于渲染的几何网格，通常使用移动立方体算法。
总结 📚
本节课中我们一起学习了基于欧拉网格的流体模拟方法。
- 网格与微分：我们使用规则网格离散空间，并利用中心差分法方便地计算导数、拉普拉斯算子等。
- 交错网格：将速度定义在网格面上的交错网格，能自然地表达流体的通量，并简化不可压缩条件的实施。
- 纳维-斯托克斯方程：我们学习了描述不可压缩粘性流体运动的核心方程组。
- 分裂解法：通过将方程分裂为外力、对流、粘性、投影四个步骤，逐步更新速度场。其中，半拉格朗日法处理对流，投影步通过求解压强泊松方程来保证流体不可压缩。
- 状态模拟：烟的模拟通过更新密度场实现；水的模拟则需要使用水平集法等技术来捕捉动态变化的界面。
这种方法在图形学中已成为模拟水和气体的标准方案，被广泛应用于电影特效和游戏开发中。下节课，我们将探讨基于拉格朗日粒子法的流体模拟。
GAMES103-基于物理的计算机动画入门—P12-Lecture-12-SPH-and-Position-Based-Fluids—GAMES-Webinar—BV12Q4y1S73g_note
概述
在本节课中，我们将学习基于粒子的流体模拟方法，特别是光滑粒子流体动力学（SPH）的基本原理。我们将了解如何用粒子表示流体，并通过光滑插值模型来计算密度、压强和粘滞力，从而模拟流体的运动。
课程背景与安排
今天是本课程的最后一节课。课程作业的最终评分、颁奖和奖品发放将安排在年后进行。如果同学们想在过年期间补交作业，仍然可以提交，但可能无法参与评奖。完成所有作业的同学将获得一份电子证书作为完成课程的证明。
模拟领域涉及面非常广，包括流体、弹性体、碰撞处理等。很少有研究组能覆盖所有方向，大家通常专注于某一特定领域。本节课将讨论流体模拟的最后一部分内容，即基于粒子的方法。
粒子模拟方法概述
上两节课我们讨论了欧拉方法，它将空间划分为固定网格，通过改变网格内的物理量来模拟动画。今天我们将讨论拉格朗日视角的模拟方法，即使用运动的粒子来进行物理模拟。
基于粒子的流体模拟有多种变形：
- SPH（光滑粒子流体动力学）：最传统的方法。
- PBF（基于位置的流体）：使用基于位置的约束进行模拟。
- Peridynamics：将SPH与弹性体模拟结合，便于模拟破碎效果。
- MPM（物质点法） 或 PIC（粒子网格法）：混合粒子与网格的方法。
不同的方法适用于不同的效果，例如MPM常用于模拟血或沙子，Peridynamics用于形变和破碎。本节课我们将从最传统的SPH方法入手。
光滑粒子（SPH）模型
核心思路
我们用大量粒子来表达流体，每个粒子附带着物理变量（如位置、速度、质量）。当粒子运动时，这些变量也随之运动，这是典型的拉格朗日视角。
在图形学中，由于GPU通常无法直接渲染粒子，我们需要将粒子转化为三角网格进行离线渲染，或在游戏中用带纹理的球体或方块来近似表示。
从简单平均到光滑插值
假设场景中有许多粒子，每个粒子都有一个物理量 \( A \)（标量或矢量）。我们想在空间中某个位置（例如图中红点）估算该物理量的值。这本质上是一个插值问题。
1. 简单平均模型
最简单的想法是在一定半径内，对所有粒子的值求平均：
\[ A_{i} = \frac{1}{n} \sum_{j} A_{j} \]
其中 \( n \) 是半径内的粒子数。但这种方法没有考虑粒子的空间分布。
2. 考虑体积权重的模型
我们为每个粒子引入一个体积 \( V \) 作为权重：
\[ A_{i} = \sum_{j} A_{j} V_{j} \]
这里假设了单位球的总体积为1。这个模型考虑了粒子所占的空间，但还不够“光滑”。
3. 引入光滑核函数
我们希望插值结果能随位置连续、平滑地变化，而不是在粒子进出插值半径时剧烈跳变。因此，我们引入一个与距离相关的权重函数——光滑核函数 \( W \)。
\[ A_{i} = \sum_{j} A_{j} V_{j} W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
其中 \( h \) 是光滑长度。核函数满足：距离近则权重大，距离远则权重小。
密度、体积与最终插值公式
然而，粒子的体积 \( V \) 并非常数，它会随着粒子的分布疏密而变化。我们需要动态计算它。
我们假设每个粒子有恒定的质量 \( m \)。这里定义的密度 \( \rho \) 是粒子分布的密度，描述粒子在空间中的拥挤程度，而非水的物理密度。
根据光滑插值公式，粒子 \( i \) 的密度可以通过其邻居估算：
\[ \rho_i = \sum_{j} m_j W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
直观理解：周围邻居的质量和越大，该处就越拥挤，密度越高。
得到密度后，我们可以计算体积：
\[ V_i = \frac{m_i}{\rho_i} \]
将体积公式代回最初的插值公式，就得到了SPH的最终插值形式：
\[ A_{i} = \sum_{j} A_{j} \frac{m_j}{\rho_j} W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
在实际计算中，通常分两步：先估算每个粒子的密度和体积，再利用此公式插值其他物理量。
微分算子的计算
使用光滑插值模型的一个巨大优势是，可以方便地计算物理量的梯度（∇）和拉普拉斯算子（∇²）。
在计算梯度时，我们近似认为邻居粒子的物理量 \( A_j \) 和体积 \( V_j \) 是常数（因为粒子 \( i \) 的微小运动对邻居影响较小）。这样，梯度算子就只作用于已知的核函数 \( W \)：
\[ \nabla A_{i} \approx \sum_{j} A_{j} V_{j} \nabla W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
由于核函数是解析定义的，其梯度有现成公式。拉普拉斯算子的计算同理：
\[ \nabla^2 A_{i} \approx \sum_{j} A_{j} V_{j} \nabla^2 W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
核函数示例
常见的核函数使用多项式形式，因其计算和求导简单。例如下面这个函数：
首先计算归一化距离 \( q = \frac{||\mathbf{x}_i - \mathbf{x}_j||}{h} \)，然后核函数定义为：
W(q, h) = \frac{315}{64\pi h^3}
\begin{cases}
(1 - q^2)^3 & 0 \leq q \leq 1 \\
0 & q > 1
\end{cases}
其梯度（一阶导）和拉普拉斯算子（二阶导）也有对应的多项式表达式。在实际代码实现中，直接套用这些公式即可。
SPH流体模拟算法
上一节我们介绍了SPH的核心插值模型，本节我们来看看如何将其应用于流体模拟。其思路与之前讲过的粒子模拟算法相似，主要考虑三种力：重力、压强力和粘滞力。
1. 重力
重力计算简单，直接作为外力施加给每个粒子：
\[ \mathbf{F}_i^{gravity} = m_i \mathbf{g} \]
2. 压强力
压强力源于压力的不平衡（压力差），而非压力本身。数学上，压强力与压强的负梯度相关：
\[ \mathbf{F}_i^{pressure} = -V_i \nabla p \]
我们假设压强场也是光滑的，利用SPH梯度公式计算：
\[ \mathbf{F}_i^{pressure} = -\sum_{j} V_j p_j \nabla W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
其中，粒子压强 \( p \) 通常由密度通过经验公式计算，例如：
\[ p_i = k ( \rho_i - \rho_0 )^7 \]
这里 \( k \) 和 \( \rho_0 \) 是常数。
3. 粘滞力
粘滞力的效果是使流体粒子的速度趋于一致，数学上可用速度场的拉普拉斯算子来描述：
\[ \mathbf{F}_i^{viscosity} = \mu m_i \nabla^2 \mathbf{v}_i \]
其中 \( \mu \) 是粘滞系数。利用SPH拉普拉斯公式计算：
\[ \mathbf{F}_i^{viscosity} = \mu \sum_{j} m_j (\mathbf{v}_j - \mathbf{v}_i) \nabla^2 W(||\mathbf{x}_i - \mathbf{x}_j||, h) \]
注意速度是矢量，计算时需对每个分量（x, y, z）分别进行。
模拟流程
以下是SPH流体模拟的基本算法步骤：
- 寻找邻居：对于每个粒子，找到其周围一定半径内的所有邻居粒子。
- 计算密度：根据邻居粒子质量，用SPH公式估算每个粒子的密度 \( \rho \)。
- 计算压强：利用密度，通过经验公式计算每个粒子的压强 \( p \)。
- 计算合力：
- 加上重力 \( \mathbf{F}^{gravity} \)。
- 加上利用SPH梯度公式计算的压强力 \( \mathbf{F}^{pressure} \)。
- 加上利用SPH拉普拉斯公式计算的粘滞力 \( \mathbf{F}^{viscosity} \)。
- 更新状态：根据合力 \( \mathbf{F}_{total} \) 更新粒子速度，再根据速度更新粒子位置。
挑战与扩展
性能挑战与优化
模拟逼真的流体需要百万甚至千万级的粒子。计算每个粒子的邻居是主要性能瓶颈，穷举法不可行。常用优化方法包括：
- 空间划分：将空间划分为均匀网格（Spatial Hashing），每个粒子存入对应网格，快速定位邻居。
- 层次结构：对于非均匀分布的粒子，可使用八叉树（Octree）等数据结构。
- 自适应粒子：在需要高细节的区域（如水面、泡沫）使用更小、更密的粒子，在内部区域使用更大、更疏的粒子以提升效率。
表面重建
粒子系统渲染前需要将点云转换为三角网格表面（三维重建）。简单方法是将每个粒子视为球体，计算其符号距离函数（SDF），然后进行等值面提取（如Marching Cubes）。通常还需对生成的网格进行平滑处理以消除噪点，同时保留细节。
研究前沿
SPH及其变体仍在不断发展，研究热点包括：
- 效率与实时性：如何在资源受限（如游戏）环境下进行高效模拟。
- 不可压缩性：精确保持流体体积不变是一个难点，MPM等方法试图更好地解决此问题。
- 边界处理：精确处理流体与固体、空气的交互边界。
- AI驱动模拟：利用机器学习加速或控制模拟过程，但目前尚难以完全替代物理引擎。
总结与展望
本节课我们一起学习了基于粒子的流体模拟方法SPH。我们从光滑插值模型出发，讲解了如何通过粒子估算密度、体积及其他物理量，并利用核函数方便地计算梯度与拉普拉斯算子，从而求解流体的压强力和粘滞力，最终实现流体运动的模拟。
图形学中的物理模拟，其核心价值在于在效果逼真与计算效率之间取得平衡，以满足实时应用（如游戏、虚拟数字人）的需求。这与追求绝对精确的科学计算模拟有显著区别。
对于希望深入该领域的同学，建议：
- 夯实基础：理解本节课及本课程系列的内容，是阅读相关论文的坚实基础。
- 专注方向：模拟领域分支众多，建议选择一个感兴趣的方向（如流体、弹性体、碰撞）深入钻研，而非泛泛而读。
- 动手实践：在阅读论文的同时，尝试实现其中的算法，理解会更深刻。
- 关注现实问题：了解工业界（游戏、影视、数字服装等）的真实需求，思考如何用技术解决实际问题，这往往能产生更有价值的研究。
物理模拟是计算机图形学中富有挑战且极具趣味的方向，希望本课程能为大家打开一扇门。祝大家新年快乐，在图形学的道路上不断进步！
GAMES103-基于物理的计算机动画入门—P2-Lecture-02-Math-Background–Vector–Matrix-and-Tensor-Calculus—GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将学习物理模拟中至关重要的数学基础。主要内容分为三部分：矢量、矩阵以及相关的微积分概念。这些知识是理解后续物理模拟算法的基石。
矢量（Vector）📐
上一节我们介绍了课程概述，本节中我们来看看矢量。矢量是描述方向和大小的一种数学对象，在二维或三维空间中尤为常见。
矢量的定义与表示
一个三维矢量 p 可以表示为 p = (p_x, p_y, p_z)，它属于实数空间 ℝ³。原点是一个特殊的矢量 0 = (0, 0, 0)。在印刷体和论文中，通常用黑体表示矢量，用斜体表示标量，以此区分，而非使用箭头符号。
坐标系：左手系与右手系
三维坐标系分为左手系和右手系。在图形学中，两者都有应用。例如，OpenGL通常使用右手系，而Unity和DirectX使用左手系。判断方法遵循“右手定则”或“左手定则”：从x轴转向y轴，大拇指指向的方向即为z轴正方向。
堆叠矢量（Stack Vector）
矢量不一定总具有直观的几何意义。例如，一个由11个顶点构成的物体，可以将所有顶点的坐标按顺序排列成一个33维的大矢量 X = [x₀, x₁, …, x₁₀]ᵀ。这个堆叠矢量用于描述物体的整体状态，其维度是顶点数与空间维度的乘积。
矢量的基本运算
以下是矢量的一些基本运算：
- 加减法：对应元素相加减。几何上，矢量加法遵循三角形法则或平行四边形法则。
- p ± q = (p_x ± q_x, p_y ± q_y, p_z ± q_z)
- 与标量的乘法：矢量与标量t相乘，表示其缩放。
- p(t) = p + tv （描述点沿方向 v 的运动）
- p(t) = (1-t)p + tq （描述点 p 和 q 之间的线性插值）
- 范数（Norm）：表示矢量的大小。
- 2-范数（欧几里得范数）：‖p‖ = √(p_x² + p_y² + p_z²)
- 1-范数（曼哈顿范数）：‖p‖₁ = |p_x| + |p_y| + |p_z|
- 无穷范数：‖p‖∞ = max(|p_x|, |p_y|, |p_z|)
- 单位矢量：范数为1的矢量。任何非零矢量 v 可通过 v / ‖v‖ 进行归一化。
矢量的乘法
我们有两种重要的矢量乘法。
点乘（内积）
点乘的结果是一个标量。
- 定义：p · q = p_x q_x + p_y q_y + p_z q_z
- 几何意义：p · q = ‖p‖ ‖q‖ cos θ，其中θ是两矢量的夹角。
- 性质：
- 满足交换律和分配律。
- p · p = ‖p‖²。
- 若 p · q = 0 且 p, q 非零，则 p 与 q 垂直。
叉乘（外积）
叉乘的结果是一个新的矢量。
- 定义：r = p × q = (p_y q_z - p_z q_y, p_z q_x - p_x q_z, p_x q_y - p_y q_x)
- 几何意义：结果矢量 r 同时垂直于 p 和 q，其大小 ‖r‖ = ‖p‖ ‖q‖ sin θ。
- 性质：
- 不满足交换律，p × q = - (q × p)。
- 满足分配律。
- 若 p × q = 0 且 p, q 非零，则 p 与 q 平行。
点乘的应用实例
点乘在图形学中有广泛的应用。
点到直线的投影
给定点 q 和过点 o、方向为 v 的直线。点 q 在直线上的投影点 s 及其距离标量s可通过点乘求得：
- s = (q - o) · (v / ‖v‖)
- s = o + s v
定义平面与有符号距离
给定平面上一点 c 和其单位法向量 n，可以定义一个平面。对于空间中任意点 p，有符号距离 d = (p - c) · n。d的符号表明了 p 相对于平面的位置（上、面上、下），其绝对值是点到平面的真实距离。
点与球的碰撞检测
一个沿方向 v 运动的点 p(t) = p + tv 与静止的球（球心 c，半径 r）发生碰撞时，满足方程 ‖p(t) - c‖² = r²。将其展开并利用点乘性质，可得到一个关于时间t的一元二次方程，通过求解该方程即可判断碰撞是否发生及发生的时间。
叉乘的应用实例
叉乘同样在图形学中扮演关键角色。
三角形的面积与法向
给定三角形顶点 x₀, x₁, x₂，定义两条边向量 e₁₀ = x₁ - x₀, e₂₀ = x₂ - x₀。
- 法向量：n = (e₁₀ × e₂₀) / ‖e₁₀ × e₂₀‖
- 面积：A = ½ ‖e₁₀ × e₂₀‖
点与三角形的位置关系（内外测试）
要判断点 p 是否在三角形 x₀, x₁, x₂ 内部，可以检查由 p 与三角形每条边构成的小三角形的法向是否与大三角形法向 n 同向。具体通过计算三个标量值来判断：
- a₀ = n · ( (x₁ - p) × (x₂ - **p`) ) / 2
- a₁ = n · ( (x₂ - p) × (x₀ - **p`) ) / 2
- a₂ = n · ( (x₀ - p) × (x₁ - **p`) ) / 2
若 a₀, a₁, a₂ 同号（通常为正），则 p 在三角形内部。这些a值实际上是带符号的面积。
重心坐标（Barycentric Coordinates）
将上述符号面积 a₀, a₁, a₂ 除以三角形总面积 A，可得到重心坐标权重 (b₀, b₁, b₂)，满足 b₀ + b₁ + b₂ = 1。三角形内任意点 p 的位置可由其顶点插值得到：p = b₀ x₀ + b₁ x₁ + b₂ x₂。这在渲染中进行颜色插值（如Gouraud着色）时非常有用。
四面体的体积
给定四面体四个顶点 x₀, x₁, x₂, x₃，其（带符号的）体积 V 为：
- V = (1/6) ( (x₁ - x₀) × (x₂ - **x₀`) ) · (x₃ - x₀)
体积的正负由顶点顺序决定，为零则表示四点共面。类似于三角形，也可以定义四面体的重心坐标。
点与三角形的碰撞检测（进阶）
结合点面相交和内外测试，可以检测运动点与三角形是否碰撞。首先，求解点 p(t) 与三角形所在平面的交点时间 t（令四面体体积公式为零）。然后，将交点代入三角形内外测试公式，检查其是否在三角形内部。
矩阵（Matrix）🔲
上一节我们深入探讨了矢量及其运算，本节中我们来看看矩阵。矩阵可以看作是一组矢量的有序排列，是描述线性变换的强大工具。
矩阵的基本定义与运算
一个3x3的实数矩阵 A 可以表示为：
A = [ a₀₀ a₀₁ a₀₂ ]
[ a₁₀ a₁₁ a₁₂ ]
[ a₂₀ a₂₁ a₂₂ ]
以下是矩阵的一些基本概念和运算：
- 转置（Transpose）：交换矩阵的行和列，记作 Aᵀ。
- 对称矩阵：满足 A = Aᵀ 的矩阵。
- 对角矩阵：非对角线元素全为零的矩阵。
- 单位矩阵：对角线元素全为1的对角矩阵，记作 I。任何矩阵或矢量与 I 相乘保持不变。
- 矩阵乘法：不满足交换律（AB ≠ BA），但满足结合律。
- 逆矩阵：若存在矩阵 B 使得 AB = BA = I，则 B 是 A 的逆矩阵，记作 A⁻¹。并非所有矩阵都可逆。
图形学中的特殊矩阵
在图形学中，以下几类矩阵尤为重要。
正交矩阵（Orthogonal Matrix）
一个矩阵 A 是正交矩阵，如果其列向量（或行向量）是两两正交的单位向量。正交矩阵有一个非常好的性质：Aᵀ = A⁻¹。在图形学中，旋转矩阵就是正交矩阵。它将物体的局部坐标系（x, y, z 轴）旋转到世界坐标系中的新方向（u, v, w），即 A = [u v w]。
缩放矩阵（Scaling Matrix）
缩放矩阵是一个对角矩阵 D，其对角线元素 (d_x, d_y, d_z) 分别表示在x, y, z轴方向上的缩放因子。
矩阵分解
矩阵分解帮助我们理解线性变换的本质。
奇异值分解（Singular Value Decomposition, SVD）
任何矩阵 A 都可以分解为 A = U D Vᵀ。其中 U 和 V 是正交矩阵，D 是对角矩阵（元素称为奇异值）。SVD的几何解释：任何线性变换都可以通过三个步骤实现：1. 旋转（Vᵀ）；2. 沿坐标轴缩放（D）；3. 再次旋转（U）。
特征值分解（Eigenvalue Decomposition）
对于对称矩阵 A，可以分解为 A = U D Uᵀ。其中 U 是正交矩阵，其列向量为特征向量；D 是对角矩阵，对角线元素为特征值。这可以看作是SVD在对称情况下的特例。
对称正定矩阵（Symmetric Positive Definite, SPD）
对称正定矩阵在求解物理模拟中的线性系统时至关重要。
- 定义：一个对称矩阵 A 是正定的，如果对于任何非零矢量 v，都有 vᵀA v > 0。若 vᵀA v ≥ 0，则为半正定。
- 直观理解：正定矩阵可以类比为正实数。一个对角元素全为正的对角矩阵显然是正定的。通过特征值分解可知，一个对称矩阵是正定的，当且仅当其所有特征值均为正。
- 性质与判定：
- 正定矩阵必然可逆。
- 若矩阵对角占优（即每一行/列上，对角线元素的绝对值大于该行/列其他元素绝对值之和），则该矩阵是正定的（这是一个充分条件，非必要条件）。
总结 📝
本节课我们一起学习了物理模拟所需的数学基础。我们从矢量的定义、运算（点乘、叉乘）及其在几何计算（投影、距离、碰撞检测）中的应用开始。然后，我们探讨了矩阵，包括其基本运算、图形学中特殊的旋转与缩放矩阵，以及重要的矩阵分解（SVD和特征值分解）。最后，我们介绍了在后续求解线性系统中非常关键的对称正定矩阵的概念。
这些数学工具是构建物理模拟算法的语言。理解它们将帮助我们更好地学习后续关于刚体模拟、有限元方法等内容。下节课我们将继续讲解张量微积分和线性系统求解，并逐渐过渡到具体的物理模拟实践。
GAMES103-基于物理的计算机动画入门—P3-Lecture-03-Rigid-Body-Dynamics—GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将学习刚体动力学的基础知识。课程内容分为两部分：首先回顾并完成上周关于线性代数和微积分的数学基础，然后深入探讨单个刚体的运动模拟，包括平移和旋转。
数学基础回顾与补充 📐
上一节我们介绍了矩阵正定性等概念，本节中我们来看看线性系统和微积分在物理模拟中的应用。
矩阵正定性示例
以下是一个证明示例：若矩阵 A 对称正定，则矩阵 B 也半正定。
B 的形式为：
B = [ A, -A;
-A, A ]
证明过程如下：
- 设任意向量可拆分为两部分：
[x; y]。 - 根据半正定定义，需证明
[x^T, y^T] * B * [x; y] >= 0。 - 展开计算可得：
(x - y)^T * A * (x - y)。 - 由于 A 正定，上式恒大于等于零（当
x = y时等于零）。 - 因此，B 半正定。
此结论在模拟中很有用，因为许多系统矩阵具有类似结构。
线性系统求解
许多数学问题最终归结为求解线性系统 Ax = b，其中 A 是矩阵，b 是边界条件向量，x 是未知向量。
直接计算 A 的逆矩阵通常不可行，因为计算量大且逆矩阵可能失去稀疏性。主要有两类方法：
1. 直接法（如LU分解）
直接法将矩阵 A 分解为下三角矩阵 L 和上三角矩阵 U 的乘积：A = LU。
求解分两步：
- 解
Ly = b（前向代入）。 - 解
Ux = y（后向代入）。
以下是前向代入的伪代码示例：
# 假设 L 是下三角矩阵，b 是已知向量
y = zero_vector(n)
for i in range(n):
sum = 0
for j in range(i):
sum += L[i][j] * y[j]
y[i] = (b[i] - sum) / L[i][i]
直接法特点：
- 矩阵稀疏性影响分解结果。
- 计算分分解和求解两步，若 A 固定可复用分解结果。
- 并行化相对困难。
- 常用库：Intel MKL Pardiso。
2. 迭代法
迭代法通过不断更新 x 来逼近解，基本形式为：
x_{k+1} = x_k + α * M^{-1} * (b - A * x_k)
其中 α 是松弛系数，M 是易于求逆的矩阵（如 A 的对角部分），(b - A*x_k) 是残差。
迭代法收敛条件与矩阵的谱半径有关。常用方法包括雅可比法（取 M 为 A 的对角线）和高斯-赛德尔法（取 M 为 A 的下三角部分）。
迭代法与直接法对比：
- 优点：易于实现、易于并行、适合不求精确解的场景。
- 缺点：存在收敛性问题、求精确解可能较慢。
微积分基础
本节回顾向量微积分，为后续内容做准备。
梯度
对于标量函数 f(x)，其梯度 ∇f 是一个向量，指向函数值增长最快的方向。
∇f = [ ∂f/∂x, ∂f/∂y, ∂f/∂z ]^T
雅可比矩阵、散度与旋度
对于向量函数 F(x)，其雅可比矩阵 J 包含所有一阶偏导。
散度 div(F) 是雅可比矩阵的迹（对角线之和）。旋度 curl(F) 描述了场的旋转特性，在流体模拟中用于漩涡计算。
泰勒展开
函数 f(x) 在 x0 处的泰勒展开为：
f(x) ≈ f(x0) + ∇f(x0)^T (x - x0) + 1/2 (x - x0)^T H(x0) (x - x0) + ...
其中 H 是海森矩阵（二阶偏导矩阵）。矩阵的正定性与海森矩阵相关，决定了函数在该点的凹凸性。
实例：弹簧模型
考虑一维弹簧，原长 L，弹性系数 k。当前位置为 x，则弹簧能量 E 和力 F 为：
E = 1/2 * k * (|x| - L)^2
F = -∇E = -k * (|x| - L) * (x / |x|)
力的方向沿弹簧轴向，大小与形变成正比。
进一步求力关于位置的导数，可得到刚度矩阵（海森矩阵）：
∂F/∂x = k * [ (x*x^T)/|x|^2 + ( (|x|-L)/|x| ) * (I - (x*x^T)/|x|^2 ) ]
对于多顶点系统，所有变量和力可组合成大向量和大矩阵，其刚度矩阵具有分块形式，与我们开头证明的矩阵形式相似。
刚体动力学 🤖
上一节我们完成了数学基础的铺垫，本节中我们来看看刚体运动的核心原理。
刚体状态描述
刚体运动包括平移和旋转，无形状变化。其状态由两部分描述：
- 位置
x：一个三维向量，描述物体中心的位置。 - 旋转
R：描述物体的朝向。可用旋转矩阵、欧拉角或四元数表示。
在Unity中，物体的 Transform 组件包含了 position 和 rotation 属性，分别对应位置和旋转（内部用四元数存储）。
运动方程与时间积分
刚体运动遵循牛顿定律。我们有两个核心变量：
- 速度
v = dx/dt - 角速度
ω（描述旋转快慢和轴）
运动方程（积分形式）为：
v_new = v_old + (1/m) * ∫ F dt
x_new = x_old + ∫ v dt
其中 F 是合力，m 是质量。
数值积分方法
求解积分需用数值方法。以一维速度积分 ∫ v dt 为例，即估算速度曲线下的面积。
-
显式欧拉法：用起始时刻速度估算。
∫ v dt ≈ v(t0) * Δt一阶精度，可能低估面积。
-
隐式欧拉法：用结束时刻速度估算。
∫ v dt ≈ v(t1) * Δt一阶精度，可能高估面积。
-
中点法：用中间时刻速度估算。
∫ v dt ≈ v(t0 + Δt/2) * Δt二阶精度，更精确。
对于刚体模拟，常采用一种“蛙跳”格式，交替更新速度和位置，本质上等价于对两者都使用中点法。
平移运动模拟
模拟平移运动的更新步骤如下：
- 计算所有顶点上的合力
F_total。 - 更新速度：
v_new = v_old + (F_total / m) * Δt。 - 更新位置：
x_new = x_old + v_new * Δt。
常见的力包括：
- 重力：
F_gravity = m * g，方向竖直向下。 - 空气阻力：可简化为对速度的衰减，如
v_new = 0.99 * v_old，简单且稳定。
旋转运动模拟
旋转的模拟更为复杂。
旋转的表示
我们推荐使用四元数 q 表示旋转。它由实部 s 和虚部向量 v 构成：q = [s, v]。一个绕单位轴 n 旋转 θ 角度的四元数为：
q = [cos(θ/2), sin(θ/2) * n]
Unity 内部即使用四元数，可通过 Transform.rotation 访问。
旋转运动方程
旋转的模拟需要以下物理量：
- 力矩
τ：力引起旋转的“推力”。对于顶点i，力矩τ_i = r_i × F_i，其中r_i是从质心到顶点的向量。总力矩为所有顶点力矩之和。 - 惯性张量
I：旋转中的“质量”，是一个3x3矩阵。它依赖于物体的形状和质量分布。惯性张量从物体参考坐标系变换到世界坐标系的公式为：I_world = R * I_ref * R^T。
旋转的更新步骤如下：
- 计算总力矩
τ_total。 - 更新角速度：
ω_new = ω_old + I^{-1} * τ_total * Δt。 - 更新旋转四元数：
q_new = q_old + (Δt/2) * [0, ω] * q_old。注意这里是四元数乘法。 - 对
q_new进行归一化，保持其为单位四元数。
刚体模拟器框架
一个完整的刚体模拟器在每个时间步需执行以下操作：
平移部分：
- 计算合力
F。 v = v + (F / m) * Δtx = x + v * Δt
旋转部分：
- 计算总力矩
τ。 ω = ω + I^{-1} * τ * Δtq = q + (Δt/2) * Quaternion(0, ω) * qq = q.Normalize()
在Unity中实现时需注意：
- 位置
x对应Transform.position。 - 旋转
q对应Transform.rotation。 - 速度
v和角速度ω需在自定义脚本中声明和更新。 - Unity 提供四元数乘法，但不直接提供四元数与标量乘法或加法，需手动对四个分量操作。
- Unity 的
Matrix4x4类可用来进行矩阵运算，但需注意其是4x4矩阵。
实现建议
- 分步实现：先实现平移运动，再实现旋转运动。测试旋转时，可先固定角速度
ω为常数，验证旋转更新正确后再加入力矩计算。 - 注意归一化：更新四元数后务必进行归一化，防止数值误差累积。
- 重力与旋转：均匀重力场中，重力不产生力矩，因此重力不会使物体自发旋转。
总结 🎯
本节课中我们一起学习了刚体动力学模拟的核心内容。
我们首先补充了相关的数学基础，包括通过示例理解了特定结构矩阵的正定性、线性系统求解的直接法与迭代法，以及向量微积分在物理公式推导中的应用。
接着，我们深入探讨了刚体模拟。刚体的状态由位置和旋转描述。我们学习了使用四元数表示旋转的优势和方法。运动模拟的核心是对运动方程进行时间积分，我们介绍了几种积分方法及其精度。对于平移运动，我们给出了基于力的更新步骤。对于更复杂的旋转运动，我们引入了力矩和惯性张量的概念，并给出了基于四元数和角速度的更新流程。
最终，我们整合出了一个刚体模拟器的基本框架。掌握这些原理，是编写一个简单刚体物理引擎的基础。在接下来的课程中，我们将探讨刚体之间的碰撞检测与响应问题。
GAMES103-基于物理的计算机动画入门—P4-Lecture-04-Rigid-Body-Contacts–Lab-1----GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将要学习刚体碰撞检测与响应的核心方法，并了解一个无需物理知识的替代方案——形状匹配。课程内容将围绕作业展开，首先介绍作业要求，然后回顾刚体模拟的基础概念，最后深入讲解碰撞处理与形状匹配技术。
作业介绍与回顾 📋
首先，我们交代本次作业。这是课程的第一次作业，目标是实现一个刚体碰撞模拟。作业分为基础任务和附加题。
基础任务是完成一个名为 RigidBody 的 Unity 脚本，实现基于冲量法的刚体碰撞。脚本中已提供计算惯性张量、叉乘矩阵等辅助功能。完成后的效果是：运行程序，按下 L 键，兔子模型会撞向墙壁并反弹落下。
附加题是使用“形状匹配”技术实现刚体模拟。这种方法不涉及任何物理公式，完全基于图形学和优化。在 Unity 中，通过切换脚本（从 RigidBody 切换到 ShapeMatching）即可看到不同的模拟效果。此方法在处理摩擦时可能会有轻微滑动。
作业周期为两周。对于不熟悉 Unity 的同学，建议先学习 Unity 的基本操作和脚本编写。
刚体运动回顾 🔄
上一节我们介绍了刚体模拟，但时间仓促。本节我们首先回顾刚体运动的核心算法，特别是针对没有大学物理基础的同学。
刚体运动分为线性（平移）运动和旋转运动。线性运动的更新基于牛顿第二定律：
位置更新公式：v_new = v_old + (F_total / m) * Δtx_new = x_old + v_new * Δt
其中，F_total 是所有顶点上力的矢量和，m 是刚体总质量。
对于旋转运动，我们需要对应的概念。造成旋转趋势的不是力，而是 力矩。
力矩计算公式：τ_i = (R * r_i) × f_i
这里，R 是描述当前刚体朝向的旋转矩阵，r_i 是顶点在参考姿态下相对于质心的向量，f_i 是作用在该顶点上的力。力矩 τ 的方向垂直于力臂和力的方向，大小与两者模长及夹角的正弦成正比。
将所有顶点的力矩求和，得到总力矩 τ_total。
接下来，我们需要更新角速度。在线性运动中，加速度是 F/m。在旋转中，角加速度是 τ 除以 惯性张量 I。但 I 不是一个标量，而是一个矩阵。
惯性张量计算公式（在参考坐标系下）：I_ref = Σ_i [ m_i * ( (r_i^T * r_i) * E - r_i * r_i^T ) ]
其中，E 是单位矩阵。这是因为物体对旋转的“抵抗”与旋转轴的方向有关。例如，质量分布远离旋转轴时，更难转动。
当刚体旋转后，其惯性张量会变化。但我们可以通过旋转矩阵快速计算当前的惯性张量，而无需重新求和：
当前惯性张量计算公式：I = R * I_ref * R^T
得到惯性张量后，即可更新角速度：
角速度更新公式：ω_new = ω_old + I^{-1} * τ_total * Δt
最后，使用新的角速度更新描述刚体朝向的四元数 q。其更新公式涉及四元数乘法，具体推导可参考课程附录。
以上就是刚体运动模拟的基本框架。接下来，我们将进入本节课的核心：碰撞处理。
单点碰撞检测与响应 🎯
本节中，我们来看看如何处理一个单独质点的碰撞，暂时不考虑刚体的旋转。我们主要介绍两种方法：惩罚法和冲量法。
首先，我们需要一个工具来判断碰撞是否发生：有符号距离函数。
有符号距离函数
定义一个函数 Φ(x)，对于空间中的任意点 x，它返回到某个物体表面的距离。符号表明点在表面的内侧（负值）还是外侧（正值）。表面就是 Φ(x) = 0 的点的集合。
以下是几个简单几何体的距离函数示例：
- 平面（平面上一点
p，法向n）：Φ(x) = (x - p) · n - 球体（球心
c，半径r）：Φ(x) = ||x - c|| - r - 圆柱体（轴上一点
p，轴向a，半径r）：d = || (x-p) - ((x-p)·a) * a ||Φ(x) = d - r
对于复杂物体，可以通过布尔操作组合简单物体的距离函数：
- 交集（物体A 与 物体B）：点在内部当所有
Φ_i(x) < 0。近似距离函数取max(Φ_i(x))。 - 并集（物体A 或 物体B）：点在内部当任一
Φ_i(x) < 0。在外侧时，距离函数可近似为min(Φ_i(x))。
利用距离函数，碰撞检测变得简单：若 Φ(x) < 0，则发生碰撞。
惩罚法
惩罚法的思路是：检测到穿透后，施加一个力将点推出去。这个力通常像弹簧力，大小与穿透深度成正比，方向是表面的法向（即距离函数的梯度方向）。
基本惩罚力公式：f = -k * Φ(x) * n
其中 n = ∇Φ(x)，k 是弹性系数。
为了避免明显的穿透，可以设置一个缓冲厚度 ε。当 Φ(x) < ε 时就施加力，力的大小与 (ε - Φ(x)) 成正比。
惩罚法的主要问题是参数 k 难以调节：太小则无法阻止穿透，太大则容易导致模拟不稳定（过冲）。一种改进是使用 对数障碍函数，让力随着点接近表面而急剧增大：
对数障碍函数力公式：f = (ρ / Φ(x)) * n
其中 ρ 是强度参数。但此法要求步长非常小，且绝对不能发生穿透。
冲量法
冲量法的思路是：在检测到碰撞的瞬间，直接更新点的速度和位置，使其立即响应，而不是等到下一时刻通过力来改变。
以下是处理一个质点碰撞的冲量法步骤：
-
位置修正：将穿透的点沿法向推回表面。
x_new = x_old - Φ(x) * n（因为Φ(x)为负） -
速度修正：修改点的速度，模拟反弹和摩擦。
- 将速度
v分解为法向分量v_n和切向分量v_t。v_n = (v·n) * nv_t = v - v_n - 计算新的法向速度（反弹）和切向速度（摩擦）。
v_n_new = -μ_n * v_n（μ_n为弹性系数，0 ≤ μ_n ≤ 1）v_t_new = a * v_t（a为摩擦衰减系数） - 根据库仑摩擦定律，
a应尽可能小，但不能使切向速度反向。其计算公式为：a = max(0, 1 - μ_t * (1 + μ_n) * |v_n| / |v_t|)
其中μ_t是摩擦系数。a=0表示静摩擦（物体停下），a>0表示动摩擦。 - 合成新速度：
v_new = v_n_new + v_t_new
- 将速度
冲量法能更精确地控制碰撞后的反弹和摩擦效果。接下来，我们将把冲量法应用到刚体上。
刚体的碰撞处理 🤖💥
对于刚体，我们不能直接修改单个顶点的速度，因为模拟的状态变量是质心的速度 v、角速度 ω、位置 x 和朝向 q。我们需要通过施加冲量 J 来间接改变顶点的速度。
假设在顶点 i 处检测到碰撞。该顶点的世界坐标和速度为：x_i = x + R * r_iv_i = v + ω × (R * r_i)
我们的目标是：找到一个冲量 J，施加在碰撞点，使得该点碰撞后的速度 v_i_new 等于用前述单点冲量法计算出的理想速度 v_i_desired。
推导表明，冲量 J 与顶点速度变化量 Δv_i = v_i_desired - v_i 存在线性关系：
冲量与速度变化的关系：Δv_i = K * J
其中矩阵 K = (1/m) * E - ( [R*r_i]× )^T * I^{-1} * [R*r_i]×，[·]× 表示叉乘矩阵。
因此，我们可以计算出所需的冲量 J = K^{-1} * Δv_i。得到 J 后，再更新刚体的整体状态：
刚体状态更新公式：v_new = v_old + J / mω_new = ω_old + I^{-1} * ( (R*r_i) × J )
以下是刚体碰撞处理的算法流程：
- 遍历刚体每个顶点，计算其世界坐标
x_i。 - 使用有符号距离函数检测
x_i是否发生碰撞（Φ(x_i) < 0）。 - 如果发生碰撞，计算该顶点的当前速度
v_i。 - 判断
v_i是否仍指向物体内部（v_i · n < 0）。若是，则继续。 - 使用单点冲量法公式，计算该顶点期望的碰撞后速度
v_i_desired。 - 计算矩阵
K和速度变化Δv_i。 - 计算冲量
J = K^{-1} * Δv_i。 - 用冲量
J更新刚体的质心速度v和角速度ω。
实现细节：
- 如果多个顶点同时碰撞，可以取这些碰撞顶点位置的平均值作为一个代表点进行处理，以避免过度反应。
- 由于重力持续作用，物体在平面上可能持续微幅抖动。可以在速度很小时施加阻尼来缓解此现象。
处理多个刚体间的碰撞更为复杂，需要求解一个线性互补问题，本课暂不深入。
形状匹配：一种无物理的替代方案 🧩
上一节我们介绍了基于物理的冲量法。本节中，我们来看看一种完全不同的、不依赖任何物理公式的方法——形状匹配。这对于物理基础薄弱的同学可能更友好。
其核心思想分为两步：
- 自由移动：将刚体视为一堆独立的质点，每个质点根据自己的受力（如重力、碰撞力）自由运动一个时间步。这会导致形状扭曲。
- 形状匹配：将扭曲后的一堆质点，通过一个最优的刚体变换（旋转
R和平移c），重新“匹配”回它原来的参考形状。
数学描述：
假设自由移动后，顶点位置为 y_i。我们想找到质心 c 和旋转矩阵 R，使得变换后的顶点 c + R * r_i 尽可能接近 y_i。这转化为一个优化问题：
优化目标：
最小化 Σ_i || (c + R * r_i) - y_i ||^2
求解过程：
- 计算最优质心：目标函数对
c求导，发现最优质心就是所有顶点y_i的平均值。c = (1/N) * Σ_i y_i - 计算最优变换：先放宽条件，令
R为一个任意矩阵A。目标函数对A求导，可解得：A = ( Σ_i (y_i - c) * r_i^T ) * ( Σ_i r_i * r_i^T )^{-1} - 提取旋转：对矩阵
A进行 极分解，将其分解为旋转矩阵R和一个对称矩阵（代表形变）的乘积：A = R * S。我们丢弃形变部分S，保留旋转部分R作为最终结果。
得到 c 和 R 后，我们就可以更新刚体的“质心”和“朝向”。同时，可以根据顶点位置的变化反推出顶点的速度，用于下一帧的“自由移动”步骤。
形状匹配的优缺点：
- 优点：实现简单，无需物理参数；易于与基于粒子的其他模拟（如流体、软体）结合。
- 缺点：难以严格处理摩擦等约束；可能需要多次迭代来满足多个约束；碰撞响应可能不够精确（如出现滑动）。
此方法适用于对碰撞精度要求不高的场合，或者刚体作为复杂系统（如衣服上的纽扣）的一部分时。
总结 🎓
本节课中，我们一起学习了刚体碰撞的核心内容。
我们首先回顾了刚体运动模拟的基础，包括力矩、惯性张量和运动更新公式。接着，我们深入探讨了碰撞处理。从单点碰撞入手，介绍了惩罚法和冲量法两种响应策略。然后，我们将冲量法扩展到刚体，详细推导了如何通过施加冲量来更新刚体的整体运动状态。最后，我们介绍了一种无需物理知识的替代方案——形状匹配，阐述了其原理和实现步骤。
本次作业要求实现基于冲量法的刚体碰撞，附加题则是实现形状匹配算法。希望本教程能帮助你理解这些概念，顺利完成作业。
GAMES103-基于物理的计算机动画入门—P5-Lecture-05-Physics-based-Cloth-Simulation—GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将要学习基于物理的布料模拟。我们将从基础的弹簧模型开始，探讨如何用它来描述和模拟布料，并分析显式积分与隐式积分两种方法的优劣。课程的核心目标是理解如何从物理能量和力的模型出发，构建一个稳定的布料模拟系统。
课程概述与安排
我们之前几周主要讨论了刚体模拟，现在将开启一个新的主题：布料模拟。布料模拟与头发模拟有一定联系，但各有特点。我们将用三周时间深入探讨布料模拟。
第一周，也就是今天，我们将重点讲解基于物理的仿真模拟，即如何根据能量和力的模型推导出布料模拟的方式。
第二周，我们将讨论约束，包括游戏开发中常用的PBD（位置动力学）方法、Projective Dynamics以及专门处理约束的Constraints Dynamics。
第三周，我们将讨论碰撞处理。布料的碰撞处理是最具挑战性的，掌握其方法后，处理其他模拟的碰撞会相对容易。
本节课将涵盖以下内容：首先讲解如何使用弹簧描述和模拟布料，包括显式和隐式积分方法。然后讨论弯曲模型的重要性，分析弹簧弯曲模型的不足及其他替代模型。最后，如果有时间，会介绍基于Shape Matching的非弹簧模拟方法。
弹簧模型基础
上一节我们介绍了课程的整体安排，本节中我们来看看如何用弹簧模型来描述一块布料。
一个理想的弹簧满足胡克定律。该定律指出，弹簧的力试图恢复其原长，且力的大小与拉伸长度成正比。
考虑一个一维弹簧，其一端固定在原点，当前位置为 x，原长为 l。那么弹簧的能量 E 为：
E = 1/2 * k * (x - l)^2
其中 k 是弹簧的弹性系数。
力 F 是能量对位置的负导数：
F = -dE/dx = -k * (x - l)
这个公式与我们中学物理所学的完全一致。
我们可以将此概念推广到二维或三维。假设一根弹簧连接了两个点 i 和 j，原长为 l。那么弹簧的能量为：
E = 1/2 * k * (||x_i - x_j|| - l)^2
其中 ||x_i - x_j|| 是两点间的当前距离。
力是能量对点位置的负梯度。对于点 i 和 j，其受力分别为：
F_i = -k * (||x_ij|| - l) * (x_ij / ||x_ij||)
F_j = -F_i
其中 x_ij = x_i - x_j。这表明作用力与反作用力大小相等，方向相反。
构建弹簧网络系统
上一节我们介绍了单根弹簧的模型，本节中我们来看看如何用多根弹簧构建一个布料网络。
根据物理原理，能量和力是可以叠加的。如果一个顶点连接了多根弹簧，那么作用在该顶点上的总力就是所有相关弹簧力的矢量和。
在实际模拟中，我们需要构造一个弹簧网络来描述布料。主要有两种方式：
结构化网格：假设布料由规整的方格构成。每个网格交点是一个顶点。我们需要以下几种弹簧：
- 结构弹簧：连接横向和纵向相邻顶点，抵抗拉伸。
- 剪切弹簧：连接对角线上的顶点（如45度和135度方向），防止布料在斜向过度拉伸。
- 弯曲弹簧：连接隔一个顶点的两个顶点（例如，跨越一个顶点的弹簧），用于抵抗布料沿边的自由弯曲，增加布料的挺括感。
一种常见的优化是采用交错的对角线弹簧布局，这样可以在减少弹簧数量的同时，兼顾各个方向的抵抗，避免模拟产生方向偏好。
非结构化三角网格：在服装设计等领域，布料的版型通常是不规则的三角网格。处理方式如下：
- 将三角网格的每一条边都视为一根结构弹簧。
- 对于网格内部的每一条边，将其所对的两个顶点（即不在这条边上的两个顶点）用一根弹簧连接起来，这根弹簧就作为弯曲弹簧。
从三角网格数据构造弹簧系统
上一节我们知道了需要哪些弹簧，本节中我们来看看如何从给定的三角网格数据中自动提取这些弹簧信息。
程序中的三角网格通常由两个列表表示：
- 顶点位置列表：存储每个顶点的3D坐标。
- 三角形索引列表：存储每个三角形所包含的三个顶点在顶点列表中的索引。
为了构造弹簧系统，我们需要从这些数据中提取出所有唯一的边（用于结构弹簧）和所有内部边（用于生成弯曲弹簧）。
以下是处理流程：
- 生成原始边列表：遍历每个三角形，将其三条边（由两个顶点索引表示）与所属三角形索引一起，存储为一个三元组
(v_a, v_b, triangle_id)。为确保每条边有唯一的表示，我们总是将顶点索引按从小到大排序（例如，边(3,0)存储为(0,3)）。 - 排序与去重：对所有原始边三元组按顶点索引进行排序。排序后，重复的内部边会相邻出现。
- 提取结构弹簧：遍历排序后的列表。如果当前边与下一条边的顶点索引完全相同，则说明是内部边，跳过其中一条（去重）。剩下的唯一边就是最终的结构弹簧边。
- 提取弯曲弹簧：在第二步去重时，被跳过的重复边就是内部边。对于每条内部边，我们知道它相邻的两个三角形。这两个三角形中，不在这条内部边上的那两个顶点，就是弯曲弹簧需要连接的两个顶点。
通过以上步骤，我们就可以从三角网格数据中构造出完整的弹簧系统。
具体实现代码将在作业示例中提供，感兴趣的同学可以深入研究。
显式积分模拟
有了弹簧系统，我们就可以开始模拟了。首先介绍显式积分方法。
一个简单的粒子系统更新步骤如下：
- 对于每个顶点
i，计算其受到的总力F_i（包括重力、弹簧力等）。 - 根据牛顿第二定律更新速度：
v_i_new = v_i + Δt * F_i / m_i - 根据速度更新位置：
x_i_new = x_i + Δt * v_i_new
对于弹簧系统，我们只需在第一步的力计算中加入所有相关弹簧的力即可。计算每根弹簧的力时，我们需要：
- 弹簧连接的顶点索引
i和j。 - 弹簧的原长
l（可预先计算）。 - 根据公式
F = -k * (||x_ij|| - l) * (x_ij / ||x_ij||)计算力，并分别加到顶点i和j的合力上。
显式积分的问题与隐式积分
上一节我们介绍了简单的显式积分，本节中我们来看看它存在的一个严重问题：数值不稳定性。
当弹簧的弹性系数 k 很大，或模拟的时间步长 Δt 较大时，显式积分容易产生“过冲”现象。这是因为过大的力会使顶点位置更新过度，越过平衡点，甚至在下一次迭代中产生更大的力，导致顶点位置发散，最终模拟崩溃。
解决思路之一是缩小时间步长 Δt，但这会显著降低模拟效率。因此，在需要稳定模拟或处理“刚性”系统时，我们更常使用隐式积分方法。
隐式积分使用未来时刻的状态来计算力，其更新公式如下：
v_new = v + Δt * M^{-1} * F(x_new, v_new)
x_new = x + Δt * v_new
这里 F(x_new, v_new) 是在新位置和新速度下的力，是未知量。这使得方程难以直接求解。
隐式积分转化为优化问题
为了求解隐式积分方程，我们首先假设力是保守力，且只与位置有关（如重力、弹簧力），即 F = F(x)。我们可以通过代入消元法，将方程组化简为只关于新位置 x_new 的方程。
更重要的是，我们可以证明，求解这个方程等价于求解以下非线性优化问题：
x_new = argmin_x { (1/(2Δt^2)) * (x - x̃)^T M (x - x̃) + E(x) }
其中 x̃ = x + Δt * v 是仅根据当前速度预测的位置，M 是质量矩阵（通常是对角阵），E(x) 是系统的势能（如弹簧势能）。
这个优化问题的含义是：新位置 x_new 使得“惯性项”（粒子试图保持预测运动趋势）和“势能项”（系统试图降低势能）的加权和最小。
通过这种转化，我们将物理模拟问题转化为一个数学优化问题，从而可以利用成熟的优化算法来求解。
牛顿法求解优化问题
上一节我们将模拟问题转化为了优化问题，本节中我们来看看如何使用牛顿法来求解它。
对于标量函数 f(x)，寻找其极小值点等价于寻找其梯度（一阶导数）为零的点：∇f(x) = 0。牛顿法的核心思想是利用泰勒展开进行迭代逼近。
假设当前迭代点为 x_k，我们希望在这一点附近对梯度函数 ∇f(x) 做线性近似：
∇f(x) ≈ ∇f(x_k) + H(x_k) * (x - x_k)
其中 H(x_k) 是函数 f 在 x_k 处的海森矩阵（二阶导数）。
令近似梯度为零，即可解出下一步的迭代点 x_{k+1}：
x_{k+1} = x_k - H(x_k)^{-1} * ∇f(x_k)
对于我们的布料模拟优化问题，f(x) 就是之前定义的目标函数。其梯度 ∇f(x) 为：
∇f(x) = (1/Δt^2) * M * (x - x̃) + ∇E(x)
而 ∇E(x) 就是负的力 -F(x)。海森矩阵 H(x) 为：
H(x) = (1/Δt^2) * M + H_E(x)
其中 H_E(x) 是势能 E(x) 的海森矩阵，即力的雅可比矩阵（刚度矩阵）。
因此，布料模拟的隐式积分牛顿法迭代步骤为：
- 初始化新位置
x（例如用x̃）。 - 计算当前梯度
∇f(x)和海森矩阵H(x)。 - 求解线性方程组
H(x) * Δx = -∇f(x)，得到位置增量Δx。 - 更新位置：
x = x + Δx。 - 重复步骤2-4，直到收敛（如
Δx足够小）。 - 最后，用更新后的位置计算新速度：
v_new = (x_new - x) / Δt。
弹簧系统的海森矩阵与数值求解
上一节我们介绍了牛顿法的框架，本节中我们深入看一下弹簧系统海森矩阵的特点和数值求解的细节。
对于弹簧系统，整个网络的海森矩阵 H 由每根弹簧的贡献叠加而成。每根弹簧对其连接的两个顶点 i 和 j 的贡献是一个 6x6 的块矩阵，放置在大矩阵的 (i,i), (i,j), (j,i), (j,j) 块位置上。
一个重要的性质是，当弹簧被拉伸时，其贡献的块矩阵是半正定的；但当弹簧被压缩时，这部分贡献可能变得非正定。这可能导致整个海森矩阵 H 非正定。
在优化中，如果海森矩阵始终正定，则目标函数是凸函数，保证存在唯一全局极小值。非正定则意味着可能存在多个局部极小值（例如，一根被挤压的弹簧可能向左或向右弯曲，形成两种不同的稳定状态）。在实际模拟中，非正定主要影响的是某些线性求解器的稳定性。一个常见的处理技巧是，当检测到弹簧处于压缩状态时，忽略其海森矩阵中可能导致非正定的项。
求解牛顿法中的线性方程组 H * Δx = -∇f 是关键步骤。主要有两类方法：
- 直接法（如LU分解）：精度高，但对于大规模矩阵（顶点数多）内存消耗大。
- 迭代法（如共轭梯度法）：内存消耗小，适合大规模问题，但收敛速度和精度依赖于矩阵条件数，通常需要预条件子。
在实际应用中，根据问题规模和可用硬件进行选择。
经典论文与课程总结
在本节课中，我们一起学习了基于物理的布料模拟基础。
我们从最简单的弹簧模型出发，解释了如何用胡克定律描述布料的拉伸行为。接着，探讨了如何将三角网格转化为弹簧网络系统，包括结构弹簧和弯曲弹簧的构造方法。
我们分析了显式积分方法的简单直观性及其固有的数值不稳定性问题，这引出了对隐式积分的需求。通过将隐式积分方程转化为等价的非线性优化问题，我们能够更稳定地求解系统状态。牛顿法为求解该优化问题提供了框架，其中涉及梯度、海森矩阵的计算以及大规模线性方程组的求解。
最后，推荐一篇该领域的经典论文《Stable but Responsive Cloth》，它最早将隐式积分引入布料模拟，尽管其推导视角与本节课的优化视角略有不同，但本质相通，是深入理解布料模拟的重要资料。
本节课的内容，特别是隐式积分部分，是许多物理模拟方法（如有限元法）的核心基础。理解这一框架，对于后续学习更复杂的模拟技术至关重要。
GAMES103-基于物理的计算机动画入门—P6-Lecture-06-Constrained-Approaches–PBD–PD-and-others–Lab-2----GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将学习布料模拟中的约束方法，特别是位置动力学（PBD）和拉伸限制（Strain Limiting）。这些方法为处理布料等可变形体的模拟提供了稳定且高效的解决方案。课程内容将涵盖基本概念、算法实现及其优缺点，帮助你理解如何在实际项目中应用这些技术。
作业延期与提交说明 📝
关于第一次作业的提交，我们注意到部分同学时间紧张。因此，作业可以适当延期提交。但请注意，每晚提交一天，成绩将扣除20%。若延迟超过五天，该次作业成绩将不计入总分统计，但你仍可提交并获得助教的反馈。我们鼓励大家尽量在五天内完成作业。
作业相关问题可以在微信群的小程序圈子中提问。如果是普遍性问题，建议在公共论坛提问；如果是具体代码问题，可以直接联系助教。实践是学习计算机课程的关键，完成作业能帮助你深入理解知识，巩固学习成果。
第二次作业介绍：布料模拟 👕
第二次作业与当前课程内容紧密相关，目标是模拟一块布料的动态效果。作业提供了一个名为“Cloth”的场景，其中包含一块布料和一个用于交互的小球。
布料有两个固定点（顶点0和顶点20），模拟时用户可以拖拽小球与布料进行交互。上图展示了使用位置动力学（PBD）技术实现的效果。布料的弹性表现可以通过调整迭代次数来控制：迭代次数越少，布料显得越有弹性；迭代次数越多，布料则显得越硬。请注意，当前模拟未处理布料的自相交，仅处理了布料与球的碰撞。
作业包含两个部分：
- 位置动力学（PBD）方法
- 隐式积分弹簧系统方法
若要使用隐式积分方法，你需要在提供的两个脚本之间进行切换。
隐式积分部分旨在帮助大家巩固上周所学的知识。在实际作业实现中，我们采用了一种近似的简化方法，使用一个对角矩阵来近似海森矩阵（Hessian），从而将解线性系统的步骤简化为一个直接的更新操作。具体的实现公式将在作业描述中给出。这种方法虽然并非标准的牛顿法，但其收敛原理我们将在后续课程中详细探讨。
布料模拟中的弯曲模型 🔄
上一节我们介绍了弹簧系统的基本模拟。本节中，我们来看看弹簧系统在模拟布料弯曲时存在的缺陷，并探讨更优的弯曲模型。
弹簧系统的一个主要问题是模拟弯曲抵抗时不够准确。例如，使用连接两个相对顶点的“弯曲弹簧”时，当两个三角形几乎共面并发生轻微弯曲时，弹簧长度的变化极小，导致产生的抵抗力非常弱，这与真实世界中纸张或布料的弯曲感觉不符。
一个更直观的思路是利用两个三角形之间的二面角（Dihedral Angle）来定义弯曲力。这种方法被称为二面角方法。
对于一个由四个顶点构成的二面角，弯曲力可以写成一个关于二面角θ的函数。每个顶点上的力 f_i 可以表示为：
f_i = k(θ) * u_i
其中，k(θ) 决定了力的大小，u_i 决定了力的方向。
推导u_i的过程基于几个合理假设：
- 顶点1和2上的力应分别沿着各自三角形的法向 n1 和 n2。
- 弯曲不应导致连接顶点3和4的边发生形变，因此作用在顶点3和4上的力之差应垂直于该边。
- 所有弯曲力作为内力，其合力应为零。
基于这些假设，可以推导出各个u_i的表达式。而力的大小函数 k(θ) 的一种常见形式为：
k(θ) = k * (e^2 / (|n1| * |n2|)) * sin((π - θ)/2)
其中，k是常数，e是边长，|n1|和|n2|是两个三角形的面积（即未归一化的法向量模长）。若要定义非平面的静止状态（初始夹角θ0），公式可修改为包含 sin((π - θ0)/2) 项。
这个基于力的模型优点是与显式积分兼容，但缺点是没有明确的能量定义，用于隐式积分时求导会比较麻烦。
另一种更简单的模型是基于余切权重的拉普拉斯弯曲模型（Discrete Shells）。该模型基于两个假设：布料静止时为平面，且模拟中拉伸形变很小。
它定义了一个能量函数 E(x)：
E(x) = (k/2) * (x^T * Q * Q^T * x)
其中，x 是所有顶点位置排成的向量，Q 是一个由余切权重（cotangent weights）构成的常数矩阵。这个能量函数本质上是在计算离散曲率的平方。当布料完全平整时，曲率为零，能量也为零；弯曲时，能量增大，产生恢复力。
该模型的优点是能量函数是二次的，其梯度（力）和海森矩阵都非常容易计算，f = -k * Q * Q^T * x，H = k * Q * Q^T，非常适合隐式积分。缺点是当布料发生较大拉伸时，其准确性会下降。
布料模拟中的锁定问题 🔒
接下来，我们讨论布料模拟中一个重要的现象：锁定问题（Locking Issue）。
理想情况下，布料的拉伸和弯曲应该是互不影响的。但在弹簧系统模型中，当弹簧非常刚硬（stiff）时，它可能会阻止布料在某些方向上的弯曲。例如，一个由两个三角形组成的布料，如果沿着对角线弯曲可能很自由，但若想沿着另一条边弯曲，中间的弹簧就会像支柱一样顶住，导致无法弯曲。
这种现象的本质是自由度丢失。对于一个流形网格，其自由度数量约为 3 + 边界边数量。当网格分辨率较低且弹簧非常刚硬时，约束（边）的数量几乎与变量（顶点坐标）数量相当，导致系统没有足够的自由度来表现弯曲，从而被“锁死”。
解决锁定问题没有完美方案，常见的工程性方法包括：
- 降低拉伸弹簧的弹性系数。
- 允许弹簧在一定长度范围内自由活动而不产生力。
- 提高网格分辨率（让顶点更密集）可以在视觉上缓解问题。
但需要注意的是，这些方法并不能从根本上解决数学模型本身的局限性。有限元方法同样存在类似的锁定问题。
约束方法：位置动力学（PBD）🎯
上一节我们探讨了模型本身的问题。本节开始，我们将学习一类强大的模拟方法——约束方法，它旨在稳定地模拟高刚度材料。
首先从单根弹簧的约束开始。我们希望弹簧保持原长 l0，即约束函数 C(xi, xj) = |xi - xj| - l0 = 0。
投影函数（Projection Function） 的目标是：找到新的顶点位置 xi’, xj’，在满足约束 C=0 的前提下，使位置变化量 Δx 尽可能小，同时保持质心不变。这可以转化为一个数学优化问题。
推导得到的投影公式为：
xi’ = xi + (wi / (wi + wj)) * (|xj - xi| - l0) * ((xj - xi) / |xj - xi|)
xj’ = xj - (wj / (wi + wj)) * (|xj - xi| - l0) * ((xj - xi) / |xj - xi|)
其中，wi, wj 是顶点的权重，通常与质量成反比。固定点可以通过设置极大质量来实现。
对于多根弹簧（多个约束），有两种主要的迭代处理方法：
- 高斯-赛德尔（Gauss-Seidel）方法：按顺序逐一处理每条边，立即更新其顶点位置，然后处理下一条边。处理完所有边后，再从头开始多次迭代，直至约束被较好地满足。这种方法简单，但更新顺序会影响结果，且难以并行化。
- 雅可比（Jacobi）方法：对每条边，计算其建议的顶点位置更新量，但不立即应用。对于每个顶点，收集所有相邻边对其的更新建议，然后取这些建议的平均值（或加权平均）来最终更新顶点位置。这种方法易于并行化，能减少顺序依赖带来的偏差。
基于投影思想，位置动力学（Position Based Dynamics, PBD） 算法流程如下：
- 对所有顶点进行常规的粒子运动更新（如显式积分）：v = v + Δt * f_ext / m, x = x + Δt * v。
- 基于当前顶点位置 x，通过多次迭代（高斯-赛德尔或雅可比方法）应用约束投影，得到满足约束的新位置 x’。
- 根据位置变化更新速度：v = (x’ - x) / Δt。
- 用 x’ 覆盖 x。
PBD的特点：
- 优点：概念简单，易于实现和并行化，通用性强（可用于布料、流体、刚体等），在低分辨率模型上效率高。
- 缺点：缺乏明确的物理参数（如杨氏模量），材料的“刚度”由迭代次数和网格密度等非物理参数控制，没有收敛到精确解的理论保证，在高分辨率下效率下降较快。
PBD非常适合游戏等实时应用中对低精度可变形体的模拟。
约束方法：拉伸限制（Strain Limiting）📏
PBD缺乏物理含义，而拉伸限制（Strain Limiting） 可以看作是其与物理模拟结合的一个改进。它的核心思想是：在正常的物理模拟步骤之后，增加一个后处理步骤，将过大的形变“拉回”到可接受的范围内，从而增强模拟的稳定性。
对于弹簧，我们定义拉伸率 σ = l / l0。我们允许拉伸率在一个范围内 [σ_min, σ_max] 内变化。投影函数的目标不再是严格保持原长，而是将当前拉伸率 σ 修正到目标范围 clamp(σ, σ_min, σ_max) 内。
修正后的目标长度 l’ = clamp(σ, σ_min, σ_max) * l0。然后使用与PBD类似的投影公式，将弹簧长度调整到 l’ 而非 l0。当 σ_min = σ_max = 1 时，就退化成了PBD。
对于三角形面积约束，思路类似：计算当前面积A，将其修正到允许范围 [A_min, A_max] 内，得到目标面积 A’。然后通过缩放三角形（保持质心不变）来使面积接近 A’，缩放因子 s = sqrt(A’ / A)。
拉伸限制的用途：
- 增强稳定性：防止物理模拟在发生大形变时出现数值不稳定、抖动甚至崩溃。
- 模拟非线性材料：例如某些布料，在小拉伸时很柔软，达到一定阈值后变得非常刚硬。可以在小形变时用物理模拟，大形变时启用拉伸限制来近似这种效果。
- 缓解锁定问题：通过允许一定范围内的形变，相当于降低了有效刚度，使弯曲更容易发生。
课程总结 🎓
本节课我们一起学习了布料模拟中的核心约束方法。
- 我们首先分析了弹簧系统在弯曲模拟上的不足，并介绍了二面角模型和基于余切权重的拉普拉斯弯曲模型。
- 然后，我们探讨了布料模拟中常见的锁定问题及其成因。
- 接着，我们深入学习了两种重要的约束求解方法：位置动力学（PBD） 和 拉伸限制（Strain Limiting）。PBD通过直接投影位置来满足约束，简单高效，适合实时应用；而拉伸限制则作为物理模拟的稳定器，将过大形变限制在合理范围内。
- 这些方法为处理高刚度材料和保证模拟稳定性提供了实用工具。
在下节课中，我们将继续探讨更高级的约束方法，如投影动力学（Projective Dynamics），它们试图在物理准确性和计算效率之间取得更好的平衡。
GAMES103-基于物理的计算机动画入门—P7-Lecture-07-Other-Constrained-Methods-and-Finite-Element-Method-I—GAMES-Webinar—BV12Q4y1S73g_note
在本节课中，我们将学习两种基于约束的高级模拟方法：Projective Dynamics 和 Constrained Dynamics，并初步了解有限元方法（FEM）的基本概念和线性三角形单元的力计算。
课程内容调整说明 📅
由于课程计划调整，本次课程将首先完成约束方法的讲解，然后开始介绍有限元方法。碰撞检测的内容将移至有限元方法之后讲解。
第一部分：Projective Dynamics 方法 🔄
上一节我们介绍了PBD和Shape Matching等约束方法。本节中，我们来看看Projective Dynamics方法，它同样基于约束，但构建方式有所不同。
Projective Dynamics的核心思想是：不直接用投影函数修改顶点位置，而是用投影结果定义一个能量函数。
方法原理
对于一个弹簧边约束，我们定义能量函数为：
E = (1/2) * k * ||(xi - xj) - (xi_new - xj_new)||^2
其中，xi_new和xj_new是根据约束投影计算出的“目标”位置，但我们不直接赋值。
经过推导可以发现，这个能量函数和力公式与普通弹簧系统完全一致：
E_spring = (1/2) * k * (||xi - xj|| - L)^2
F_spring = -k * (||xi - xj|| - L) * ( (xi - xj) / ||xi - xj|| )
那么，Projective Dynamics的特殊之处在哪里？关键在于其Hessian矩阵（刚度矩阵）。
常数Hessian矩阵的优势
在Projective Dynamics中，我们做了一个关键假设：投影得到的目标位置 x_new 是常数（与当前顶点位置 x 无关）。这使得能量函数 E(x) 成为顶点位置 x 的纯二次函数。
因此，其Hessian矩阵 H = d²E/dx² 是一个常数矩阵。对于一个由弹簧边构成的网格，其Hessian矩阵非常容易构造：
- 对角元：顶点度数 * k
- 非对角元（若顶点i和j有边相连）：-k
常数Hessian矩阵带来了巨大优势：在隐式积分求解线性系统 (M - Δt² H) Δv = Δt F 时，我们只需要对矩阵 (M - Δt² H) 做一次LU分解或Cholesky分解。在后续所有时间步中，都可以复用这个分解结果，只需进行高效的回代求解，从而大幅提升CPU上的模拟效率。
方法总结
Projective Dynamics通过引入基于约束投影的中间变量，构造了一个具有常数Hessian矩阵的系统。其优点是在CPU上效率很高，且具有明确的物理含义（等同于弹簧系统）。缺点是收敛速度后期可能变慢，且在GPU上直接法求解不占优，约束改变时矩阵也需要更新。
第二部分：Constrained Dynamics 方法 🔗
接下来，我们探讨另一种处理强约束的方法——Constrained Dynamics。此方法常用于模拟多刚体系统（如人体关节）。
核心思路
该方法旨在处理刚度（Stiffness）极大甚至无穷大的约束。我们引入一个辅助变量——拉格朗日乘子 λ。
首先，将弹簧能量改写为关于约束 φ(x) = ||xi - xj|| - L 的形式：
E = (1/2) * φᵀ * C⁻¹ * φ
其中 C 是柔度矩阵（C = (1/k) I）。
通过计算力的梯度并引入 λ = -C⁻¹ φ，我们可以得到力的新形式：
F = -Jᵀ * λ
这里 J 是雅可比矩阵 J = dφ/dx。
构建求解系统
结合动量守恒方程和约束方程（在新状态下 φ_new = 0 的泰勒展开），我们得到一个包含两个变量（速度 v 和拉格朗日乘子 λ）的线性系统：
[ M -Δt Jᵀ ] [ v_new ] = [ M * v_old ]
[ Δt J C ] [ λ_new ] [ -φ ]
有两种求解方式：
- 联合求解（Primal-Dual）：直接求解上述完整系统得到
v_new和λ_new。 - 消元求解：先将
λ消去，得到一个只关于v的系统（该形式与隐式积分系统等价），求解后再计算λ。
方法特点与应用
Constrained Dynamics 的优点是可以完美处理刚性约束（当 k → ∞ 时，C → 0，系统依然可解）。它广泛应用于多刚体动力学模拟，例如游戏中的布娃娃（Ragdoll）动画。
第三部分：有限元方法（FEM）初步 📐
现在，我们开始进入有限元方法的世界。有限元是模拟连续体（如弹性体、布料）形变的重要工具。
基本概念：线性三角形单元
我们以二维三角形单元为例。假设三角形在参考（未形变）状态的顶点为 X0, X1, X2，在当前（形变后）状态的顶点为 x0, x1, x2。
有限元的核心假设是：三角形内部的形变是均匀的，即内部任意点的位移映射可以通过一个线性函数描述：
x = F * X + c
其中 F 是一个2x2矩阵，称为形变梯度（Deformation Gradient）；c 是平移向量。
计算形变梯度 F
利用三角形两条边的向量变化，可以计算出 F：
F = D_s * D_m⁻¹
其中：
D_s = [x1-x0, x2-x0]是当前状态的边向量矩阵。D_m = [X1-X0, X2-X0]是参考状态的边向量矩阵。
格林应变张量（Green Strain）
形变梯度 F 包含了旋转和形变。为了剔除旋转分量，只度量纯形变，我们使用格林应变张量 G：
G = (1/2) * (FᵀF - I)
G 是一个对称矩阵，其大小与形变成正比，且在物体纯旋转时为零。
能量与应力
我们定义单位面积的能量密度 W 为应变 G 的函数。对于三角形，总弹性势能为：
E = A_ref * W(G)
其中 A_ref 是三角形在参考状态下的面积。
一种常用的简单模型是Saint Venant-Kirchhoff模型：
W(G) = μ * trace(GᵀG) + (λ/2) * [trace(G)]²
其中 μ 和 λ 是拉梅参数，控制材料的剪切和体积变形阻力。
能量密度对应变的导数称为第二皮奥拉-基尔霍夫应力（Second Piola-Kirchhoff Stress） S：
S = dW/dG = 2μG + λ * trace(G) * I
计算顶点力
最终目标是计算作用在顶点上的力 F_i = -dE/dx_i。通过链式法则和一系列推导（此处省略繁琐的中间步骤），可以得到每个顶点力的简洁计算公式：
对于顶点1和顶点2（以边 X1-X0 和 X2-X0 为基准）：
f1 = -A_ref * F * S * (D_m⁻¹)的列1
f2 = -A_ref * F * S * (D_m⁻¹)的列2
顶点0的力由力平衡得出：
f0 = -f1 - f2
这个矩阵形式公式极大简化了实现难度，避免了直接对复杂能量函数求导的繁琐过程。
总结 🎯
本节课我们一起学习了：
- Projective Dynamics：通过约束投影构造能量，得到常数Hessian矩阵以加速CPU求解。
- Constrained Dynamics：引入拉格朗日乘子处理强约束，常用于多刚体系统。
- 有限元方法初步：介绍了线性三角形单元的基本概念，包括形变梯度
F、格林应变G、能量密度W、应力S，并给出了计算顶点力的核心矩阵公式。
这些方法为模拟复杂物体提供了更多工具。下一节课，我们将继续深入有限元方法，并探讨另一种实现方式——有限体积法（FVM）。
GAMES103-基于物理的计算机动画入门—P8-Lecture-08-Finite-Element-Method-II–Lab-3----GAMES-Webinar—BV12Q4y1S73g_note
概述
在本节课中，我们将继续学习有限元方法，并重点介绍一种与之等价但推导更清晰的有限体积法。我们将学习如何通过计算表面牵引力的积分来得到力，并理解不同应力定义（如柯西应力、第一/第二皮奥拉-基尔霍夫应力）之间的关系与转换。最后，我们会介绍超弹性材料模型，并简要探讨非线性优化在物理模拟中的应用。
作业说明：Lab 3
上一节我们介绍了有限元方法的基本概念，本节中我们来看看如何将其应用于实际模拟。课程作业（Lab 3）已经发布。
作业提供了一个Unity脚本和代码框架。你需要完成的核心任务是：使用有限元方法计算四面体网格形变产生的力，并更新顶点位置和速度，以模拟一个弹性小房子掉落到地面并弹起的动画。
以下是作业的关键点：
- 模型：使用名为“house”的四面体网格。
- 材料模型：采用St. Venant-Kirchhoff (STVK) 模型，因其公式相对简单。
- 积分方法：使用显式积分进行模拟。
- 注意事项：
- STVK模型无法处理四面体的反转（即顶点穿透到另一侧）。如果形变过大导致四面体反转，它将无法恢复。
- 显式积分存在稳定性问题。如果形变过于剧烈或时间步长不当，模拟可能会“爆炸”。
- 编程代码量可能是所有作业中最短的，但理解背后的原理需要仔细阅读课件。
核心挑战在于如何根据有限元理论计算每个四面体作用于顶点的力。理解这一点是完成作业的关键。
从有限元到有限体积法 🔄
上节课我们通过定义形变能，并对位置求导来得到力。本节我们将从另一个角度——有限体积法——来推导力，这对于四面体、三角形这类线性单元而言，与有限元法是等价的，且推导过程更清晰直观。
有限体积法的核心思想是：力来源于物体内部应力在边界表面的作用。
假设一个弹性体被一个界面分割为两部分。该界面上单位面积（或单位长度）所受的力称为牵引力。那么，整个界面上的合力就是牵引力在界面上的积分。
如何得到牵引力呢？这需要引入应力的概念。应力本质上是一个矩阵（二阶张量），它将表面的法向量映射为该点的牵引力。具体关系由以下公式描述：traction = stress * normal
其中，traction是牵引力向量，stress是应力张量，normal是表面单位法向量。
因此，一个区域所受的合力，就是其边界上 stress * normal 的积分。
计算三角形对顶点的力贡献 📐
让我们具体到一个三角形单元。考虑顶点 x0，我们假设它代表周围的一个区域。那么，作用在这个区域上的合力，就是其边界上牵引力的积分。
对于三角形而言，x0 所代表的区域边界由两条从 x0 出发的边的中点连线构成。我们假设三角形内部的应力和形变是均匀的（线性单元假设），因此应力 σ 在三角形内为常数。
根据散度定理，对一个闭合曲线的法向量积分结果为零。因此，计算 x0 所受的力，可以转化为计算两条边（x0x1 和 x0x2 的中点连线）上 -σ * normal 的积分。
由于边上法向量恒定，积分简化为应力乘以法向量再乘以边长的一半。由此，我们得到三角形对顶点 x0 的力贡献公式：f0 = - (σ * n10) * |x10|/2 - (σ * n20) * |x20|/2
其中，n10 和 n20 是两条边的单位法向量，|x10| 和 |x20| 是边长。
为什么取边的中点？ 这是为了将三角形产生的力平均地分配到三个顶点上，满足动量守恒。
扩展到四面体情况 🔶
在三维情况下，我们需要处理四面体。一个顶点（如 x0）与三个面相邻。此时，力是应力在三个相邻面上的面积分。
与三角形类似，为了将每个面的贡献均匀分配到其三个顶点上，我们引入系数 1/3。顶点 x0 所受的力来自三个相邻面（如 012, 023, 031）的贡献之和。
经过推导和简化（利用叉积计算面积和法向量），我们可以得到四面体对顶点 x0 的力贡献的紧凑公式：f0 = - (σ * (x20 × x30 + x30 × x10 + x10 × x20)) / 6
其中 × 表示叉积。类似地，可以计算 f1, f2, f3。
目前，公式中的应力 σ 尚未定义。接下来我们将解决这个问题。
理解不同的应力定义 🧮
这里的关键在于，上节课有限元法推导出的应力与本节有限体积法直接用来计算力的应力，并不是同一种应力。
它们的区别在于所作用的法向量和牵引力所处的空间不同：
- 参考配置：物体未形变时的状态。
- 当前配置：物体形变后的状态。
根据法向量和牵引力所处的配置不同，定义了多种应力：
- 第二皮奥拉-基尔霍夫应力：法向量和牵引力都定义在参考配置。上节课通过能量密度对格林应变求导得到的就是它。
- 第一皮奥拉-基尔霍夫应力：法向量定义在参考配置，牵引力定义在当前配置。
- 柯西应力：法向量和牵引力都定义在当前配置。这正是我们有限体积法公式中直接需要的应力。
我们需要找到一种方法，从上节课可计算的第二皮奥拉-基尔霍夫应力，得到本节需要的柯西应力。
应力之间的转换关系
不同应力之间可以通过形变梯度 F 进行转换。形变梯度描述了从参考配置到当前配置的局部线性变换。
以下是核心转换关系：
- 第一皮奥拉-基尔霍夫应力 P 与第二皮奥拉-基尔霍夫应力 S 的关系：
P = F * S - 柯西应力 σ 与第一皮奥拉-基尔霍夫应力 P 的关系：
σ = (1 / det(F)) * P * F^T
因此，结合两者，我们可以从已知的 S 计算出所需的 σ：σ = (1 / det(F)) * F * S * F^T
应用于力计算公式
现在我们将柯西应力 σ 的表达式代入四面体的力计算公式中。但直接使用上述转换公式会引入复杂的计算。我们可以采用一个技巧：
观察四面体力公式 f0 = - (σ * (x20 × x30 + ...)) / 6，括号内的叉积项实际上与当前配置下的法向量有关。如果我们将其替换为参考配置下顶点位置 X 的叉积，那么对应的应力就应该使用第一皮奥拉-基尔霍夫应力 P。
这样做的好处是，参考配置下的顶点位置 X 是固定的，因此由它们计算的叉积项是常数，可以预先计算并存储，大大减少了实时模拟的计算量。
经过替换和整理（P = F * S），我们得到最终用于计算的实用公式：f0 = - (F * S * B0) / 6
其中 B0 是一个基于参考配置顶点位置预先计算的常向量：B0 = X20 × X30 + X30 × X10 + X10 × X20。F 是形变梯度，S 是第二皮奥拉-基尔霍夫应力。
类似地，f1 = - (F * S * B1) / 6, f2 = - (F * S * B2) / 6。根据动量守恒，f3 = - (f0 + f1 + f2)。
算法伪代码
以下是模拟循环中计算一个四面体内力的核心伪代码：
// 预先计算（在初始化阶段）
Dm = matrix(x10, x20, x30) // 参考配置下的边矩阵
inv_Dm = inverse(Dm) // 其逆矩阵
B0 = X20 × X30 + X30 × X10 + X10 × X20 // 预计算常向量
B1 = ... // 类似计算 B1, B2
B2 = ...
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/5bc90e468d9b6b2f277c79542d262c4f_10.png>
// 每帧计算（对于每个四面体）
// 1. 计算形变梯度 F
Ds = matrix(x1-x0, x2-x0, x3-x0) // 当前配置下的边矩阵
F = Ds * inv_Dm
// 2. 计算格林应变 G 和第二皮奥拉-基尔霍夫应力 S
G = (F^T * F - I) / 2
S = dΨ/dG // 根据选定的能量密度函数 Ψ 计算，例如 STVK 模型
// 3. 计算第一皮奥拉-基尔霍夫应力 P
P = F * S
<https://github.com/OpenDocCN/cs-notes-pt3-zh/raw/master/docs/games/img/5bc90e468d9b6b2f277c79542d262c4f_12.png>
// 4. 计算各顶点力（以 f1 为例）
f1 = - (P * B1) / 6
// 类似计算 f2, f3
f2 = - (P * B2) / 6
f3 = - (P * B3) / 6
// 由动量守恒得到 f0
f0 = - (f1 + f2 + f3)
// 5. 将 f0, f1, f2, f3 分别加到对应顶点的合力上
超弹性材料模型 🧱
之前我们使用的 STVK 模型虽然简单，但有明显缺陷：当四面体被严重压缩或反转时，其产生的恢复力会变小甚至方向错误，导致无法恢复原状，模拟失真。
在图形学和力学中，更常用的是超弹性材料模型。这类模型通过定义能量密度函数 Ψ 来描述材料特性，应力由能量密度导出（S = dΨ/dG）。
对于各向同性材料（如橡胶，各个方向性质相同），其能量密度可以表示为形变梯度 F 的主拉伸 λ1, λ2, λ3 的函数。主拉伸是 F 经过奇异值分解后得到的对角矩阵元素，代表了三个主方向上的拉伸比率。
常见的各向同性超弹性模型有：
- 新胡克模型：
Ψ = μ/2 * (I1 - 3) + λ/2 * (J - 1)^2，其中I1是主拉伸平方和，J是体积变化率。它能更好地阻止体积压缩，避免 STVK 的反转问题。 - 穆尼-里夫林模型：新胡克模型的增强版。
- 冯·米塞斯模型：常用于模拟生物软组织。
在实现上，使用这些模型需要计算形变梯度 F 的奇异值分解以获得主拉伸，然后根据具体的 Ψ 公式计算应力 S。这比 STVK 复杂，但物理上更准确、更稳健。
非线性优化简介 ⚙️
物理模拟中的隐式积分等方法，最终需要求解一个非线性优化问题（即最小化系统总能量）。这里简要介绍其核心思想。
目标是找到变量 x 使函数 f(x) 最小化。基本思路是迭代：从初始点 x_k 出发，寻找一个下降方向 d_k 和合适的步长 α_k，使得 x_{k+1} = x_k + α_k * d_k 的函数值减小。
如何选择下降方向 d_k？
- 最速下降法：
d_k = -∇f(x_k)，即负梯度方向。简单但收敛可能慢。 - 牛顿法：
d_k = -H^{-1} ∇f(x_k)，其中H是海森矩阵（二阶导数）。利用了曲率信息，收敛快，但需计算并求解海森矩阵，代价高。 - 拟牛顿法/投影法：用近似矩阵
M代替海森矩阵H，d_k = -M^{-1} ∇f(x_k)。在收敛速度和单步计算成本之间寻求平衡，在图形学中应用广泛。
如何确定步长 α_k？常用回溯线搜索：从一个初始步长开始，不断按比例缩小，直到满足函数值充分下降的条件。
在物理模拟中，我们追求的是总计算时间 = 每次迭代耗时 × 迭代次数 的最小化。图形学的研究重点常在于为特定硬件和问题设计高效的优化方法，以在速度和精度间取得最佳平衡。
总结
本节课我们一起深入探讨了有限元方法的另一种视角——有限体积法。
- 我们学习了如何通过积分表面牵引力来计算顶点力，并推导了三角形和四面体的具体公式。
- 我们厘清了柯西应力、第一/第二皮奥拉-基尔霍夫应力等不同应力定义的区别与联系，并掌握了它们之间的转换方法，从而得到了实用的力计算公式。
- 我们认识了 STVK 模型的局限性，并介绍了更通用、更稳健的超弹性材料模型（如新胡克模型）。
- 最后，我们简要了解了非线性优化的基本概念，它是许多高级模拟方法（如隐式积分）背后的数学工具。
希望本讲内容能帮助你更好地完成 Lab 3 作业，并加深对基于物理的动画中弹性体模拟的理解。下节课我们将进入一个新的重要主题：碰撞检测与处理。
GAMES103-基于物理的计算机动画入门—P9-Lecture-09-Collision-Handling—GAMES-Webinar—BV12Q4y1S73g_note
概述
在本节课中，我们将学习计算机动画中一个至关重要且实践性极强的主题：碰撞处理。我们将从如何检测碰撞开始，逐步深入到如何处理碰撞响应，涵盖离散与连续碰撞检测、内点法与冲击优化等核心方法，并探讨在复杂情况（如布料自碰撞）下的解决方案。
作业技巧与有限元模拟调试 🛠️
上一节我们结束了关于弹性体有限元的讨论。在进入新主题前，先针对上周布置的作业提供一些调试技巧。
作业提供了一个包含约1000个四面体的房屋模型。直接模拟整个模型难以调试，因为一旦结果出错，很难定位问题所在。
以下是调试建议：
- 从单个四面体开始测试：不要直接运行房屋模型。首先尝试使用一个单独的四面体进行模拟。
- 测试静止状态：关闭重力，让四面体处于静止状态。此时，形变梯度
F应为单位矩阵I。由此计算出的格林应变ε应为零矩阵0，进而应力σ和力f也应为零。这可以验证你的力计算在平衡状态下是否正确。- 公式：
F = I->ε = 0->σ = 0->f = 0
- 公式：
- 测试拉伸状态：将四面体沿某个方向（如X轴）拉伸。弹性体的力总是阻碍形变的，因此每个顶点上的力理论上应试图将四面体拉回原状。检查力的方向是否符合这个趋势。
- 操作方法：通过缩放四面体顶点的位置来实现拉伸。
- 关于性能：在Unity中运行此类模拟帧率可能较低，这受脚本执行和硬件配置影响。纯粹的CPU实现通常效率更高。四面体有限元模拟在游戏开发中通常只适用于分辨率很低的情况。
关于作业的其他问题：
- 速度场平滑：可对速度场进行拉普拉斯平滑，即每个点的速度与其邻居点的速度取加权平均。
- SBD的力公式：建议自行推导验证，确保公式正确性。
- 模型中的三角网格：用于从四面体网格构造渲染表面，若无直接用途可忽略。
碰撞处理概述 🤔
本节我们将开始讨论碰撞处理这个新话题。碰撞处理的概念相对简单，但要想实现得鲁棒高效，算法复杂度非常高，需要处理各种边界情况。
本节课的目标是为大家建立基础概念，以便未来需要时能更快地上手。课程将分为以下几个部分：
- 碰撞检测：如何发现碰撞。
- 碰撞响应方法：介绍两种主流方法——内点法和冲击优化法。
- 离散碰撞处理：介绍相交解除方法。
布料碰撞是所有碰撞处理中最困难的。若能处理好布料碰撞，其他碰撞问题便可视为其简化版。
碰撞检测 🔍
处理碰撞的第一步是知道碰撞是否发生。碰撞检测通常分为两个阶段，以提高效率。
广义阶段碰撞剔除
当场景中有大量元素（如数万个三角形）时，对所有元素进行两两碰撞检测的计算量是平方级的，无法承受。因此，第一阶段（Broad-phase Collision Culling）的目的是剔除那些绝不可能发生碰撞的元素对，输出一组候选碰撞对。
以下是两种常见的广义阶段剔除方法：
1. 空间划分法
核心思想是将整个3D空间划分为均匀网格。每个三角形根据其位置（或运动轨迹所占用的空间）被存储到与之相交的网格单元格中。
检测流程：
- 对于每个三角形，找出它占据的所有网格单元格。
- 要判断某三角形可能与哪些其他三角形相交，只需检查它所在单元格内存储的所有其他三角形。
- 将这些三角形对作为候选对输出。
优点：实现相对简单，对GPU友好。
缺点：需要为整个空间预分配网格内存，可能造成浪费；物体运动后需要更新网格信息。
优化方案——空间哈希与排序：
为了节省内存，不直接分配网格数组，而是为每个“三角形-单元格”关联创建一个记录（如 (cell_id, tri_id)）。然后对所有记录按 cell_id 排序。排序后，连续且 cell_id 相同的记录就代表了该单元格内所有的三角形。这避免了为空单元格分配内存。
进一步优化——莫顿码：
在三维空间中，按行、列、层顺序访问网格会导致内存访问不连续（尤其在访问空间邻居时）。莫顿码是一种空间填充曲线，它将多维空间索引编码成一维，使得空间上相邻的单元格在内存地址上也更接近，从而优化缓存访问效率。
2. 包围体层次结构
核心思想是对物体本身进行层次划分，而非对空间进行划分。例如，一个机器人模型可以按肢体划分成多个部分，每个部分用一个简单的包围体（如轴对齐包围盒AABB、球体）包裹。
检测流程（以BVH树为例）：
- 从根节点开始，检查两个包围体是否相交。
- 若不相交，则其下所有子元素都不可能相交，整条分支可剔除。
- 若相交，则递归地检查它们的子节点。
- 对于自碰撞检测，需要额外检查同一节点下不同子节点之间的包围体是否相交。
优点：更新效率高（只需更新包围体），适合结构化的、运动连贯的物体。
缺点：实现较复杂；树形结构对GPU不友好；对于密集近邻的剔除效率不如空间划分法。
其他方法：有研究尝试用形变能量的大小来预测碰撞可能性（形变越大，碰撞可能越高），但对布料等常有大形变却无碰撞的情况效果有限。
方法对比总结：
| 特性 | 空间划分法 | 包围体层次法 |
|---|---|---|
| 实现难度 | 较简单 | 较复杂 |
| GPU友好度 | 友好 | 不友好 |
| 更新成本 | 高（需重建） | 低（更新包围体即可） |
| 适用场景 | 元素分布均匀/随机 | 物体结构性强、层次清晰 |
狭义阶段碰撞检测 🎯
经过广义阶段剔除，我们得到了一批候选碰撞对。在狭义阶段（Narrow-phase Collision Detection），我们需要对这些候选对进行精确的检测，判断它们是否真的发生了碰撞或相交。根据后续处理方式的不同，检测也分为两类。
离散碰撞检测
DCD检测的是在某个特定时刻，几何体之间是否相交（即穿透），而不关心运动过程。
基本操作：对于三角形网格，最基本的是检测边与三角形是否相交。
- 计算边与三角形所在平面的交点参数
t。 - 判断该交点是否位于三角形内部。
- 若是，则报告相交。
优点：计算简单、快速、稳定。
缺点：可能产生隧道效应。如果物体运动速度过快，可能在两帧之间直接穿过另一个物体，而在这两个离散时刻都未检测到相交。这对于薄物体（如布料、兔耳朵）尤为严重。
连续碰撞检测
CCD检测的是在两个时刻之间的运动过程中，几何体是否发生了碰撞。它关注的是第一个接触点，而非穿透状态。
基本元素：对于运动的三角形网格，基本的连续碰撞检测有两种：
- 点-三角形检测：一个运动的点与一个运动的三角形。
- 边-边检测：两条运动的边。
以点-三角形检测为例：
假设四个点都在做线性运动。碰撞发生时，这四个点共面，即它们构成的四面体体积为零。根据体积公式，可以推导出一个关于时间 t 的一元三次方程。
- 公式推导：
Volume = ( (x1(t)-x0(t)) × (x2(t)-x0(t)) ) · (x3(t)-x0(t)) = 0，展开后得到a*t^3 + b*t^2 + c*t + d = 0。
解此方程得到可能的碰撞时间，再验证该时刻点是否落在三角形内部。
为什么需要边-边检测？因为两个三角形可能只在边上发生碰撞，而点-三角形检测无法捕捉这种情况。
CCD的挑战：
- 数值稳定性：求解一元三次方程时，应避免使用求根公式（浮点误差大），推荐使用数值方法（如二分法、牛顿法），并只关注时间区间
[0, 1]内的根。 - 计算成本：比DCD高，因为需求解方程。
- 实现复杂度：高，容易遇到各种边界情况。
尽管CCD计算更复杂，但在广义阶段剔除大部分元素后，其成本通常可以接受。CCD是避免隧道效应、实现高质量碰撞的必要手段。
碰撞响应：内点法与冲击优化法 ⚖️
当我们检测到碰撞（或相交）后，下一步就是处理它，将系统修正到一个无碰撞的状态。假设模拟后得到了一个有碰撞的状态，而我们的目标是安全区域（无碰撞状态）。有两种主流的数学思路来解决这个问题。
内点法
核心思想：从安全区域内（上一帧的无碰撞状态）出发，采取小步长迭代，确保每一步迭代后的中间状态都保持在安全区域内，最终逼近模拟目标附近的无碰撞解。
类比：小心翼翼地从安全区走向目的地，绝不踏出安全区边界。
优点：
- 必然成功：即使未完全收敛，只要在安全区内停止，结果就是可接受的（无碰撞）。
- 数值稳定。
缺点：
- 速度慢：需要小步长保证安全，且可能迭代步数很多。
- 计算量大：每步都需要进行碰撞检测以确保安全。
实现示例——对数障碍函数：
在能量函数中添加一个关于距离的对数障碍项。当距离趋近于零时，该项产生的排斥力趋近于无穷大，从而防止穿透。图形学中常会进行截断以限制影响范围。
- 优化目标：
min ||x - x_target||^2 + λ * Σ -log(d_i(x))，其中d_i(x) > 0为距离约束。
冲击优化法
核心思想：胆子更大。直接从模拟得到的不安全状态出发，通过优化或投影操作，试图一步或几步将其拉回安全区域，不关心中间过程是否安全。
类比：从目的地（不安全）直接跳回或拉回安全区。
优点：
- 速度快：通常碰撞只发生在局部，可以集中计算资源进行优化，步长可以较大。
- 高效：只处理出问题的区域。
缺点：
- 可能失败：如果初始状态离安全区太远，可能无法收敛到无碰撞解。
- 不保证严格无碰撞：结果可能仍有轻微穿透。
实现方式：通常构建为一个约束优化问题，使用拉格朗日乘子法、增广拉格朗日法等求解。例如，检测到穿透时，施加一个冲量将物体推开。
刚性冲击法
当上述方法都失败或过于复杂时，一种“最后一搏”的方案是刚性冲击法。它将发生复杂碰撞的整个区域锁定为一个刚体，禁止其内部相对运动，从而从根本上避免继续碰撞。
优点：简单、安全、绝对有效。
缺点：会引入明显的视觉瑕疵，失去局部形变细节。
实际工作流建议：
- 首先尝试快速的冲击优化法。
- 若失败，且计算资源允许，则回退到鲁棒但较慢的内点法。
- 若仍无法解决或需要实时性能，则对碰撞区域启用刚性冲击法作为保底。
碰撞响应：相交解除 ✂️
对于某些应用（如游戏、CAD软件），目标可能不是物理上精确的碰撞响应，而是快速消除视觉上的穿插。这就是“相交解除”的思路。
核心思想：承认当前帧可能存在相交，并采用一种后处理技术，直接修改顶点位置来消除穿插，而不严格考虑物理正确性或运动过程。
对于有体积的物体
这相对简单。常用符号距离场（SDF）来表示物体体积。对于每个发现穿透的点，计算其到SDF表面的最短方向，并将其沿该方向推出表面即可。
- 之前作业中的球与布料碰撞即采用此思路。
对于布料自碰撞
布料是开放的曲面，无法定义内外，因此SDF方法失效。需要更复杂的算法。
思路一（2003年）：将相交的布料区域分段，假设面积较小的片段是穿插产生的“多余”部分，并设法将其移除或平滑。但该方法难以处理边界处的穿插。
思路二（2006年）：将布料间的相交视为一条空间曲线。通过优化，最小化这条相交曲线的长度。当曲线长度为零时，穿插即被解除。该方法能更好地处理边界情况，但也不是万能的。
相交解除是一种实用的、面向视觉效果的策略，常作为物理碰撞响应失败后的补充或替代方案。
总结与课程预告 📚
本节课我们一起学习了碰撞处理的完整流程：
-
碰撞检测：
- 广义阶段：通过空间划分或包围体层次结构剔除不可能碰撞的元素对，提高效率。
- 狭义阶段：
- 离散碰撞检测：检测特定时刻的相交，快速但可能有隧道效应。
- 连续碰撞检测：检测运动过程中的碰撞，精确但计算复杂，需解一元三次方程。
-
碰撞响应：
- 内点法：从安全区出发，步步为营，保证中间状态安全，鲁棒但慢。
- 冲击优化法：直接修正不安全状态，高效但可能失败。
- 刚性冲击法：将碰撞区域冻结为刚体，简单粗暴的保底方案。
- 相交解除：以消除视觉穿插为目标的实用后处理技术，尤其适用于布料等复杂情况。
碰撞处理是一个概念易懂但实现艰难的领域，充满了算法细节和工程挑战。理解这些基本方法为应对实际开发中的碰撞问题奠定了重要基础。
下周预告：我们将开启全新的主题——流体模拟。课程将分为三部分：基础介绍、基于网格的欧拉方法、以及基于粒子的SPH方法。
GAMES105-计算机角色动画基础—P1-Lecture01-Introduction-to-Character-Animation—GAMES-Webinar—BV1GG4y1p7fF_note
在本节课中，我们将要学习计算机角色动画的基础概念、主要技术分类以及该领域的发展脉络。课程旨在为初学者提供一个清晰的概览，理解如何让虚拟角色“动”起来。
课程与讲师介绍
大家好，欢迎参加GAMES105课程。本课程名为“计算机角色动画基础”。与其他GAMES系列课程类似，我们主要关注三维计算机角色动画，而非二维动画。课程会涉及二维内容，但核心研究对象是三维角色动画。
我是刘立宾，现任北京大学智能学院助理教授。我的研究方向主要是基于物理的角色动画，特别是如何控制角色运动。在加入北大之前，我曾在迪士尼研究院工作。右下角的二维码是我们北京大学可视计算实验室的公众号，欢迎大家关注，我们会通过公众号发布相关活动信息。
课程主页位于GitHub。上课时间为每周一晚上八点到九点，持续约12周，从今天开始到12月底或1月初结束。
前置要求：学习本课程需要了解线性代数、微积分和基本的编程技巧。课程作业将主要使用Python完成，我们会提供简单的代码框架。作业提交将通过专门的网站进行，课程邀请码已显示在屏幕上。此外，我们有一个课程交流群，二维码长期有效。
什么是计算机角色动画？
首先，我们需要明确什么是计算机角色动画，以及我们关注其中的哪些内容。
提到角色动画，大家可能会想到动画电影，如《疯狂动物城》、《寻梦环游记》等。虽然动画电影是重要的应用领域，但当前电影制作中使用的计算机角色动画技术比例其实很小。本课程讲解的许多技术尚未在电影中广泛应用，但随着AI生成技术的进步，未来可能会有更多应用。
实际上，角色动画技术应用最广泛的领域是游戏。因为游戏需要根据用户输入实时生成可交互的内容，这与一次性渲染的电影不同。游戏的需求直接推动了角色动画技术的进步与发展。
近年来，我们也看到许多新应用，例如虚拟偶像、虚拟主播、VR社交、VR游戏，以及大场景人群仿真等。这些都属于计算机图形学的组成部分。
角色动画的核心组成部分
如果想让电影中的一个角色运动起来，至少涉及三个领域：建模、动画和渲染。
- 建模 关注静态外观或某一帧的表现。
- 渲染 负责将外观展示出来。
- 动画 则解决相邻帧之间如何生成的问题，它建模的是时间序列上的规律。
这与物理仿真有显著不同。物理仿真（如刚体、软体、流体仿真）的目标是模拟客观物理现象的时间演变。而角色动画更关注行为建模，例如人或动物的动作。
两者之间存在联系。在我看来，物理仿真加上控制，就构成了动画。纯粹的物理仿真是被动的，而控制引入了主观意愿。物理仿真通常有精确的数学描述（如微分方程），而动作行为很难用单一数学公式完全描述，往往需要通过大量观察进行统计建模。
角色动画关注的对象是“角色”，通常具有数十到上百个关节参数。理论上，几乎所有动作都可以由动画师手动“K帧”（逐帧调整姿态）完成，从而得到高质量动画。但这使得角色动画在需要交互的场景（如游戏）中成为一种劳动密集型技术。
因此，计算机角色动画的研究目标，是理解人和动物产生动作的规律并建立模型。一方面，这能提供更智能的动画编辑工具；另一方面，可以生成全新的动作。简而言之，我们的目标是将角色动画从劳动密集型转变为计算密集型。
角色动画的生成流程
如何生成一段动画或让一个虚拟形象动起来？通常流程如下：
- 几何建模：获得虚拟形象的模型（通过雕刻、扫描等方式）。
- 骨骼绑定：将模型绑定到骨骼上。这涉及骨骼设置和计算蒙皮权重。
- 运动生成：驱动骨骼运动，从而带动外部模型形变。如何生成骨骼的运动，是本课程主要关注的内容。
在真实世界中，人做一个动作（如伸手）的过程是：大脑想法 -> 神经信号 -> 肌肉脉冲 -> 肌肉收缩 -> 施加力/力矩于骨骼系统 -> 在物理规律下产生运动。
计算机角色动画技术可以根据是否对物理过程进行建模来区分。
两大技术分类：运动学 vs. 动力学
基于运动学（关键帧）的方法
这类方法不模拟物理过程，而是隐含了从意图到姿态的转换。它直接更新角色的状态（如姿势、速度）。角色可以实现瞬移等不符合物理规律的动作。这是工业界应用非常广泛的一类方法。
基于动力学（物理）的方法
这类方法希望复现真实的物理过程。它会对物理仿真进行建模，通过仿真生成最终动作。在此过程中，不能直接干预角色姿态，因此瞬移并非其目标。
总结：角色动画方法大致分为两类：使用物理仿真的基于物理的角色动画，和不使用物理仿真、直接改变姿态的基于运动学（关键帧）的方法。
控制维度：低级与高级
控制角色动作还可以从不同维度进行：
- 低级控制：像动画师一样，逐帧调整每个关节的旋转。优点是可以精确控制每个细节；缺点是效率低、成本高。
- 高级控制：给角色一个高级目标（如“去拿杯水”），角色自动生成完整动作序列。优点是用很少的信息生成复杂动画；缺点是无法精确控制细节。
角色动画的发展趋势，正是从低级控制逐渐向高级控制过渡。
接下来，我们将简要介绍角色动画各个方向的基本技术。
基于运动学的方法
上一节我们介绍了角色动画的两大分类，本节中我们来看看基于运动学（关键帧）方法的具体技术。
关键帧动画与迪士尼12准则
关键帧动画是最基础的动画技术。早在计算机图形学出现之前，动画师就通过逐帧绘制来制作动画（如1917年的《猫和老鼠》）。他们总结出了著名的迪士尼动画12准则，部分准则是为了模仿真实物理过程，部分则与艺术夸张和叙事相关。
在三维计算机图形学时代，同样的技术可以应用于三维领域。动画师使用Maya、Blender等软件，一帧一帧地调整角色姿态，连接起来形成动画。制作高质量的关键帧动画需要相当的训练和技术。
从算法角度看，有两个非常重要的技术：
前向运动学与逆向运动学
- 前向运动学：给定每个关节的旋转角度，计算末端（如手）的位置。
- 公式：
末端位置 = FK(关节1角度, 关节2角度, ...)
- 公式：
- 逆向运动学：给定末端（如手）的目标位置，反推需要改变每个关节多少旋转角度。
- 公式：
(关节1角度, 关节2角度, ...) = IK(末端目标位置)
IK是角色动画中的关键技术，它让姿态调整变得直观。
- 公式：
补间动画
制作关键帧动画时，通常只制作几个关键姿态，然后通过算法自动生成中间的过渡帧，即补间动画。基本思路是使用不同的插值方法来生成平滑的动画。
然而，关键帧动画需要逐帧控制，是一种非常低效的低级控制方法。
动作捕捉与数据重用
为了解决关键帧动画效率低的问题，动作捕捉技术被广泛应用。
动作捕捉技术
- 光学动捕：使用专用设备（如Vicon），在身体上贴反光标记点，通过多视角相机捕捉。质量高，用于电影制作。
- 惯性动捕：使用便携的惯性传感器。价格更亲民，常用于虚拟偶像直播。
- 视频动捕：从普通视频中估计动作。使用范围广，但目前质量难以与前两者相比。
无论哪种动捕，都会面临一个问题：如何将捕捉到的动作映射到不同的虚拟角色上？这需要动作重定向技术。
动作重定向
动作重定向是将一个动作适配到骨骼结构、尺寸可能完全不同的新虚拟角色上的技术。
动作捕捉本质上只是记录和回放动作，无法生成新动作。为了在交互式应用（如游戏）中重用数据，产生了以下技术：
状态机与运动图
一个简单的想法是使用状态机。例如，角色有“奔跑”和“攻击”两个动作，用户按下攻击键时，从奔跑状态切换到攻击状态，再切回奔跑。大部分游戏引擎都支持状态机模型。
学术上，2002年提出的运动图技术与此类似。它从运动数据中自动构建状态机，允许在不同动作片段间切换。后续工作对其进行了改进，例如在状态节点内进行动作插值以实现更精确控制，或在其上训练AI来完成更高级任务。
运动图的优点是能重用现有数据生成新运动。缺点是构造复杂，随着动作数量增加，图会变得非常庞大且容易出错。
运动匹配
运动匹配是育碧公司在2016年左右提出的技术。其核心思想不是播放完整动作片段，而是在更细粒度（如每一帧）进行控制。每一帧结束时，通过最近邻搜索找到一个新姿态，该姿态既要满足控制目标（如移动方向），又要与当前状态连贯。
运动匹配很大程度上是工程实现，需要精心设计损失函数和动作库。它简单实用，现已集成到一些游戏引擎中。它同样解决了运动图的难点，并能提供高级控制。
生成模型与跨模态生成
我们提到，角色动画的目标是理解动作的内在统计规律。生成模型为计算机角色动画带来了巨大进步。
生成模型
生成模型可以从大量动作数据中学习规律，根据用户输入自动生成下一帧姿态。例如，可以学习控制指令（如摇杆输入）与动作之间的关系。
常见的生成模型包括GAN、VAE、标准化流，以及最新的扩散模型。这些模型只需收集原始动作数据（如让人随意走动两小时）进行训练即可生成动作，而传统方法（如运动图）需要手动分割和构建状态机。
跨模态生成
当前的研究前沿是实现跨模态生成，即用高级的、其他模态的指令来控制角色动作。
- 音乐生成舞蹈：输入一段音乐，生成对应的舞蹈动作。
- 语音生成动作：输入一段语音（如演讲），生成对应的肢体语言和动作。
- 文本生成动作：输入一句话，生成与之匹配的动作。
这非常接近我们“通过高级控制生成复杂动作”的原始目标。
基于运动学方法的局限性
尽管基于运动学的方法取得了很大进展，但它仍有相当大的局限性：
- 物理准确性不足：容易产生穿模、脚底打滑等问题。
- 动作范围受限：本质上是现有动作的重放或组合。当环境变化需要全新的交互时（如被绊倒、完成危险特技），如果动作库中没有类似数据，则很难生成。
为了解决这些问题，我们需要回到基于物理仿真的方法。
基于物理仿真的方法
上一节我们探讨了基于运动学方法的优势与局限，本节我们将深入基于物理仿真的角色动画。
基于物理的方法旨在在虚拟世界中复现真实世界产生动作的完整过程。回想一下，基于运动学的方法是：输入信号 + 当前姿态 -> 运动模型 -> 下一帧姿态。而基于物理的方法多了一步：运动模型输出的不是直接姿态，而是控制量（如力），这些力通过物理仿真最终改变姿态。
布娃娃仿真
最基本的物理仿真应用是布娃娃仿真，即关闭角色的控制，只让其在物理规律下运动。这常用于表现角色死亡、失去意识或被绊倒瞬间的反应。此时生成的动作是无意识的。
物理仿真的应用场景
- VR/AR交互：在VR环境中，用户可能从任意方向击打角色，只有物理仿真能生成真实的受力反馈。
- VR社交：在只有头显和手柄的情况下，通过物理仿真可以合理估计出下半身的运动，实现全身虚拟形象。
- 精细操作：如使用筷子，物理仿真能逼真地生成手指的细微动作。
物理控制的基本技术：PD控制
在物理仿真中，我们通常将角色简化为由关节力矩驱动的多刚体系统，而非复杂的肌肉模型。
一类常用方法是PD控制。例如，我想把手从A高度举到B高度，需要施加多少力？一个粗略的计算是：力 = k_p * (目标位置 - 当前位置) + k_d * (目标速度 - 当前速度)。通过给出目标轨迹并使用PD控制，可以生成关节力矩，驱动角色运动。
这项技术很早（1995年）就已出现，但至今未在游戏中大规模应用，主要原因在于控制非常困难。
轨迹优化与简化模型
为了设计出合理的控制轨迹，产生了其他方法：
- 轨迹优化/时空优化：通过数学优化自动计算出合适的控制轨迹。它允许高级输入（如指定脚部落点），但求解高维非线性优化问题非常困难且缓慢。
- 简化模型：借鉴机器人学思路，使用高度简化的模型（如倒立摆）来描述和指导动作。这种方法非常稳定，能抵抗外力干扰，但生成的动作缺乏细节，且通常只针对特定动作（如行走）设计，难以实现复杂动作（如后空翻）。
强化学习
强化学习为物理角色动画带来了突破。智能体通过与环境不断交互（试错）来更新运动策略，最终学会某种动作。2015年DeepMind在雅达利游戏上的突破鼓舞了该方向。
通过强化学习，可以实现运动跟踪控制，即在物理仿真环境下复现给定的运动数据。在此基础上，可以进一步结合状态机、运动匹配等思想，完成更复杂的技能组合。
控制生成模型
与基于运动学的生成模型类似，在物理仿真领域，我们可以学习控制生成模型。它学习的是角色在产生大量数据时，其控制策略的规律。我们可以在这个生成模型的状态空间中进行采样，每个采样都会产生不同的动作，同时也能响应高级指令和环境干扰。
基于生成模型的物理控制，正在将基于物理的角色动画推向“通过高级控制生成动作”的目标。未来，跨模态生成（如语言、音乐驱动）也有望在物理仿真中实现。
课程总结与展望
本节课我们一起回顾了计算机角色动画领域过去30年的主要研究方向。
基于运动学的方法（如运动匹配、生成模型）已逐渐在真实场景中落地。而基于物理的方法，虽然30年来一直被寄予厚望，但直到最近，随着控制生成模型等技术的出现，才真正接近可用的状态。我们可能正处在一个转折点，基于物理的方法有望在虚拟世界中更真实地还原动作产生过程。
课程安排与信息
最后，重申一下本课程的相关信息。
课程定位：本课程不会教授具体软件（如Maya、Unreal Engine）的使用，也不会培养动画师。我们将从科学和技术的角度，探讨角色动画背后的方法、理论和算法，特别是角色运动学、动力学、物理仿真与控制。最终目标是让大家能自己实现一个可交互的虚拟角色。
课程大纲：共12次课。今天为第一次。后续将从基础数学开始，先介绍基于运动学的方法，中间涉及蒙皮绑定基础知识，后半部分专注于物理仿真与基于物理的角色动画。大纲可能根据反馈调整。
课程项目：共有5个小项目，其中4个关于关节角色动画（每2-3周一个），1个关于模型绑定。
实现方式：所有项目主要使用Python实现。我们会使用一些Python下的物理引擎，并在其上实现仿真、控制和动画效果。
前置课程：本课程内容自成体系，但建议大家有时间可以学习前面的GAMES系列课程（101-104），它们提供了计算机图形学、几何、建模、仿真和游戏引擎方面的良好基础。本课程会在特定点上进行更深入的探讨。
课程资料：本课程没有指定教材，主要资料是课件和授课PPT。对物理仿真或机器人学相关背景的了解会对学习有所帮助。
好的，我们今天的课程就到这里。下周我们将从数学基础和前向运动学相关内容开始讲起。谢谢大家，我们下周再见！
GAMES105-计算机角色动画基础—P10-Lecture09-Actuating-Simulated-Character—GAMES-Webinar—BV1GG4y1p7fF_note
在本节课中，我们将学习如何驱动一个物理仿真中的虚拟角色，使其产生我们期望的动画。我们将从回顾刚体仿真的基本概念开始，然后深入探讨如何通过施加力和力矩来控制角色，并重点介绍一种基础且重要的控制方法——比例微分（PD）控制。
回顾：刚体仿真基础
上一节我们介绍了如何对一个虚拟角色进行物理仿真。在计算机角色动画中，我们处理的虚拟角色通常由不会变形的刚性肢体（刚体）通过关节连接而成。
为了计算这样一个系统的运动，我们需要求解其运动方程。对于一个质点，牛顿第二定律 F = m * a 描述了力、质量和加速度的关系。通过积分加速度，我们可以更新质点的速度和位置。
然而，刚体比质点更复杂，因为它具有形状和朝向。朝向（旋转）是一个非线性量，这带来了计算上的挑战。刚体的运动状态由位置、朝向、线速度和角速度共同描述。
为了改变刚体的运动，我们需要施加外力（F）和外力矩（τ）。线力改变线速度，力矩改变角速度。刚体的“惯性”由质量（m，抵抗线运动变化）和惯性张量（I，抵抗旋转运动变化）描述。惯性张量与物体的质量和形状分布有关。
多个刚体通过关节连接时，关节会施加约束力，确保连接点不会分离。整个多刚体系统的运动方程可以写成一个紧凑的矩阵形式：M(q) * v̇ + C(q, v) = τ_ext + J_c^T * λ
其中 M 是质量矩阵，q 和 v 是状态和速度，C 包含科里奥利力和离心力，τ_ext 是外部广义力（包括我们施加的控制力），J_c^T * λ 代表约束力。
在物理引擎中定义角色时，我们需要为每个刚体提供质量、惯性张量、初始状态（位置、朝向、速度、角速度）以及用于碰撞检测的简化几何形状（如胶囊体、立方体）。同时，需要定义连接刚体的关节类型和初始位置。
仿真器的核心流程是一个前向动力学过程：根据当前状态和所有外力（包括控制力和约束力），求解运动方程得到加速度，然后通过积分更新速度和状态，从而得到下一时刻的动作。
如何驱动角色：施加力与力矩
上一节我们介绍了仿真系统如何运作，本节中我们来看看如何主动地驱动角色运动。如果不对角色施加任何控制力，它只会在重力作用下瘫倒在地，这种状态称为“布娃娃”（Ragdoll）仿真。
为了让角色按照我们的意愿运动，我们需要施加控制力。对于单个刚体，我们可以在其质心上施加合力和合力矩。对于由关节连接的多刚体系统（如人形角色），有两种主要的驱动方式：
- 直接对每个刚体施加力和力矩：这种方式自由度高，但控制复杂。
- 施加关节力矩：这是更常见和直观的方式。关节力矩模拟了生物关节处肌肉或机器人关节处电机产生的扭矩效果。
关节力矩的本质：当我们在一个关节上施加力矩 τ 时，其物理效果等价于在该关节连接的两个刚体上分别施加一个大小相等、方向相反的力矩。通常，在子刚体上施加 +τ，在父刚体上施加 -τ。这样，两个刚体会产生相对的旋转趋势，而整个系统的总力矩和为零，不会影响系统整体的质心运动。
控制器的作用：控制器是一个模块，它根据角色当前的状态（位置、速度等）以及我们期望的运动目标（如目标姿势或路径），实时计算出每个关节需要施加的力矩 τ，然后交给仿真器执行。这个过程可以看作是一个逆向动力学问题：已知期望的运动（或加速度），求解需要施加的力。
核心控制策略：比例微分（PD）控制 🎯
为了计算具体的关节力矩，我们需要一个控制策略。比例微分控制是一种基础、强大且广泛应用的反饋控制方法。
什么是比例微分控制？
我们可以通过一个简单例子来理解：控制一个沿竖直杆滑动的物块，使其到达并停留在目标高度 x_d。
- 比例（P）控制：施加的力与当前位置和目标位置的误差成正比。
F_p = k_p * (x_d - x)。k_p称为刚度系数。误差越大，施加的力越大。 - 问题：仅使用比例控制，物块会在目标位置附近持续振荡，无法稳定停下。
- 微分（D）控制：引入与速度成正比的阻尼力。
F_d = -k_d * v。k_d称为阻尼系数。速度越大，反向的阻尼力越大，用于消耗能量、抑制振荡。 - 比例微分（PD）控制：结合两者。
F = k_p * (x_d - x) - k_d * v。
对于关节控制，公式是类似的。对于某个关节，我们希望其角度 θ 跟踪目标角度 θ_d：τ = k_p * (θ_d - θ) - k_d * θ̇
其中 τ 是需要施加的关节力矩，θ̇ 是当前角速度。
参数的影响
以下是 k_p（刚度）和 k_d（阻尼）两个参数对控制效果的影响：
- 刚度
k_p过小：角色显得“柔软无力”，难以到达目标姿势。 - 刚度
k_p过大：角色显得“非常僵硬”，能快速跟踪目标，但可能引发数值不稳定，动作不自然。 - 阻尼
k_d过小：角色在跟踪目标时会产生明显的、持续的振荡。 - 阻尼
k_d过大：角色运动“缓慢迟滞”，需要很长时间才能到达目标姿势，但最终能够到达。
全身控制与轨迹跟踪
对于全身角色，我们对每一个关节都独立地使用PD控制器。我们为每个关节设计一条随时间变化的目标角度轨迹 θ_d(t)。在每一仿真时刻，控制器读取当前所有关节的状态 [θ, θ̇] 和对应的目标 [θ_d(t), θ̇_d(t)]，通过PD公式计算出所有关节的力矩 τ，驱动角色运动。
欠驱动系统与“上帝之手” ✋
然而，直接使用PD控制跟踪运动捕获数据或关键帧动画时，效果往往不理想，角色很容易失去平衡摔倒。其根本原因在于：人形角色是一个欠驱动系统。
- 完全驱动系统（如固定基座的机械臂）：控制力的自由度 ≥ 系统状态自由度。理论上可以精确控制每一个状态。
- 欠驱动系统（如自由站立的人形角色）：控制力的自由度 < 系统状态自由度。关节力矩的总和为零，因此无法直接控制整体的质心位置和全身朝向这些全局状态。
在欠驱动系统中，微小的跟踪误差（如质心稍微偏离）会因为没有直接的控制手段去修正而不断累积，最终导致失控摔倒。
为了解决这个问题，一个常见但取巧的方法是引入根节点力（Root Force），或称“上帝之手”。即在角色的根刚体（如骨盆）上直接施加一个额外的、虚拟的PD控制力，用于跟踪根节点的目标轨迹。
openEuler 是由开放原子开源基金会孵化的全场景开源操作系统项目，面向数字基础设施四大核心场景（服务器、云计算、边缘计算、嵌入式），全面支持 ARM、x86、RISC-V、loongArch、PowerPC、SW-64 等多样性计算架构
更多推荐
0
0- 0
已为社区贡献32条内容
所有评论(0)