New particle collisions code:
authorJanne Karhu <jhkarh@gmail.com>
Fri, 18 Mar 2011 15:31:32 +0000 (15:31 +0000)
committerJanne Karhu <jhkarh@gmail.com>
Fri, 18 Mar 2011 15:31:32 +0000 (15:31 +0000)
* The old collisions code detected particle collisions by calculating the
  collision times analytically from the collision mesh faces. This was
  pretty accurate, but didn't support rotating/deforming faces at all, as
  the equations for these quickly become quite nasty.
* The new code uses a simple "distance to plane/edge/vert" function and
  iterates this with the Newton-Rhapson method to find the closest particle
  distance during a simulation step.
* The advantage in this is that the collision object can now move, rotate,
  scale or even deform freely and collisions are still detected reliably.
* For some extreme movements the calculation errors could stack up so much
  that the detection fails, but this can be easily fixed by increasing the
  particle size or simulation substeps.
* As a side note the algorithm doesn't really do point particles anymore,
  but uses a very small radius as the particle size when "size deflect" isn't
  selected.
* I've also updated the collision response code a bit, so now the particles
  shouldn't leak even from tight corners.

All in all the collisions code is now much cleaner and more robust than before!

source/blender/blenkernel/BKE_particle.h
source/blender/blenkernel/intern/boids.c
source/blender/blenkernel/intern/particle_system.c
source/blender/blenloader/intern/readfile.c
source/blender/editors/physics/particle_edit.c
source/blender/makesdna/DNA_modifier_types.h
source/blender/modifiers/intern/MOD_collision.c

index e4fbd5603f2b0eb27ab18d21cec1bb8b144ddaff..9b1651b12cd0d071d8969c6a5ef2eea3cc677d17 100644 (file)
@@ -155,19 +155,57 @@ typedef struct ParticleBillboardData
        short align, uv_split, anim, split_offset;
 } ParticleBillboardData;
 
