简介
隐式积分法
显示积分简单而言是通过, 过去的求解未来. 而隐式积分, 简单而言是我要求解现在, 但是我的未知量中也有现在的未知量. 简单而言就是需要通过方程组的思想来进行求解.
参考文献
代码参考师弟 ~~
对于cloth问题, 简而言之, 有两个变量需要我们求解. 即速度v和位置x.
\[\left\{\begin{array}{l}\mathbf{x}^{[1]}=\mathbf{x}^{[0]}+\Delta t \mathbf{v}^{[0]}+\Delta t^{2} \mathbf{M}^{-1} \mathbf{f}^{[1]} \\ \mathbf{v}^{[1]}=\left(\mathbf{x}^{[1]}-\mathbf{x}^{[0]}\right) / \Delta t\end{array}\right.
\]
上述为了求解\(x^{[1]}\), 我们用到了\(\mathbf{f}^{[1]}\), 两个都是未来的变量, 需要通过方程组来进行求解.
数学家将隐式积分的问题转换成求解
\[\mathbf{x}^{[1]}=\operatorname{argmin} F(\mathbf{x}) \quad$ for $\quad F(\mathbf{x})=\frac{1}{2 \Delta t^{2}}\left\|\mathbf{x}-\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\right\|_{\mathbf{M}}^{2}+E(\mathbf{x})
\]
碰撞检测
使用基于脉冲法. 其实在lab2.pdf中也有讲解
\[\mathbf{v}_{i} \leftarrow \mathbf{v}_{i}+\frac{1}{\Delta t}\left(\mathbf{c}+r \frac{\mathbf{x}_{i}-\mathbf{c}}{\left\|\mathbf{x}_{i}-\mathbf{c}\right\|}-\mathbf{x}_{i}\right), \quad \mathbf{x}_{i} \leftarrow \mathbf{c}+r \frac{\mathbf{x}_{i}-\mathbf{c}}{\left\|\mathbf{x}_{i}-\mathbf{c}\right\|}
\]
求解函数
\[\left(\frac{1}{\Delta t^{2}} \mathbf{M}+\mathbf{H}\left(\mathbf{x}^{(k)}\right)\right) \Delta \mathbf{x}=-\frac{1}{\Delta t^{2}} \mathbf{M}\left(\mathbf{x}^{(k)}-\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\right)+\mathbf{f}\left(\mathbf{x}^{(k)}\right)
\]
中的\(\Delta \mathbf{x}\)来计算新的坐标和位移
A
A 就是 \(\left(\frac{1}{\Delta t^{2}} \mathbf{M}+\mathbf{H}\left(\mathbf{x}^{(k)}\right)\right)\)
其中的 hessian 矩阵比较难求. 作者通过简化A
A =
\[\frac{1}{\Delta t^{2}} m_{i}+4 k
\]
G 梯度
G 其实是 -b
为什么梯度就是 -b 呢??
因为作者使用的是牛顿迭代法
牛顿迭代法有一个特性
\[0=F^{\prime}(x) \approx F^{\prime}\left(x^{(k)}\right)+F^{\prime \prime}\left(x^{(k)}\right)\left(x-x^{(k)}\right)
\]
一个函数的一阶导数等于其一阶导数+二阶导数×偏差.
也就是 \(-F^{\prime}(x) = F^{\prime \prime}\left(x^{(k)}\right)\left(x-x^{(k)}\right)\)
其中\(-F^{\prime}(x)\)
就是
\[\nabla F\left(\mathbf{x}^{(k)}\right)=\frac{1}{\Delta t^{2}} \mathbf{M}\left(\mathbf{x}^{(k)}-\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\right)-\mathbf{f}\left(\mathbf{x}^{(k)}\right)
\]
而 b 就是 \(-F^{\prime}(x)\)
G 中包含了 \(\mathbf{f}\left(\mathbf{x}^{(k)}\right)\) 包含两个力, 一个是重力, 另一个是弹簧的弹力.
image
核心公式
code
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class implicit_model : MonoBehaviour
{
float t = 0.0333f;
float mass = 1;
float damping = 0.99f;
float rho = 0.995f;
float spring_k = 8000;
int[] E; // 边对
float[] L; // 初始边长度
Vector3[] V;
Vector3 g = new Vector3(0, -9.8f, 0); // 重力
float r = 2.7f; // 球半径
// Start is called before the first frame update
void Start()
{
Mesh mesh = GetComponent<MeshFilter> ().mesh;
//Resize the mesh.
int n=21;
Vector3[] X = new Vector3[n*n];
Vector2[] UV = new Vector2[n*n];
int[] triangles = new int[(n-1)*(n-1)*6];
for(int j=0; j<n; j++)
for(int i=0; i<n; i++)
{
X[j*n+i] =new Vector3(5-10.0f*i/(n-1), 0, 5-10.0f*j/(n-1));
UV[j*n+i]=new Vector3(i/(n-1.0f), j/(n-1.0f));
}
int t=0;
for(int j=0; j<n-1; j++)
for(int i=0; i<n-1; i++)
{
triangles[t*6+0]=j*n+i;
triangles[t*6+1]=j*n+i+1;
triangles[t*6+2]=(j+1)*n+i+1;
triangles[t*6+3]=j*n+i;
triangles[t*6+4]=(j+1)*n+i+1;
triangles[t*6+5]=(j+1)*n+i;
t++;
}
mesh.vertices=X;
mesh.triangles=triangles;
mesh.uv = UV;
mesh.RecalculateNormals ();
//Construct the original E
int[] _E = new int[triangles.Length*2];
for (int i=0; i<triangles.Length; i+=3)
{
_E[i*2+0]=triangles[i+0];
_E[i*2+1]=triangles[i+1];
_E[i*2+2]=triangles[i+1];
_E[i*2+3]=triangles[i+2];
_E[i*2+4]=triangles[i+2];
_E[i*2+5]=triangles[i+0];
}
//Reorder the original edge list
for (int i=0; i<_E.Length; i+=2)
if(_E[i] > _E[i + 1])
Swap(ref _E[i], ref _E[i+1]);
//Sort the original edge list using quicksort
Quick_Sort (ref _E, 0, _E.Length/2-1);
int e_number = 0;
for (int i=0; i<_E.Length; i+=2)
if (i == 0 || _E [i + 0] != _E [i - 2] || _E [i + 1] != _E [i - 1])
e_number++;
E = new int[e_number * 2];
for (int i=0, e=0; i<_E.Length; i+=2)
if (i == 0 || _E [i + 0] != _E [i - 2] || _E [i + 1] != _E [i - 1])
{
E[e*2+0]=_E [i + 0];
E[e*2+1]=_E [i + 1];
e++;
}
L = new float[E.Length/2];
for (int e=0; e<E.Length/2; e++)
{
int v0 = E[e*2+0];
int v1 = E[e*2+1];
L[e]=(X[v0]-X[v1]).magnitude;
}
V = new Vector3[X.Length];
for (int i=0; i<V.Length; i++)
V[i] = new Vector3 (0, 0, 0);
}
void Quick_Sort(ref int[] a, int l, int r)
{
int j;
if(l<r)
{
j=Quick_Sort_Partition(ref a, l, r);
Quick_Sort (ref a, l, j-1);
Quick_Sort (ref a, j+1, r);
}
}
int Quick_Sort_Partition(ref int[] a, int l, int r)
{
int pivot_0, pivot_1, i, j;
pivot_0 = a [l * 2 + 0];
pivot_1 = a [l * 2 + 1];
i = l;
j = r + 1;
while (true)
{
do ++i; while( i<=r && (a[i*2]<pivot_0 || a[i*2]==pivot_0 && a[i*2+1]<=pivot_1));
do --j; while( a[j*2]>pivot_0 || a[j*2]==pivot_0 && a[j*2+1]> pivot_1);
if(i>=j) break;
Swap(ref a[i*2], ref a[j*2]);
Swap(ref a[i*2+1], ref a[j*2+1]);
}
Swap (ref a [l * 2 + 0], ref a [j * 2 + 0]);
Swap (ref a [l * 2 + 1], ref a [j * 2 + 1]);
return j;
}
void Swap(ref int a, ref int b)
{
int temp = a;
a = b;
b = temp;
}
void Collision_Handling()
{
Mesh mesh = GetComponent<MeshFilter> ().mesh;
Vector3[] X = mesh.vertices;
GameObject sphere = GameObject.Find("Sphere");
Vector3 c = sphere.transform.position;
//Handle colllision.
for(int i=0; i<X.Length; i++)
{
if (i == 0 || i == 20) continue;
float distance = (X[i] - c).magnitude;
if(distance < r)
{
V[i] += (c + r * (X[i] - c) / distance - X[i]) / t; // 需要一个很小的速度偏差
X[i] = c + r * (X[i] - c) / distance;
}
}
mesh.vertices = X;
}
void Get_Gradient(Vector3[] X, Vector3[] X_hat, float t, Vector3[] G)
{
//Momentum and Gravity.
for(int i=0; i<G.Length; i++)
{
G[i] = (X[i] - X_hat[i]) * mass / (t * t) - g * mass;
}
//Spring Force.
for(int ei = 0; ei < L.Length; ei++)
{
int i = E[ei * 2];
int j = E[ei * 2 + 1];
Vector3 xi_xj = X[i] - X[j];
float len_ij = L[ei] / xi_xj.magnitude;
Vector3 t_ij = spring_k * (1 - len_ij) * xi_xj;
G[i] += t_ij;
G[j] -= t_ij;
}
}
// Update is called once per frame
void Update ()
{
Mesh mesh = GetComponent<MeshFilter> ().mesh;
Vector3[] X = mesh.vertices;
Vector3[] last_X = new Vector3[X.Length];
Vector3[] X_hat = new Vector3[X.Length]; // x = x + vi * t
Vector3[] G = new Vector3[X.Length]; // Gradient
Vector3[] delta_x = new Vector3[X.Length];
Vector3[] last_delta_x = new Vector3[X.Length];
Vector3[] old_delta_x = new Vector3[X.Length];
//Initial Setup.
// 更新速度和位置
for (int i=0; i<X.Length; i++)
{
V[i] *= damping;
last_X[i] = X[i];
X_hat[i] = X[i] + t * V[i];
X[i] = X_hat[i];
}
for(int k=0; k<32; k++)
{
Get_Gradient(X, X_hat, t, G);
// chebyshev
float A = mass / (t * t) + 4 * spring_k;
float w = 0;
float rest = 0;
for(int i = 0; i<X.Length; i++)
{
delta_x[i] = new Vector3(0, 0, 0);
last_delta_x[i] = new Vector3(0, 0, 0);
}
// iter
for(int j=0; j<50; j++)
{
rest = 0;
for(int i=0; i<X.Length; i++)
{
rest += (-delta_x[i] * A - G[i]).sqrMagnitude;
}
if(Mathf.Sqrt(rest) < 0.00001)
{
Debug.Log("Iter :" + k);
break;
}
if(j == 0)
{
w = 1;
}else if(j == 1)
{
w = 2 / (2 - rho * rho);
}else
{
w = 4 / (4 - rho * rho * w);
}
for(int i=0; i<X.Length; i++)
{
old_delta_x[i] = delta_x[i];
delta_x[i] += -delta_x[i] - G[i] / A;
delta_x[i] = w * delta_x[i] + (1 - w) * last_delta_x[i];
last_delta_x[i] = old_delta_x[i];
}
}
Debug.Log("rest:" + rest);
//Update X by gradient.
for(int i=0; i<X.Length; i++)
{
if (i == 0 || i == 20) continue;
X[i] += delta_x[i];
}
}
//Finishing.
for(int i=0; i<V.Length; i++)
{
V[i] += (X[i] - X_hat[i]) / t;
}
mesh.vertices = X;
Collision_Handling ();
mesh.RecalculateNormals ();
}
}
---------------------------我的天空里没有太阳，总是黑夜，但并不暗，因为有东西代替了太阳。虽然没有太阳那么明亮，但对我来说已经足够。凭借着这份光，我便能把黑夜当成白天。我从来就没有太阳，所以不怕失去。
--------《白夜行》
