其实这个去年就做好了，后面研究流体各种原因就咕咕咕了现在就来写一下大概思路，个人觉得比刚体模拟要稍微简单一些。该布料解算的模型是弹簧质点系统，物理计算是用隐式积分进行迭代，而碰撞检测和碰撞响应就和上一篇基于物理冲量模拟的刚体碰撞（图形学Games103）​

是一样的，不过简化了摩擦系数和弹性系数。

一，无序弹簧质点系统的构造：

     一个几何网格体，我们就把组成他的三角面作为主要研究对象。为了代码方面使用一个n*n的面片来做说明。

用2个数据结构去存它：一个float3列表去存顶点的坐标，这样我们就可以通过顶点id去获取顶点的坐标；一个int列表去存三角面的边的数据，每9个int为一个三角形，然后3个int为一条边（前2个int为边的两端的顶点id，第三个为三角形的id）。

然后我们这样得到的边的数据是有重合的，为了得到非重合的边的数据，需要把他去排序

排序完后，检查前后2个float3，如果其x和y2个元素一致，则说明共线重合。 我们把无重合的边收集起来，他就是我们需要的边的存在，存边的数据时，我们只记录边的两端的顶点id，一条边用2个int去记录。

然后我们还需求记录初始状态的每条边的边长，其上的代码如下。

count = n * n;       //21 *21的一个布料       int[] triangles = new int[(n-1)*(n-1)*6];//三角形       X=new Vector3[count];//函数外面申请内存地址       ResizeMesh(X,triangles,n);//调整模型的形状      //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];      }      //根据大小进行排序      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 原始的边的顶点坐标            E = new int[e_number * 2];//E 非重合的边的顶点id //边的数量为E.Length/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存边的长度      //X 为顶点的坐标      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;      }

为了更好处理边的数据，我定义了个结构体。

