From b555992fa4a4d91b1e27fa7e9be9fd7fc4023fe8 Mon Sep 17 00:00:00 2001
From: Tobias Frisch <tfrisch@uni-koblenz.de>
Date: Sat, 5 Feb 2022 13:53:05 +0100
Subject: [PATCH] Corrected epsilon to have stable elastic impact
Signed-off-by: Tobias Frisch <tfrisch@uni-koblenz.de>
---
projects/wobble_bobble/shaders/particle.inc | 2 +-
.../shaders/transform_particles_to_grid.comp | 2 +-
.../shaders/update_particle_velocities.comp | 36 ++++++++++++-------
3 files changed, 25 insertions(+), 15 deletions(-)
diff --git a/projects/wobble_bobble/shaders/particle.inc b/projects/wobble_bobble/shaders/particle.inc
index db022be7..72791f9a 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 1bcea69d..b85ed31e 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 2d773d85..cba521a5 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;
--
GitLab