Fix for [#34099] Particles leaking from moving meshes

This was caused by a floating point precision error. During collision detection, Newton-Raphson iteration is used to find the exact time of the collision. But when using subframes, the initial Newton step was too small. Now the initial step is given in absolute units. When subframes = 0, this should behave almost the same as before.

Thanks to Janne Karhu, Lukas Toenne and Ton Roosendaal for their help with this patch, and to AutoCRC for funding.
This commit is contained in:
Alex Fraser
2013-02-19 02:24:52 +00:00
parent 442d16b468
commit 4e73ff0d63
2 changed files with 42 additions and 15 deletions

View File

@@ -220,7 +220,10 @@ typedef struct ParticleCollision {
ParticleCollisionElement pce;
float total_time, inv_timestep;
/* total_time is the amount of time in this subframe
* inv_total_time is the opposite
* inv_timestep is the inverse of the amount of time in this frame */
float total_time, inv_total_time, inv_timestep;
float radius;
float co1[3], co2[3];

View File

@@ -3028,13 +3028,20 @@ static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, f
normalize_qt(pa->state.rot);
}
/************************************************/
/* Collisions */
/************************************************/
/************************************************
* Collisions
*
* The algorithm is roughly:
* 1. Use a BVH tree to search for faces that a particle may collide with.
* 2. Use Newton's method to find the exact time at which the collision occurs.
* http://en.wikipedia.org/wiki/Newton's_method
*
************************************************/
#define COLLISION_MAX_COLLISIONS 10
#define COLLISION_MIN_RADIUS 0.001f
#define COLLISION_MIN_DISTANCE 0.0001f
#define COLLISION_ZERO 0.00001f
#define COLLISION_INIT_STEP 0.00008f
typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
{
@@ -3189,16 +3196,20 @@ static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce
/* find first root in range [0-1] starting from 0 */
static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
{
float t0, t1, d0, d1, dd, n[3];
float t0, t1, dt_init, d0, d1, dd, n[3];
int iter;
pce->inv_nor = -1;
/* Initial step size should be small, but not too small or floating point
* precision errors will appear. - z0r */
dt_init = COLLISION_INIT_STEP * col->inv_total_time;
/* start from the beginning */
t0 = 0.f;
collision_interpolate_element(pce, t0, col->f, col);
d0 = distance_func(col->co1, radius, pce, n);
t1 = 0.001f;
t1 = dt_init;
d1 = 0.f;
for (iter=0; iter<10; iter++) {//, itersum++) {
@@ -3208,11 +3219,6 @@ static float collision_newton_rhapson(ParticleCollision *col, float radius, Part
d1 = distance_func(pce->p, radius, pce, n);
/* no movement, so no collision */
if (d1 == d0) {
return -1.f;
}
/* particle already inside face, so report collision */
if (iter == 0 && d0 < 0.f && d0 > -radius) {
copy_v3_v3(pce->p, col->co1);
@@ -3220,7 +3226,24 @@ static float collision_newton_rhapson(ParticleCollision *col, float radius, Part
pce->inside = 1;
return 0.f;
}
/* Zero gradient (no movement relative to element). Can't step from
* here. */
if (d1 == d0) {
/* If first iteration, try from other end where the gradient may be
* greater. Note: code duplicated below. */
if (iter == 0) {
t0 = 1.f;
collision_interpolate_element(pce, t0, col->f, col);
d0 = distance_func(col->co2, radius, pce, n);
t1 = 1.0f - dt_init;
d1 = 0.f;
continue;
}
else
return -1.f;
}
dd = (t1-t0)/(d1-d0);
t0 = t1;
@@ -3228,14 +3251,14 @@ static float collision_newton_rhapson(ParticleCollision *col, float radius, Part
t1 -= d1*dd;
/* particle movin away from plane could also mean a strangely rotating face, so check from end */
/* Particle moving away from plane could also mean a strangely rotating
* face, so check from end. Note: code duplicated above. */
if (iter == 0 && t1 < 0.f) {
t0 = 1.f;
collision_interpolate_element(pce, t0, col->f, col);
d0 = distance_func(col->co2, radius, pce, n);
t1 = 0.999f;
t1 = 1.0f - dt_init;
d1 = 0.f;
continue;
}
else if (iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
@@ -3683,6 +3706,7 @@ static void collision_check(ParticleSimulationData *sim, int p, float dfra, floa
memset(&col, 0, sizeof(ParticleCollision));
col.total_time = timestep * dfra;
col.inv_total_time = 1.0f/col.total_time;
col.inv_timestep = 1.0f/timestep;
col.cfra = cfra;