记录流体模拟的学习历程：（二）基于质点弹簧系统的布料模拟

上次做了质点弹簧系统，基于opengl的质点弹簧系统，这次在那个的基础上改了改，做了一个弹簧质点模型的布料模拟，在制作的图中出了很多问题，现在展现的结果也有很多需要改进的地方。

先展示一下成果。（我再问一下怎么制作没有水印的动图？）

视觉效果可能比较差，和真实的布料也会有些差距，主要还是先理解模拟过程和公式等。以后会继续改进。

依旧是使用opengl绘制点和线。

而且是2D视角的

下次会改进一下，使用三角形绘制并进行渲染，然后改成3D视角

下面来介绍一下基本的概念和公式方法。

把构成网格的顶点看作是质点，把顶点之间的连线看作是弹簧，弹簧可以分为 ３ 种：结构弹簧、剪切弹簧、弯曲弹簧。结构弹簧的作用是为了保持质点之间的初始距离，用于防止在经纬二个方向有过度变形；剪切弹簧是为了模拟织物倾斜方向的作用力。在这里我只加了前两种弹簧，并且点分布的比较均匀，没有像图c一样有一些偏移量。

接下来是受力分析

分析的对象是构成织物的所有质点，其中主要分为内力和外力，内力有弹簧的力和阻尼力，这个模拟中外力只添加了重力（下次改成3D版本会添加风力和缓冲力等）。

针对每个质点，我没要累加其周围八个质点带给他的力（中间的点所连接的八个点）

弹簧力根据胡克定律计算，公式如下

其中Lt是弹簧的初始长度，Pnv是当前研究的点，Pij是周围的八个点。

下面是阻尼力的计算公式

阻尼力是用来缓冲弹簧力较大时引起的过度变形，具体表现为织物的硬度。

外力暂时只研究了重力。

暂时就加这么多

由于本人没有可靠的参数（胡克定律系数，阻尼系数等），在模拟时出现了很多状况，最终效果也不是很好，所以在这里求一套可用的参数。。。

最后 上代码。

#include "ClothSimulation.h"
//这个可以不看，主要是用代码生成了各个点的位置和线的位置
void ClothSimulation::InitPositions()
	{
		int raw = 0;
		for (int i = 0; i < length; i++) {
			std::vector<glm::vec3> rawFactor;
			int colum = 0;
			for (int j = 0; j < length; j++)
			{
				if (i == length - 1 && j == length - 1){
					glm::vec3 position(76, 76, 0);
					colum += initLenth;
					rawFactor.push_back(position);
				}
				else if(i == length - 1 && j == 0)
				{
					glm::vec3 position(21, 76, 0);
					colum += initLenth;
					rawFactor.push_back(position);
				}
				else
				{
					glm::vec3 position(20 + colum, 20.0f + raw, 0);
					colum += initLenth;
					rawFactor.push_back(position);
				}
				
			}
			ParticlePositions.push_back(rawFactor);
			raw += initLenth;
		}
	}
//这个也可以不看，就是生成了所有的边
void ClothSimulation::InitEdges() {
	//for (size_t i = 0; i < length; i++) {
	//	for (int j = 0; j < length; j++) {
	//		//static_cast:将i的数据类型转化为float
	//		ParticlePositions[i][j].x = static_cast<float>(i);
	//		ParticlePositions[i][j].y = static_cast<float>(i);
	//	}
	//}
	for (int i = 0; i < length; i++)
	{
		for (int j = 0; j < length; j++) 
		{
			if (j + 1 == length || i + 1 == length)
				continue;
			//横线
			edges.push_back(ClothEdge{ ClothPoint {i, j},ClothPoint{i,j + 1} });
			edges.push_back(ClothEdge{ ClothPoint {j, i},ClothPoint{j+1, i} });
			edges.push_back(ClothEdge{ ClothPoint {i, j},ClothPoint{i + 1, j + 1} });
			//edges.push_back(ClothEdge{ ClothPoint {j + 1, i},ClothPoint{i, j + 1} });
		}
	}
	for (int i = 0; i < length - 1; i++) {
		//放入边界
		edges.push_back(ClothEdge{ ClothPoint {i, length - 1},ClothPoint{i + 1, length - 1} });
		edges.push_back(ClothEdge{ ClothPoint {length - 1, i},ClothPoint{length - 1, i + 1} });
	}
}
//初始化数组
void ClothSimulation::InitForcesArray() {
	forces.resize(length);
	velocities.resize(length);
	for (int i = 0; i < length; i++) {
		forces[i].resize(length);
		velocities[i].resize(length);
	}
}
//每帧计算每个质点收到的各种力
void ClothSimulation::CalculateForces() {
	//计算重力
	glm::vec3 gravity(0.0f, -2.0f, 0.0f);
	for (int i = 0; i < length; i++) {
		for (int j = 0; j < length; j++) {
			forces[i][j] = mass * gravity;
		}
	}
	for (int i = 0; i < length; i++) 
	{
		for (int j = 0; j < length; j++) 
		{
			//每个点要计算周围八个点带来的弹簧拉力
			for (int k = i - 1; k < i + 2; k++) 
			{
				for (int n = j - 1; n < j + 2; n++) 
				{
					//自身点，跳过
					if(k == i && n == j)
						continue;
					//边界点，跳过
					if (k < 0 || n < 0 || k > length - 1 || n > length - 1)
						continue;
					//计算弹簧力
					glm::vec3 pos0 = (ParticlePositions[i][j]);
					glm::vec3 pos1 = (ParticlePositions[k][n]);
					glm::vec3 r = pos0 - pos1;
				    float distance = glm::distance(pos0, pos1);
					glm::vec3 elasticForce(-stiffness * (distance - initLenth) * glm::normalize(r));
					forces[i][j] += elasticForce;
					//计算阻尼力
					glm::vec3 velocity1 = (velocities[i][j]);
					glm::vec3 velocity2 = (velocities[k][n]);
					glm::vec3 v = velocity1 - velocity2;
					glm::vec3 dampForce(-dampCoefficient * v);
					forces[i][j] += dampForce;
				}
			}
		}
	}
}
//更新质点的状态（速度、加速度、位置等）
void ClothSimulation::updateStates(float timeIntervalInSeconds) {
	for (int i = 0; i < length; i++)
	{
		for (int j = 0; j < length; j++)
		{
			glm::vec3 newAcceleration(forces[i][j] / mass);
			glm::vec3 newVelocity(velocities[i][j] +
				(timeIntervalInSeconds * newAcceleration));
			//速度到达一定程度让他停止
			if (glm::length(newVelocity) <= 0.3f)
				continue;
			//设置两个固定的点
			if (i == length - 1 && j == length - 1)
				continue;
			if (i == length - 1 && j == 0)
				continue;
			glm::vec3 pos0 = (ParticlePositions[i][j]);
			
			glm::vec3 newPosition(ParticlePositions[i][j] +
				(timeIntervalInSeconds * newVelocity));
			
			velocities[i][j] = newVelocity;
			ParticlePositions[i][j] = newPosition;
			
		}
	}
}

主要代码就这么多，绘制部分就是OpenGL画点和线，下次把3D模式的做出来，然后用三角形表示的时候再放绘制的代码。