games103作业2

PBD方法

首先是每个质点的力的分析，不考虑碰撞和弹簧弹力的情况下，每个质点受重力的影响，所以需要对每个质点进行速度和位置的重力影响更新。

float 		t= 0.0333f;
float		damping= 0.99f;
int[] 		E;
float[] 	L;
Vector3[] 	V;
Vector3 gravity = new Vector3(0.0f, -9.8f, 0.0f);
......
for(int i=0; i<X.Length; i++)
{
	if(i==0 || i==20)	continue;
	V[i] = V[i] + gravity * t;
	V[i] *= damping;
	X[i] = X[i] + V[i] * t;
	//Initial Setup
	//...
}

思考一个问题，现实生活中布料的每个质点在拉扯变大以后，会越来越难以拉扯，基于胡可定律的弹簧模型中需要增大弹性系数k来模拟这种现象，但这会造成显式积分和隐式积分都出现问题，增大了模拟计算量。基于约束的方法被提出的动机就是想要解决这个问题。也就是PBD

使用Jacobi的方式对质点进行约束位置更新。然后通过PBD的算法流程对位置和速度进行更新。

	void Strain_Limiting()
	{
		Mesh mesh = GetComponent<MeshFilter> ().mesh;
		Vector3[] vertices = mesh.vertices;
		Vector3[] vertices_new = new Vector3[vertices.Length];
		int[] n = new int[vertices.Length];
		for(int i = 0; i < vertices.Length; i++)
		{
			vertices_new[i] = new Vector3(0.0f, 0.0f, 0.0f);
			n[i] = 0;
		}
		
		for(int e = 0; e < L.Length; e++)//注意是消重的
		{
			int a = E[e * 2 + 0];
			int b = E[e * 2 + 1];
			Vector3 a_b = vertices[a] - vertices[b];
			float halfDistance = (a_b.magnitude - L[e])*0.5f;
			Vector3 pointMove = halfDistance * a_b.normalized;

			vertices_new[a] = vertices_new[a] + vertices[a] - pointMove;
			vertices_new[b] = vertices_new[b] + vertices[b] + pointMove;
			n[a]++;
			n[b]++;
			
		}
		for(int i = 0; i < vertices.Length; i++)
		{
			if (i == 0 || i == 20) continue;
			V[i] = V[i] + ((vertices_new[i] + 0.2f * vertices[i]) / (n[i] + 0.2f) - vertices[i]) / t;
			vertices[i] = (vertices_new[i] + 0.2f * vertices[i]) / (n[i] + 0.2f);
		}
		//Apply PBD here.
		//...
		mesh.vertices = vertices;
	}

布料效果

球的撞击

在之前的课程中，求的是刚体对碰撞体进行撞击，所以最后要进行约束回来，但是这里不需要，这是流体。

这里计算碰撞位移是在PBD以后，才计算是否发生碰撞以及碰撞后的速度和位移变换。（这些都算在一个帧内进行更新）

	void Collision_Handling()
	{
		Mesh mesh = GetComponent<MeshFilter> ().mesh;
		Vector3[] X = mesh.vertices;
		Vector3 spherePosition = sphere.transform.position;
		for(int i = 0; i < X.Length; i++)
		{
			if (i == 0 || i == 20) continue;
			if ((X[i] - spherePosition).magnitude > r)
			{
				continue;
			}

			//发生碰撞，得到碰撞点
			Vector3 collosionPoint = r * (X[i] - spherePosition).normalized + spherePosition;
			Vector3 normal = (collosionPoint - spherePosition).normalized;
			float jud = Vector3.Dot(V[i], normal);
			/*V[i] = V[i]+ (collosionPoint - X[i]) / t;
            X[i] = collosionPoint;*/
			
			Vector3 v_N = jud * normal;
			Vector3 v_T = V[i] - v_N;
			//作业这里的意思是碰撞以后，位移到球体表面
			v_N = v_N + (collosionPoint - X[i]) / t;
			X[i] = collosionPoint;
			V[i] = v_N + v_T;//忽略摩擦
		}
		//For every vertex, detect collision and apply impulse if needed.
		//...
		mesh.vertices = X;
	}

