Cloth simulation in C++ using OpenGL

A cloth simulation implemented with a compute shader, based on the Position Based Dynamics algorithm.

Position Based dynamics

It was introduced by Matthias Müller in 2006. Instead of working with forces and velocities, it works directly with positions and the constraints applied to them. It is popular in the game industry because it is simple, robust, and efficient. It is worth noting that this is not a physically accurate method, but it produces a realistic result. Over several iterations, a set of constraints is applied, and a weight parameter moves the positions to satisfy them.
I use four constraints to generate the final result:

• Gravity constraint that pulls each point towards the gravity vector.
• Collision constraint that pushes each point away from colliders.
• Distance constraint that moves each point to keep a fixed distance with neighbors.
• Bending constraint that prevents critical angles between points and their neighbors.

This is for example the distance constraint:

Distance constraint for cloth

I find the distance between p1 and p2 vectors, and also the vector pointing from p1 to p2 (or from p2 to p1, depending which we want to fix) and push it that way depending on the distance and the weight parameters.
|p1 - p2| - d = 0 equation should be satisfied eventually.

Compute shaders

A compute shader is a shader stage that is not tied to rendering; it can be used more broadly (e.g. simulations, particle systems). Here I use it for the cloth simulation. The number of work groups can be defined in up to three dimensions, and any of them can be 1 to reduce the dimensionality. In this example I generate an N x N grid of points, where N is the resolution, then bind that data to the compute shader using 2D workgroups.
I am going to demonstrate an example compute shader usage in C++ OpenGL using the distance constraint.
Firstly, I load the compute shader and compile it using the shader library:

vertexDistanceShader.loadShader(GL_COMPUTE_SHADER, "../shaders/distance.comp");
vertexDistanceShader.compile();

Secondly, I generate our buffers I will use for the compute shaders:

glGenBuffers(1, &positionBuffer);
glBindBuffer(GL_SHADER_STORAGE_BUFFER, positionBuffer);
glBufferData(GL_SHADER_STORAGE_BUFFER, vNum * sizeof(xyzw), NULL, GL_STATIC_DRAW);
xyzw* pos = (xyzw*)glMapBufferRange(GL_SHADER_STORAGE_BUFFER, 0, vNum * sizeof(xyzw), GL_MAP_WRITE_BIT | GL_MAP_INVALIDATE_BUFFER_BIT);
for (unsigned int i = 0; i < N; ++i)
	for (unsigned int j = 0; j < N; ++j)
	{
		pos[j*N+i].x = (double)i / ((double)N - 1.0) - 0.5;
		pos[j*N+i].y = 0;
		pos[j*N+i].z = (double)j / ((double)N - 1.0) - 0.5;
		pos[j*N+i].w = 1.0f;
	}
glUnmapBuffer(GL_SHADER_STORAGE_BUFFER);

I also initialize them with positions on a plane, here I generate the grid.
I now bind the points for the compute shader and each buffer for their appropriate location:

vertexDistanceShader.enable();
glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, positionBuffer);
glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, velocityBuffer);
glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, positionBufferTmp);
vertexDistanceShader.setUniform1f("d", d);
glDispatchCompute(N / Nwg, N / Nwg, 1);
vertexDistanceShader.disable();

Note that the initial d distance parameter is set as well. Now I dispatch the buffer for the compute shader.
A piece of the core of the compute shader looks like this:

for (int x = -1; x <= 1; ++x){
	for (int y = -1; y <= 1; ++y){
		if (gl_GlobalInvocationID.y + y < 0 || gl_GlobalInvocationID.y + y > 2 * gl_WorkGroupSize.y - 1 ||
			gl_GlobalInvocationID.x + x < 0 || gl_GlobalInvocationID.x + x > 2 * gl_WorkGroupSize.x - 1)
			continue;
		uint _gid = (gl_GlobalInvocationID.y + y) * gl_NumWorkGroups.x * gl_WorkGroupSize.x + (gl_GlobalInvocationID.x + x);		
		if (_gid == gid) continue;
		float dis = d;
		if (x != 0 && y != 0) dis *= sqrt(2);
		vec3 dir = positionT[_gid].xyz - positionT[gid].xyz;
		dp += dir / length(dir) * (length(dir) - dis) * 0.5 * fix;
	}
}

Every neighbor is collected around each point and the vector to move it to the correct position is calculated, then eventually added to the position.