From 50c6820cf9182d66fb7eb8e794138b0e7ad27c84 Mon Sep 17 00:00:00 2001
From: Tobias Frisch <tfrisch@uni-koblenz.de>
Date: Sat, 5 Feb 2022 17:05:15 +0100
Subject: [PATCH] Added types of simulations

Signed-off-by: Tobias Frisch <tfrisch@uni-koblenz.de>
---
 projects/wobble_bobble/shaders/grid.vert      |  16 ++-
 projects/wobble_bobble/shaders/particle.inc   |  37 ++++++
 .../shaders/transform_particles_to_grid.comp  |   6 +-
 .../shaders/update_particle_velocities.comp   |  36 ++++--
 projects/wobble_bobble/src/main.cpp           | 116 +++++++++++++-----
 5 files changed, 160 insertions(+), 51 deletions(-)

diff --git a/projects/wobble_bobble/shaders/grid.vert b/projects/wobble_bobble/shaders/grid.vert
index 3636d87b..54de3f3e 100644
--- a/projects/wobble_bobble/shaders/grid.vert
+++ b/projects/wobble_bobble/shaders/grid.vert
@@ -6,6 +6,10 @@
 layout(set=0, binding=0) uniform texture3D gridImage;
 layout(set=0, binding=1) uniform sampler gridSampler;
 
+layout(set=0, binding=2) uniform simulationBlock {
+    Simulation simulation;
+};
+
 layout(location = 0) in vec2 vertexPos;
 
 layout(location = 0) out vec2 passPos;
