[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [32073] trunk/blender/source/blender/ blenkernel/intern/particle_system.c: Fix for [#23872] particle deflection in conjunction with SPH particles is apparently buggy

Janne Karhu jhkarh at gmail.com
Thu Sep 23 11:31:13 CEST 2010


Revision: 32073
          http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=32073
Author:   jhk
Date:     2010-09-23 11:31:13 +0200 (Thu, 23 Sep 2010)

Log Message:
-----------
Fix for [#23872] particle deflection in conjunction with SPH particles is apparently buggy
* Fix turned into a thorough cleanup and reorganization of particle collision response code.
* Collisions are now much more accurate, stable and even a bit more in agreement with real world physics.
* Only still remaining problem is rotating/deforming deflector objects, but that's something for the future.
* Visible changes should only be positive, i.e. no leaking particles, no strange instabilities etc.

Modified Paths:
--------------
    trunk/blender/source/blender/blenkernel/intern/particle_system.c

Modified: trunk/blender/source/blender/blenkernel/intern/particle_system.c
===================================================================
--- trunk/blender/source/blender/blenkernel/intern/particle_system.c	2010-09-23 08:15:53 UTC (rev 32072)
+++ trunk/blender/source/blender/blenkernel/intern/particle_system.c	2010-09-23 09:31:13 UTC (rev 32073)
@@ -2837,11 +2837,18 @@
 
 	} while(t2);
 }
-/* particle - mesh collision code */
-/* in addition to basic point to surface collisions handles friction & damping,*/
-/* angular momentum <-> linear momentum and swept sphere - mesh collisions */
-/* 1. check for all possible deflectors for closest intersection on particle path */
-/* 2. if deflection was found kill the particle or calculate new coordinates */
+/* 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;
@@ -2849,19 +2856,21 @@
 	ParticleCollision col;
 	ColliderCache *coll;
 	BVHTreeRayHit hit;
-	float ray_dir[3], zerovec[3]={0.0,0.0,0.0};
+	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);
+	float timestep = psys_get_timestep(sim) * dfra;
+	float inv_timestep = 1.0f/timestep;
 	int deflections=0, max_deflections=10;
 
-	VECCOPY(col.co1, pa->prev_state.co);
-	VECCOPY(col.co2, pa->state.co);
-	
-	VECCOPY(col.ve1, pa->prev_state.vel);
-	VECCOPY(col.ve2, pa->state.vel);
-	mul_v3_fl(col.ve1, timestep * dfra);
-	mul_v3_fl(col.ve2, timestep * dfra);
-	
+	/* 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);
+
+	/* 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.t = 0.0f;
 
 	/* override for boids */
@@ -2876,7 +2885,7 @@
 	if(sim->colliders) while(deflections < max_deflections){
 		/* 1. */
 
-		VECSUB(ray_dir, col.co2, col.co1);
+		sub_v3_v3v3(ray_dir, col.co2, col.co1);
 		hit.index = -1;
 		hit.dist = col.ray_len = len_v3(ray_dir);
 
@@ -2904,46 +2913,50 @@
 		/* 2. */
 		if(hit.index>=0) {
 			PartDeflect *pd = col.hit_ob->pd;
-			int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0;
 			float co[3]; /* point of collision */
-			float vec[3]; /* movement through collision */
-			float acc[3]; /* acceleration */
+			float x = hit.dist/col.ray_len; /* location of collision between this iteration */
+			float df = col.t + x * (1.0f - col.t); /* time of collision between frame change*/
+			float dt1 = (df - col.t) * timestep; /* iteration time of collision (in seconds) */
+			float dt2 = (1.0f - df) * timestep; /* time left after collision (in seconds) */
+			int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
 
-			float x = hit.dist/col.ray_len; /* location of collision between this iteration */
-			float le = len_v3(col.ve1)/col.ray_len;
-			float ac = len_v3(col.ve2)/col.ray_len - le; /* (taking acceleration into account) */
-			float t = (-le + sqrt(le*le + 2*ac*x))/ac; /* time of collision between this iteration */
-			float dt = col.t + x * (1.0f - col.t); /* time of collision between frame change*/
-			float it = 1.0 - t;
+			deflections++;
 			
 			interp_v3_v3v3(co, col.co1, col.co2, x);
-			VECSUB(vec, col.co2, col.co1);
-
-			VECSUB(acc, col.ve2, col.ve1);
 			
-			mul_v3_fl(col.vel, 1.0f-col.t);
+			/* 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));
 
 			/* particle dies in collision */
 			if(through == 0 && (part->flag & PART_DIE_ON_COL || pd->flag & PDEFLE_KILL_PART)) {
 				pa->alive = PARS_DYING;
-				pa->dietime = pa->state.time + (cfra - pa->state.time) * dt;
-				
-				/* we have to add this for dying particles too so that reactors work correctly */
-				VECADDFAC(co, co, col.nor, (through ? -0.0001f : 0.0001f));
+				pa->dietime = pa->state.time + (cfra - pa->state.time) * df;
 
-				VECCOPY(pa->state.co, co);
-				interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, dt);
-				interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, dt);
-				interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, dt);
+				copy_v3_v3(pa->state.co, co);
+				interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, df);
+				interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, df);
+				interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, df);
 
 				/* particle is dead so we don't need to calculate further */
