约束想解决的问题:
在现实世界中,很多的布料对于拉伸有很强的约束作用,但在模拟中如果将弹簧系数K调的很大,会造成问题
1.显式积分中,模型会不稳定
2.隐式积分中,会造成ill_conditioned
问题可以解决,比如说更小的时间步长或者更多次的迭代,但这些解决方法都会带来很大的计算量,花更多的时间
因此使用约束是想要在弹簧系数很大的时候,如何减少模拟所花费的计算量和时间
约束1.拉伸约束
对于一根弹簧可以有弹簧的长度等于原长,可以得到一个约束
因此我们可以构造一个函数,使得弹簧长度始终为原长,拉伸状态下,需要减少长度,压缩转态下,需要扩大长度,可以得到两个顶点新的位置
这个问题可以是一个投影问题,若是在高纬空间里面,我们需要做的就是保证x的位置在合理区域里面,挪动的距离还要尽可能少
有
具体推导函数来自论文,作为结论使用(这个推导较为容易,质心不变为前提,质量和挪动的量相关,质量一般在模拟的时候都设置为1,若是固定点,就设置为无穷大(当然可以直接不更新))
多弹簧(高斯—赛德尔方法),就是一根一根的处理这些弹簧,每处理一根,另外的也会变化,也就是需要不断迭代运算(方法很好理解,但名字很高大上。。。)
这个方法虽然容易,但是变形是和边的顺序是有关的,会有偏向性而且会影响收敛的速度,而且不管怎么做,总会有不收敛的边虽然这个方法叫做高斯-赛德尔方法,但事实上更像机器学习的梯度下降,只不过有方向
多弹簧(雅各布方法 jacobi)
减少边的顺序造成的影响,可以进行一些并行的计算,并行计算比高斯法更容易一些
但我每次计算更新时,先保存下来,有很多的更新值,我最后做一个类似取平均值的加权计算,得到一个总的更新,收敛的速度上会慢一些
基于上面的方法,可以做PBD模拟
第一步:先做简单的自由活动的位置更新,更新速度和位置
第二步:开始施加约束,
1.对顶点位置做约束
2.根据新位置,我们可以对速度做新的计算
3.最后,将原来的位置覆盖掉
这个方法的问题在于,这个模拟过程中只注重了视觉的表达,但实际上并未涉及到什么物理的方法,并没有什么弹性系数等等东西
在这种模拟方法中,弹性效果的体现仅仅在于通过迭代数量和网格来展现,迭代数量越少,网格数越多(因为需要更多次数的迭代才会收敛),也就越有弹性
其余的约束还有体积约束,能量约束等
缺点是没有什么物理含义,而且在高分辨率网格上效率下降很快,实时低精度模拟套路
Muller. 2008. Hierarchical Position Based Dynamics. VRIPHYS.
为了使得PBD具有物理含义,我们可以这样做,前面的位置速度更新,我们使用正常的物理方法,只是使用后面的约束提高其稳定性
此时,后面的约束只是用来保证模拟的稳定,解决抖动问题
例如说对于一阶弹簧,我们可以定义拉伸比例
我希望靠近满足比例,意思就是将约束放宽的PBD
三角形面积约束:
也可以约束三角形面积的大小,实际操作起来的话计算一个缩放比就可以了
有稳定性的问题(大形变),模拟非线性的表现,可以在比较小的阶段使用物理模拟,在形变比较大的阶段使用约束限制,,解决锁死问题
PBD代码中的拉伸约束实现
void Strain_Limiting()
{
Mesh mesh = GetComponent<MeshFilter>().mesh;
Vector3[] vertices = mesh.vertices;
Vector3[] sum_x = new Vector3[vertices.Length];
int[] sum_n = new int[vertices.Length];
//Apply PBD here.
for (int i = 0; i < vertices.Length; i++)
{
sum_x[i] = new Vector3(0, 0, 0);
sum_n[i] = 0;
}
for (int e = 0; e < L.Length; e++)
{
int i = E[e * 2];
int j = E[e * 2 + 1];
Vector3 xij = vertices[i] - vertices[j];
sum_x[i] = sum_x[i] + 0.5f * (vertices[i] + vertices[j] + L[e] * xij * (1.0f / xij.magnitude));
sum_x[j] = sum_x[j] + 0.5f * (vertices[i] + vertices[j] - L[e] * xij * (1.0f / xij.magnitude));
sum_n[i]++;
sum_n[j]++;
}
for (int i = 0; i < vertices.Length; i++)
{
if (i == 0 || i == 20) continue;
V[i] = V[i] + (1.0f / t) * ((0.2f * vertices[i] + sum_x[i]) / (0.2f + (float)sum_n[i]) - vertices[i]);
vertices[i] = (0.2f * vertices[i] + sum_x[i]) / (0.2f + (float)sum_n[i]);
}
mesh.vertices = vertices;
}
代码逻辑为
比较直观,根据伪代码一步一步执行即可
完整代码如下所示
public class PBD_model : MonoBehaviour
{
float t = 0.0333f;
float damping = 0.99f;
int[] E;
float[] L;
Vector3[] V;
Vector3 gravity = new Vector3(0, -9.8f, 0);
// Use this for initialization
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[] T = 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++)
{
T[t * 6 + 0] = j * n + i;
T[t * 6 + 1] = j * n + i + 1;
T[t * 6 + 2] = (j + 1) * n + i + 1;
T[t * 6 + 3] = j * n + i;
T[t * 6 + 4] = (j + 1) * n + i + 1;
T[t * 6 + 5] = (j + 1) * n + i;
t++;
}
mesh.vertices = X;
mesh.triangles = T;
mesh.uv = UV;
mesh.RecalculateNormals();
//Construct the original edge list
int[] _E = new int[T.Length * 2];
for (int i = 0; i < T.Length; i += 3)
{
_E[i * 2 + 0] = T[i + 0];
_E[i * 2 + 1] = T[i + 1];
_E[i * 2 + 2] = T[i + 1];
_E[i * 2 + 3] = T[i + 2];
_E[i * 2 + 4] = T[i + 2];
_E[i * 2 + 5] = T[i + 0];
}
//Reorder the original edge list
//Guarantee smaller index in front of bigger one
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 i = E[e * 2 + 0];
int j = E[e * 2 + 1];
L[e] = (X[i] - X[j]).magnitude;
}
V = new Vector3[X.Length];
for (int i = 0; i < X.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 Strain_Limiting()
{
Mesh mesh = GetComponent<MeshFilter>().mesh;
Vector3[] vertices = mesh.vertices;
Vector3[] sum_x = new Vector3[vertices.Length];
int[] sum_n = new int[vertices.Length];
//Apply PBD here.
for (int i = 0; i < vertices.Length; i++)
{
sum_x[i] = new Vector3(0, 0, 0);
sum_n[i] = 0;
}
for (int e = 0; e < L.Length; e++)
{
int i = E[e * 2];
int j = E[e * 2 + 1];
Vector3 xij = vertices[i] - vertices[j];
sum_x[i] = sum_x[i] + 0.5f * (vertices[i] + vertices[j] + L[e] * xij * (1.0f / xij.magnitude));
sum_x[j] = sum_x[j] + 0.5f * (vertices[i] + vertices[j] - L[e] * xij * (1.0f / xij.magnitude));
sum_n[i]++;
sum_n[j]++;
}
for (int i = 0; i < vertices.Length; i++)
{
if (i == 0 || i == 20) continue;
V[i] = V[i] + (1.0f / t) * ((0.2f * vertices[i] + sum_x[i]) / (0.2f + (float)sum_n[i]) - vertices[i]);
vertices[i] = (0.2f * vertices[i] + sum_x[i]) / (0.2f + (float)sum_n[i]);
}
mesh.vertices = vertices;
}
void Collision_Handling()
{
Mesh mesh = GetComponent<MeshFilter>().mesh;
Vector3[] X = mesh.vertices;
//For every vertex, detect collision and apply impulse if needed.
float radius = 2.7f;
GameObject sphere = GameObject.Find("Sphere");
Vector3 center = sphere.transform.position;
for (int i = 0; i < X.Length; i++)
{
if (i == 0 || i == 20)
continue;
Vector3 d = X[i] - center;
if (d.magnitude < radius)
{
Vector3 A = center + radius * d.normalized;
V[i] = V[i] + (A - X[i]) / t;
X[i] = A;
}
}
mesh.vertices = X;
}
// Update is called once per frame
void Update()
{
Mesh mesh = GetComponent<MeshFilter>().mesh;
Vector3[] X = mesh.vertices;
for (int i = 0; i < X.Length; i++)
{
if (i == 0 || i == 20) continue;
//Initial Setup
V[i] = V[i] * damping;
V[i] = V[i] + gravity * t;
X[i] = X[i] + V[i] * t;
}
mesh.vertices = X;
for (int l = 0; l < 32; l++)
Strain_Limiting();
Collision_Handling();
mesh.RecalculateNormals();
}
}