Added Example 40: Cloth simulation using compute shader

Example40/.gitignore

@@ -0,0 +1 @@
+/x86__Windows__MinGW_Debug/

Example40/shader/cloth.comp.glsl

@@ -1,30 +1,232 @@
 #version 430 core
 
-layout(std430, binding=1) buffer VertexIn {
-    vec4 b_vertexIn[];
+// Number of iterations the spring constraints and collisions are updated.
+#define ITERATIONS 4
+
+// Unit: kg * m * s*s
+#define FORCE vec4(0.0, -0.2, 0.0, 0.0)
+
+// Unit: kg
+#define MASS 0.5
+
+// Offset, that a particle is not inside the sphere surface.
+#define RADIUS_TOLERANCE 0.01
+
+uniform int u_verticesPerRow;
+
+uniform float u_deltaTime;
+
+// Distance at rest between vertical and horizontal particle.
+uniform float u_distanceRest;
+
+// One coliding sphere.
+uniform vec4 u_sphereCenter;
+uniform float u_sphereRadius;
+
+layout(std430, binding=1) buffer VertexInPrevious {
+    vec4 b_vertexInPrevious[];
+};
+
+layout(std430, binding=2) buffer VertexInCurrent {
+    vec4 b_vertexInCurrent[];
 };
 
-layout(std430, binding=2) buffer VertexOut {
+layout(std430, binding=3) buffer VertexOut {
     vec4 b_vertexOut[];
 };
 
-layout(std430, binding=3) buffer NormalOut {
+layout(std430, binding=4) buffer NormalOut {
     vec3 b_normalOut[];
     // Padding[]
 };
 
-// If ROWS and COLUMS change in the main file, update local size as well.
-// (ROWS+1)*(COLUMNS+1)
-layout (local_size_x = 121) in;
+// (ROWS+1)^2
+layout (local_size_x = 1024) in;
+
+void calculateDeltaVector(out vec4 deltaVector, vec4 a, vec4 b)
+{
+	vec4 vector = a - b;
+
+	float deltaLength = length(vector);
+
+	// 0.5 because this happens on each side.
+	deltaVector = vector * 0.5 * (u_distanceRest - deltaLength) / deltaLength;
+}
 
+// see http://web.archive.org/web/20070610223835/http://www.teknikus.dk/tj/gdc2001.htm
+//     http://developer.download.nvidia.com/SDK/10/direct3d/Source/Cloth/doc/Cloth.pdf
 void main(void)
 {
-    uint i = gl_GlobalInvocationID.x;
+    int currentIndex = int(gl_GlobalInvocationID.x) + int(gl_GlobalInvocationID.y) * u_verticesPerRow;
+
+    int currentIndexRow = currentIndex % u_verticesPerRow;
+    int currentIndexColumn = currentIndex / u_verticesPerRow;
+
+	//
+	// Summing all forces, to get the acceleration. In this case, only one force directing downwards is used to simulate gravity.
+	// 
+
+	vec4 a = FORCE / MASS;
+
+	//
+	// Verlet integration.
+	//
+
+	float squareDeltaTime = u_deltaTime * u_deltaTime;
+
+	vec4 x = b_vertexInCurrent[currentIndex];
+
+	vec4 xPrev = b_vertexInPrevious[currentIndex];
 
-	// TODO Do cloth simulation, as seen:
-	// http://web.archive.org/web/20070610223835/http://www.teknikus.dk/tj/gdc2001.htm
-	// http://developer.download.nvidia.com/SDK/10/direct3d/Source/Cloth/doc/Cloth.pdf
+    b_vertexOut[currentIndex] = 2.0 * x - xPrev + a * squareDeltaTime;
+
+    memoryBarrierBuffer();
+
+	//
+	// Relaxation iterations.
+	//
+	for (int iteration = 0; iteration < ITERATIONS; iteration++)
+	{
+	    //
+	    // Process spring constraints.
+	    //
+
+		// Alternating "even" and "odd" cells to avoid dependencies.  
+		for (int evenOdd = 0; evenOdd < 2; evenOdd++)
+		{
+			int rightIndex;
+			int bellowIndex;
+			int diagonalIndex;
+
+			vec4 deltaVector;
+
+			if (currentIndexRow % 2 == evenOdd && currentIndexColumn % 2 == evenOdd)
+			{    
+				rightIndex = currentIndex + 1;
+				if (rightIndex % u_verticesPerRow > currentIndexRow)
+				{
+					calculateDeltaVector(deltaVector, b_vertexOut[currentIndex], b_vertexOut[rightIndex]);
+
+					b_vertexOut[currentIndex] += deltaVector;
+					b_vertexOut[rightIndex] -= deltaVector;
+
+					memoryBarrierBuffer();  
+				}    			
+
+
+				bellowIndex = currentIndex + u_verticesPerRow;
+				if (bellowIndex < u_verticesPerRow * u_verticesPerRow)			
+				{
+					calculateDeltaVector(deltaVector, b_vertexOut[currentIndex], b_vertexOut[bellowIndex]);
+
+					b_vertexOut[currentIndex] += deltaVector;
+					b_vertexOut[bellowIndex] -= deltaVector;
+
+					memoryBarrierBuffer();  
+				}
+
+
+				diagonalIndex = currentIndex + 1 + u_verticesPerRow;
+				if (diagonalIndex < u_verticesPerRow * u_verticesPerRow && diagonalIndex % u_verticesPerRow > currentIndexRow)
+				{
+					calculateDeltaVector(deltaVector, b_vertexOut[rightIndex], b_vertexOut[diagonalIndex]);
+
+					b_vertexOut[rightIndex] += deltaVector;
+					b_vertexOut[diagonalIndex] -= deltaVector;  
+
+					memoryBarrierBuffer();
+
+					calculateDeltaVector(deltaVector, b_vertexOut[bellowIndex], b_vertexOut[diagonalIndex]);
+
+					b_vertexOut[bellowIndex] += deltaVector;
+					b_vertexOut[diagonalIndex] -= deltaVector;
+
+					memoryBarrierBuffer();  
+				}
+			}
+			else if (evenOdd == 1 && currentIndexRow % 2 == 1 && currentIndexColumn == 0)
+			{
+				// First row spring, when processing odd pass. 
+
+				rightIndex = currentIndex + 1;
+				if (rightIndex < u_verticesPerRow)
+				{
+					calculateDeltaVector(deltaVector, b_vertexOut[currentIndex], b_vertexOut[rightIndex]);
+
+					b_vertexOut[currentIndex] += deltaVector;
+					b_vertexOut[rightIndex] -= deltaVector;
+
+					memoryBarrierBuffer();  
+				}
+			}
+
+			barrier();
+		}
+
+		//
+		// Process collision constraints.
+		//
+
+		vec4 sphereVector = b_vertexOut[currentIndex] - u_sphereCenter;
+
+		// If particle is inside sphere ...
+		if (length(sphereVector) < u_sphereRadius + RADIUS_TOLERANCE)
+		{
+			// ... move it outside the sphere.
+			b_vertexOut[currentIndex] = u_sphereCenter + normalize(sphereVector) * (u_sphereRadius + RADIUS_TOLERANCE); 
+
+			memoryBarrierBuffer();
+		}
+
+		barrier();		
+	}	
+
+    //
+    // Calculate normal.
+    //
+
+    vec3 normal = vec3(0.0, 0.0, 0.0);
+
+	vec3 tangent;
+	vec3 bitangent;
 
-    b_vertexOut[i] = b_vertexIn[i];
-    b_normalOut[i] = vec3(0.0, 1.0, 0.0);
+	// Taking all neighbour particles, if available, into account.
+    if (currentIndexRow + 1 < u_verticesPerRow)
+    {
+    	tangent = normalize((b_vertexOut[currentIndex + 1] - b_vertexOut[currentIndex]).xyz);
+
+	    if (currentIndexColumn + 1 < u_verticesPerRow)
+	    {
+	    	bitangent = normalize((b_vertexOut[currentIndex + u_verticesPerRow] - b_vertexOut[currentIndex]).xyz);
+
+	    	normal += normalize(cross(bitangent, tangent));
+	    }    	
+	    if (currentIndexColumn - 1 >= 0)
+	    {
+	    	bitangent = normalize((b_vertexOut[currentIndex] - b_vertexOut[currentIndex - u_verticesPerRow]).xyz);
+
+	    	normal += normalize(cross(bitangent, tangent));
+	    }    	
+    }
+    if (currentIndexRow - 1 >= 0)
+    {
+    	tangent = normalize((b_vertexOut[currentIndex] - b_vertexOut[currentIndex - 1]).xyz);
+
+	    if (currentIndexColumn + 1 < u_verticesPerRow)
+	    {
+	    	bitangent = normalize((b_vertexOut[currentIndex + u_verticesPerRow] - b_vertexOut[currentIndex]).xyz);
+
+	    	normal += normalize(cross(bitangent, tangent));
+	    }    	
+	    if (currentIndexColumn - 1 >= 0)
+	    {
+	    	bitangent = normalize((b_vertexOut[currentIndex] - b_vertexOut[currentIndex - u_verticesPerRow]).xyz);
+
+	    	normal += normalize(cross(bitangent, tangent));
+	    }    	
+    }
+
+    //
+
+    b_normalOut[currentIndex] = normalize(normal);
 }

Example40/shader/cloth.frag.glsl

@@ -1,10 +1,14 @@
 #version 430 core
 
-in vec4 v_color;
+uniform vec3 u_lightDirection;
+uniform vec4 u_color;
+
+in vec3 v_normal;
 
 out vec4 fragColor;
 
 void main(void)
 {
-	fragColor = v_color;
+	// Two sided rendering, so use absolute value.
+	fragColor = u_color * abs(dot(normalize(v_normal), u_lightDirection));
 }

Example40/shader/cloth.vert.glsl

@@ -1,26 +1,16 @@
 #version 430 core
 
-layout(std430, binding=1) buffer Vertex {
-    vec4 b_vertex[];
-};
-
-layout(std430, binding=2) buffer Normal {
-    vec3 b_normal[];
-    // Padding[]
-};
-
 uniform mat4 u_modelViewProjectionMatrix;
 uniform mat3 u_normalMatrix;
-uniform vec3 u_lightDirection;
-uniform vec4 u_color;
 
-out vec4 v_color;
+in vec4 a_vertex;
+in vec4 a_normal;
+
+out vec3 v_normal;
 
 void main(void)
 {
-	vec3 normal = u_normalMatrix * b_normal[gl_VertexID];
-
- 	v_color = u_color * max(dot(normal, u_lightDirection), 0.0);
+	v_normal = u_normalMatrix * a_normal.xyz;
 
-	gl_Position = u_modelViewProjectionMatrix * b_vertex[gl_VertexID];
+	gl_Position = u_modelViewProjectionMatrix * a_vertex;
 }

Example40/shader/sphere.frag.glsl

@@ -0,0 +1,13 @@
+#version 430 core
+
+uniform vec3 u_lightDirection;
+uniform vec4 u_color;
+
+in vec3 v_normal;
+
+out vec4 fragColor;
+
+void main(void)
+{
+	fragColor = u_color * max(dot(normalize(v_normal), u_lightDirection), 0.0);
+}

Example40/shader/sphere.vert.glsl

@@ -0,0 +1,16 @@
+#version 430 core
+
+uniform mat4 u_modelViewProjectionMatrix;
+uniform mat3 u_normalMatrix;
+
+in vec4 a_vertex;
+in vec3 a_normal;
+
+out vec3 v_normal;
+
+void main(void)
+{
+	v_normal = u_normalMatrix * a_normal;
+
+	gl_Position = u_modelViewProjectionMatrix * a_vertex;
+}