基于质点弹簧系统的布料模拟修改版

最新推荐文章于 2023-05-03 17:23:06 发布

-岁月如歌

最新推荐文章于 2023-05-03 17:23:06 发布

阅读量628

点赞数

分类专栏：流体仿真 opengl 文章标签： opengl 图形学

本文链接：https://blog.csdn.net/qq_40744606/article/details/109004802

版权

流体仿真同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

opengl

3 篇文章 0 订阅

订阅专栏

修改了上一篇博客的一些算法，质点的存储方式，弹簧的存储方式等，是一次大改。

不再直接画点和线，而是使用三角形进行绘制。

增加了简单的光照，渲染，一些最简单的碰撞检测。

添加了风力。

这次的效果比之前好了很多。

gif效果图（有点大所以不能上传）

首先是质点和弹簧的存储方式，不再用二维数组存储。而是用指针存储数组的方式。下面是具体实现代码。

上班时候养成了良好的注释习惯，尽管看代码就好。

首先是质点类和弹簧类。

class Particle {
public:
	//质量
	float mass = 0.01f;
	//是否固定
	bool fixed = false;
	//位置
	glm::vec3 position;
	//速度
	glm::vec3 velocity;
	//法向量
	glm::vec3 normal;
	//受到的力
	glm::vec3 forces;
};

class Spring {
public:
	//弹簧连接的第一个和第二个质点下标
	int firstParticle;
	int endParticle;
	//运动时的弹力
	float elasticForce;
	//胡克定律系数
	float stiffness = 15.0f;
	//初始长度
	float initLength;
};

然后对质点位置和弹簧关系进行初始化。

void initCloth() {
	//计算质点的数量
	particlNumber = gridSize * gridSize;

	//计算弹簧数量
	springNumber = (gridSize - 1) * gridSize * 2;
	springNumber += (gridSize - 1) * (gridSize - 1) * 2;
	springNumber += (gridSize - 2) * gridSize * 2;

	//创建数组空间
	particleArray = new Particle[particlNumber];
	particleArray2 = new Particle[particlNumber];
	springArray = new Spring[springNumber];
	//初始化质点位置
	for (int i = 0; i < gridSize; ++i)
	{
		for (int j = 0; j < gridSize; ++j)
		{
			
			particleArray[i * gridSize + j].position = glm::vec3
				(2.0f,
				float(j) - float(gridSize - 1) / 2,
				float(i) - float(gridSize - 1) / 2);
			particleArray[i * gridSize + j].velocity = glm::vec3(0.0f, 0.0f, 0.0f);
			particleArray[i * gridSize + j].forces = glm::vec3(0.0f, 0.0f, 0.0f);
			particleArray[i * gridSize + j].normal = glm::vec3(0.0f, 0.0f, 0.0f);
		}
	}
	particleArray[gridSize - 1].fixed = true;
	particleArray[gridSize * gridSize - 1].fixed = true;
	//Copy the balls into the other array
	for (int i = 0; i < particlNumber; ++i)
		particleArray2[i] = particleArray[i];
	//Set the currentParticle and nextParticle pointers
	currentParticle = particleArray;
	nextParticle = particleArray2;


	//Initialise the springArray
	Spring* currentSpring = &springArray[0];

	//创造从左到右的弹簧
	for (int i = 0; i < gridSize; ++i)
	{
		for (int j = 0; j < gridSize - 1; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = i * gridSize + j + 1;
			currentSpring->initLength = initLength;
			++currentSpring;
		}
	}

	//创造从上到下的弹簧
	for (int i = 0; i < gridSize - 1; ++i)
	{
		for (int j = 0; j < gridSize; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = (i + 1) * gridSize + j;
			currentSpring->initLength = initLength;
			++currentSpring;
		}
	}

	//斜向下的弹簧
	for (int i = 0; i < gridSize - 1; ++i)
	{
		for (int j = 0; j < gridSize - 1; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = (i + 1) * gridSize + j + 1;
            //斜着的弹簧长度是横竖弹簧的1.414倍
			currentSpring->initLength = initLength * sqrt(2.0f);
			++currentSpring;
		}
	}

	//斜向上的弹簧
	for (int i = 0; i < gridSize - 1; ++i)
	{
		for (int j = 1; j < gridSize; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = (i + 1) * gridSize + j - 1;
			currentSpring->initLength = initLength * sqrt(2.0f);

			++currentSpring;
		}
	}

	//隔一个的弹簧
	for (int i = 0; i < gridSize; ++i)
	{
		for (int j = 0; j < gridSize - 2; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = i * gridSize + j + 2;
			currentSpring->initLength = initLength * 2;

			++currentSpring;
		}
	}

	//隔一个的弹簧
	for (int i = 0; i < gridSize - 2; ++i)
	{
		for (int j = 0; j < gridSize; ++j)
		{
			currentSpring->firstParticle = i * gridSize + j;
			currentSpring->endParticle = (i + 2) * gridSize + j;
			currentSpring->initLength = initLength * 2;

			++currentSpring;
		}
	}
}

