Cloth simulation in C++ using OpenGL

Cloth simulation created with compute shader based on position based dynamics algorithm.

Position Based dynamics

This method was introduced by Matthias Müller in 2006. The basic idea is instead of calculating with forces and velocities, it directly calculates with positions and constraints attached to them. This is a popular technique used in game industry, because it is simple, robust and efficient. However, it must be kept in mind this is not a physically correct method, but it produces a realistic result. Repeated a number of times, in each iteration a number of constraints are given, and by a weight parameter, the positions are arranged to satisfy these constraints.
I use four constraints to generate the final result:

• Gravity constraint that pulls each point to 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

Compute shader is a shader stage that is not necessarily used for rendering, it can be used in a wider range (e.g. simulations, particle systems). In the current example I am using it for cloth simulation. The number of work groups that execute the computation can be defined with three numbers (three dimensions) and any of these can be 1, reducing the dimension. In the current example I generate an N * N point grid, where N is a constant that describes the resolution and then these data are binded to the compute shader with 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.