@@ -32,13 +36,17 @@ void main()	{
     gridID = actual_mod(gridID, gridResolution);
 
     vec3 position = (vec3(gridID) + vec3(0.5f)) / gridResolution;
-    float size = 1.0f / length(vec3(gridResolution));
+
+    vec3 size = vec3(1.0f) / vec3(gridResolution);
+    float volume = size.x * size.y * size.z;
+    float radius = cube_radius(volume);
 
     vec4 gridData = texture(sampler3D(gridImage, gridSampler), position);
 
     float mass = gridData.w;
-    float density = mass / sphere_volume(size);
-    float alpha = clamp(density / 10000.0f, 0.0f, 1.0f);
+    float density = mass / volume;
+
+    float alpha = clamp(density / simulation.density, 0.0f, 1.0f);
 
     passPos = vertexPos;
     passVelocity = gridData.xyz;
@@ -46,5 +54,5 @@ void main()	{
 
     // align voxel to face camera
     gl_Position = mvp * vec4(position, 1);      // transform position into projected view space
-    gl_Position.xy += vertexPos * mix(0.0f, size * 2.0f, alpha);  // move position directly in view space
+    gl_Position.xy += vertexPos * (radius * 2.0f) * alpha;  // move position directly in view space
 }
\ No newline at end of file
diff --git a/projects/wobble_bobble/shaders/particle.inc b/projects/wobble_bobble/shaders/particle.inc
index 72791f9a..266e5e01 100644
--- a/projects/wobble_bobble/shaders/particle.inc
+++ b/projects/wobble_bobble/shaders/particle.inc
@@ -19,6 +19,24 @@ struct Particle {
 	mat4 mls;
 };
 
+#define SIM_FORM_SPHERE 0
+#define SIM_FORM_CUBE 1
+
+#define SIM_TYPE_HYPERELASTIC 0
+#define SIM_TYPE_FLUID 1
+
+struct Simulation {
+	float density;
+	float size;
+	float lame1;
+	float lame2;
+	int form;
+	int type;
+	float K;
+	float E;
+	float gamma;
+};
+
 const float PI = 3.1415926535897932384626433832795;
 
 float sphere_volume(float radius) {
@@ -29,6 +47,25 @@ float sphere_radius(float volume) {
 	return pow(volume * 3.0f / 4.0f / PI, 1.0f / 3.0f);
 }
 
+float cube_volume(float radius) {
+	return 8.0f * (radius * radius * radius);
+}
+
+float cube_radius(float volume) {
+	return pow(volume / 8.0f, 1.0f / 3.0f);
+}
+
+float simulationMass(float density, float size, int form) {
+	switch (form) {
+		case SIM_FORM_SPHERE:
+			return density * sphere_volume(size);
+		case SIM_FORM_CUBE:
+			return density * cube_volume(size);
+		default:
+			return 0.0f;
+	}
+}
+
 float weight_A(float x) {
 	return max(1.0f - x, 0.0f);
 }
diff --git a/projects/wobble_bobble/shaders/transform_particles_to_grid.comp b/projects/wobble_bobble/shaders/transform_particles_to_grid.comp
index 1bcea69d..236ff763 100644
--- a/projects/wobble_bobble/shaders/transform_particles_to_grid.comp
+++ b/projects/wobble_bobble/shaders/transform_particles_to_grid.comp
@@ -12,8 +12,6 @@ layout(set=0, binding=0, std430) readonly buffer particleBuffer {
 layout(set=0, binding=1, rgba32f) restrict writeonly uniform image3D gridImage;
 
 layout( push_constant ) uniform constants {
-    float lame1;
-    float lame2;
     float t;
     float dt;
     float speedfactor;
@@ -71,6 +69,10 @@ void main()	{
         memoryBarrierShared();
     }
 
+    if (any(isnan(gridValue.xyz)) || any(isinf(gridValue.xyz))) {
+        gridValue.xyz = vec3(0.0f);
+    }
+
     gridValue.xyz += vec3(0.0f, -9.81f * dts * gridValue.w, 0.0f);
 
     bvec3 lowerID = lessThanEqual(gl_GlobalInvocationID, ivec3(0));
diff --git a/projects/wobble_bobble/shaders/update_particle_velocities.comp b/projects/wobble_bobble/shaders/update_particle_velocities.comp
index cba521a5..60b5a72e 100644
--- a/projects/wobble_bobble/shaders/update_particle_velocities.comp
+++ b/projects/wobble_bobble/shaders/update_particle_velocities.comp
@@ -10,12 +10,14 @@ layout(set=0, binding=0, std430) restrict buffer particleBuffer {
     Particle particles [];
 };
 
-layout(set=0, binding=1) uniform texture3D gridImage;
-layout(set=0, binding=2) uniform sampler gridSampler;
+layout(set=0, binding=1) uniform simulationBlock {
+    Simulation simulation;
+};
+
+layout(set=0, binding=2) uniform texture3D gridImage;
+layout(set=0, binding=3) uniform sampler gridSampler;
 
 layout( push_constant ) uniform constants {
-    float lame1;
-    float lame2;
     float t;
     float dt;
     float speedfactor;
@@ -86,18 +88,30 @@ void main()	{
         }
 
         float J = max(determinant(F), EPSILON);
-
         float volume = sphere_volume(size);
 
-        mat3 F_T = transpose(F);
-        mat3 F_T_inv = inverse(F_T);
+        mat3 stress = mat3(0.0f);
 
-        mat3 P_term_0 = lame2 * (F - F_T_inv);
-        mat3 P_term_1 = lame1 * log(J) * F_T_inv;
+        switch (simulation.type) {
+            case SIM_TYPE_HYPERELASTIC:
+                mat3 F_T = transpose(F);
+                mat3 F_T_inv = inverse(F_T);
 
-        mat3 P = P_term_0 + P_term_1;
+                mat3 P_term_0 = simulation.lame2 * (F - F_T_inv);
+                mat3 P_term_1 = simulation.lame1 * log(J) * F_T_inv;
 
-        mat3 stress = P * F_T;
+                mat3 P = P_term_0 + P_term_1;
+
+                stress = P * F_T;
+                break;
+            case SIM_TYPE_FLUID:
+                float pressure = simulation.K * (1.0f / pow(J, simulation.gamma) - 1.0f);
+
+                stress = mat3(-pressure * J);
+                break;
+            default:
+                break;
+        }
 
         mls_Q -= dts * volume * stress * D_inv;
 
diff --git a/projects/wobble_bobble/src/main.cpp b/projects/wobble_bobble/src/main.cpp
index 5fabe497..117a7cf9 100644
--- a/projects/wobble_bobble/src/main.cpp
+++ b/projects/wobble_bobble/src/main.cpp
@@ -19,9 +19,25 @@ struct Particle {
 	glm::mat4 mls;
 };
 
-struct Physics {
+#define SIM_FORM_SPHERE 0
+#define SIM_FORM_CUBE 1
+
+#define SIM_TYPE_HYPERELASTIC 0
+#define SIM_TYPE_FLUID 1
+
+struct Simulation {
+	float density;
+	float size;
 	float lame1;
-    float lame2;
+	float lame2;
+	int form;
+	int type;
+	float K;
+	float E;
+	float gamma;
+};
+
+struct Physics {
 	float t;
 	float dt;
 	float speedfactor;
@@ -170,17 +186,33 @@ vkcv::ComputePipelineHandle createComputePipeline(vkcv::Core& core, vkcv::shader
 }
 
 void resetParticles(vkcv::Buffer<Particle>& particles, const glm::vec3& velocity,
-					float density, float size) {
+					float density, float size, int form) {
 	std::vector<Particle> particles_vec (particles.getCount());
 	
-	distributeParticlesCube(
-			particles_vec.data(),
-			particles_vec.size(),
-			glm::vec3(0.5f),
-			size,
-			density * sphere_volume(size),
-			velocity
-	);
+	switch (form) {
+		case SIM_FORM_SPHERE:
+			distributeParticlesSphere(
+					particles_vec.data(),
+					particles_vec.size(),
+					glm::vec3(0.5f),
+					size,
+					density * sphere_volume(size),
+					velocity
+			);
+			break;
+		case SIM_FORM_CUBE:
+			distributeParticlesCube(
+					particles_vec.data(),
+					particles_vec.size(),
+					glm::vec3(0.5f),
+					size,
+					density * sphere_volume(size),
+					velocity
+			);
+			break;
+		default:
+			break;
+	}
 	
 	particles.fill(particles_vec);
 }
@@ -217,18 +249,12 @@ int main(int argc, const char **argv) {
 			swapchainExtent.width,
 			swapchainExtent.height
 	).getHandle();
-
-	glm::vec3 initialVelocity (0.0f, 0.1f, 0.0f);
-	float density = 2500.0f;
-	float radius = 0.1f;
 	
 	vkcv::Buffer<Particle> particles = core.createBuffer<Particle>(
 			vkcv::BufferType::STORAGE,
 			1024
 	);
 	
-	resetParticles(particles, initialVelocity, density, radius);
-	
 	vkcv::Image grid = core.createImage(
 			vk::Format::eR32G32B32A32Sfloat,
 			64,
@@ -255,6 +281,26 @@ int main(int argc, const char **argv) {
 			vkcv::SamplerBorderColor::FLOAT_ZERO_TRANSPARENT
 	);
 	
+	vkcv::Buffer<Simulation> simulation = core.createBuffer<Simulation>(
+			vkcv::BufferType::UNIFORM, 1, vkcv::BufferMemoryType::HOST_VISIBLE
+	);
+	
+	Simulation* sim = simulation.map();
+	
+	glm::vec3 initialVelocity (0.0f, 0.1f, 0.0f);
+	
+	sim->density = 2500.0f;
+	sim->size = 0.1f;
+	sim->lame1 = 10.0f;
+	sim->lame2 = 20.0f;
+	sim->form = SIM_FORM_SPHERE;
+	sim->type = SIM_TYPE_HYPERELASTIC;
+	sim->K = 2.2f;
+	sim->E = 35.0f;
+	sim->gamma = 1.330f;
+	
+	resetParticles(particles, initialVelocity, sim->density, sim->size, sim->form);
+	
 	vkcv::shader::GLSLCompiler compiler;
 	
 	std::vector<vkcv::DescriptorSetHandle> initParticleWeightsSets;
@@ -296,8 +342,9 @@ int main(int argc, const char **argv) {
 	{
 		vkcv::DescriptorWrites writes;
 		writes.storageBufferWrites.push_back(vkcv::BufferDescriptorWrite(0, particles.getHandle()));
-		writes.sampledImageWrites.push_back(vkcv::SampledImageDescriptorWrite(1, grid.getHandle()));
-		writes.samplerWrites.push_back(vkcv::SamplerDescriptorWrite(2, gridSampler));
+		writes.uniformBufferWrites.push_back(vkcv::BufferDescriptorWrite(1, simulation.getHandle()));
+		writes.sampledImageWrites.push_back(vkcv::SampledImageDescriptorWrite(2, grid.getHandle()));
+		writes.samplerWrites.push_back(vkcv::SamplerDescriptorWrite(3, gridSampler));
 		core.writeDescriptorSet(updateParticleVelocitiesSets[0], writes);
 	}
 	
@@ -404,6 +451,7 @@ int main(int argc, const char **argv) {
 		vkcv::DescriptorWrites writes;
 		writes.sampledImageWrites.push_back(vkcv::SampledImageDescriptorWrite(0, grid.getHandle()));
 		writes.samplerWrites.push_back(vkcv::SamplerDescriptorWrite(1, gridSampler));
+		writes.uniformBufferWrites.push_back(vkcv::BufferDescriptorWrite(2, simulation.getHandle()));
 		core.writeDescriptorSet(gfxSetGrid, writes);
 	}
 	
@@ -435,7 +483,6 @@ int main(int argc, const char **argv) {
 	gfxPipelineConfigGrid.m_VertexLayout = vertexLayoutGrid;
 	gfxPipelineConfigGrid.m_DescriptorLayouts = { gfxSetLayoutGrid };
 	gfxPipelineConfigGrid.m_UseDynamicViewport = true;
-	gfxPipelineConfigGrid.m_blendMode = vkcv::BlendMode::Additive;
 	
 	vkcv::VertexLayout vertexLayoutParticles ({
 		vkcv::VertexBinding(0, gfxProgramParticles.getVertexAttachments())
@@ -548,8 +595,6 @@ int main(int argc, const char **argv) {
 	
 	bool renderGrid = true;
 	
-	float lame1 = 10.0f;
-	float lame2 = 20.0f;
 	float speed_factor = 1.0f;
 	
 	auto start = std::chrono::system_clock::now();
@@ -585,8 +630,6 @@ int main(int argc, const char **argv) {
 		current = next;
 		
 		Physics physics;
-		physics.lame1 = lame1;
-		physics.lame2 = lame2;
 		physics.t = static_cast<float>(0.000001 * static_cast<double>(time.count()));
 		physics.dt = static_cast<float>(0.000001 * static_cast<double>(deltatime.count()));
 		physics.speedfactor = speed_factor;
@@ -642,6 +685,7 @@ int main(int argc, const char **argv) {
 			
 			core.recordBeginDebugLabel(cmdStream, "UPDATE PARTICLE VELOCITIES", {0.78f, 0.89f, 0.94f, 1.0f});
 			core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
+			core.recordBufferMemoryBarrier(cmdStream, simulation.getHandle());
 			core.prepareImageForSampling(cmdStream, grid.getHandle());
 			
 			core.recordComputeDispatchToCmdStream(
@@ -665,6 +709,7 @@ int main(int argc, const char **argv) {
 		
 		if (renderGrid) {
 			core.recordBeginDebugLabel(cmdStream, "RENDER GRID", { 0.13f, 0.20f, 0.22f, 1.0f });
+			core.recordBufferMemoryBarrier(cmdStream, simulation.getHandle());
 			core.prepareImageForSampling(cmdStream, grid.getHandle());
 			
 			core.recordDrawcallsToCmdStream(
@@ -714,25 +759,27 @@ int main(int argc, const char **argv) {
 		
 		gui.beginGUI();
 		ImGui::Begin("Settings");
+		
+		ImGui::Combo("Type", &(sim->type), "Hyperelastic\0Fluid", 2);
 
-		ImGui::SliderFloat("Density", &density, std::numeric_limits<float>::epsilon(), 5000.0f);
+		ImGui::SliderFloat("Density", &(sim->density), std::numeric_limits<float>::epsilon(), 5000.0f);
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::SmallButton("Reset##density")) {
-			density = 2500.0f;
+			sim->density = 2500.0f;
 		}
 		
-		ImGui::SliderFloat("Radius", &radius, 0.0f, 0.5f);
+		ImGui::SliderFloat("Radius", &(sim->size), 0.0f, 0.5f);
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::SmallButton("Reset##radius")) {
-			radius = 0.1f;
+			sim->size = 0.1f;
 		}
-
-		ImGui::BeginGroup();
-		ImGui::SliderFloat("Compression Modulus", &lame1, 0.0f, 100.0f);
-		ImGui::EndGroup();
 		
 		ImGui::BeginGroup();
-		ImGui::SliderFloat("Elasticity Modulus", &lame2, 0.0f, 100.0f);
+		ImGui::SliderFloat("Bulk Modulus", &(sim->K), 0.0f, 100.0f);
+		ImGui::SliderFloat("Young's Modulus", &(sim->E), 0.0f, 100.0f);
+		ImGui::SliderFloat("Heat Capacity Ratio", &(sim->gamma), 1.0f, 2.0f);
+		ImGui::SliderFloat("Lame1", &(sim->lame1), 0.0f, 100.0f);
+		ImGui::SliderFloat("Lame2", &(sim->lame2), 0.0f, 100.0f);
 		ImGui::EndGroup();
 
 		ImGui::Spacing();
@@ -745,7 +792,7 @@ int main(int argc, const char **argv) {
 		ImGui::DragFloat3("Initial Velocity", reinterpret_cast<float*>(&initialVelocity), 0.001f);
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::Button("Reset##particle_velocity")) {
-			resetParticles(particles, initialVelocity, density, radius);
+			resetParticles(particles, initialVelocity, sim->density, sim->size, sim->form);
 		}
 		
 		ImGui::End();
@@ -754,5 +801,6 @@ int main(int argc, const char **argv) {
 		core.endFrame(windowHandle);
 	}
 	
+	simulation.unmap();
 	return 0;
 }
-- 
GitLab