Added modulus for material properties and explosions

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

projects/wobble_bobble/shaders/update_grid_forces.comp

@@ -14,6 +14,8 @@ layout(set=0, binding=2, std430) buffer particleBuffer {
 };
 
 layout( push_constant ) uniform constants {
+    float K;
+    float E;
     float t;
     float dt;
 };
@@ -23,8 +25,24 @@ layout( push_constant ) uniform constants {
 shared Particle shared_particles [SHARED_PARTICLES_BATCH_SIZE];
 
 void main()	{
+    const float K9_E = (9.0f * K - E);
+
+    float lame1 = 3.0f * K * (3.0f * K - E);
+    float lame2 = 3.0f * K * E;
+    float poisson = (3.0f * K - E);
+    float M = 3.0f * K * (3.0f * K + E);
+
+    if (K9_E > 0.0f) {
+        lame1 /= K9_E;
+        lame2 /= K9_E;
+        M /= K9_E;
+    }
+
+    if (K > 0.0f) {
+        poisson /= (6.0f * K);
+    }
+
     const vec3 position = (vec3(gl_GlobalInvocationID) + vec3(0.5f)) / imageSize(gridImage);
-    const float elasticity_module = 12.5f * 1000000.0f; // 10..15 = Wood
 
     memoryBarrierImage();
 
@@ -59,14 +77,17 @@ void main()	{
         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;
-            mat3 deformation = mat3(shared_particles[gl_LocalInvocationIndex].deformation);
+            mat3 F = mat3(shared_particles[gl_LocalInvocationIndex].deformation);
+
+            float J = determinant(F);
+
+            if ((J > 0.0f) && (volume > 0.0f)) {
+                mat3 F_invT = inverse(transpose(F));
 
-            float J = determinant(deformation);
+                vec3 v = shared_particles[gl_LocalInvocationIndex].minimal.velocity;
+                float m = shared_particles[gl_LocalInvocationIndex].minimal.mass;
 
-            if ((abs(J) > 0.0f) && (volume > 0.0f)) {
-                mat3 epsilon = deformation - mat3(1.0f);
-                mat3 delta = elasticity_module * epsilon;
-                mat3 delta_cauchy = J * delta * inverse(transpose(deformation));
+                mat3 delta = lame2 * (F - F_invT) + lame1 * log(J) * F_invT;
 
                 vec3 weight_gradient = voxel_particle_grad_weight(
                     position,
@@ -75,7 +96,7 @@ void main()	{
 
                 force -= (
                     volume *
-                    delta_cauchy *
+                    delta * transpose(F) *
                     weight_gradient
                 );
             }

projects/wobble_bobble/shaders/update_particle_deformation.comp

@@ -36,11 +36,10 @@ void main()	{
                     if (length(offset) < minimal.size * 2.0f) {
                         vec4 gridSample = texture(sampler3D(gridImage, gridSampler), voxel);
 
-                        // TODO: needs to be activated once everything else is stable!
-                        /*velocity_gradient += outerProduct(
+                        velocity_gradient += outerProduct(
                             gridSample.xyz,
                             voxel_particle_grad_weight(voxel, minimal)
-                        );*/
+                        );
                     }
                 }
             }

projects/wobble_bobble/shaders/update_particle_positions.comp

@@ -10,6 +10,8 @@ layout(set=0, binding=0, std430) buffer particleBuffer {
 };
 
 layout( push_constant ) uniform constants {
+    float K;
+    float E;
     float t;
     float dt;
 };

projects/wobble_bobble/src/main.cpp

@@ -14,7 +14,9 @@ struct Particle {
 	glm::mat4 deformation;
 };
 
-struct Time {
+struct Physics {
+	float K;
+	float E;
 	float t;
 	float dt;
 };
@@ -537,6 +539,11 @@ int main(int argc, const char **argv) {
 	
 	bool initializedParticleVolumes = false;
 	bool renderGrid = true;
+	
+	// Glass is glass and glass breaks...
+	float compression_modulus = 65.0f;
+	float elasticity_modulus = 45.0f;
+	
 	float alpha = 0.5f;
 	float beta = 0.75f;
 	
@@ -572,12 +579,14 @@ int main(int argc, const char **argv) {
 		
 		current = next;
 		
-		Time timeConstants;
-		timeConstants.t = static_cast<float>(0.000001 * static_cast<double>(time.count()));
-		timeConstants.dt = static_cast<float>(0.000001 * static_cast<double>(deltatime.count()));
+		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.t = static_cast<float>(0.000001 * static_cast<double>(time.count()));
+		physics.dt = static_cast<float>(0.000001 * static_cast<double>(deltatime.count()));
 		
-		vkcv::PushConstants timePushConstants (sizeof(timeConstants));
-		timePushConstants.appendDrawcall(timeConstants);
+		vkcv::PushConstants physicsPushConstants (sizeof(physics));
+		physicsPushConstants.appendDrawcall(physics);
 		
 		Tweaks tweaks;
 		tweaks.alpha = alpha;
@@ -586,7 +595,7 @@ int main(int argc, const char **argv) {
 		vkcv::PushConstants tweakPushConstants (sizeof(tweaks));
 		tweakPushConstants.appendDrawcall(tweaks);
 		
-		cameraManager.update(timeConstants.dt);
+		cameraManager.update(physics.dt);
 
 		glm::mat4 mvp = cameraManager.getActiveCamera().getMVP();
 		vkcv::PushConstants cameraPushConstants (sizeof(glm::mat4));
@@ -640,7 +649,7 @@ int main(int argc, const char **argv) {
 				{ vkcv::DescriptorSetUsage(
 						0, core.getDescriptorSet(updateGridForcesSets[0]).vulkanHandle
 				) },
-				timePushConstants
+				physicsPushConstants
 		);
 		
 		core.recordImageMemoryBarrier(cmdStream, grid.getHandle());
@@ -686,7 +695,7 @@ int main(int argc, const char **argv) {
 				{ vkcv::DescriptorSetUsage(
 						0, core.getDescriptorSet(updateParticlePositionsSets[0]).vulkanHandle
 				) },
-				timePushConstants
+				physicsPushConstants
 		);
 		
 		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
@@ -748,6 +757,9 @@ int main(int argc, const char **argv) {
 		
 		ImGui::Begin("Settings");
 		
+		ImGui::SliderFloat("Compression Modulus", &compression_modulus, 0.0f, 500.0f);
+		ImGui::SliderFloat("Elasticity Modulus", &elasticity_modulus, 0.0f, 1000.0f);
+		
 		ImGui::Checkbox("Render Grid", &renderGrid);
 		ImGui::SliderFloat("Alpha (PIC -> FLIP)", &alpha, 0.0f, 1.0f);
 		ImGui::SameLine(0.0f, 10.0f);
@@ -760,7 +772,7 @@ int main(int argc, const char **argv) {
 		    beta = 0.75f;
 		}
 		
-		ImGui::DragFloat3("Initial velocity", reinterpret_cast<float*>(&initialVelocity));
+		ImGui::DragFloat3("Initial Velocity", reinterpret_cast<float*>(&initialVelocity));
 		ImGui::SameLine(0.0f, 10.0f);
 		if (ImGui::Button("Reset##particle_velocity")) {
 			resetParticles(particles, initialVelocity);