diff --git a/projects/wobble_bobble/shaders/particle.inc b/projects/wobble_bobble/shaders/particle.inc
index db022be78f120c37c886edc3c740e6a391ff97b1..72791f9a15d18874603e87a2fde3a48a81a3829f 100644
--- a/projects/wobble_bobble/shaders/particle.inc
+++ b/projects/wobble_bobble/shaders/particle.inc
@@ -1,7 +1,7 @@
 #ifndef PARTICLE_INC
 #define PARTICLE_INC
 
-#define EPSILON 0.00001f
+#define EPSILON 0.00000001f
 
 struct ParticleMinimal {
     vec3 position;
diff --git a/projects/wobble_bobble/shaders/transform_particles_to_grid.comp b/projects/wobble_bobble/shaders/transform_particles_to_grid.comp
index 1bcea69def425ce20f282b038c5b941206c8f60d..b85ed31ef61a97f919840244dead5b0aa15dad3e 100644
--- a/projects/wobble_bobble/shaders/transform_particles_to_grid.comp
+++ b/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));
diff --git a/projects/wobble_bobble/shaders/update_particle_velocities.comp b/projects/wobble_bobble/shaders/update_particle_velocities.comp
index 2d773d85b7d4bd2e483b189dc2225cb42dc84e6f..cba521a5ba540abd730803a857efeab336987741 100644
--- a/projects/wobble_bobble/shaders/update_particle_velocities.comp
+++ b/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;