代码中所需要的变量都放在这里了。


//球体半径
float sphereRadius = 4.0f;
//背景颜色
COLOR backgroundColor(0.75f, 0.75f, 0.85f, 1.0f);
//网格中每行列的质点数量
const int gridSize = 15;
//阻尼系数
float dampFactor = 1.8f;
//弹簧初始长度
float initLength = 1.0f;
//质点数量
int particlNumber;
//布料拖拽系数
float dragCoefficient = 0.15f;
//质点存储数组
Particle* particleArray = NULL;
Particle* particleArray2 = NULL;

Particle* currentParticle = NULL;
Particle* nextParticle = NULL;
//重力
glm::vec3 gravity(0.0f, -0.98f, 0.0f);
//风力
glm::vec3 wind(0.7f, 0.0f, 0.0f);
//空气密度
float airDensity = 1.2f;
//弹簧数量
int springNumber;
//存储弹簧的数组
Spring* springArray = NULL;

接下来是各个三角形的绘制。

glBegin(GL_TRIANGLES);
	{
		for (int i = 0; i < gridSize - 1; ++i)
		{
			for (int j = 0; j < gridSize - 1; ++j)
			{
				//绘制两种三角形（合并成一个正方形）
				
				glNormal3fv(glm::value_ptr(currentParticle[i * gridSize + j].normal));
				glVertex3fv(glm::value_ptr(currentParticle[i * gridSize + j].position));
				glNormal3fv(glm::value_ptr(currentParticle[i * gridSize + j + 1].normal));
				glVertex3fv(glm::value_ptr(currentParticle[i * gridSize + j + 1].position));
				glNormal3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j].normal));
				glVertex3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j].position));

				glNormal3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j].normal));
				glVertex3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j].position));
				glNormal3fv(glm::value_ptr(currentParticle[i * gridSize + j + 1].normal));
				glVertex3fv(glm::value_ptr(currentParticle[i * gridSize + j + 1].position));
				glNormal3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j + 1].normal));
				glVertex3fv(glm::value_ptr(currentParticle[(i + 1) * gridSize + j + 1].position));
			}
		}
	}
	glEnd();
	glDisable(GL_LIGHTING);

这样就绘制除了初始的布料。

然后用glu库绘制一个球体。球体的中心在坐标原点，这样进行碰撞处理就变得十分简单。

