soft tetrahedrons

include/vulkan/vertex.hpp

@@ -9,7 +9,7 @@ struct Vertex {
 	alignas(16) glm::vec3 color;
 	alignas(16) glm::vec3 normal;
 	alignas(16) glm::vec3 velocity;
-	glm::vec3 prevPosition;
+	alignas(16) glm::vec3 prevPosition;
 	float inverseMass;
 
 	static VkVertexInputBindingDescription getBindingDescription();

shaders/pbd.comp

@@ -17,7 +17,13 @@ struct Edge {
 	float restLength;
 };
 
-
+struct Tetrahedron {
+	uint a;
+	uint b;
+	uint c;
+	uint d;
+	float restVolume;
+};
 
 
 layout (std430, set = 0, binding = 0) buffer VertexBuffer {
@@ -28,6 +34,10 @@ layout (std430, set = 0, binding = 2) buffer EdgeBuffer {
 	Edge edges[];
 };
 
+layout (std430, set = 0, binding = 4) buffer TetrahedronBuffer {
+	Tetrahedron tetrahedra[];
+};
+
 layout (std140, set = 0, binding = 5) uniform Sizes {
 	uint vertexCount;
 	uint faceCount;
@@ -57,30 +67,68 @@ void preSolve(uint vID){
 }
 
 void solveEdge(uint eID){
-	Edge e = edges[eID];
+	Edge edge = edges[eID];
 	
-	Vertex v1 = vertices[e.a];
-	Vertex v2 = vertices[e.b];
+	Vertex v1 = vertices[edge.a];
+	Vertex v2 = vertices[edge.b];
 
 	vec3 diff = v1.position - v2.position;
 	float currentLength = length(diff);
 
-	float alpha = 0.0001 / (dt * dt);
+	float alpha = 0.3 / dt / dt;
 
-	float s = (currentLength - e.restLength) / (1 + 1 + alpha);
+	float s = -(currentLength - edge.restLength) / (v1.w + v2.w + alpha);
 
 	vec3 g1 = normalize(diff);
 	vec3 g2 = -g1;
 
-	vec3 delta1 = -s * 1 * g1;
-	vec3 delta2 = -s * 1 * g2;
+	vec3 delta1 = s * v1.w * g1;
+	vec3 delta2 = s * v2.w * g2;
 
-	vertices[e.a].position += delta1;
-	vertices[e.b].position += delta2;
+	vertices[edge.a].position += delta1;
+	vertices[edge.b].position += delta2;
 }
 
 void solveTetrahedron(uint tetID){
+	Tetrahedron tetrahedron = tetrahedra[tetID];
+
+	Vertex va = vertices[tetrahedron.a];
+	Vertex vb = vertices[tetrahedron.b];
+	Vertex vc = vertices[tetrahedron.c];
+	Vertex vd = vertices[tetrahedron.d];
+
+	vec3 a = va.position;
+	vec3 b = vb.position;
+	vec3 c = vc.position;
+	vec3 d = vd.position;
+
+	float volumeError = dot(d - a, cross(b - a, c - a)) / 6 - tetrahedron.restVolume;
+
+    vec3 ga = cross(d - b, c - b) / 6;
+    vec3 gb = cross(c - a, d - a) / 6;
+    vec3 gc = cross(d - a, b - a) / 6;
+    vec3 gd = cross(b - a, c - a) / 6;
+
+    float w =
+    va.w * dot(ga, ga) +
+    vb.w * dot(gb, gb) +
+    vc.w * dot(gc, gc) +
+    vd.w * dot(gd, gd);
+
+    if (w == 0) return;
+
+    float alpha = 0 / dt / dt;
+    float s = -volumeError / (w + alpha);
+
+    vec3 addA = s * ga * va.w;
+    vec3 addB = s * gb * vb.w;
+    vec3 addC = s * gc * vc.w;
+    vec3 addD = s * gd * vd.w;
 
+	vertices[tetrahedron.a].position += addA;
+	vertices[tetrahedron.b].position += addB;
+	vertices[tetrahedron.c].position += addC;
+	vertices[tetrahedron.d].position += addD;
 }
 
 void postSolve(uint vID){

src/simulation.cpp

@@ -23,7 +23,7 @@ Simulation::Simulation() {
 
 	properties.gravity = {0, -9.81, 0};
 	properties.k = 10;
-	properties.dt = 1.f / 144.f / static_cast<float>(properties.k);
+	properties.dt = 1.f / 60.f / static_cast<float>(properties.k);
 
 	propertiesBuffer = make_unique<Buffer>(
 			sizeof(Properties), &properties, sizeof(properties),
@@ -148,8 +148,6 @@ void Simulation::createMeshBuffers() {
 				: Buffer(size, data, size, VK_BUFFER_USAGE_STORAGE_BUFFER_BIT | additionalUsageFlags, VMA_MEMORY_USAGE_AUTO_PREFER_DEVICE, 0) {}
 	};
 
-	vertices[0].position += glm::vec3(0, 0.2, 0);
-
 	vertexBuffer = make_unique<SimulationBuffer>(vertices.data(), vertices.size() * sizeof(Vertex), VK_BUFFER_USAGE_VERTEX_BUFFER_BIT);
 	faceBuffer = make_unique<SimulationBuffer>(faces.data(), faces.size() * sizeof(Face), VK_BUFFER_USAGE_INDEX_BUFFER_BIT);
 	edgeBuffer = make_unique<SimulationBuffer>(constraintData.edges.data(), constraintData.edges.size() * sizeof(Edge));

src/soft_body.cpp

@@ -9,7 +9,7 @@
 
 SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
 	tetgenbehavior behavior;
-	behavior.parse_commandline(std::string("pYa0.001Qfez").data());
+	behavior.parse_commandline(std::string("pYQfeza0.1").data());
 
 	tetgenio in;
 	in.numberofpoints = static_cast<int>(mesh->vertices.size());
@@ -79,13 +79,36 @@ SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
 		constraintData.edges.emplace_back(Edge(a, b, length));
 	}
 
+	float density = 1;
+	vector<float> masses(vertices.size());
+
 	constraintData.tetrahedra.reserve(out.numberoftetrahedra);
 	for (size_t i = 0; i < out.numberoftetrahedra; i++){
 		uint32_t a = out.tetrahedronlist[i * 4 + 0];
 		uint32_t b = out.tetrahedronlist[i * 4 + 1];
 		uint32_t c = out.tetrahedronlist[i * 4 + 2];
 		uint32_t d = out.tetrahedronlist[i * 4 + 3];
-		constraintData.tetrahedra.emplace_back(Tetrahedron(a, b, c, d, 0));
+
+		glm::vec3 p1 = vertices[a].position;
+		glm::vec3 p2 = vertices[b].position;
+		glm::vec3 p3 = vertices[c].position;
+		glm::vec3 p4 = vertices[d].position;
+
+		float volume = glm::dot(p4 - p1, glm::cross(p2 - p1, p3 - p1)) / 6;
+		assert(volume > 0);
+
+		float pointMass = volume * density / 4;
+
+		masses[a] += pointMass;
+		masses[b] += pointMass;
+		masses[c] += pointMass;
+		masses[d] += pointMass;
+
+		constraintData.tetrahedra.emplace_back(Tetrahedron(a, b, c, d, volume));
+	}
+
+	for (size_t i = 0; i < masses.size(); i++){
+		vertices[i].inverseMass = 1.0f / masses[i];
 	}
 
 	splitConstraints();