Merge branch 'USTC-CG:main' into main

Framework3D/source/nodes/nodes/geometry/mass_spring/MassSpring.h

@@ -1,7 +1,7 @@
 #pragma once 
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
-#include <unordered_set>
+#include <set>
 #include "utils.h"
 #include <chrono>
 
@@ -13,16 +13,9 @@
 
 namespace USTC_CG::node_mass_spring {
 
-using Edge = std::pair<int, int>;
-struct hashEdge {
-    size_t operator()(const Edge &x) const
-    {
-        return x.first ^ x.second;
-    }
-};
-
 using namespace Eigen;
-using EdgeSet = std::unordered_set<Edge, hashEdge>;
+using Edge = std::pair<int, int>;
+using EdgeSet = std::set<Edge>;
 using MatrixXd = Eigen::MatrixXd;
 using SparseMatrix_d = Eigen::SparseMatrix<double>;
 using Trip_d = Eigen::Triplet<double>;

Framework3D/source/nodes/nodes/geometry/node_mass_spring.cpp

@@ -18,8 +18,6 @@ namespace USTC_CG::node_mass_spring {
 
 // -------------------- helper functions (No need to modify) --------------------
 
-using EdgeSet = std::unordered_set<Edge, hashEdge>;
-
 Eigen::MatrixXi usd_faces_to_eigen(
     const pxr::VtArray<int>& faceVertexCount,
     const pxr::VtArray<int>& faceVertexIndices)

Homeworks/8_mass_spring/README.md

@@ -52,6 +52,9 @@ ID_姓名_Homework*/
 |   └── ...
 ├── node_mass_spring.cpp
 ├── report.pdf                    // 实验报告
+├── CompositionGraph.json               // 节点连接信息
+├── GeoNodeSystem.json
+├── RenderGraph.json
 └── ...                           // 其他补充文件
 ```
 

Homeworks/8_mass_spring/documents/README-part2.md

@@ -16,7 +16,7 @@ $$
 
 对于(1)式中的能量，我们之前看待的角度为：给定原长为 $L_i$ 的弹簧，让 $\mathbf{x}_i$ 变化，不同的 $\mathbf{x}_i$ 带来了不同的能量。 $\mathbf{x}_i$ 具有旋转的自由度，只要其范数 $= L_i$ 即可最小化能量。
 
-但是Liu等人提出，我们可以从另一个角度看待这个能量：如果固定 $\mathbf{x}_ i$ ，让弹簧原始的长度与方向变化，也就是把弹簧表示为一个向量 $\mathbf{d}$ ，那么 $(\|\mathbf{x} _ {i}\| - L)^2$ 可以被看作是一个优化问题的解： $\min_{\| \mathbf{d}\| =r}\| \mathbf{x}_{i} - \mathbf{d}  \|^2$ 。 通过“三角形两边之差一定小于第三边”的原理我们可以容易地证明这一点，如下图所示。
+但是Liu等人提出，我们可以从另一个角度看待这个能量：如果固定 $\mathbf{x}_ i$ ，让弹簧原始的长度与方向变化，也就是把弹簧表示为一个向量 $\mathbf{d}$ ，那么 $(\|\mathbf{x} _ {i}\| - L)^2$ 可以被看作是一个优化问题的解： $\min_{\| \mathbf{d}\| =L} \| \mathbf{x}_{i} - \mathbf{d}  \|^2$ 。 通过“三角形两边之差一定小于第三边”的原理我们可以容易地证明这一点，如下图所示。
 
 <div  align="center">    
  <img src="../images/liu13-1.png" style="zoom:40%" />
@@ -106,7 +106,7 @@ $$
 
 需要注意的是：求解的时候因为 $\mathbf{y}$ 的定义为  $\mathbf{y} := \mathbf{x}^n + h \mathbf{v}^n + h^2 \mathbf{M}^{-1} \mathbf{f}_{\text{ext}}$，而每次迭代求解的 $\mathbf{x}$ 为 $\mathbf{x}^{n+1}$，所以 $\mathbf{y}$ 不需要更新，仍然用上一步的 $\mathbf{x}^n$ 和 $\mathbf{v}^n$。
 
-整个方法的流程示意图如下：
+整个方法的流程示意图如下，鼓励大家在通过本文档了解基本思想后，自己阅读原论文获得更全面的理解。
 
 <div  align="center">    
  <img src="../images/liu13-pipeline.png" style="zoom:100%" />

Homeworks/8_mass_spring/documents/README.md

@@ -63,7 +63,7 @@ $$
 E_i = \frac{k}{2} (\|\mathbf{x} _ {i1} - \mathbf{x}_{i2}\| - L)^2 
 $$
 
-我们记录 $\mathbf{x}_{i} := \mathbf{x} _ {i1} - \mathbf{x} _ {i2}$
+我们记 $\mathbf{x}_{i} := \mathbf{x} _ {i1} - \mathbf{x} _ {i2}$
 那么总体的能量为：
 
 $$
@@ -230,7 +230,7 @@ $$
 % \end{align}
 $$
 
-那么总体的Hessian $\mathbf{H} \in \mathbf{R}^{3n\times 3n}$ 为单根弹簧Hessian $\mathbf{H} \in \mathbf{R}^{3 \times 3}$ 按照顶点索引组装起来：
+那么总体的Hessian $\mathbf{H} \in \mathbf{R}^{3n\times 3n}$ 为单根弹簧Hessian $\mathbf{H} \in \mathbf{R}^{3 \times 3}$ 按照顶点索引组装起来（比如第i个顶点通过弹簧和第j、k个顶点相连，那么3nx3n的Hessian中第3i行3i列、第3i行3j列，第3j行3i列，第3j行3j列对应的3x3 block会有值，同理第3i行3i列、第3i行3k列，第3k行3i列，第3k行3k列对应的3x3 block也会有值，那么第3i行3i列的3x3 Hessian block就有多个弹簧的贡献加在一起）：
 
 <div  align="center">    
  <img src="../images/hessian_assemble.png" style="zoom:40%"/>
@@ -256,6 +256,8 @@ $$
 \nabla^2 g \Delta \mathbf{x} = -\nabla g \tag{7}
 $$ 
 
+在最优化中，求解一个优化问题可能需要迭代多次，直到收敛（ $\|\nabla g\|$ 小于阈值 ），我们这里出于效率考虑，在一个时间步内只进行一次牛顿法的迭代，你也可以尝试在一个时间步内让牛顿法迭代多次直到收敛，并比较两种做法的仿真结果。
+
 > 这里其实还涉及到一个Line Search的部分，一般基于线搜索方法优化问题的流程：1. 先确定搜索方向 $\mathbf{p}$（我们这里 $\mathbf{p} = (\nabla^2 g)^{-1}\mathbf{\nabla}g$ ），2. 然后确定要前进的步长 $\alpha$（该步骤称为Line Search），3. 最后更新 $\mathbf{x}^{n+1} = \mathbf{x}^n - \alpha \mathbf{p}$  。
 >由于牛顿法的推荐步长是1，这里我们就不额外进行Line Search
 
@@ -310,7 +312,7 @@ Eigen::SparseMatrix<double> MassSpring::computeHessianSparse(double stiffness)
 > 能量 $g$ 的Hessian的正定性问题：牛顿法并不是无条件收敛，也就是牛顿法给出的下降方向不一定能够使得能量真的下降！即不满足 $((\nabla^2 g)^{-1}\nabla g)^{\top} \nabla g > 0$ . 只有 $\nabla^2 g$ 正定的时候才能保证收敛。
 > 
 > 你可以首先不管这个问题，看看仿真结果如何。如果出现问题，为了让 $\nabla^2 g$ 正定，你可以尝试：
-> 1. 在计算弹簧能量的Hessian时，**当$L_i > \|\mathbf{x}_i \|$ 时**，令第 $i$ 根弹簧 $\mathbf{H}_i$ 近似为 $\mathbf{H}_i \approx k \frac{\mathbf{x}_i {\mathbf{x}_i}^\top}{\|\mathbf{x}_i\|^2}$ . 
+> 1. 在计算弹簧能量的Hessian时，**当 $L_i > \|\mathbf{x}_i \|$ 时**，令第 $i$ 根弹簧 $\mathbf{H}_i$ 近似为 $\mathbf{H}_i \approx k \frac{\mathbf{x}_i {\mathbf{x}_i}^\top}{\|\mathbf{x}_i\|^2}$ . 
 > 2. 为 $\nabla^2 g$ 对角线加上 $\epsilon \mathbf{I}$， $\epsilon$ 为可调参数，来让Hessian最小的特征值大于0. 
 > 3. 对 $\nabla^2 g$ 做SVD分解，然后精确地获取其最小特征值，令其大于0，再重新用SVD得到新的Hessian（速度预期会很慢）
 > 
@@ -339,7 +341,7 @@ $$
 $$  
 
 
-如果实现正确，将1. 劲度系数`stiffness`和2.时间步长`h`设置为合理的值（隐式时间积分不需要调阻尼系数），并考虑了Hessian的正定性：就可以看到下面的仿真结果（gif经过加速），可以实现比半隐式时间积分大20倍甚至更多的时间步长（但由于需要组装Hessian并求解线性方程组，隐式时间积分每一步的时间会比半隐式时间积分长）：
+如果实现正确，将1. 劲度系数`stiffness`(如100-1000)和2.时间步长`h`（如0.01）设置为合理的值（隐式时间积分不需要调阻尼系数），并考虑了Hessian的正定性：就可以看到下面的仿真结果（gif经过加速），可以实现比半隐式时间积分大20倍甚至更多的时间步长（但由于需要组装Hessian并求解线性方程组，隐式时间积分每一步的时间会比半隐式时间积分长）：
 
 <div  align="center">    
  <img src="../images/implicit-euler.gif" style="zoom:80%" />
@@ -435,7 +437,7 @@ Eigen::MatrixXd MassSpring::getSphereCollisionForce(Eigen::Vector3d center, doub
 
 在实践中，我们发现牛顿法求解弹簧质点仿真还是有点慢。（我们提供了`TIC`和`TOC`宏来打印程序运行时间）
 
-如何加速？我们将在Part2进行介绍，这也是本次作业的选做内容。
+如何加速？我们将在[Part2](./README-part2.md)进行介绍，这也是本次作业的选做内容。
 
 ## 参考资料
 1. GAMES 103 Lecture 2 & 5