//使用glu库绘制球体。这种绘制一定要在while中进行  不然绘制不出来
void drawSphere() {
    //绘制中心点在原点的球体。
    static GLUquadricObj* sphere = gluNewQuadric();
    glEnable(GL_LIGHTING);
    /*glMaterialfv函数的使用
    指定用于光照计算的当前材质属性。参数face的取值可以是GL_FRONT、GL_BACK或GL_FRONT_AND_BACK，指出材质属性将应用于物体的哪面。
    最后的参数是一个指向数组的指针
    第一个参数        默认值           作用
    GL_AMBIENT（0.2，0.2，0.2，1.0）材质的环境颜色
    GL_DIFFUSE（0.8，0.8，0.8，1.0）材质的散射颜色
    GL_AMBIENT_AND_DIFFUSE 材质的环境颜色和散射颜色
    GL_SPECULAR（0.0，0.0，0.0，1.0）材质的镜面反射颜色
    GL_SHININESS 0.0 镜面反射指数  指数越大，高光部分越小
    GL_EMISSION（0.0，0.0，0.1，1.0）材质的发射光颜色
    GL_COLOR_INDEXES（0， 1， 1） 环境颜色索引、散射颜色索引和镜面反射颜色索引*/
    glMaterialfv(GL_FRONT, GL_AMBIENT, COLOR(1.0f, 0.0f, 0.0f, 0.0f));
    glMaterialfv(GL_FRONT, GL_DIFFUSE, COLOR(1.0f, 0.0f, 0.0f, 0.0f));
    glMaterialfv(GL_FRONT, GL_SPECULAR, blue);
    glMaterialf(GL_FRONT, GL_SHININESS, 32.0f);
    //开启或禁用多边形正面或者背面上的光照、阴影和颜色计算及操作，消除不必要的渲染计算。
    glEnable(GL_CULL_FACE);
    //三四个参数指定围绕z轴的细分数（类似于经度线）。  不是很懂
    gluSphere(sphere, sphereRadius, 48, 24);

	glDisable(GL_LIGHTING);
	glDisable(GL_CULL_FACE);
	glMaterialfv(GL_FRONT, GL_SPECULAR, black);
}

以上是基本的准备工作，接下来才是这个系统的核心，就是力的计算和速度、位置的处理。

具体公式等，参考的是《一种基于OpenGL顶点处理机制的织物交互仿真方法》，公式大多参考这个文献上写出的公式。具体公式就看上一个博客就好。
在这里插入图片描述

下边先看风力的计算代码，风力存储在了质点上。

	for (int i = 0; i < gridSize - 1; i++) {
		for (int j = 0; j < gridSize - 1; j++) {
			//求出速度平均值
			glm::vec3 aveVelocity = 
			(currentParticle[i * gridSize + j].velocity +
			currentParticle[i * gridSize + j + 1].velocity +
			currentParticle[(i + 1) * gridSize + j].velocity) / 3.0f;
			//计算相对速度
			glm::vec3 relVelocity = aveVelocity - wind;
			//计算三角形面积
			float a = glm::length(currentParticle[i * gridSize + j].position - currentParticle[i * gridSize + j + 1].position);
			float b = glm::length(currentParticle[i * gridSize + j].position - currentParticle[(i + 1) * gridSize + j].position);
			float c = glm::length(currentParticle[i * gridSize + j + 1].position - currentParticle[(i + 1) * gridSize + j].position);
			float p = (a + b + c) * (a + b - c) * (a + c - b) * (b + c - a)/4;
			float s = sqrt(p);
			//求解风力方程。
			currentParticle[i].forces 
				= currentParticle[i * gridSize + j + 1].forces 
				= currentParticle[(i + 1) * gridSize + j].forces
				= -0.5f * airDensity * glm::abs(relVelocity) * dragCoefficient * s
				* (relVelocity * currentParticle[i].normal) * currentParticle[i].normal;
		}
	}

然后是弹力的计算，弹力赋值到了弹簧的类上。

//计算弹力
	for (int i = 0; i < springNumber; ++i)
	{
		//运动过程中弹簧的长度
		float springLength = glm::length(currentParticle[springArray[i].firstParticle].position -
			currentParticle[springArray[i].endParticle].position);
		//长度的变化量
		float extension = springLength - springArray[i].initLength;
		//胡克定律 F = k * （L - L0）   不过后边为什么除以初始长度？
		//将计算的弹力赋值到弹簧上
		springArray[i].elasticForce = springArray[i].stiffness * extension / springArray[i].initLength;
	}

