简介

隐式积分法
显示积分简单而言是通过, 过去的求解未来. 而隐式积分, 简单而言是我要求解现在, 但是我的未知量中也有现在的未知量. 简单而言就是需要通过方程组的思想来进行求解.

参考文献

代码参考师弟 ~~

对于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 ().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 _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; epivot_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 ().mesh;
		Vector3[] X = mesh.vertices;
		GameObject sphere = GameObject.Find("Sphere");
		Vector3 c = sphere.transform.position;
		//Handle colllision.
		for(int i=0; i ().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