-				deflections=max_deflections;
+				return;
 			}
+			/* figure out velocity and other data after collision */
 			else {
-				float nor_vec[3], tan_vec[3], tan_vel[3];
+				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;
-				float inp, inp_v;
+
+				/* 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);
@@ -2952,119 +2965,118 @@
 				CLAMP(frict,0.0,1.0);
 
 				/* treat normal & tangent components separately */
-				inp = dot_v3v3(col.nor, vec);
-				inp_v = dot_v3v3(col.nor, col.vel);
+				v0_dot = dot_v3v3(col.nor, v0);
+				madd_v3_v3v3fl(v0_tan, v0, col.nor, -v0_dot);
 
-				VECADDFAC(tan_vec, vec, col.nor, -inp);
-				VECADDFAC(tan_vel, col.vel, col.nor, -inp_v);
-				if((part->flag & PART_ROT_DYN)==0)
-					interp_v3_v3v3(tan_vec, tan_vec, tan_vel, frict);
+				vc_dot = dot_v3v3(col.nor, col.vel);
+				madd_v3_v3v3fl(vc_tan, col.vel, col.nor, -vc_dot);
 
-				VECCOPY(nor_vec, col.nor);
-				inp *= 1.0f - damp;
+				/* 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);
 
-				if(through)
-					inp_v *= damp;
+						/* 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) */
 
-				/* special case for object hitting the particle from behind */
-				if(through==0 && ((inp_v>0 && inp>0 && inp_v>inp) || (inp_v<0 && inp<0 && inp_v<inp)))
-					mul_v3_fl(nor_vec, inp_v);
-				else
-					mul_v3_fl(nor_vec, inp_v + (through ? 1.0f : -1.0f) * inp);
+						/* rolling friction is around 0.01 of sliding friction (could be made a parameter) */
+						mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
 
-				/* angular <-> linear velocity - slightly more physical and looks even nicer than before */
-				if(part->flag & PART_ROT_DYN) {
-					float surface_vel[3], rot_vel[3], friction[3], dave[3], dvel[3];
+						/* surface_velocity is opposite to cm velocity */
+						mul_v3_v3fl(vr_tan, v1_tan, -1.0f);
 
-					/* apparent velocity along collision surface */
-					VECSUB(surface_vel, tan_vec, tan_vel);
+						/* get back to global coordinates */
+						add_v3_v3(v1_tan, vc_tan);
 
-					/* direction of rolling friction */
-					cross_v3_v3v3(rot_vel, pa->state.ave, col.nor);
-					/* convert to current dt */
-					mul_v3_fl(rot_vel, (timestep*dfra) * (1.0f - col.t));
-					mul_v3_fl(rot_vel, pa->size);
+						/* convert to angular velocity*/
+						cross_v3_v3v3(ave, vr_tan, col.nor);
+						mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001));
 
-					/* apply sliding friction */
-					VECSUB(surface_vel, surface_vel, rot_vel);
-					VECCOPY(friction, surface_vel);
+						/* 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);
+					}
+				}
 
-					mul_v3_fl(surface_vel, 1.0 - frict);
-					mul_v3_fl(friction, frict);
+				/* 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;
+				}
 
-					/* sliding changes angular velocity */
-					cross_v3_v3v3(dave, col.nor, friction);
-					mul_v3_fl(dave, 1.0f/MAX2(pa->size, 0.001));
+				/* damping and flipping of velocity around normal */
+				v0_dot *= 1.0f - damp;
+				vc_dot *= through ? damp : 1.0f;
 
-					/* we assume rolling friction is around 0.01 of sliding friction */
-					mul_v3_fl(rot_vel, 1.0 - frict*0.01);
+				/* special case for object hitting the particle from behind */

@@ Diff output truncated at 10240 characters. @@




More information about the Bf-blender-cvs mailing list