diff options
author | Janne Karhu <jhkarh@gmail.com> | 2011-03-18 18:31:32 +0300 |
---|---|---|
committer | Janne Karhu <jhkarh@gmail.com> | 2011-03-18 18:31:32 +0300 |
commit | 60ce95f5622d3947e30fe960eda22ea305660619 (patch) | |
tree | 20a322ff834993f9794ebb5c110437da41273812 /source | |
parent | 7e53769d09863b59beae4155ea0ad13c6ab2fac2 (diff) |
New particle collisions code:
* 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!
Diffstat (limited to 'source')
-rw-r--r-- | source/blender/blenkernel/BKE_particle.h | 60 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/boids.c | 45 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/particle_system.c | 987 | ||||
-rw-r--r-- | source/blender/blenloader/intern/readfile.c | 2 | ||||
-rw-r--r-- | source/blender/editors/physics/particle_edit.c | 145 | ||||
-rw-r--r-- | source/blender/makesdna/DNA_modifier_types.h | 2 | ||||
-rw-r--r-- | source/blender/modifiers/intern/MOD_collision.c | 17 |
7 files changed, 829 insertions, 429 deletions
diff --git a/source/blender/blenkernel/BKE_particle.h b/source/blender/blenkernel/BKE_particle.h index e4fbd5603f2..9b1651b12cd 100644 --- a/source/blender/blenkernel/BKE_particle.h +++ b/source/blender/blenkernel/BKE_particle.h @@ -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 */ diff --git a/source/blender/blenkernel/intern/boids.c b/source/blender/blenkernel/intern/boids.c index 9fb7f6408ac..286cd3b2499 100644 --- a/source/blender/blenkernel/intern/boids.c +++ b/source/blender/blenkernel/intern/boids.c @@ -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 */ diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c index 5f7c62042ec..b5eb9b6ab29 100644 --- a/source/blender/blenkernel/intern/particle_system.c +++ b/source/blender/blenkernel/intern/particle_system.c @@ -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); diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c index 4c04db5f56f..35aeb756274 100644 --- a/source/blender/blenloader/intern/readfile.c +++ b/source/blender/blenloader/intern/readfile.c @@ -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; diff --git a/source/blender/editors/physics/particle_edit.c b/source/blender/editors/physics/particle_edit.c index 8a304765a7a..d3e61e785e1 100644 --- a/source/blender/editors/physics/particle_edit.c +++ b/source/blender/editors/physics/particle_edit.c @@ -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++; } diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h index 5496c460b6f..de84f69c78d 100644 --- a/source/blender/makesdna/DNA_modifier_types.h +++ b/source/blender/makesdna/DNA_modifier_types.h @@ -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; diff --git a/source/blender/modifiers/intern/MOD_collision.c b/source/blender/modifiers/intern/MOD_collision.c index 2ed66f2374b..83ba8a12163 100644 --- a/source/blender/modifiers/intern/MOD_collision.c +++ b/source/blender/modifiers/intern/MOD_collision.c @@ -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); } |