Corrected epsilon to have stable elastic impact

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

projects/wobble_bobble/shaders/particle.inc

 #ifndef PARTICLE_INC
 #define PARTICLE_INC
 
-#define EPSILON 0.00001f
+#define EPSILON 0.00000001f
 
 struct ParticleMinimal {
     vec3 position;

projects/wobble_bobble/shaders/transform_particles_to_grid.comp

@@ -71,7 +71,7 @@ void main()	{
         memoryBarrierShared();
     }
 
-    gridValue.xyz += vec3(0.0f, -9.81f * dts * gridValue.w, 0.0f);
+    gridValue.xyz += vec3(0.0f, -9.81f * dts * gridValue.w * 0.0f, 0.0f);
 
     bvec3 lowerID = lessThanEqual(gl_GlobalInvocationID, ivec3(0));
     bvec3 negativeVelocity = lessThan(gridValue.xyz, vec3(0.0f));

projects/wobble_bobble/shaders/update_particle_velocities.comp

@@ -76,25 +76,35 @@ void main()	{
 
         mat3 F = mat3(particle.deformation);
 
-        if (abs(determinant(affine_D)) >= EPSILON) {
-            mat3 D_inv = inverse(affine_D);
-            float J = determinant(F);
+        mat3 D_inv = inverse(affine_D);
+        float D_det = determinant(D_inv);
 
-            if (J > 0.0f) {
-                float volume = sphere_volume(size);
+        if ((isnan(D_det)) || (isinf(D_det))) {
+            D_inv = mat3(0.0f);
+        } else {
+            D_inv *= min(abs(D_det), 1.0f / EPSILON) / abs(D_det);
+        }
 
-                mat3 F_T = transpose(F);
+        float J = max(determinant(F), EPSILON);
 
-                mat3 stress = lame2 * (F * F_T - mat3(1.0f)) + lame1 * log(J);
+        float volume = sphere_volume(size);
 
-                mls_Q -= dts * volume * stress * D_inv;
-            }
+        mat3 F_T = transpose(F);
+        mat3 F_T_inv = inverse(F_T);
 
-            affine_C = affine_B * D_inv;
-            mls_Q += affine_C * mass;
+        mat3 P_term_0 = lame2 * (F - F_T_inv);
+        mat3 P_term_1 = lame1 * log(J) * F_T_inv;
 
-            F = (mat3(1.0f) + dts * affine_C) * F;
-        }
+        mat3 P = P_term_0 + P_term_1;
+
+        mat3 stress = P * F_T;
+
+        mls_Q -= dts * volume * stress * D_inv;
+
+        affine_C = affine_B * D_inv;
+        mls_Q += affine_C * mass;
+
+        F = (mat3(1.0f) + dts * affine_C) * F;
 
         position = position + velocity_pic * dts;