[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [33164] trunk/blender/source/blender/ blenkernel: Algorithm fix for fluid particles:

Janne Karhu jhkarh at gmail.com
Thu Nov 18 20:12:36 CET 2010


Revision: 33164
          http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=33164
Author:   jhk
Date:     2010-11-18 20:12:36 +0100 (Thu, 18 Nov 2010)

Log Message:
-----------
Algorithm fix for fluid particles:
* The SPH fluid particle algorithm was implemented a bit wrong. This problem could for example result in the fluid moving sideways after being dropped straight to a horizontal collision surface, a very big no-no as far as real world physics are concerned!
* After some extensive code shuffling the algorithm is now much more true to the paper it was implemented from, and more importantly now the physics should be correct too!
* The main thing was that fluids calculations can effect many particles simultaneously, so just a single loop through all particles can't work properly. As a side note this also means that the actual fluid algorithm can't be made threaded :(
* To make things work I also had to reshuffle some general particle physics code, but there should be no functional changes what so ever to other physics types, so poke me immediately if something strange happens.

Note to users: these changes will most probably effect the way previously done sph fluid simulations look, so some parameter tweaking will be needed to get things back looking the way they were.

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

Modified: trunk/blender/source/blender/blenkernel/BKE_particle.h
===================================================================
--- trunk/blender/source/blender/blenkernel/BKE_particle.h	2010-11-18 19:11:05 UTC (rev 33163)
+++ trunk/blender/source/blender/blenkernel/BKE_particle.h	2010-11-18 19:12:36 UTC (rev 33164)
@@ -63,6 +63,7 @@
 #define LOOP_PARTICLES	for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
 #define LOOP_EXISTING_PARTICLES for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++) if(!(pa->flag & PARS_UNEXIST))
 #define LOOP_SHOWN_PARTICLES for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++) if(!(pa->flag & (PARS_UNEXIST|PARS_NO_DISP)))
+#define LOOP_DYNAMIC_PARTICLES for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++) if(pa->state.time > 0.f)
 
 #define PSYS_FRAND_COUNT	1024
 #define PSYS_FRAND(seed)	psys->frand[(seed) % PSYS_FRAND_COUNT]

Modified: trunk/blender/source/blender/blenkernel/intern/particle_system.c
===================================================================
--- trunk/blender/source/blender/blenkernel/intern/particle_system.c	2010-11-18 19:11:05 UTC (rev 33163)
+++ trunk/blender/source/blender/blenkernel/intern/particle_system.c	2010-11-18 19:12:36 UTC (rev 33164)
@@ -2291,132 +2291,121 @@
 	precalc_guides(sim, sim->psys->effectors);
 }
 
-/*************************************************
+/*********************************************************************************************************
                     SPH fluid physics 
 
- In theory, there could be unlimited implementation
-                    of SPH simulators
-**************************************************/
-void particle_fluidsim(ParticleSystem *psys, ParticleData *pa, ParticleSettings *part, ParticleSimulationData *sim, float dfra, float UNUSED(cfra), float mass){
-/****************************************************************************************************************
-* 	This code uses in some parts adapted algorithms from the pseduo code as outlined in the Research paper
-*	Titled: Particle-based Viscoelastic Fluid Simulation.
-* 	Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
-*
-*	Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
-*	Presented at Siggraph, (2005)
-*
-*****************************************************************************************************************/
-	KDTree *tree = psys->tree;
+ In theory, there could be unlimited implementation of SPH simulators
+
+ This code uses in some parts adapted algorithms from the pseduo code as outlined in the Research paper:
+
+ Titled: Particle-based Viscoelastic Fluid Simulation.
+ Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
+ Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
+
+ Presented at Siggraph, (2005)
+
+***********************************************************************************************************/
+static void particle_fluidsim(ParticleSystem *psys, int own_psys, ParticleData *pa, ParticleSettings *part, float dtime, float mass, float *gravity)
+{
+	SPHFluidSettings *fluid = psys->part->fluid;
 	KDTreeNearest *ptn = NULL;
-	
-	SPHFluidSettings *fluid = part->fluid;
-	ParticleData *second_particle;
+	ParticleData *npa;
 
-	float start[3], end[3], v[3];
 	float temp[3];
-	float q, radius, D;
-	float p, pnear, pressure_near, pressure;
-	float dtime = dfra * psys_get_timestep(sim);
+	float q, q1, u, I, D;
+	float pressure_near, pressure;
+	float p=0, pnear=0;
+
+	float radius = fluid->radius;
 	float omega = fluid->viscosity_omega;
-	float beta = fluid->viscosity_omega;
+	float beta = fluid->viscosity_beta;
 	float massfactor = 1.0f/mass;
-	int n, neighbours;
+	float spring_k = fluid->spring_k;
+	float L = fluid->rest_length;
 
-		
-	radius 	= fluid->radius;
+	int n, neighbours = BLI_kdtree_range_search(psys->tree, radius, pa->prev_state.co, NULL, &ptn);
+	int index = own_psys ? pa - psys->particles : -1;
 
-	VECCOPY(start, pa->prev_state.co);
-	VECCOPY(end, pa->state.co);
+	/* pressure and near pressure */
+	for(n=own_psys?1:0; n<neighbours; n++) {
+		sub_v3_v3(ptn[n].co, pa->prev_state.co);
+		mul_v3_fl(ptn[n].co, 1.f/ptn[n].dist);
+		q = ptn[n].dist/radius;
 
-	VECCOPY(v, pa->state.vel);
+		if(q < 1.f) {
+			q1 = 1.f - q;
 
-	neighbours = BLI_kdtree_range_search(tree, radius, start, NULL, &ptn);
-
-	/* use ptn[n].co to store relative direction */
-	for(n=1; n<neighbours; n++) {
-		sub_v3_v3(ptn[n].co, start);
-		normalize_v3(ptn[n].co);
-	}
-        
-	/* Viscosity - Algorithm 5  */
-	if (omega > 0.f || beta > 0.f) {
-		float u, I;
-
-		for(n=1; n<neighbours; n++) {
-			second_particle = psys->particles + ptn[n].index;
-			q = ptn[n].dist/radius;
-			
-			sub_v3_v3v3(temp, v, second_particle->prev_state.vel);
-			
-			u = dot_v3v3(ptn[n].co, temp);
-
-			if (u > 0){
-				I = dtime * ((1-q) * (omega * u + beta * u*u)) * 0.5f;
-				madd_v3_v3fl(v, ptn[n].co, -I * massfactor);
-			} 
-		}	
-	}
-
-	/* Hooke's spring force  */
-	if (fluid->spring_k > 0.f) {
-		float D, L = fluid->rest_length;
-		for(n=1; n<neighbours; n++) {
-			/* L is a factor of radius */
-			D = dtime * 10.f * fluid->spring_k * (1.f - L) * (L - ptn[n].dist/radius);
-			madd_v3_v3fl(v, ptn[n].co, -D * massfactor);
+			p += q1*q1;
+			pnear += q1*q1*q1;
 		}
 	}
-	/* Update particle position */	
-	VECADDFAC(end, start, v, dtime);
 
-	/* Double Density Relaxation - Algorithm 2 */
-	p = 0;
-	pnear = 0;
-	for(n=1; n<neighbours; n++) {
-		q = ptn[n].dist/radius;
-		p += ((1-q)*(1-q));
-		pnear += ((1-q)*(1-q)*(1-q));
-	}
-	p *= part->mass;
-	pnear *= part->mass;
+	p *= mass;
+	pnear *= mass;
 	pressure =  fluid->stiffness_k * (p - fluid->rest_density);
 	pressure_near = fluid->stiffness_knear * pnear;
 
-	for(n=1; n<neighbours; n++) {
+	/* main calculations */
+	for(n=own_psys?1:0; n<neighbours; n++) {
+		npa = psys->particles + ptn[n].index;
+
 		q = ptn[n].dist/radius;
+		q1 = 1.f-q;
 
-		D =  dtime * dtime * (pressure*(1-q) + pressure_near*(1-q)*(1-q))* 0.5f;
-		madd_v3_v3fl(end, ptn[n].co, -D * massfactor);
-	} 	
+		/* Double Density Relaxation - Algorithm 2 (can't be thread safe!)*/
+		D =  dtime * dtime * (pressure + pressure_near*q1)*q1 * 0.5f;
+		madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
+		if(own_psys)
+			madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
 
+		if(index < ptn[n].index) {
+			/* Viscosity - Algorithm 5 */
+			if(omega > 0.f	|| beta > 0.f) {		
+				sub_v3_v3v3(temp, pa->state.vel, npa->state.vel);
+				u = dot_v3v3(ptn[n].co, temp);
+
+				if (u > 0){
+					I = dtime * (q1 * (omega * u + beta * u*u)) * 0.5f;
+					madd_v3_v3fl(pa->state.vel, ptn[n].co, -I * massfactor);
+
+					if(own_psys)
+						madd_v3_v3fl(npa->state.vel, ptn[n].co, I * massfactor);
+				}
+			}
+
+			/* Hooke's spring force */
+			if(spring_k > 0.f) {
+				/* L is a factor of radius */
+				D = 0.5 * dtime * dtime * 10.f * fluid->spring_k * (1.f - L) * (L - q);
+
+				madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
+				if(own_psys)
+					madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
+			}
+		}
+	} 
+
 	/* Artificial buoyancy force in negative gravity direction  */
-	if (fluid->buoyancy >= 0.f && psys_uses_gravity(sim)) {
+	if (fluid->buoyancy >= 0.f && gravity) {
 		float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
-		madd_v3_v3fl(end, sim->scene->physics_settings.gravity, -B * massfactor);
+		madd_v3_v3fl(pa->state.co, gravity, -B * massfactor);
 	}
 
-	/* apply final result and recalculate velocity */
-	VECCOPY(pa->state.co, end);
-	sub_v3_v3v3(pa->state.vel, end, start);
-	mul_v3_fl(pa->state.vel, 1.f/dtime);
-
-	if(ptn){ MEM_freeN(ptn); ptn=NULL;}
+	if(ptn)
+		MEM_freeN(ptn);
 }
 
-static void apply_particle_fluidsim(ParticleSystem *psys, ParticleData *pa, ParticleSettings *part, ParticleSimulationData *sim, float dfra, float cfra){
+static void apply_particle_fluidsim(Object *ob, ParticleSystem *psys, ParticleData *pa, float dtime, float *gravity){
 	ParticleTarget *pt;
-//	float dtime = dfra*psys_get_timestep(sim);
-	float particle_mass = part->mass;
 
-	particle_fluidsim(psys, pa, part, sim, dfra, cfra, particle_mass);
+	particle_fluidsim(psys, 1, pa, psys->part, dtime, psys->part->mass, gravity);
 	
 	/*----check other SPH systems (Multifluids) , each fluid has its own parameters---*/
-	for(pt=sim->psys->targets.first; pt; pt=pt->next) {
-		ParticleSystem *epsys = psys_get_target_system(sim->ob, pt);
+	for(pt=psys->targets.first; pt; pt=pt->next) {
+		ParticleSystem *epsys = psys_get_target_system(ob, pt);
 
 		if(epsys)
-			particle_fluidsim(epsys, pa, epsys->part, sim, dfra, cfra, particle_mass);
+			particle_fluidsim(epsys, 0, pa, epsys->part, dtime, psys->part->mass, gravity);
 	}
 	/*----------------------------------------------------------------*/	 	 
 }
@@ -3372,17 +3361,15 @@
 	/* current time */
 	float ctime;
 	/* frame & time changes */
-	float dfra, dtime, pa_dtime, pa_dfra=0.0;
+	float dfra, dtime;
 	float birthtime, dietime;
-
-	int invalidParticles=0;
 	
 	/* where have we gone in time since last time */
 	dfra= cfra - psys->cfra;
 
 	timestep = psys_get_timestep(sim);
+	ctime= cfra*timestep;
 	dtime= dfra*timestep;
-	ctime= cfra*timestep;
 
 	if(dfra<0.0){
 		LOOP_EXISTING_PARTICLES {
@@ -3402,33 +3389,40 @@
 	if(part->type != PART_HAIR)
 		sim->colliders = get_collider_cache(sim->scene, NULL, NULL);
 
-	if(part->phystype==PART_PHYS_BOIDS){
-		ParticleTarget *pt = psys->targets.first;
-		bbd.sim = sim;
-		bbd.part = part;
-		bbd.cfra = cfra;
-		bbd.dfra = dfra;
-		bbd.timestep = timestep;
+	/* initialize physics type specific stuff */
+	switch(part->phystype) {
+		case PART_PHYS_BOIDS:
+		{
+			ParticleTarget *pt = psys->targets.first;
+			bbd.sim = sim;
+			bbd.part = part;
+			bbd.cfra = cfra;
+			bbd.dfra = dfra;
+			bbd.timestep = timestep;
 
-		psys_update_particle_tree(psys, cfra);
+			psys_update_particle_tree(psys, cfra);
 
-		boids_precalc_rules(part, cfra);
+			boids_precalc_rules(part, cfra);
 
-		for(; pt; pt=pt->next) {
-			if(pt->ob)
-				psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
+			for(; pt; pt=pt->next) {
+				if(pt->ob)
+					psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
+			}
+			break;
 		}
-	}
-	else if(part->phystype==PART_PHYS_FLUID){
-		ParticleTarget *pt = psys->targets.first;
-		psys_update_particle_tree(psys, cfra);
-		

@@ Diff output truncated at 10240 characters. @@




More information about the Bf-blender-cvs mailing list