下边写的有点乱，计算弹力和阻尼力和计算新的位置、速度等代码放在一起了。

//将当前的质点信息赋值给下一时刻的质点（挺绕的，看代码理解。）
	for (int i = 0; i < particlNumber; ++i)
	{
		//若果质点固定，则不更新位置，速度设置为0  
		//不固定，则计算速度和位置
		if (currentParticle[i].fixed)
		{
			nextParticle[i].position = currentParticle[i].position;
			nextParticle[i].velocity = glm::vec3(0.0f, 0.0f, 0.0f);
		}
		else
		{
			//计算质点受到的所有力。
			glm::vec3 force = gravity + particleArray[i].forces;
			/*遍历所有弹簧，如果当前遍历的弹簧的第一个端点是当前for循环的质点，则通过
			当前质点的下标算出相邻质点的下标，并求出质点之间的力的作用。
			至于为什么要算如下两次，请搞清楚创建弹簧时的具体信息
			*/
			for (int j = 0; j < springNumber; ++j)
			{
				//当弹簧的起点是当前点
				if (springArray[j].firstParticle == i)
				{
					//弹力方向是从当前点指向下一个点，并且对向量进行单位化处理
					glm::vec3 tensionDirection = glm::normalize(currentParticle[springArray[j].endParticle].position -
						currentParticle[i].position);
					//质点所受到的力    加上弹力
					force += springArray[j].elasticForce * tensionDirection;
					//计算阻尼力  f = 阻尼系数 * 速度差
					glm::vec3 FVD = currentParticle[springArray[j].endParticle].velocity -
						currentParticle[i].velocity;
					force += dragCoefficient * FVD;
				}
				//当弹簧的终点是当前点。操作和上一个类似
				if (springArray[j].endParticle == i)
				{
					glm::vec3 tensionDirection = glm::normalize(currentParticle[springArray[j].firstParticle].position -
						currentParticle[i].position);
					force += springArray[j].elasticForce * tensionDirection;
					//计算阻尼力  f = 阻尼系数 * 速度差
					glm::vec3 FVD = currentParticle[springArray[j].endParticle].velocity -
						currentParticle[i].velocity;
					force += dragCoefficient * FVD;
				}
			}

			
			//计算加速度
			glm::vec3 acceleration = force / currentParticle[i].mass;
			//更新质点速度
			nextParticle[i].velocity = currentParticle[i].velocity + acceleration * (float)deltaTime;

			//计算新的位置
			nextParticle[i].position = currentParticle[i].position +
				(nextParticle[i].velocity + currentParticle[i].velocity) * (float)deltaTime / 2.0f;

			/* 球体的碰撞检测
			球体的中心店在坐标圆点，所以只需要检测质点到坐标原点的距离是否小于圆的半径
			为了避免穿模，所以半径乘以一定值。
			*/
			float squareLength = nextParticle[i].position.x * nextParticle[i].position.x
				+ nextParticle[i].position.y * nextParticle[i].position.y
				+ nextParticle[i].position.z * nextParticle[i].position.z;
			if (squareLength < sphereRadius * 1.08f * sphereRadius * 1.08f)
				nextParticle[i].position = glm::normalize(nextParticle[i].position)* sphereRadius * 1.08f;

			//检测质点是否和地板发生碰撞。
			if (nextParticle[i].position.y < -8.5f)
				nextParticle[i].position.y = -8.5f;
		}
	}
	currentParticle = nextParticle;

基本就这些，剩下的不重要的代码就不放在这里了。

布料的仿真暂时告一段落。

明天开始研究SPH算法，研究流体的仿真，希望一个月内可以成功进行流体的仿真。

编辑于 4 分钟前