效果

隐式积分的方式

求解X的步骤为一个优化问题，不断的优化X得到最终的解。在这个过程中实际上就是F(x)为最小的时候X的值为多少。当F(x)的Hessian矩阵正定的时候，我们知道，F(x)存在唯一的最小值。

对于上述式子的求解方式为

作业中的求解方式

首先我们需要对F(x)求解gradient和Hessian矩阵。对于gradient矩阵作业给出的方式是忽略速度的影响的。于是我们可以根据重力和弹簧的作用力对每个质点求解gradient。

求解gradient

void Get_Gradient(Vector3[] X, Vector3[] X_hat, float t, Vector3[] G)
{
	//Momentum and Gravity.
	float inv_t = 1.0f / (t * t);
	for (int i = 0; i < X.Length; i++)
	{
		G[i] = inv_t * mass * (X[i] - X_hat[i]) - mass * gravity;
	}
	//Spring Force.
	for(int e = 0; e < L.Length; e++)
	{
		int i = E[e * 2 + 0];
		int j = E[e * 2 + 1];
		Vector3 transform = (X[i] - X[j]) - L[e] * (X[i] - X[j]).normalized;
		G[i] += spring_k * transform;
		G[i] -= spring_k * transform;
	}
}

求解质点的新位置

在求解的时候实际是F(x)的Hessian矩阵的逆矩阵右乘F(x)的gradient矩阵得到，然后对进行更新。作业中对于Hessian矩阵，给定了一个magic matrix(一个对角矩阵)。所以只需要对gradient矩阵每个值都乘一个相同的数就可以了，在每次迭代的时候对X的值进行更新

作业中不管是否为较小值，硬性迭代32次

for(int k=0; k<32; k++)
{
	Get_Gradient(X, X_hat, t, ref G);
	for(int i = 0; i < X.Length; i++)
	{
		if (i == 0 || i == 20) continue;
		float hessian = 1.0f / (mass * (1.0f / (t * t)) + 4 * spring_k);
		X[i] -= hessian * G[i];
	}
	//Update X by gradient.
	
}
for(int i = 0; i < V.Length; i++)
{
	V[i] = (X[i] - X_hat[i]) / t;
}

此时效果会非常慢。

我们可以利用下面的切比雪夫版Jacobi方式进行加速

不等式中，在每次对进行值的更新以后判断残差值是否在合理的范围内。也就是优化的思想，换个思路，这里我们迭代更新的是X，是否也可以使用类似的思想对X进行更新，让他趋近收敛。（因为这里的w就是让方程更加趋近收敛）

for(int k=0; k<32; k++)
{
	Get_Gradient(X, X_hat, t, ref G);
	jud = true;
    if (k == 0) w = 1.0f;
    else if (k == 1) w = 2.0f / (2.0f - rho * rho);
	else
	{
        w = 4.0f / (4.0f - rho * rho * w);
    }
    for (int i = 0; i < X.Length; i++)
	{
		if (i == 0 || i == 20) continue;
		float hessian = 1.0f / (mass * (1.0f / (t * t)) + 4 * spring_k);
		Vector3 old_X = X[i];
		X[i] -= hessian * G[i];
		X[i] = w * X[i] + (1 - w) * last_X[i];
		last_X[i] = old_X;
		if ((X[i]-old_X).magnitude > 0.0001f)
		{
			jud = false;
		}
		
	}
	//Update X by gradient.
	if (jud) break;
}

撞击

同PBD中介绍，就不多阐述了。

void Collision_Handling()
{
	Mesh mesh = GetComponent<MeshFilter> ().mesh;
	Vector3[] X = mesh.vertices;
	Vector3 spherePosition = sphere.transform.position;
	for(int i = 0; i < X.Length; i++)
	{
		float d = (X[i] - spherePosition).magnitude;
		if (d > r) continue;
		V[i] = V[i] + (spherePosition + r * (X[i] - spherePosition).normalized - X[i])/t;
		X[i] = spherePosition + r * (X[i] - spherePosition).normalized;
	}
	//Handle colllision.

	mesh.vertices = X;
}

最终效果