Added several physical properties to match materials

Signed-off-by: Tobias Frisch <tfrisch@uni-koblenz.de>

projects/wobble_bobble/shaders/grid.vert

 #version 450
+#extension GL_GOOGLE_include_directive : enable
+
+#include "particle.inc"
 
 layout(set=0, binding=0) uniform texture3D gridImage;
 layout(set=0, binding=1) uniform sampler gridSampler;
@@ -34,6 +37,8 @@ void main()	{
     vec4 gridData = texture(sampler3D(gridImage, gridSampler), position);
 
     float mass = gridData.w;
+    float density = mass / sphere_volume(size);
+    float alpha = clamp(density / 100000000.0f, 0.0f, 1.0f);
 
     passPos = vertexPos;
     passVelocity = gridData.xyz;
@@ -41,5 +46,5 @@ void main()	{
 
     // align voxel to face camera
     gl_Position = mvp * vec4(position, 1);      // transform position into projected view space
-    gl_Position.xy += vertexPos * (size * 2.0f) * (mass * 100.0f);  // move position directly in view space
+    gl_Position.xy += vertexPos * mix(0.0f, size * 2.0f, alpha);  // move position directly in view space
 }
\ No newline at end of file

projects/wobble_bobble/shaders/init_particle_volumes.comp

@@ -9,11 +9,44 @@ layout(set=0, binding=0, std430) buffer particleBuffer {
     Particle particles [];
 };
 
+layout(set=0, binding=1) uniform texture3D gridImage;
+layout(set=0, binding=2) uniform sampler gridSampler;
+
 void main()	{
     memoryBarrierBuffer();
+    memoryBarrierImage();
 
     if (gl_GlobalInvocationID.x < particles.length()) {
-        // nothing
+        ParticleMinimal minimal = particles[gl_GlobalInvocationID.x].minimal;
+
+        ivec3 gridResolution = textureSize(sampler3D(gridImage, gridSampler), 0);
+        ivec3 gridWindow = ivec3(minimal.size * 2.0f * gridResolution);
+
+        float volume = sphere_volume(minimal.size);
+
+        float mass = 0.0f;
+        int i, j, k;
+
+        for (i = -gridWindow.x; i <= gridWindow.x; i++) {
+            for (j = -gridWindow.y; j <= gridWindow.y; j++) {
+                for (k = -gridWindow.z; k <= gridWindow.z; k++) {
+                    vec3 offset = vec3(i, j, k) / gridResolution;
+                    vec3 voxel = minimal.position + offset;
+
+                    if (length(offset) < minimal.size * 2.0f) {
+                        vec4 gridSample = texture(sampler3D(gridImage, gridSampler), voxel);
+
+                        mass += gridSample.w * voxel_particle_weight(voxel, minimal);
+                    }
+                }
+            }
+        }
+
+        if (volume > 0.0f) {
+            //volume = minimal.mass / (mass / volume);
+        }
+
+        particles[gl_GlobalInvocationID.x].minimal.size = sphere_radius(volume);
     }
 
     memoryBarrierBuffer();

projects/wobble_bobble/shaders/particle.inc

@@ -13,6 +13,16 @@ struct Particle {
     mat4 deformation;
 };
 
+const float PI = 3.1415926535897932384626433832795;
+
+float sphere_volume(float radius) {
+	return 4.0f * (radius * radius * radius) * PI / 3.0f;
+}
+
+float sphere_radius(float volume) {
+	return pow(volume * 3.0f / 4.0f / PI, 1.0f / 3.0f);
+}
+
 float weight_A(float x) {
 	return max(1.0f - x, 0.0f);
 }

projects/wobble_bobble/shaders/update_grid_forces.comp

 #version 450
 #extension GL_GOOGLE_include_directive : enable
 
-const float PI = 3.1415926535897932384626433832795;
-
 layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
 
 #include "particle.inc"
@@ -42,7 +40,9 @@ void main()	{
         poisson /= (6.0f * K);
     }
 
-    const vec3 position = (vec3(gl_GlobalInvocationID) + vec3(0.5f)) / imageSize(gridImage);
+    const vec3 gridResolution = vec3(imageSize(gridImage));
+    const vec3 position = (vec3(gl_GlobalInvocationID) + vec3(0.5f)) / gridResolution;
+    const float h3 = 1.0f / gridResolution.x / gridResolution.y / gridResolution.z;
 
     memoryBarrierImage();
 
@@ -75,8 +75,7 @@ void main()	{
         memoryBarrierShared();
 
         for (uint i = 0; i < SHARED_PARTICLES_BATCH_SIZE; i++) {
-            float size = shared_particles[gl_LocalInvocationIndex].minimal.size;
-            float volume = 4.0f / 3.0f * PI * size * size * size;
+            float volume = sphere_volume(shared_particles[gl_LocalInvocationIndex].minimal.size);
             mat3 F = mat3(shared_particles[gl_LocalInvocationIndex].deformation);
 
             float J = determinant(F);

projects/wobble_bobble/src/main.cpp

@@ -26,6 +26,14 @@ struct Tweaks {
 	float beta;
 };
 
+float sphere_volume(float radius) {
+	return 4.0f * (radius * radius * radius) * M_PI / 3.0f;
+}
+
+float sphere_radius(float volume) {
+	return pow(volume * 3.0f / 4.0f / M_PI, 1.0f / 3.0f);
+}
+
 std::random_device random_dev;
 std::uniform_int_distribution<int> dist(0, RAND_MAX);
 
@@ -55,11 +63,14 @@ void distributeParticles(Particle *particles, size_t count, const glm::vec3& cen
 		particles[i].size = size;
 		particles[i].velocity = velocity;
 		
-		volume += size;
+		volume += sphere_volume(size);
 	}
 	
+	// Keep the same densitiy as planned!
+	mass *= (volume / sphere_volume(radius));
+	
 	for (size_t i = 0; i < count; i++) {
-		particles[i].mass = (mass * particles[i].size / volume);
+		particles[i].mass = (mass * sphere_volume(particles[i].size) / volume);
 		particles[i].deformation = glm::mat4(1.0f);
 	}
 }
@@ -105,15 +116,16 @@ vkcv::ComputePipelineHandle createComputePipeline(vkcv::Core& core, vkcv::shader
 	return core.createComputePipeline(config);
 }
 
-void resetParticles(vkcv::Buffer<Particle>& particles, const glm::vec3& velocity) {
+void resetParticles(vkcv::Buffer<Particle>& particles, const glm::vec3& velocity,
+					float density, float size) {
 	std::vector<Particle> particles_vec (particles.getCount());
 	
 	distributeParticles(
 			particles_vec.data(),
 			particles_vec.size(),
 			glm::vec3(0.5f),
-			0.05f,
-			0.27f,
+			size,
+			density * sphere_volume(size),
 			velocity
 	);
 	
@@ -154,13 +166,15 @@ int main(int argc, const char **argv) {
 	).getHandle();
 	
 	glm::vec3 initialVelocity (0.0f, 0.0f, 0.0f);
+	float density = 2500.0f;
+	float radius = 0.1f;
 	
 	vkcv::Buffer<Particle> particles = core.createBuffer<Particle>(
 			vkcv::BufferType::STORAGE,
 			1024
 	);
 	
-	resetParticles(particles, initialVelocity);
+	resetParticles(particles, initialVelocity, density, radius);
 	
 	vkcv::Image grid = core.createImage(
 			vk::Format::eR32G32B32A32Sfloat,
@@ -230,6 +244,8 @@ 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));
 		core.writeDescriptorSet(initParticleVolumesSets[0], writes);
 	}
 	
@@ -542,7 +558,10 @@ int main(int argc, const char **argv) {
 	
 	// Glass is glass and glass breaks...
 	float compression_modulus = 65.0f;
+	int compression_exponent = 9;
+	
 	float elasticity_modulus = 45.0f;
+	int elasticity_exponent = 9;
 	
 	float alpha = 0.5f;
 	float beta = 0.75f;
@@ -580,8 +599,8 @@ int main(int argc, const char **argv) {
 		current = next;
 		
 		Physics physics;
-		physics.K = static_cast<float>(compression_modulus * std::pow(10.0, 9.0));
-		physics.E = static_cast<float>(elasticity_modulus * std::pow(10.0, 9.0));
+		physics.K = static_cast<float>(compression_modulus * std::pow(10.0, compression_exponent));
+		physics.E = static_cast<float>(elasticity_modulus * std::pow(10.0, elasticity_exponent));
 		physics.t = static_cast<float>(0.000001 * static_cast<double>(time.count()));
 		physics.dt = static_cast<float>(0.000001 * static_cast<double>(deltatime.count()));
 		
@@ -626,6 +645,8 @@ int main(int argc, const char **argv) {
 		
 		if (!initializedParticleVolumes) {
 			core.recordBeginDebugLabel(cmdStream, "INIT PARTICLE VOLUMES", { 0.78f, 0.89f, 0.94f, 1.0f });
+			core.prepareImageForSampling(cmdStream, grid.getHandle());
+			
 			core.recordComputeDispatchToCmdStream(
 					cmdStream,
 					initParticleVolumesPipeline,
@@ -642,6 +663,8 @@ int main(int argc, const char **argv) {
 		}
 		
 		core.recordBeginDebugLabel(cmdStream, "UPDATE GRID FORCES", { 0.47f, 0.77f, 0.85f, 1.0f });
+		core.prepareImageForStorage(cmdStream, grid.getHandle());
+		
 		core.recordComputeDispatchToCmdStream(
 				cmdStream,
 				updateGridForcesPipeline,
@@ -754,18 +777,49 @@ int main(int argc, const char **argv) {
 		core.submitCommandStream(cmdStream);
 		
 		gui.beginGUI();
-		
 		ImGui::Begin("Settings");
 		
+		ImGui::SliderFloat("Density", &density, std::numeric_limits<float>::epsilon(), 5000.0f);
+		ImGui::SameLine(0.0f, 10.0f);
+		if (ImGui::SmallButton("Reset##density")) {
+			density = 2500.0f;
+		}
+		
+		ImGui::SliderFloat("Radius", &radius, 0.0f, 0.5f);
+		ImGui::SameLine(0.0f, 10.0f);
+		if (ImGui::SmallButton("Reset##radius")) {
+			radius = 0.1f;
+		}
+		
+		ImGui::BeginGroup();
 		ImGui::SliderFloat("Compression Modulus", &compression_modulus, 0.0f, 500.0f);
+		ImGui::SliderInt("##compression_exponent", &compression_exponent, 1, 9);
+		ImGui::SameLine(0.0f, 10.0f);
+		if (ImGui::SmallButton("Reset##compression")) {
+			compression_modulus = 65.0f;
+			compression_exponent = 9;
+		}
+		ImGui::EndGroup();
+		
+		ImGui::BeginGroup();
 		ImGui::SliderFloat("Elasticity Modulus", &elasticity_modulus, 0.0f, 1000.0f);
+		ImGui::SliderInt("##elasticity_exponent", &elasticity_exponent, 1, 9);
+		ImGui::SameLine(0.0f, 10.0f);
+		if (ImGui::SmallButton("Reset##elasticity")) {
+			elasticity_modulus = 45.0f;
+			elasticity_exponent = 9;
+		}
+		ImGui::EndGroup();
 		
+		ImGui::Spacing();
 		ImGui::Checkbox("Render Grid", &renderGrid);
+		
 		ImGui::SliderFloat("Alpha (PIC -> FLIP)", &alpha, 0.0f, 1.0f);
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::SmallButton("Reset##alpha")) {
 		    alpha = 0.5f;
 		}
+		
 		ImGui::SliderFloat("Beta (Alpha -> APIC)", &beta, 0.0f, 1.0f);
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::SmallButton("Reset##beta")) {
@@ -775,11 +829,11 @@ int main(int argc, const char **argv) {
 		ImGui::DragFloat3("Initial Velocity", reinterpret_cast<float*>(&initialVelocity));
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::Button("Reset##particle_velocity")) {
-			resetParticles(particles, initialVelocity);
+			resetParticles(particles, initialVelocity, density, radius);
+			initializedParticleVolumes = false;
 		}
 		
 		ImGui::End();
-		
 		gui.endGUI();
 		
 		core.endFrame(windowHandle);