计算机图形学入门Games103——约束，布料模拟(PBD) lab2

约束想解决的问题：

在现实世界中，很多的布料对于拉伸有很强的约束作用，但在模拟中如果将弹簧系数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();

    }


}