public struct Edge {          public float length;//该边的长度    public int p0, p1;//对应的顶点的id }

最终，我们需要在start函数里进行所有内容的初始化，包括之前没有提及的初始的速度，这里全部给初速度为0。

void Start()    {       count = n * n;       //21 *21的一个布料       int[] triangles = new int[(n-1)*(n-1)*6];//三角形       X=new Vector3[count];//函数外面申请内存地址       ResizeMesh(X,triangles,n);//调整模型的形状      //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];      }      //根据大小进行排序      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 原始的边的顶点坐标            E = new int[e_number * 2];//E 非重合的边的顶点id //边的数量为E.Length/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存边的长度      //X 为顶点的坐标      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;      } //初始化顶点速度为0      V = new float3[count];      for (int i=0; i<count; i++)         V[i] = new Vector3 (0, 0, 0);      InitEdges();          }

二：隐式积分

先计算梯度。

梯度的推导详情大家请看王老师的课程，王老师把非线性的数字求解问题变成了优化问题的思路是非常棒。其中的力分为弹簧力和重力2部分。

弹簧力f(x)=k(1-\frac{L}{||x_{i}-x_{j}||})(x_{i}-x_{j})

重力就很简单了f(x)=mg,其对应的代码如下：

void Get_Gradient(Vector3[] X, Vector3[] X_hat, float t, Vector3[] g) {     //Momentum and Gravity.   //逐顶点计算 // Vector3[] g = new Vector3[count];   for (int i = 0; i < count; i++)   {            Vector3 xi = X[i];      Vector3 xi_hat = X_hat[i];      g[i] = (xi - xi_hat) * mi / (dt * dt);//未考虑弹簧               }   //逐边计算   for (int k = 0; k < edges.Length; k++)   {      Edge edge = edges[k];      float Le = edge.length;      int i = edge.p0;      int j = edge.p1;      g[i] += spring_k * (1 - Le / (X[i] - X[j]).magnitude) * (X[i] - X[j]);      g[j] -= spring_k * (1 - Le / (X[i] - X[j]).magnitude) * (X[i] - X[j]);   }   //考虑重力   for (int i = 0; i < count; i++)   {      g[i]-=Vector3.down*9.8f*mi;   } }

三：迭代

     迭代的方法类似牛顿法，牛顿法也是我们常用的模拟方法之一，他是用一阶导数去近视逼近求解，它有点rayMatching的那味了，关于他的具体介绍请移步王老师的课。

     在这里，由于直接构造海森矩阵，或者在unity力使用线性解算器，都是一件麻烦的事情，故简化了计算方法。我们这里直接把海森矩阵当成对角阵进行处理的话就会简单很多。故我们的顶点更新就简化成了

步骤如下：

1.计算梯度。

2.更新顶点坐标。

3.反复迭代1-2步骤，直到有不错的效果，默认给的迭代32次。这个可以随着自己的效果喜好进行自定义调整。

4.更新速度。

代码如下：

   void Update ()   {            X     = _mesh.vertices;   //apply damping              Vector3[] X_hat    = new Vector3[count];      Vector3[] G       = new Vector3[count];//梯度      //初始化      for (int i = 0; i < count; i++)      {         V[i] *= damping;         X_hat[i] = dt * (Vector3)V[i]+ X[i];      }      for (int i = 0; i < count; i++)      {         X[i] =X_hat[i] ;      }       //迭代        for(int k=0; k<interate; k++)      {         Get_Gradient(X, X_hat, dt, G);                 //Update X by gradient.                 for (int i = 0; i < count; i++)         {            if (i == 0 || i == 20 || i == 7 || i==49) continue;//固定钉住                        X[i]  += 1/( mi/(dt*dt)+4*spring_k ) * -G[i];         }      }      //计算速度      for (int i = 0; i < count; i++)      {         V[i] +=(float3) (X[i] - X_hat[i]) / dt;      }                  _mesh.vertices = X;      Collision_Handling ();      _mesh.RecalculateNormals ();   }

下面是迭代2次和迭代32次的效果差异。

32次

32次迭代

2次

迭代2次

四：碰撞检测和处理

       碰撞的检测使用的是SDF，这个老生常谈的东西了。视频中我计算了3个模型的碰撞，最简单的球体，其次的胶囊体，还有稍微麻烦点的圆锥台。碰撞响应这里做得很简单粗暴，在检测到碰撞后，把进入sdf的顶点直接向SDF的梯度方向推出去即可，由于这里是并未把速度方向分解成垂直SDF和平行SDF方向，所以就没有摩擦力和弹性力。当然读者要加也简单，如同之前的刚体碰撞处理一样。其代码如下。

void Collision_Handling()   {      Mesh mesh = GetComponent<MeshFilter> ().mesh;      X = mesh.vertices; //处理球体      Vector3 c = sphere.position;            for (int i = 0; i < count; i++)      {         Vector3 xi = X[i];         float sdf = SDF.Sphere(xi,sphere.position,r);         if (sdf<=0)         {            V[i] +=  -sdf* SDF. Sphere_SDFGradient(xi , c)  / dt;            X[i] +=  -sdf* (Vector3)SDF. Sphere_SDFGradient(xi , c) ;         }        } //处理胶囊体      Vector3 pivot = _capsuleCollider.transform.position;      float scale = _capsuleCollider.transform.localScale.x;      float h = _capsuleCollider.height*scale*1.1F;      float radius = _capsuleCollider.radius*scale*1.1F;//乘1.1是为了避免穿插（因为是基于顶点的  还是要大一点比较好            for (int i = 0; i < count; i++)      {         Vector3 xi = X[i];         float sdf = SDF.Caapsule(xi,pivot,h,dir,_capsuleCollider.transform.rotation,radius);         if (sdf<=0)         {            V[i] += -sdf* SDF. Caapsule_SDFGradient(xi,pivot,h,dir,radius) / dt;            X[i] += -sdf* (Vector3)SDF. Caapsule_SDFGradient(xi,pivot,h,dir,radius) ;         }        } //处理圆锥台      for (int i = 0; i < count; i++)      {         Vector3 xi = X[i];         float sdf = SDF.CappedCone(xi, cone_a.position, cone_b.position, ra, rb);         bool isfa;         if (sdf<0)         {            V[i] += -sdf* SDF. CappedCone_SDFGradient(xi, cone_a.position, cone_b.position, ra, rb,out isfa) / dt;            X[i] += -sdf* (Vector3)SDF. CappedCone_SDFGradient(xi, cone_a.position, cone_b.position, ra, rb,out isfa) ;         }      }      mesh.vertices = X;   }

做完后，你也可以给它加个风场玩玩。我这里就做了一个简单的直线风场，只能控制大小和方向和爆发风力（可能不太明显哈）。做一些3d空间的扰动也是可以，这里只提供思路。

简易风场

五：总结

      最核心的处理，就已经完成了，其中最复杂的还是隐式积分这块，我只着重讲了作业的一些要点和思路。以上的内容也不保证完全正确，有疑惑欢迎楼下讨论。

附王老师的课的链接

https://games-cn.org/games103/

知乎版本有更好的代码体验

https://zhuanlan.zhihu.com/p/503495059