compliance

include/constraints.hpp

@@ -10,6 +10,7 @@ struct Edge {
 	uint32_t a;
 	uint32_t b;
 	float length;
+	float compliance;
 };
 
 struct Face {
@@ -20,6 +21,7 @@ struct Face {
 
 struct Triangle : Face {
 	float area;
+	float compliance;
 };
 
 struct Tetrahedron {
@@ -28,6 +30,7 @@ struct Tetrahedron {
 	uint32_t c;
 	uint32_t d;
 	float volume;
+	float compliance;
 };
 
 struct ConstraintData {

include/soft_body.hpp

@@ -20,12 +20,10 @@ class Mesh;
 
 class SoftBody {
 public:
-	explicit SoftBody(Mesh* mesh, float compliance);
+	explicit SoftBody(Mesh* mesh, float edgeCompliance, float triangleCompliance=0, float tetrahedronCompliance=0);
 	uint32_t firstIndex = 0;
 	int32_t vertexOffset = 0;
 
-	float compliance;
-
 	vector<Vertex> vertices;
 	vector<Face> faces;
 	ConstraintData constraintData;

shaders/pbd.comp

@@ -15,6 +15,7 @@ struct Edge {
 	uint a;
 	uint b;
 	float restLength;
+	float compliance;
 };
 
 struct Tetrahedron {
@@ -23,6 +24,7 @@ struct Tetrahedron {
 	uint c;
 	uint d;
 	float restVolume;
+	float compliance;
 };
 
 
@@ -62,8 +64,32 @@ layout (push_constant, std430) uniform PushConstants {
 
 void preSolve(uint vID){
 	vertices[vID].prevPosition = vertices[vID].position;
-	// vertices[vID].velocity += dt * gravity;
+	if (vertices[vID].w == 0){
+		return;
+	}
+	vertices[vID].velocity += dt * gravity;
 	vertices[vID].position += dt * vertices[vID].velocity;
+
+	float dist = vertices[vID].position.y + 5;
+	if (dist < 0){
+		vec3 normal = vec3(0, 1, 0);
+		vertices[vID].position -= dist * normal;
+
+		// contact friction
+		{
+			float mStatic = 1;
+			float mKinetic = 1;
+
+			vec3 a = vertices[vID].position - vertices[vID].prevPosition;
+			vec3 D = a - dot(a, normal) / dot(normal, normal) * normal;
+
+			vec3 add = -D;
+			if (length(D) >= mStatic * -dist){
+				add *= min(mKinetic * -dist / length(D), 1);
+			}
+			vertices[vID].position += add;
+		}
+	}
 }
 
 void solveEdge(uint eID){
@@ -75,7 +101,7 @@ void solveEdge(uint eID){
 	vec3 diff = v1.position - v2.position;
 	float currentLength = length(diff);
 
-	float alpha = 0.3 / dt / dt;
+	float alpha = edge.compliance / dt / dt;
 
 	float s = -(currentLength - edge.restLength) / (v1.w + v2.w + alpha);
 
@@ -104,26 +130,27 @@ void solveTetrahedron(uint tetID){
 
 	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;
+	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);
 
-    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;
 
-    if (w == 0) return;
+	float alpha = tetrahedron.compliance / dt / dt;
 
-    float alpha = 0 / dt / dt;
-    float s = -volumeError / (w + alpha);
+	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;
+	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;
@@ -132,6 +159,9 @@ void solveTetrahedron(uint tetID){
 }
 
 void postSolve(uint vID){
+	if (vertices[vID].w == 0){
+		return;
+	}
 	vertices[vID].velocity = (vertices[vID].position - vertices[vID].prevPosition) / dt;
 }
 
@@ -146,9 +176,8 @@ void main() {
 		case 1:
 			if (id < edgePartition.size){
 				solveEdge(id + edgePartition.offset);
-			}
-			if (id < tetrahedronPartition.size){
-				solveTetrahedron(id + tetrahedronPartition.offset);
+			} else if (id - edgePartition.size < tetrahedronPartition.size){
+				solveTetrahedron(id - edgePartition.size + tetrahedronPartition.offset);
 			}
 			break;
 		case 2:

src/constraints.cpp

@@ -11,13 +11,13 @@ void ConstraintData::writePartitionInformation() {
 }
 
 void DistanceConstraint::writeData(ConstraintData &dataLists) const {
-	dataLists.edges.push_back(Edge(a, b, length));
+	dataLists.edges.push_back(Edge(a, b, length, compliance));
 }
 
 void AreaConstraint::writeData(ConstraintData &dataLists) const {
-	dataLists.triangles.push_back(Triangle(Face(a, b, c), area));
+	dataLists.triangles.push_back(Triangle(Face(a, b, c), area, compliance));
 }
 
 void VolumeConstraint::writeData(ConstraintData &dataLists) const {
-	dataLists.tetrahedra.push_back(Tetrahedron(a, b, c, d, volume));
+	dataLists.tetrahedra.push_back(Tetrahedron(a, b, c, d, volume, compliance));
 }

src/simulation.cpp

@@ -48,21 +48,21 @@ Simulation::Simulation() {
 }
 
 void Simulation::createMeshBuffers() {
-	Mesh sphere("models/sphere_medium.ply");
+	Mesh sphere("models/icosphere.ply");
 	Mesh bunny("models/bunny_medium.ply");
 
-	auto body = std::make_unique<SoftBody>(&sphere, 0.3f);
+	auto body = std::make_unique<SoftBody>(&sphere, 1.f / 50);
 
-	for (size_t i = 0; i < 40; i++){
+	for (size_t i = 0; i < 10; i++){
 		auto copy = std::make_unique<SoftBody>(*body.get());
 		copy->applyVertexOffset({i * 2, 0, 0});
 		softBodies.push_back(std::move(copy));
 	}
 
-	body = std::make_unique<SoftBody>(&bunny, 0.3f);
-	for (size_t i = 0; i < 80; i++){
+	body = std::make_unique<SoftBody>(&bunny, 1.f / 3);
+	for (size_t i = 0; i < 5; i++){
 		auto copy = std::make_unique<SoftBody>(*body.get());
-		copy->applyVertexOffset({i * 2, 2, 0});
+		copy->applyVertexOffset({i * 2, 0, 2});
 		softBodies.push_back(std::move(copy));
 	}
 
@@ -209,7 +209,7 @@ void Simulation::recordComputeCommands(VkCommandBuffer cmdBuffer) {
 			auto edgePartition = constraintData.edgePartitions[partition];
 			auto tetrahedronPartition = constraintData.tetrahedronPartitions[partition];
 
-			if (edgePartition.size == 0 && tetrahedronPartition.size == 0)
+			if (edgePartition.size + tetrahedronPartition.size == 0)
 				continue;
 
 			ConstraintData::Partition partitions[2] = {edgePartition, tetrahedronPartition};
@@ -218,7 +218,7 @@ void Simulation::recordComputeCommands(VkCommandBuffer cmdBuffer) {
 							   offsetof(PBDPushData, edgePartition),
 							   sizeof(partitions), partitions);
 
-			uint32_t invocations = getGroupCount(std::max(edgePartition.size, tetrahedronPartition.size), BlOCK_SIZE);
+			uint32_t invocations = getGroupCount(edgePartition.size + tetrahedronPartition.size, BlOCK_SIZE);
 			vkCmdDispatch(cmdBuffer, invocations, 1, 1);
 
 			vkCmdPipelineBarrier(cmdBuffer, VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT, VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT, 0, 1, &barrier, 0, nullptr, 0, nullptr);

src/soft_body.cpp

@@ -7,9 +7,9 @@
 #include <omp.h>
 #include <iostream>
 
-SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
+SoftBody::SoftBody(Mesh* mesh, float edgeCompliance, float triangleCompliance, float tetrahedronCompliance) {
 	tetgenbehavior behavior;
-	behavior.parse_commandline(std::string("pYQfeza0.1").data());
+	behavior.parse_commandline(std::string("pYQfeza1").data());
 
 	tetgenio in;
 	in.numberofpoints = static_cast<int>(mesh->vertices.size());
@@ -67,7 +67,7 @@ SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
 		uint32_t c = out.trifacelist[i * 3 + 2];
 		if (out.trifacemarkerlist[i] != 0)
 			faces.emplace_back(Face(a, b, c));
-		constraintData.triangles.emplace_back(Triangle(Face(a, b, c), 0));
+		constraintData.triangles.emplace_back(Triangle(Face(a, b, c), 0, triangleCompliance));
 	}
 	faces.shrink_to_fit();
 
@@ -76,7 +76,7 @@ SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
 		uint32_t a = out.edgelist[i * 2 + 0];
 		uint32_t b = out.edgelist[i * 2 + 1];
 		float length = glm::length(vertices[a].position - vertices[b].position);
-		constraintData.edges.emplace_back(Edge(a, b, length));
+		constraintData.edges.emplace_back(Edge(a, b, length, edgeCompliance));
 	}
 
 	float density = 1;
@@ -104,7 +104,7 @@ SoftBody::SoftBody(Mesh* mesh, float compliance) : compliance(compliance) {
 		masses[c] += pointMass;
 		masses[d] += pointMass;
 
-		constraintData.tetrahedra.emplace_back(Tetrahedron(a, b, c, d, volume));
+		constraintData.tetrahedra.emplace_back(Tetrahedron(a, b, c, d, volume, tetrahedronCompliance));
 	}
 
 	for (size_t i = 0; i < masses.size(); i++){
@@ -120,14 +120,14 @@ void SoftBody::applyVertexOffset(const glm::vec3 &offset) {
 	}
 }
 
-vector<unordered_set<const Constraint *>> createPartitions(const Graph &graph){
+vector<unordered_set<const Constraint *>> createPartitions(const Graph &graph, const vector<const Constraint *> &ordering){
 	unordered_map<const Constraint *, uint32_t> constraintToColor;
 	vector<unordered_set<const Constraint *>> colorToConstraintList;
 
-	for (const auto& [constraint, adjacentList] : graph){
+	for (const Constraint * constraint : ordering){
 
 		unordered_set<uint32_t> forbiddenColors;
-		for (auto adjacent : adjacentList)
+		for (auto adjacent : graph.at(constraint))
 			if (constraintToColor.contains(adjacent))
 				forbiddenColors.insert(constraintToColor.at(adjacent));
 
@@ -154,6 +154,41 @@ vector<unordered_set<const Constraint *>> createPartitions(const Graph &graph){
 	return colorToConstraintList;
 }
 
+vector<const Constraint *> smallestLastOrdering(const Graph &graph) {
+	vector<const Constraint *> ordering(graph.size());
+
+	unordered_map<size_t, unordered_set<const Constraint *>> degreesToConstraints;
+	unordered_map<const Constraint *, size_t> constraintsToDegrees;
+	for (const auto &[constraint, adjacent] : graph) {
+		degreesToConstraints[adjacent.size()].insert(constraint);
+		constraintsToDegrees[constraint] = adjacent.size();
+	}
+
+
+	for (auto iter = ordering.rbegin(); iter != ordering.rend(); iter++) {
+		size_t minDegree = 0;
+		while (degreesToConstraints[minDegree].empty())
+			minDegree++;
+		const Constraint * minConstraint = *degreesToConstraints[minDegree].begin();
+		degreesToConstraints[minDegree].erase(minConstraint);
+		constraintsToDegrees.erase(minConstraint);
+
+		for (const auto &neighbor: graph.at(minConstraint)) {
+			if (constraintsToDegrees.contains(neighbor)) {
+				size_t oldDegree = constraintsToDegrees[neighbor];
+				size_t newDegree = oldDegree - 1;
+
+				constraintsToDegrees[neighbor] = newDegree;
+				degreesToConstraints[oldDegree].erase(neighbor);
+				degreesToConstraints[newDegree].insert(neighbor);
+			}
+		}
+
+		*iter = minConstraint;
+	}
+	return ordering;
+}
+
 void SoftBody::splitConstraints() {
 	omp_set_num_threads(omp_get_num_procs() - 2);
 
@@ -210,7 +245,9 @@ void SoftBody::splitConstraints() {
 		});
 	}
 
-	vector<unordered_set<const Constraint *>> partitions = createPartitions(graph);
+	auto ordering = smallestLastOrdering(graph);
+
+	vector<unordered_set<const Constraint *>> partitions = createPartitions(graph, ordering);
 	constraintData.partitionCount = partitions.size();
 
 	reorderConstraintIndices(partitions);
@@ -228,8 +265,4 @@ void SoftBody::reorderConstraintIndices(const vector<unordered_set<const Constra
 	std::swap(reordered.tetrahedronPartitions, constraintData.tetrahedronPartitions);
 	std::swap(reordered.edges, constraintData.edges);
 	std::swap(reordered.tetrahedra, constraintData.tetrahedra);
-	//std::swap(reordered, constraintData);
-	//reordered.triangles.swap(constraintData.triangles);
-	//reordered.tetrahedra.swap(constraintData.tetrahedra);
-	//std::swap(reordered.partitionCount, constraintData.partitionCount);
 }