+typedef struct ParticleCollisionElement
+{
+       Object *ob;
+
+       /* pointers to original data */
+       float *x[4], *v[4];
+
+       /* values interpolated from original data*/
+       float x0[3], x1[3], x2[3], p[3];
+       
+       /* results for found intersection point */
+       float nor[3], vel[3], uv[2];
+
+       /* count of original data (1-4) */
+       int tot;
+
+       /* flags for inversed normal / particle already inside element at start */
+       short inv_nor, inside;
+} ParticleCollisionElement;
+
 /* container for moving data between deflet_particle and particle_intersect_face */
 typedef struct ParticleCollision
 {
-       struct Object *ob, *hit_ob; // collided and current objects
-       struct CollisionModifierData *md, *hit_md; // collision modifiers for current and hit object;
-       float nor[3]; // normal at collision point
-       float vel[3]; // velocity of collision point
-       float co1[3], co2[3]; // ray start and end points
-       float ve1[3], ve2[3]; // particle velocities
-       float ray_len; // original length of co2-co1, needed for collision time evaluation
+       struct Object *current;
+       struct Object *hit;
+       struct Object *prev;
+       struct Object *skip;
+       struct Object *emitter;
+
+       struct CollisionModifierData *md; // collision modifier for current object;
+
        float f;        // time factor of previous collision, needed for substracting face velocity
-       float cfra; // start of the timestep (during frame change, since previous integer frame)
-       float dfra; // duration of timestep in frames
+       float fac1, fac2;
+
+       float cfra, old_cfra;
+
+       float original_ray_length; //original length of co2-co1, needed for collision time evaluation
+
+       int prev_index;
+
+       ParticleCollisionElement pce;
+
+       float total_time, inv_timestep;
+
+       float radius;
+       float co1[3], co2[3];
+       float ve1[3], ve2[3];
+
+       float acc[3], boid_z;
+
+       int boid;
 } ParticleCollision;
 
 typedef struct ParticleDrawData {
@@ -289,10 +327,8 @@ void psys_interpolate_face(struct MVert *mvert, struct MFace *mface, struct MTFa
 float psys_particle_value_from_verts(struct DerivedMesh *dm, short from, struct ParticleData *pa, float *values);
 void psys_get_from_key(struct ParticleKey *key, float *loc, float *vel, float *rot, float *time);
 
-/* only in edisparticle.c*/
-int psys_intersect_dm(struct Scene *scene, struct Object *ob, struct DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_uv, float *face_minmax, float *pa_minmax, float radius, float *ipoint);
 /* BLI_bvhtree_ray_cast callback */
-void particle_intersect_face(void *userdata, int index, const struct BVHTreeRay *ray, struct BVHTreeRayHit *hit);
+void BKE_psys_collision_neartest_cb(void *userdata, int index, const struct BVHTreeRay *ray, struct BVHTreeRayHit *hit);
 void psys_particle_on_dm(struct DerivedMesh *dm, int from, int index, int index_dmcache, float *fw, float foffset, float *vec, float *nor, float *utan, float *vtan, float *orco, float *ornor);
 
 /* particle_system.c */
index 9fb7f6408ac1d10661c65872e1cf33630c9e38a0..286cd3b249964d6366f6d472d9c7be2d20693fa3 100644 (file)
@@ -213,10 +213,8 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
                sub_v3_v3v3(ray_dir, col.co2, col.co1);
                mul_v3_fl(ray_dir, acbr->look_ahead);
                col.f = 0.0f;
-               col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-               col.dfra = bbd->dfra;
                hit.index = -1;
-               hit.dist = col.ray_len = len_v3(ray_dir);
+               hit.dist = col.original_ray_length = len_v3(ray_dir);
 
                /* find out closest deflector object */
                for(coll = bbd->sim->colliders->first; coll; coll=coll->next) {
@@ -224,18 +222,18 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
                        if(coll->ob == bpa->ground)
                                continue;
 
-                       col.ob = coll->ob;
+                       col.current = coll->ob;
                        col.md = coll->collmd;
 
                        if(col.md && col.md->bvhtree)
-                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
                }
                /* then avoid that object */
                if(hit.index>=0) {
-                       t = hit.dist/col.ray_len;
+                       t = hit.dist/col.original_ray_length;
 
                        /* avoid head-on collision */
-                       if(dot_v3v3(col.nor, pa->prev_state.ave) < -0.99) {
+                       if(dot_v3v3(col.pce.nor, pa->prev_state.ave) < -0.99) {
                                /* don't know why, but uneven range [0.0,1.0] */
                                /* works much better than even [-1.0,1.0] */
                                bbd->wanted_co[0] = BLI_frand();
@@ -243,7 +241,7 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
                                bbd->wanted_co[2] = BLI_frand();
                        }
                        else {
-                               VECCOPY(bbd->wanted_co, col.nor);
+                               copy_v3_v3(bbd->wanted_co, col.pce.nor);
                        }
 
                        mul_v3_fl(bbd->wanted_co, (1.0f - t) * val->personal_space * pa->size);
@@ -779,27 +777,26 @@ static Object *boid_find_ground(BoidBrainData *bbd, ParticleData *pa, float *gro
                /* first try to find below boid */
                copy_v3_v3(col.co1, pa->state.co);
                sub_v3_v3v3(col.co2, pa->state.co, zvec);
-               sub_v3_v3(col.co2, zvec);
                sub_v3_v3v3(ray_dir, col.co2, col.co1);
                col.f = 0.0f;
-               col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-               col.dfra = bbd->dfra;
                hit.index = -1;
-               hit.dist = col.ray_len = len_v3(ray_dir);
+               hit.dist = col.original_ray_length = len_v3(ray_dir);
+               col.pce.inside = 0;
 
                for(coll = bbd->sim->colliders->first; coll; coll = coll->next){
-                       col.ob = coll->ob;
+                       col.current = coll->ob;
                        col.md = coll->collmd;
+                       col.fac1 = col.fac2 = 0.f;
 
                        if(col.md && col.md->bvhtree)
-                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
                }
                /* then use that object */
                if(hit.index>=0) {
-                       t = hit.dist/col.ray_len;
+                       t = hit.dist/col.original_ray_length;
                        interp_v3_v3v3(ground_co, col.co1, col.co2, t);
-                       normalize_v3_v3(ground_nor, col.nor);
-                       return col.hit_ob;
+                       normalize_v3_v3(ground_nor, col.pce.nor);
+                       return col.hit;
                }
 
                /* couldn't find below, so find upmost deflector object */
@@ -808,24 +805,22 @@ static Object *boid_find_ground(BoidBrainData *bbd, ParticleData *pa, float *gro
                sub_v3_v3(col.co2, zvec);
                sub_v3_v3v3(ray_dir, col.co2, col.co1);
                col.f = 0.0f;
-               col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-               col.dfra = bbd->dfra;
                hit.index = -1;
-               hit.dist = col.ray_len = len_v3(ray_dir);
+               hit.dist = col.original_ray_length = len_v3(ray_dir);
 
                for(coll = bbd->sim->colliders->first; coll; coll = coll->next){
-                       col.ob = coll->ob;
+                       col.current = coll->ob;
                        col.md = coll->collmd;
 
                        if(col.md && col.md->bvhtree)
-                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
                }
                /* then use that object */
                if(hit.index>=0) {
-                       t = hit.dist/col.ray_len;
+                       t = hit.dist/col.original_ray_length;
                        interp_v3_v3v3(ground_co, col.co1, col.co2, t);
-                       normalize_v3_v3(ground_nor, col.nor);
-                       return col.hit_ob;
+                       normalize_v3_v3(ground_nor, col.pce.nor);
+                       return col.hit;
                }
 
                /* default to z=0 */
index 5f7c62042ec6c0cd00bbc5a46e1d96bcd4db0b1d..b5eb9b6ab29125001de01634a85740763cfd0fec 100644 (file)
@@ -2607,466 +2607,691 @@ static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, f
 /************************************************/
 /*                     Collisions                                                      */
 /************************************************/
-/* convert from triangle barycentric weights to quad mean value weights */
-static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
+#define COLLISION_MAX_COLLISIONS       10
+#define COLLISION_MIN_RADIUS 0.001f
+#define COLLISION_MIN_DISTANCE 0.0001f
+#define COLLISION_ZERO 0.00001f
+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)
 {
-       float co[3], vert[4][3];
+       float p0[3], e1[3], e2[3], d;
 
-       VECCOPY(vert[0], v1);
-       VECCOPY(vert[1], v2);
-       VECCOPY(vert[2], v3);
-       VECCOPY(vert[3], v4);
+       sub_v3_v3v3(e1, pce->x1, pce->x0);
+       sub_v3_v3v3(e2, pce->x2, pce->x0);
+       sub_v3_v3v3(p0, p, pce->x0);
 
-       co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
-       co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
-       co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];
+       cross_v3_v3v3(nor, e1, e2);
+       normalize_v3(nor);
 
-       interp_weights_poly_v3( w,vert, 4, co);
-}
+       d = dot_v3v3(p0, nor);
+
+       if(pce->inv_nor == -1) {
+               if(d < 0.f)
+                       pce->inv_nor = 1;
+               else
+                       pce->inv_nor = 0;
+       }
+
+       if(pce->inv_nor == 1) {
+               mul_v3_fl(nor, -1.f);
+               d = -d;
+       }
 
-/* check intersection with a derivedmesh */
-int psys_intersect_dm(Scene *scene, Object *ob, DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_w,
-                                                 float *face_minmax, float *pa_minmax, float radius, float *ipoint)
+       return d - radius;
+}
+static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *nor)
 {
-       MFace *mface=0;
-       MVert *mvert=0;
-       int i, totface, intersect=0;
-       float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
-       float cur_ipoint[3];
-       
-       if(dm==0){
-               psys_disable_all(ob);
+       float v0[3], v1[3], v2[3], c[3];
 
-               dm=mesh_get_derived_final(scene, ob, 0);
-               if(dm==0)
-                       dm=mesh_get_derived_deform(scene, ob, 0);
+       sub_v3_v3v3(v0, pce->x1, pce->x0);
+       sub_v3_v3v3(v1, p, pce->x0);
+       sub_v3_v3v3(v2, p, pce->x1);
 
-               psys_enable_all(ob);
+       cross_v3_v3v3(c, v1, v2);
 
-               if(dm==0)
-                       return 0;
+       return fabs(len_v3(c)/len_v3(v0)) - radius;
+}
+static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *nor)
+{
+       return len_v3v3(p, pce->x0) - radius;
+}
+static void collision_interpolate_element(ParticleCollisionElement *pce, float t, float fac, ParticleCollision *col)
+{
+       /* t is the current time for newton rhapson */
+       /* fac is the starting factor for current collision iteration */
+       /* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
+       float f = fac + t*(1.f-fac);
+       float mul = col->fac1 + f * (col->fac2-col->fac1);
+       if(pce->tot > 0) {
+               madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
+
+               if(pce->tot > 1) {
+                       madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
+
+                       if(pce->tot > 2)
+                               madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
+               }
        }
+}
+static void collision_point_velocity(ParticleCollisionElement *pce)
+{
+       float v[3];
 
-       
+       copy_v3_v3(pce->vel, pce->v[0]);
 
-       if(pa_minmax==0){
-               INIT_MINMAX(p_min,p_max);
-               DO_MINMAX(co1,p_min,p_max);
-               DO_MINMAX(co2,p_min,p_max);
+       if(pce->tot > 1) {
+               sub_v3_v3v3(v, pce->v[1], pce->v[0]);
+               madd_v3_v3fl(pce->vel, v, pce->uv[0]);
+
+               if(pce->tot > 2) {
+                       sub_v3_v3v3(v, pce->v[2], pce->v[0]);
+                       madd_v3_v3fl(pce->vel, v, pce->uv[1]);
+               }
        }
-       else{
-               VECCOPY(p_min,pa_minmax);
-               VECCOPY(p_max,pa_minmax+3);
+}
+static float collision_point_distance_with_normal(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
+{
+       if(fac >= 0.f)
+               collision_interpolate_element(pce, 0.f, fac, col);
+
+       switch(pce->tot) {
+               case 1:
+               {
+                       sub_v3_v3v3(nor, p, pce->x0);
+                       return normalize_v3(nor);
+               }
+               case 2:
+               {
+                       float u, e[3], vec[3];
+                       sub_v3_v3v3(e, pce->x1, pce->x0);
+                       sub_v3_v3v3(vec, p, pce->x0);
+                       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+                       madd_v3_v3v3fl(nor, vec, e, -u);
+                       return normalize_v3(nor);
+               }
+               case 3:
+                       return nr_signed_distance_to_plane(p, 0.f, pce, nor);
        }
+       return 0;
+}
+static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
+{
+       collision_interpolate_element(pce, 0.f, fac, col);
 
-       totface=dm->getNumFaces(dm);
-       mface=dm->getFaceDataArray(dm,CD_MFACE);
-       mvert=dm->getVertDataArray(dm,CD_MVERT);
-       
-       /* lets intersect the faces */
-       for(i=0; i<totface; i++,mface++){
-               if(vert_cos){
-                       VECCOPY(v1,vert_cos+3*mface->v1);
-                       VECCOPY(v2,vert_cos+3*mface->v2);
-                       VECCOPY(v3,vert_cos+3*mface->v3);
-                       if(mface->v4)
-                               VECCOPY(v4,vert_cos+3*mface->v4)
+       switch(pce->tot) {
+               case 1:
+               {
+                       sub_v3_v3v3(co, p, pce->x0);
+                       normalize_v3(co);
+                       madd_v3_v3v3fl(co, pce->x0, co, col->radius);
+                       break;
                }
-               else{
-                       VECCOPY(v1,mvert[mface->v1].co);
-                       VECCOPY(v2,mvert[mface->v2].co);
-                       VECCOPY(v3,mvert[mface->v3].co);
-                       if(mface->v4)
-                               VECCOPY(v4,mvert[mface->v4].co)
+               case 2:
+               {
+                       float u, e[3], vec[3], nor[3];
+                       sub_v3_v3v3(e, pce->x1, pce->x0);
+                       sub_v3_v3v3(vec, p, pce->x0);
+                       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+                       madd_v3_v3v3fl(nor, vec, e, -u);
+                       normalize_v3(nor);
+
+                       madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
+                       madd_v3_v3fl(co, nor, col->radius);
+                       break;
                }
+               case 3:
+               {
+                               float p0[3], e1[3], e2[3], nor[3];
 
-               if(face_minmax==0){
-                       INIT_MINMAX(min,max);
-                       DO_MINMAX(v1,min,max);
-                       DO_MINMAX(v2,min,max);
-                       DO_MINMAX(v3,min,max);
-                       if(mface->v4)
-                               DO_MINMAX(v4,min,max)
-                       if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
-                               continue;
+                               sub_v3_v3v3(e1, pce->x1, pce->x0);
+                               sub_v3_v3v3(e2, pce->x2, pce->x0);
+                               sub_v3_v3v3(p0, p, pce->x0);
+
+                               cross_v3_v3v3(nor, e1, e2);
+                               normalize_v3(nor);
+
+                               if(pce->inv_nor == 1)
+                                       mul_v3_fl(nor, -1.f);
+
+                               madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
+                               madd_v3_v3fl(co, e1, pce->uv[0]);
+                               madd_v3_v3fl(co, e2, pce->uv[1]);
+                       break;
                }
-               else{
-                       VECCOPY(min, face_minmax+6*i);
-                       VECCOPY(max, face_minmax+6*i+3);
-                       if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
-                               continue;
+       }
+}
+/* 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];
+       int iter;
+
+       pce->inv_nor = -1;
+
+       /* 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;
+       d1 = 0.f;
+
+       for(iter=0; iter<10; iter++) {//, itersum++) {
+               /* get current location */
+               collision_interpolate_element(pce, t1, col->f, col);
+               interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
+
+               d1 = distance_func(pce->p, radius, pce, n);
+
+               /* no movement, so no collision */
+               if(d1 == d0) {
+                       return -1.f;
                }
 
-               if(radius>0.0f){
-                       if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
-                               if(cur_d<*min_d){
-                                       *min_d=cur_d;
-                                       VECCOPY(ipoint,cur_ipoint);
-                                       *min_face=i;
-                                       intersect=1;
-                               }
-                       }
-                       if(mface->v4){
-                               if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
-                                       if(cur_d<*min_d){
-                                               *min_d=cur_d;
-                                               VECCOPY(ipoint,cur_ipoint);
-                                               *min_face=i;
-                                               intersect=1;
-                                       }
-                               }
-                       }
+               /* particle already inside face, so report collision */
+               if(iter == 0 && d0 < 0.f && d0 > -radius) {
+                       copy_v3_v3(pce->p, col->co1);
+                       copy_v3_v3(pce->nor, n);
+                       pce->inside = 1;
+                       return 0.f;
                }
-               else{
-                       if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
-                               if(cur_d<*min_d){
-                                       *min_d=cur_d;
-                                       min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
-                                       min_w[1]= cur_uv[0];
-                                       min_w[2]= cur_uv[1];
-                                       min_w[3]= 0.0f;
-                                       if(mface->v4)
-                                               intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
-                                       *min_face=i;
-                                       intersect=1;
-                               }
-                       }
-                       if(mface->v4){
-                               if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
-                                       if(cur_d<*min_d){
-                                               *min_d=cur_d;
-                                               min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
-                                               min_w[1]= 0.0f;
-                                               min_w[2]= cur_uv[0];
-                                               min_w[3]= cur_uv[1];
-                                               intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
-                                               *min_face=i;
-                                               intersect=1;
-                                       }
-                               }
+               
+               dd = (t1-t0)/(d1-d0);
+
+               t0 = t1;
+               d0 = d1;
+
+               t1 -= d1*dd;
+
+               /* particle movin away from plane could also mean a strangely rotating face, so check from end */
+               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;
+                       d1 = 0.f;
+
+                       continue;
+               }
+               else if(iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
+                       return -1.f;
+
+               if(d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
+                       if(t1 >= -COLLISION_ZERO && t1 <= 1.f) {
+                               if(distance_func == nr_signed_distance_to_plane)
+                                       copy_v3_v3(pce->nor, n);
+
+                               CLAMP(t1, 0.f, 1.f);
+
+                               return t1;
                        }
+                       else
+                               return -1.f;
                }
        }
-       return intersect;
+       return -1.0;
+}
+static int collision_sphere_to_tri(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
+{
+       ParticleCollisionElement *result = &col->pce;
+       float ct, u, v;
+
+       pce->inv_nor = -1;
+       pce->inside = 0;
+
+       ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
+
+       if(ct >= 0.f && ct < *t && (result->inside==0 || pce->inside==1) ) {
+               float e1[3], e2[3], p0[3];
+               float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
+
+               sub_v3_v3v3(e1, pce->x1, pce->x0);
+               sub_v3_v3v3(e2, pce->x2, pce->x0);
+               /* XXX: add radius correction here? */
+               sub_v3_v3v3(p0, pce->p, pce->x0);
+
+               e1e1 = dot_v3v3(e1, e1);
+               e1e2 = dot_v3v3(e1, e2);
+               e1p0 = dot_v3v3(e1, p0);
+               e2e2 = dot_v3v3(e2, e2);
+               e2p0 = dot_v3v3(e2, p0);
+
+               inv = 1.f/(e1e1 * e2e2 - e1e2 * e1e2);
+               u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
+               v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
+
+               if(u>=0.f && u<=1.f && v>=0.f && u+v<=1.f) {
+                       *result = *pce;
+
+                       /* normal already calculated in pce */
+
+                       result->uv[0] = u;
+                       result->uv[1] = v;
+
+                       *t = ct;
+                       return 1;
+               }
+       }
+       return 0;
+}
+static int collision_sphere_to_edges(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t, int quad)
+{
+       ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
+       ParticleCollisionElement *result = &col->pce;
+
+       float ct;
+       int i;
+
+       for(i=0; i<3; i++) {
+               /* in case of a quad, no need to check "edge" that goes through face */
+               if((pce->x[3] && i==2) || (quad && i==0))
+                       continue;
+
+               cur = edge+i;
+               cur->x[0] = pce->x[i]; cur->x[1] = pce->x[(i+1)%3];
+               cur->v[0] = pce->v[i]; cur->v[1] = pce->v[(i+1)%3];
+               cur->tot = 2;
+               cur->inside = 0;
+
+               ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
+
+               if(ct >= 0.f && ct < *t) {
+                       float u, e[3], vec[3];
+
+                       sub_v3_v3v3(e, cur->x1, cur->x0);
+                       sub_v3_v3v3(vec, cur->p, cur->x0);
+                       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+                       if(u < 0.f || u > 1.f)
+                               break;
+
+                       *result = *cur;
+
+                       madd_v3_v3v3fl(result->nor, vec, e, -u);
+                       normalize_v3(result->nor);
+
+                       result->uv[0] = u;
+
+                       
+                       hit = cur;
+                       *t = ct;
+               }
+
+       }
+
+       return hit != NULL;
 }
+static int collision_sphere_to_verts(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
+{
+       ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
+       ParticleCollisionElement *result = &col->pce;
+
+       float ct;
+       int i;
+
+       for(i=0; i<3; i++) {
+               /* in case of quad, only check one vert the first time */
+               if(pce->x[3] && i != 1)
+                       continue;
+
+               cur = vert+i;
+               cur->x[0] = pce->x[i];
+               cur->v[0] = pce->v[i];
+               cur->tot = 1;
+               cur->inside = 0;
+
+               ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
+               
+               if(ct >= 0.f && ct < *t) {
+                       *result = *cur;
+
+                       sub_v3_v3v3(result->nor, cur->p, cur->x0);
+                       normalize_v3(result->nor);
+
+                       hit = cur;
+                       *t = ct;
+               }
 
-void particle_intersect_face(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
+       }
+
+       return hit != NULL;
+}
+/* Callback for BVHTree near test */
+void BKE_psys_collision_neartest_cb(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
 {
        ParticleCollision *col = (ParticleCollision *) userdata;
+       ParticleCollisionElement pce;
        MFace *face = col->md->mfaces + index;
        MVert *x = col->md->x;
        MVert *v = col->md->current_v;
-       float vel[3], co1[3], co2[3], uv[2], ipoint[3], temp[3], t;
-       float x0[3], x1[3], x2[3], x3[3];
-       float *t0=x0, *t1=x1, *t2=x2, *t3=(face->v4 ? x3 : NULL);
-
-       /* move collision face to start of timestep */
-       madd_v3_v3v3fl(t0, x[face->v1].co, v[face->v1].co, col->cfra);
-       madd_v3_v3v3fl(t1, x[face->v2].co, v[face->v2].co, col->cfra);
-       madd_v3_v3v3fl(t2, x[face->v3].co, v[face->v3].co, col->cfra);
-       if(t3)
-               madd_v3_v3v3fl(t3, x[face->v4].co, v[face->v4].co, col->cfra);
-
-       /* calculate average velocity of face */
-       copy_v3_v3(vel, v[ face->v1 ].co);
-       add_v3_v3(vel, v[ face->v2 ].co);
-       add_v3_v3(vel, v[ face->v3 ].co);
-       mul_v3_fl(vel, 0.33334f*col->dfra);
-
-       /* substract face velocity, in other words convert to 
-          a coordinate system where only the particle moves */
-       madd_v3_v3v3fl(co1, col->co1, vel, -col->f);
-       sub_v3_v3v3(co2, col->co2, vel);
+       float t = hit->dist/col->original_ray_length;
+       float uv[2] = {0.f,0.f};
+       int collision = 0, quad = 0;
+
+       pce.x[0] = x[face->v1].co;
+       pce.x[1] = x[face->v2].co;
+       pce.x[2] = x[face->v3].co;
+       pce.x[3] = face->v4 ? x[face->v4].co : NULL;
+
+       pce.v[0] = v[face->v1].co;
+       pce.v[1] = v[face->v2].co;
+       pce.v[2] = v[face->v3].co;
+       pce.v[3] = face->v4 ? v[face->v4].co : NULL;
+
+       pce.tot = 3;
+       pce.inside = 0;
 
        do
        {       
-               if(ray->radius == 0.0f) {
-                       if(isect_line_tri_v3(co1, co2, t0, t1, t2, &t, uv)) {
-                               if(t >= 0.0f && t < hit->dist/col->ray_len) {
-                                       hit->dist = col->ray_len * t;
-                                       hit->index = index;
-
-                                       /* calculate normal that's facing the particle */
-                                       normal_tri_v3( col->nor,t0, t1, t2);
-                                       VECSUB(temp, co2, co1);
-                                       if(dot_v3v3(col->nor, temp) > 0.0f)
-                                               negate_v3(col->nor);
-
-                                       VECCOPY(col->vel,vel);
-
-                                       col->hit_ob = col->ob;
-                                       col->hit_md = col->md;
-                               }
-                       }
+               collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
+               if(col->pce.inside == 0) {
+                       collision += collision_sphere_to_edges(col, ray->radius, &pce, &t, quad);
+                       collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
                }
-               else {
-                       if(isect_sweeping_sphere_tri_v3(co1, co2, ray->radius, t0, t1, t2, &t, ipoint)) {
-                               if(t >=0.0f && t < hit->dist/col->ray_len) {
-                                       hit->dist = col->ray_len * t;
-                                       hit->index = index;
-
-                                       interp_v3_v3v3(temp, co1, co2, t);
-                                       
-                                       VECSUB(col->nor, temp, ipoint);
-                                       normalize_v3(col->nor);
 
-                                       VECCOPY(col->vel,vel);
+               if(collision) {
+                       hit->dist = col->original_ray_length * t;
+                       hit->index = index;
+                               
+                       collision_point_velocity(&col->pce);
 
-                                       col->hit_ob = col->ob;
-                                       col->hit_md = col->md;
-                               }
-                       }
+                       col->hit = col->current;
                }
 
-               t1 = t2;
-               t2 = t3;
-               t3 = NULL;
+               pce.x[1] = pce.x[2];
+               pce.x[2] = pce.x[3];
+               pce.x[3] = NULL;
 
-       } while(t2);
+               pce.v[1] = pce.v[2];
+               pce.v[2] = pce.v[3];
+               pce.v[3] = NULL;
+               quad++;
+
+       } while(pce.x[2]);
 }
-/* Particle - Mesh collision code
- * Features:
- * - point and swept sphere to mesh surface collisions
- * - moving colliders (but not yet rotating or deforming colliders)
- * - friction & damping
- * - angular momentum <-> linear momentum
- * - high accuracy by re-applying particle acceleration after collision
- * - behaves relatively well even if limit of 10 collisions per simulation step is exceeded
- * Main parts:
- * 1. check for all possible deflectors for closest intersection on particle path
- * 2. if deflection was found calculate new coordinates or kill the particle
- */
-static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, float cfra){
-       Object *ground_ob = NULL;
-       ParticleSettings *part = sim->psys->part;
-       ParticleData *pa = sim->psys->particles + p;
-       ParticleCollision col;
+static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, ListBase *colliders)
+{
        ColliderCache *coll;
-       BVHTreeRayHit hit;
-       float ray_dir[3], acc[3];
-       float radius = ((part->flag & PART_SIZE_DEFL)?pa->size:0.0f), boid_z = 0.0f;
-       float timestep = psys_get_timestep(sim) * dfra;
-       float inv_timestep = 1.0f/timestep;
-       int deflections=0, max_deflections=10;
+       float ray_dir[3];
 
-       /* get acceleration (from gravity, forcefields etc. to be re-applied after collision) */
-       sub_v3_v3v3(acc, pa->state.vel, pa->prev_state.vel);
-       mul_v3_fl(acc, inv_timestep);
+       if(colliders->first == NULL)
+               return 0;
 
-       /* set values for first iteration */
-       copy_v3_v3(col.co1, pa->prev_state.co);
-       copy_v3_v3(col.co2, pa->state.co);
-       copy_v3_v3(col.ve1, pa->prev_state.vel);
-       copy_v3_v3(col.ve2, pa->state.vel);
-       col.f = 0.0f;
+       sub_v3_v3v3(ray_dir, col->co2, col->co1);
+       hit->index = -1;
+       hit->dist = col->original_ray_length = len_v3(ray_dir);
+       col->pce.inside = 0;
 
-       /* override for boids */
-       if(part->phystype == PART_PHYS_BOIDS) {
-               BoidParticle *bpa = pa->boid;
-               radius = pa->size;
-               boid_z = pa->state.co[2];
-               ground_ob = bpa->ground;
+       /* even if particle is stationary we want to check for moving colliders */
+       /* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
+       if(hit->dist == 0.0f)
+               hit->dist = col->original_ray_length = 0.000001f;
+
+       for(coll = colliders->first; coll; coll=coll->next){
+               /* for boids: don't check with current ground object */
+               if(coll->ob == col->skip)
+                       continue;
+
+               /* particles should not collide with emitter at birth */
+               if(coll->ob == col->emitter && pa->time < col->cfra && pa->time >= col->old_cfra)
+                       continue;
+
+               col->current = coll->ob;
+               col->md = coll->collmd;
+               col->fac1 = (col->old_cfra - coll->collmd->time_x) / (coll->collmd->time_xnew - coll->collmd->time_x);
+               col->fac2 = (col->cfra - coll->collmd->time_x) / (coll->collmd->time_xnew - coll->collmd->time_x);
+
+               if(col->md && col->md->bvhtree)
+                       BLI_bvhtree_ray_cast(col->md->bvhtree, col->co1, ray_dir, col->radius, hit, BKE_psys_collision_neartest_cb, col);
        }
 
-       /* 10 iterations to catch multiple deflections */
-       if(sim->colliders) while(deflections < max_deflections){
-               /* 1. */
+       return hit->index >= 0;
+}
+static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, int kill, int dynamic_rotation)
+{
+       ParticleCollisionElement *pce = &col->pce;
+       PartDeflect *pd = col->hit->pd;
+       float co[3];                                                                            /* point of collision */
+       float x = hit->dist/col->original_ray_length;           /* location factor of collision between this iteration */
+       float f = col->f + x * (1.0f - col->f);                         /* time factor of collision between timestep */
+       float dt1 = (f - col->f) * col->total_time;                     /* time since previous collision (in seconds) */
+       float dt2 = (1.0f - f) * col->total_time;                       /* time left after collision (in seconds) */
+       int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
+
+       /* calculate exact collision location */
+       interp_v3_v3v3(co, col->co1, col->co2, x);
+
+       /* particle dies in collision */
+       if(through == 0 && (kill || pd->flag & PDEFLE_KILL_PART)) {
+               pa->alive = PARS_DYING;
+               pa->dietime = col->old_cfra + (col->cfra - col->old_cfra) * f;
+
+               copy_v3_v3(pa->state.co, co);
+               interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
+               interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
+               interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
+
+               /* particle is dead so we don't need to calculate further */
+               return 0;
+       }
+       /* figure out velocity and other data after collision */
+       else {
+               float v0[3];    /* velocity directly before collision to be modified into velocity directly after collision */
+               float v0_nor[3];/* normal component of v0 */
+               float v0_tan[3];/* tangential component of v0 */
+               float vc_tan[3];/* tangential component of collision surface velocity */
+               float v0_dot, vc_dot;
+               float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
+               float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
+               float distance, nor[3], dot;
+
+               CLAMP(damp,0.0,1.0);
+               CLAMP(frict,0.0,1.0);
+
+               /* get exact velocity right before collision */
+               madd_v3_v3v3fl(v0, col->ve1, col->acc, dt1);
+                               
+               /* convert collider velocity from 1/framestep to 1/s TODO: here we assume 1 frame step for collision modifier */
+               mul_v3_fl(pce->vel, col->inv_timestep);
+
+               /* calculate tangential particle velocity */
+               v0_dot = dot_v3v3(pce->nor, v0);
+               madd_v3_v3v3fl(v0_tan, v0, pce->nor, -v0_dot);
+
+               /* calculate tangential collider velocity */
+               vc_dot = dot_v3v3(pce->nor, pce->vel);
+               madd_v3_v3v3fl(vc_tan, pce->vel, pce->nor, -vc_dot);
+
+               /* handle friction effects (tangential and angular velocity) */
+               if(frict > 0.0f) {
+                       /* angular <-> linear velocity */
+                       if(dynamic_rotation) {
+                               float vr_tan[3], v1_tan[3], ave[3];
+                                       
+                               /* linear velocity of particle surface */
+                               cross_v3_v3v3(vr_tan, pce->nor, pa->state.ave);
+                               mul_v3_fl(vr_tan, pa->size);
 
-               sub_v3_v3v3(ray_dir, col.co2, col.co1);
-               hit.index = -1;
-               hit.dist = col.ray_len = len_v3(ray_dir);
+                               /* change to coordinates that move with the collision plane */
+                               sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
+                                               
+                               /* The resulting velocity is a weighted average of particle cm & surface
+                                       * velocity. This weight (related to particle's moment of inertia) could
+                                       * be made a parameter for angular <-> linear conversion.
+                                       */
+                               madd_v3_v3fl(v1_tan, vr_tan, -0.4);
+                               mul_v3_fl(v1_tan, 1.0f/1.4f); /* 1/(1+0.4) */
 
-               col.cfra = fmod(cfra-dfra, 1.0f);
-               col.dfra = dfra;
+                               /* rolling friction is around 0.01 of sliding friction (could be made a parameter) */
+                               mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
 
-               /* even if particle is stationary we want to check for moving colliders */
-               /* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
-               if(hit.dist == 0.0f)
-                       hit.dist = col.ray_len = 0.000001f;
+                               /* surface_velocity is opposite to cm velocity */
+                               mul_v3_v3fl(vr_tan, v1_tan, -1.0f);
 
-               for(coll = sim->colliders->first; coll; coll=coll->next){
-                       /* for boids: don't check with current ground object */
-                       if(coll->ob == ground_ob)
-                               continue;
+                               /* get back to global coordinates */
+                               add_v3_v3(v1_tan, vc_tan);
 
-                       /* particles should not collide with emitter at birth */
-                       if(coll->ob == sim->ob && pa->time < cfra && pa->time >= sim->psys->cfra)
-                               continue;
+                               /* convert to angular velocity*/
+                               cross_v3_v3v3(ave, vr_tan, pce->nor);
+                               mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001));
 
-                       col.ob = coll->ob;
-                       col.md = coll->collmd;
+                               /* only friction will cause change in linear & angular velocity */
+                               interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
+                               interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
+                       }
+                       else {
+                               /* just basic friction (unphysical due to the friction model used in Blender) */
+                               interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
+                       }
+               }
 
-                       if(col.md && col.md->bvhtree)
-                               BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+               /* stickness was possibly added before, so cancel that before calculating new normal velocity */
+               /* otherwise particles go flying out of the surface because of high reversed sticky velocity */
+               if(v0_dot < 0.0f) {
+                       v0_dot += pd->pdef_stickness;
+                       if(v0_dot > 0.0f)
+                               v0_dot = 0.0f;
                }
 
-               /* 2. */
-               if(hit.index>=0) {
-                       PartDeflect *pd = col.hit_ob->pd;
-                       float co[3];                                                    /* point of collision */
-                       float x = hit.dist/col.ray_len;                 /* location factor of collision between this iteration */
-                       float f = col.f + x * (1.0f - col.f);   /* time factor of collision between timestep */
-                       float dt1 = (f - col.f) * timestep;             /* time since previous collision (in seconds) */
-                       float dt2 = (1.0f - f) * timestep;              /* time left after collision (in seconds) */
-                       int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
-
-                       deflections++;
-                       
-                       interp_v3_v3v3(co, col.co1, col.co2, x);
-                       
-                       /* make sure we don't hit the current face again */
-                       /* TODO: could/should this be proportional to pa->size? */
-                       madd_v3_v3fl(co, col.nor, (through ? -0.0001f : 0.0001f));
+               /* damping and flipping of velocity around normal */
+               v0_dot *= 1.0f - damp;
+               vc_dot *= through ? damp : 1.0f;
 
-                       /* particle dies in collision */
-                       if(through == 0 && (part->flag & PART_DIE_ON_COL || pd->flag & PDEFLE_KILL_PART)) {
-                               pa->alive = PARS_DYING;
-                               pa->dietime = sim->psys->cfra + (cfra - sim->psys->cfra) * f;
+               /* calculate normal particle velocity */
+               /* special case for object hitting the particle from behind */
+               if(through==0 && ((vc_dot>0.0f && v0_dot>0.0f && vc_dot>v0_dot) || (vc_dot<0.0f && v0_dot<0.0f && vc_dot<v0_dot)))
+                       mul_v3_v3fl(v0_nor, pce->nor, vc_dot);
+               else if(v0_dot > 0.f)
+                       mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? -1.0f : 1.0f) * v0_dot);
+               else
+                       mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
 
-                               copy_v3_v3(pa->state.co, co);
-                               interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
-                               interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
-                               interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
+               /* combine components together again */
+               add_v3_v3v3(v0, v0_nor, v0_tan);
 
-                               /* particle is dead so we don't need to calculate further */
-                               return;
+               if(col->boid) {
+                       /* keep boids above ground */
+                       BoidParticle *bpa = pa->boid;
+                       if(bpa->data.mode == eBoidMode_OnLand || co[2] <= col->boid_z) {
+                               co[2] = col->boid_z;
+                               v0[2] = 0.0f;
                        }
-                       /* figure out velocity and other data after collision */
-                       else {
-                               float v0[3];    /* velocity directly before collision to be modified into velocity directly after collision */
-                               float v0_nor[3];/* normal component of v0 */
-                               float v0_tan[3];/* tangential component of v0 */
-                               float vc_tan[3];/* tangential component of collision surface velocity */
-                               float check[3];
-                               float v0_dot, vc_dot, check_dot;
-                               float damp, frict;
-
-                               /* get exact velocity right before collision */
-                               madd_v3_v3v3fl(v0, col.ve1, acc, dt1);
-                               
-                               /* convert collider velocity from 1/framestep to 1/s */
-                               mul_v3_fl(col.vel, inv_timestep);
-                               
-                               /* get damping & friction factors */
-                               damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
-                               CLAMP(damp,0.0,1.0);
+               }
+               
+               /* re-apply acceleration to final location and velocity */
+               madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
+               madd_v3_v3fl(pa->state.co, col->acc, 0.5f*dt2*dt2);
+               madd_v3_v3v3fl(pa->state.vel, v0, col->acc, dt2);
+
+               /* make sure particle stays on the right side of the surface */
+               if(!through) {
+                       distance = collision_point_distance_with_normal(co, pce, -1.f, col, nor);
+                       
+                       if(distance < col->radius + COLLISION_MIN_DISTANCE)
+                               madd_v3_v3fl(co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
 
-                               frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
-                               CLAMP(frict,0.0,1.0);
+                       dot = dot_v3v3(nor, v0);
+                       if(dot < 0.f)
+                               madd_v3_v3fl(v0, nor, -dot);
 
-                               /* treat normal & tangent components separately */
-                               v0_dot = dot_v3v3(col.nor, v0);
-                               madd_v3_v3v3fl(v0_tan, v0, col.nor, -v0_dot);
+                       distance = collision_point_distance_with_normal(pa->state.co, pce, 1.f, col, nor);
 
-                               vc_dot = dot_v3v3(col.nor, col.vel);
-                               madd_v3_v3v3fl(vc_tan, col.vel, col.nor, -vc_dot);
+                       if(distance < col->radius + COLLISION_MIN_DISTANCE)
+                               madd_v3_v3fl(pa->state.co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
 
-                               /* handle friction effects (tangential and angular velocity) */
-                               if(frict > 0.0f) {
-                                       /* angular <-> linear velocity */
-                                       if(part->flag & PART_ROT_DYN) {
-                                               float vr_tan[3], v1_tan[3], ave[3];
-                                       
-                                               /* linear velocity of particle surface */
-                                               cross_v3_v3v3(vr_tan, col.nor, pa->state.ave);
-                                               mul_v3_fl(vr_tan, pa->size);
+                       dot = dot_v3v3(nor, pa->state.vel);
+                       if(dot < 0.f)
+                               madd_v3_v3fl(pa->state.vel, nor, -dot);
+               }
 
-                                               /* change to coordinates that move with the collision plane */
-                                               sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
-                                               
-                                               /* The resulting velocity is a weighted average of particle cm & surface
-                                                * velocity. This weight (related to particle's moment of inertia) could
-                                                * be made a parameter for angular <-> linear conversion.
-                                                */
-                                               madd_v3_v3fl(v1_tan, vr_tan, -0.4);
-                                               mul_v3_fl(v1_tan, 1.0f/1.4f); /* 1/(1+0.4) */
+               /* add stickness to surface */
+               madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);
 
-                                               /* rolling friction is around 0.01 of sliding friction (could be made a parameter) */
-                                               mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
+               /* set coordinates for next iteration */
+               copy_v3_v3(col->co1, co);
+               copy_v3_v3(col->co2, pa->state.co);
 
-                                               /* surface_velocity is opposite to cm velocity */
-                                               mul_v3_v3fl(vr_tan, v1_tan, -1.0f);
+               copy_v3_v3(col->ve1, v0);
+               copy_v3_v3(col->ve2, pa->state.vel);
 
-                                               /* get back to global coordinates */
-                                               add_v3_v3(v1_tan, vc_tan);
+               col->f = f;
+       }
 
-                                               /* convert to angular velocity*/
-                                               cross_v3_v3v3(ave, vr_tan, col.nor);
-                                               mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001));
+       col->prev = col->hit;
+       col->prev_index = hit->index;
 
-                                               /* only friction will cause change in linear & angular velocity */
-                                               interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
-                                               interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
-                                       }
-                                       else {
-                                               /* just basic friction (unphysical due to the friction model used in Blender) */
-                                               interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
-                                       }
-                               }
+       return 1;
+}
+static void collision_fail(ParticleData *pa, ParticleCollision *col)
+{
+       /* final chance to prevent total failure, so stick to the surface and hope for the best */
+       collision_point_on_surface(col->co1, &col->pce, 1.f, col, pa->state.co);
 
-                               /* stickness was possibly added before, so cancel that before calculating new normal velocity */
-                               /* otherwise particles go flying out of the surface because of high reversed sticky velocity */
-                               if(v0_dot < 0.0f) {
-                                       v0_dot += pd->pdef_stickness;
-                                       if(v0_dot > 0.0f)
-                                               v0_dot = 0.0f;
-                               }
+       copy_v3_v3(pa->state.vel, col->pce.vel);
+       mul_v3_fl(pa->state.vel, col->inv_timestep);
 
-                               /* damping and flipping of velocity around normal */
-                               v0_dot *= 1.0f - damp;
-                               vc_dot *= through ? damp : 1.0f;
 
-                               /* special case for object hitting the particle from behind */
-                               if(through==0 && ((vc_dot>0.0f && v0_dot>0.0f && vc_dot>v0_dot) || (vc_dot<0.0f && v0_dot<0.0f && vc_dot<v0_dot)))
-                                       mul_v3_v3fl(v0_nor, col.nor, vc_dot);
-                               else
-                                       mul_v3_v3fl(v0_nor, col.nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
+       /* printf("max iterations\n"); */
+}
 
-                               /* combine components together again */
-                               add_v3_v3v3(v0, v0_nor, v0_tan);
+/* Particle - Mesh collision detection and response
+ * Features:
+ * -friction and damping
+ * -angular momentum <-> linear momentum
+ * -high accuracy by re-applying particle acceleration after collision
+ * -handles moving, rotating and deforming meshes
+ * -uses Newton-Rhapson iteration to find the collisions
+ * -handles spherical particles and (nearly) point like particles
+ */
+static void collision_check(ParticleSimulationData *sim, int p, float dfra, float cfra){
+       Object *ground_ob = NULL;
+       ParticleSettings *part = sim->psys->part;
+       ParticleData *pa = sim->psys->particles + p;
+       ParticleCollision col;
+       BVHTreeRayHit hit;
+       int collision_count=0;
 
-                               /* keep boids above ground */
-                               if(part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
-                                       BoidParticle *bpa = pa->boid;
-                                       if(bpa->data.mode == eBoidMode_OnLand || co[2] <= boid_z) {
-                                               co[2] = boid_z;
-                                               v0[2] = 0.0f;
-                                       }
-                               }
-                               
-                               if(deflections < max_deflections) {
-                                       /* re-apply acceleration to final velocity and location */
-                                       madd_v3_v3v3fl(pa->state.vel, v0, acc, dt2);
-                                       madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
-                                       madd_v3_v3fl(pa->state.co, acc, 0.5f*dt2*dt2);
+       float timestep = psys_get_timestep(sim);
 
-                                       /* make sure particle stays on the right side of the surface */
-                                       sub_v3_v3v3(check, pa->state.co, co);
-                                       /* (collision surface has moved during the time too) */
-                                       madd_v3_v3fl(check, col.vel, -dt2);
+       memset(&col, 0, sizeof(ParticleCollision));
 
-                                       check_dot = dot_v3v3(check, col.nor);
-                                       if((!through && check_dot < 0.0f) || (through && check_dot > 0.0f))
-                                               madd_v3_v3fl(pa->state.co, col.nor, (through ? -0.0001f : 0.0001f) - check_dot);
+       col.total_time = timestep * dfra;
+       col.inv_timestep = 1.0f/timestep;
 
-                                       /* Stickness to surface */
-                                       madd_v3_v3fl(pa->state.vel, col.nor, -pd->pdef_stickness);
+       col.cfra = cfra;
+       col.old_cfra = sim->psys->cfra;
 
-                                       /* set coordinates for next iteration */
-                                       copy_v3_v3(col.co1, co);
-                                       copy_v3_v3(col.co2, pa->state.co);
+       /* get acceleration (from gravity, forcefields etc. to be re-applied in collision response) */
+       sub_v3_v3v3(col.acc, pa->state.vel, pa->prev_state.vel);
+       mul_v3_fl(col.acc, 1.f/col.total_time);
 
-                                       copy_v3_v3(col.ve1, v0);
-                                       copy_v3_v3(col.ve2, pa->state.vel);
+       /* set values for first iteration */
+       copy_v3_v3(col.co1, pa->prev_state.co);
+       copy_v3_v3(col.co2, pa->state.co);
+       copy_v3_v3(col.ve1, pa->prev_state.vel);
+       copy_v3_v3(col.ve2, pa->state.vel);
+       col.f = 0.0f;
 
-                                       col.f = f;
-                               }
-                               else {
-                                       /* final chance to prevent failure, so stick to the surface and hope for the best */
-                                       madd_v3_v3v3fl(pa->state.co, co, col.vel, dt2);
-                                       copy_v3_v3(pa->state.vel, v0);
-                               }
-                       }
+       col.radius = ((part->flag & PART_SIZE_DEFL) || (part->phystype == PART_PHYS_BOIDS)) ? pa->size : COLLISION_MIN_RADIUS;
+
+       /* override for boids */
+       if(part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
+               col.boid = 1;
+               col.boid_z = pa->state.co[2];
+               col.skip = pa->boid->ground;
+       }
+
+       /* 10 iterations to catch multiple collisions */
+       while(collision_count < COLLISION_MAX_COLLISIONS){
+               if(collision_detect(pa, &col, &hit, sim->colliders)) {
+                       
+                       collision_count++;
+
+                       if(collision_count == COLLISION_MAX_COLLISIONS)
+                               collision_fail(pa, &col);
+                       else if(collision_response(pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN)==0)
+                               return;
                }
                else
                        return;
@@ -3454,7 +3679,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
        
                                /* deflection */
                                if(sim->colliders)
-                                       deflect_particle(sim, p, pa->state.time, cfra);
+                                       collision_check(sim, p, pa->state.time, cfra);
 
                                /* rotations */
                                basic_rotate(part, pa, pa->state.time, timestep);
@@ -3473,7 +3698,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 
                                        /* deflection */
                                        if(sim->colliders)
-                                               deflect_particle(sim, p, pa->state.time, cfra);
+                                               collision_check(sim, p, pa->state.time, cfra);
                                }
                        }
                        break;
@@ -3494,7 +3719,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
                                sph_integrate(sim, pa, pa->state.time, gravity, springhash);
 
                                if(sim->colliders)
-                                       deflect_particle(sim, p, pa->state.time, cfra);
+                                       collision_check(sim, p, pa->state.time, cfra);
                                
                                /* SPH particles are not physical particles, just interpolation particles,  thus rotation has not a direct sense for them */    
                                basic_rotate(part, pa, pa->state.time, timestep);  
index 4c04db5f56fa9c6b10d8c90fe828cef0d3c82ebc..35aeb756274ddd1338229296ea5810f20f44345c 100644 (file)
@@ -4035,7 +4035,7 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
                        collmd->current_x = NULL;
                        collmd->current_xnew = NULL;
                        collmd->current_v = NULL;
-                       collmd->time = -1000;
+                       collmd->time_x = collmd->time_xnew = -1000;
                        collmd->numverts = 0;
                        collmd->bvhtree = NULL;
                        collmd->mfaces = NULL;
index 8a304765a7a10702d7931bf6caca49a2f3c0cc8d..d3e61e785e10c3dbf94591fe6d5dafbd379bf01f 100644 (file)
@@ -3131,6 +3131,149 @@ static void brush_smooth_do(PEData *data, float UNUSED(mat[][4]), float imat[][4
        (data->edit->points + point_index)->flag |= PEP_EDIT_RECALC;
 }
 
+/* convert from triangle barycentric weights to quad mean value weights */
+static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
+{
+       float co[3], vert[4][3];
+
+       VECCOPY(vert[0], v1);
+       VECCOPY(vert[1], v2);
+       VECCOPY(vert[2], v3);
+       VECCOPY(vert[3], v4);
+
+       co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
+       co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
+       co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];
+
+       interp_weights_poly_v3( w,vert, 4, co);
+}
+
+/* check intersection with a derivedmesh */
+static int particle_intersect_dm(Scene *scene, Object *ob, DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_w,
+                                                 float *face_minmax, float *pa_minmax, float radius, float *ipoint)
+{
+       MFace *mface=0;
+       MVert *mvert=0;
+       int i, totface, intersect=0;
+       float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
+       float cur_ipoint[3];
+       
+       if(dm==0){
+               psys_disable_all(ob);
+
+               dm=mesh_get_derived_final(scene, ob, 0);
+               if(dm==0)
+                       dm=mesh_get_derived_deform(scene, ob, 0);
+
+               psys_enable_all(ob);
+
+               if(dm==0)
+                       return 0;
+       }
+
+       
+
+       if(pa_minmax==0){
+               INIT_MINMAX(p_min,p_max);
+               DO_MINMAX(co1,p_min,p_max);
+               DO_MINMAX(co2,p_min,p_max);
+       }
+       else{
+               VECCOPY(p_min,pa_minmax);
+               VECCOPY(p_max,pa_minmax+3);
+       }
+
+       totface=dm->getNumFaces(dm);
+       mface=dm->getFaceDataArray(dm,CD_MFACE);
+       mvert=dm->getVertDataArray(dm,CD_MVERT);
+       
+       /* lets intersect the faces */
+       for(i=0; i<totface; i++,mface++){
+               if(vert_cos){
+                       VECCOPY(v1,vert_cos+3*mface->v1);
+                       VECCOPY(v2,vert_cos+3*mface->v2);
+                       VECCOPY(v3,vert_cos+3*mface->v3);
+                       if(mface->v4)
+                               VECCOPY(v4,vert_cos+3*mface->v4)
+               }
+               else{
+                       VECCOPY(v1,mvert[mface->v1].co);
+                       VECCOPY(v2,mvert[mface->v2].co);
+                       VECCOPY(v3,mvert[mface->v3].co);
+                       if(mface->v4)
+                               VECCOPY(v4,mvert[mface->v4].co)
+               }
+
+               if(face_minmax==0){
+                       INIT_MINMAX(min,max);
+                       DO_MINMAX(v1,min,max);
+                       DO_MINMAX(v2,min,max);
+                       DO_MINMAX(v3,min,max);
+                       if(mface->v4)
+                               DO_MINMAX(v4,min,max)
+                       if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
+                               continue;
+               }
+               else{
+                       VECCOPY(min, face_minmax+6*i);
+                       VECCOPY(max, face_minmax+6*i+3);
+                       if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
+                               continue;
+               }
+
+               if(radius>0.0f){
+                       if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
+                               if(cur_d<*min_d){
+                                       *min_d=cur_d;
+                                       VECCOPY(ipoint,cur_ipoint);
+                                       *min_face=i;
+                                       intersect=1;
+                               }
+                       }
+                       if(mface->v4){
+                               if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
+                                       if(cur_d<*min_d){
+                                               *min_d=cur_d;
+                                               VECCOPY(ipoint,cur_ipoint);
+                                               *min_face=i;
+                                               intersect=1;
+                                       }
+                               }
+                       }
+               }
+               else{
+                       if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
+                               if(cur_d<*min_d){
+                                       *min_d=cur_d;
+                                       min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
+                                       min_w[1]= cur_uv[0];
+                                       min_w[2]= cur_uv[1];
+                                       min_w[3]= 0.0f;
+                                       if(mface->v4)
+                                               intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
+                                       *min_face=i;
+                                       intersect=1;
+                               }
+                       }
+                       if(mface->v4){
+                               if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
+                                       if(cur_d<*min_d){
+                                               *min_d=cur_d;
+                                               min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
+                                               min_w[1]= 0.0f;
+                                               min_w[2]= cur_uv[0];
+                                               min_w[3]= cur_uv[1];
+                                               intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
+                                               *min_face=i;
+                                               intersect=1;
+                                       }
+                               }
+                       }
+               }
+       }
+       return intersect;
+}
+
 static int brush_add(PEData *data, short number)
 {
        Scene *scene= data->scene;
@@ -3187,7 +3330,7 @@ static int brush_add(PEData *data, short number)
                min_d=2.0;
                
                /* warning, returns the derived mesh face */
-               if(psys_intersect_dm(scene, ob,dm,0,co1,co2,&min_d,&add_pars[n].num,add_pars[n].fuv,0,0,0,0)) {
+               if(particle_intersect_dm(scene, ob,dm,0,co1,co2,&min_d,&add_pars[n].num,add_pars[n].fuv,0,0,0,0)) {
                        add_pars[n].num_dmcache= psys_particle_dm_face_lookup(ob,psmd->dm,add_pars[n].num,add_pars[n].fuv,NULL);
                        n++;
                }
index 5496c460b6f9a781a4a1db2a264030e2e763ffc1..de84f69c78d404b8ffb3b4c0bba9e755d434f94b 100644 (file)
@@ -474,7 +474,7 @@ typedef struct CollisionModifierData {
        
        unsigned int numverts;
        unsigned int numfaces;
-       float time, pad;                /* cfra time of modifier */
+       float time_x, time_xnew;                /* cfra time of modifier */
        struct BVHTree *bvhtree; /* bounding volume hierarchy for this cloth object */
 } CollisionModifierData;
 
index 2ed66f2374b2091001c1d8be611f98dc54b5ae93..83ba8a12163c63d89b90e251eb9dc6598917e2d3 100644 (file)
@@ -64,7 +64,7 @@ static void initData(ModifierData *md)
        collmd->current_x = NULL;
        collmd->current_xnew = NULL;
        collmd->current_v = NULL;
-       collmd->time = -1000;
+       collmd->time_x = collmd->time_xnew = -1000;
        collmd->numverts = 0;
        collmd->bvhtree = NULL;
 }
@@ -95,7 +95,7 @@ static void freeData(ModifierData *md)
                collmd->current_x = NULL;
                collmd->current_xnew = NULL;
                collmd->current_v = NULL;
-               collmd->time = -1000;
+               collmd->time_x = collmd->time_xnew = -1000;
                collmd->numverts = 0;
                collmd->bvhtree = NULL;
                collmd->mfaces = NULL;
@@ -139,11 +139,11 @@ static void deformVerts(ModifierData *md, Object *ob,
                current_time = BKE_curframe(md->scene);
                
                if(G.rt > 0)
-                       printf("current_time %f, collmd->time %f\n", current_time, collmd->time);
+                       printf("current_time %f, collmd->time_xnew %f\n", current_time, collmd->time_xnew);
                
                numverts = dm->getNumVerts ( dm );
                
-               if((current_time > collmd->time)|| (BKE_ptcache_get_continue_physics()))
+               if((current_time > collmd->time_xnew)|| (BKE_ptcache_get_continue_physics()))
                {
                        unsigned int i;
 
@@ -151,7 +151,7 @@ static void deformVerts(ModifierData *md, Object *ob,
                        if(collmd->x && (numverts != collmd->numverts))
                                freeData((ModifierData *)collmd);
                        
-                       if(collmd->time == -1000) // first time
+                       if(collmd->time_xnew == -1000) // first time
                        {
                                collmd->x = dm->dupVertArray(dm); // frame start position
                                
@@ -174,7 +174,7 @@ static void deformVerts(ModifierData *md, Object *ob,
                                // create bounding box hierarchy
                                collmd->bvhtree = bvhtree_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->x, numverts, ob->pd->pdef_sboft);
                                
-                               collmd->time = current_time;
+                               collmd->time_x = collmd->time_xnew = current_time;
                        }
                        else if(numverts == collmd->numverts)
                        {
@@ -182,6 +182,7 @@ static void deformVerts(ModifierData *md, Object *ob,
                                tempVert = collmd->x;
                                collmd->x = collmd->xnew;
                                collmd->xnew = tempVert;
+                               collmd->time_x = collmd->time_xnew;
                                
                                memcpy(collmd->xnew, dm->getVertArray(dm), numverts*sizeof(MVert));
                                
@@ -216,7 +217,7 @@ static void deformVerts(ModifierData *md, Object *ob,
                                        bvhtree_update_from_mvert ( collmd->bvhtree, collmd->mfaces, collmd->numfaces, collmd->current_x, collmd->current_xnew, collmd->numverts, 1 );
                                }
                                
-                               collmd->time = current_time;
+                               collmd->time_xnew = current_time;
                        }
                        else if(numverts != collmd->numverts)
                        {
@@ -224,7 +225,7 @@ static void deformVerts(ModifierData *md, Object *ob,
                        }
                        
                }
-               else if(current_time < collmd->time)
+               else if(current_time < collmd->time_xnew)
                {       
                        freeData((ModifierData *)collmd);
                }