欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。
本文模拟的是布料与球碰撞
基于物理的弹簧模拟
布料三角化&内部应力
把布料看成是三角形组成的平面，各三角形的点可独立活动。
在相邻的点之间加装弹簧来模拟布料的弹性。
每个点会受到重力，以及弹簧的弹力。
弹簧力的计算
模拟过程
每步行计算出点的受力，根据受力计算速度，从而计算出位移。
以上是显示积分过程。该算法会产生overshooting的缺陷。
隐式积分&转为优化问题
左式利用作用后的力进行计算，问题是作用后的力不好求。先转化为右式。
x
[
1
]
中
的
v
[
1
]
代
换
就
可
以
得
到
右
上
边
式
子
。
x^{[1]} 中的v^{[1]}代换就可以得到右上边式子。
x[1]中的v[1]代换就可以得到右上边式子。
利用导数为0时取得最值，可以得到优化F(x)是等价于求解上述方程的。
从而转化成了一个优化问题。
以上就是布料模拟的大致原理，具体实现是利用了预测移动后的x来得到力的大小。
算法实现
关键步骤
正常过程：
更新力
x
−
是
预
测
后
的
位
移
，
f
(
x
)
是
弹
簧
的
力
和
重
力
\overset {-} {x} 是预测后的位移，f(x)是弹簧的力和重力
x−是预测后的位移，f(x)是弹簧的力和重力
根据力计算出新的位移
碰撞过程：
计算各点与球的距离来判断是否有碰撞。
如果有碰撞则计算出点的位移与速度。
代码
代码库：https://github.com/LightningBilly/ACMAlgorithms/tree/master/图形学算法/模拟算法/布料模拟/
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;
private bool m_EnableChebyshevAcceleration = true;
// 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]; // 441
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));
// X[j*n+i] =new Vector3(5-10.0f*i/(n-1), -10.0f*j/(n-1),0);
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;
//Handle colllision.
Vector3 c = GameObject.Find("Sphere").transform.position;
// 碰撞检测，并将碰撞的点沿球心与点的方向移动
for (int i = 0; i < V.Length; i++)
{
float distance = (X[i] - c).magnitude;
if (distance < 2.7f)
{
Vector3 direction = (X[i] - c) / distance;
V[i] += 1.0f / t * (c + 2.7f * direction - X[i]);
X[i] = c + 2.7f * direction;
}
}
mesh.vertices = X;
}
void Get_Gradient(Vector3[] X, Vector3[] X_hat, float t, Vector3[] G)
{
//Momentum and Gravity.
Vector3 gravity = new Vector3(0, -9.8f, 0);
Vector3 wind = Vector3.zero;
if (Time.frameCount % 300 < 100 + 100 * Random.value)
{
wind = new Vector3(0, 0, -1) * Random.value * 0.5f;
}
// 根据点的移动，模拟力的大小，包括力与重力。
for (int i = 0; i < G.Length; i++)
{
G[i] = (1.0f / t / t) * mass * (X[i] - X_hat[i]) - (gravity + wind) * mass;
}
// Spring Force.
// 模拟弹簧的力
for(int e = 0; e < E.Length / 2; e++)
{
int i = E[e * 2];
int j = E[e * 2 + 1];
Vector3 springForce = spring_k * (1 - L[e] / (X[i] - X[j]).magnitude) * (X[i] - X[j]);
G[i] += springForce;
G[j] -= springForce;
}
}
// 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];
Vector3[] G = new Vector3[X.Length];
//Initial Setup.
for (int i = 0; i < V.Length; ++i)
{
V[i] *= damping;
X_hat[i] = X[i] + t * V[i];
// X[i] = X_hat[i];
}
float w = 0;
for(int k=0; k<32; k++)
{
Get_Gradient(X, X_hat, t, G);
//Update X by gradient.
if(m_EnableChebyshevAcceleration)
{
if (k == 0) w = 1;
else if(k == 1) w = 2 / (2 - rho * rho);
else w = 4 / (4 - rho * rho * w);
for (int i = 0; i < X.Length; i++)
{
if (i == 0 || i == 20) continue;
Vector3 xOld = X[i];
X[i] -= 1.0f / ((1.0f / t / t) * mass + 4.0f * spring_k) * G[i];
X[i] = w * X[i] + (1 - w) * last_X[i];
last_X[i] = xOld;
}
}
else
{
// 根据力的大小更新位移
for (int i = 0; i < X.Length; i++)
{
if (i == 0 || i == 20) continue;
X[i] -= 1.0f / ((1.0f / t / t) * mass + 4.0f * spring_k) * G[i];
}
}
}
// Finishing.
// 更新速度
for (int i = 0; i < V.Length; i++)
{
if (i == 0 || i == 20) continue;
V[i] += 1.0f / t * (X[i] - X_hat[i]);
}
mesh.vertices = X;
Collision_Handling ();
mesh.RecalculateNormals ();
}
}
算法效果
本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。
1515
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
