Fluid physics for particles by Raul Fernandez Hernandez (Farsthary) and Stephen Swhitehorn:

This patch add SPH (Smoothed Particle Hydrodynamics)fluid dynamics to the
blender particle system. SPH is an boundless Lagrangian interpolation
technique to solve the fluid motion equations.
 
From liquids to sand, goo and gases could be simulated using the particle
system.
 
It features internal viscosity, a double density relaxation that accounts
for surface tension effects, static internal springs for plastic fluids,
and buoyancy for gases.

---------------------------------------

This is a commit of the core fluid physics. Raul will work on proper
documentation soon and more features such as surface extraction from
the particle point cloud and increasing stability by sub-frame calculations
later.

release/scripts/ui/properties_particle.py

@@ -395,6 +395,43 @@ class PARTICLE_PT_physics(ParticleButtonsPanel):
             sub.prop(part, "integrator")
             sub.prop(part, "time_tweak")
 
+        elif part.physics_type == 'FLUID':
+            fluid = part.fluid
+            split = layout.split()
+            sub = split.column()
+
+            sub.label(text="Forces:")
+            sub.prop(part, "brownian_factor")
+            sub.prop(part, "drag_factor", slider=True)
+            sub.prop(part, "damp_factor", slider=True)
+            sub = split.column()
+            sub.prop(part, "size_deflect")
+            sub.prop(part, "die_on_collision")
+            sub.prop(part, "integrator")
+            sub.prop(part, "time_tweak")
+
+            split = layout.split()
+            sub = split.column()
+            sub.label(text="Fluid Interaction:")
+            sub.prop(fluid, "fluid_radius", slider=True)
+            sub.prop(fluid, "stiffness_k")
+            sub.prop(fluid, "stiffness_knear")
+            sub.prop(fluid, "rest_density")
+
+            sub.label(text="Viscosity:")
+            sub.prop(fluid, "viscosity_omega", text="Linear")
+            sub.prop(fluid, "viscosity_beta", text="Square")
+
+            sub = split.column()
+
+            sub.label(text="Springs:")
+            sub.prop(fluid, "spring_k", text="Force", slider=True)
+            sub.prop(fluid, "rest_length", slider=True)
+            layout.label(text="Multiple fluids interactions:")
+
+            sub.label(text="Buoyancy:")
+            sub.prop(fluid, "buoyancy", slider=True)
+
         elif part.physics_type == 'KEYED':
             split = layout.split()
             sub = split.column()
@@ -454,7 +491,7 @@ class PARTICLE_PT_physics(ParticleButtonsPanel):
             col.prop(boids, "banking", slider=True)
             col.prop(boids, "height", slider=True)
 
-        if part.physics_type == 'KEYED' or part.physics_type == 'BOIDS':
+        if part.physics_type == 'KEYED' or part.physics_type == 'BOIDS' or part.physics_type == 'FLUID':
             if part.physics_type == 'BOIDS':
                 layout.label(text="Relations:")
 
@@ -484,7 +521,7 @@ class PARTICLE_PT_physics(ParticleButtonsPanel):
                     col.active = psys.keyed_timing
                     col.prop(key, "time")
                     col.prop(key, "duration")
-                else:
+                elif part.physics_type == 'BOIDS':
                     sub = row.row()
                     #doesn't work yet
                     #sub.red_alert = key.valid
@@ -492,6 +529,12 @@ class PARTICLE_PT_physics(ParticleButtonsPanel):
                     sub.prop(key, "system", text="System")
 
                     layout.prop(key, "mode", expand=True)
+                elif part.physics_type == 'FLUID':
+                    sub = row.row()
+                    #doesn't work yet
+                    #sub.red_alert = key.valid
+                    sub.prop(key, "object", text="")
+                    sub.prop(key, "system", text="System")
 
 
 class PARTICLE_PT_boidbrain(ParticleButtonsPanel):

source/blender/blenkernel/intern/particle.c

@@ -355,6 +355,12 @@ int psys_uses_gravity(ParticleSimulationData *sim)
 /************************************************/
 /*			Freeing stuff						*/
 /************************************************/
+void fluid_free_settings(SPHFluidSettings *fluid)
+{
+	if(fluid)
+		MEM_freeN(fluid); 
+}
+
 void psys_free_settings(ParticleSettings *part)
 {
 	BKE_free_animdata(&part->id);
@@ -367,6 +373,7 @@ void psys_free_settings(ParticleSettings *part)
 	BLI_freelistN(&part->dupliweights);
 
 	boid_free_settings(part->boids);
+	fluid_free_settings(part->fluid);
 }
 
 void free_hair(Object *ob, ParticleSystem *psys, int dynamics)

source/blender/blenkernel/intern/particle_system.c

@@ -24,7 +24,7 @@
  *
  * The Original Code is: all of this file.
  *
- * Contributor(s): none yet.
+ * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
  *
  * ***** END GPL LICENSE BLOCK *****
  */
@@ -2270,6 +2270,137 @@ static void psys_update_effectors(ParticleSimulationData *sim)
 	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 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;
+	KDTreeNearest *ptn = NULL;
+	
+	SPHFluidSettings *fluid = part->fluid;
+	ParticleData *second_particle;
+
+	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 omega = fluid->viscosity_omega;
+	float beta = fluid->viscosity_omega;
+	float massfactor = 1.0f/mass;
+	int n, neighbours;
+
+		
+	radius 	= fluid->radius;
+
+	VECCOPY(start, pa->prev_state.co);
+	VECCOPY(end, pa->state.co);
+
+	sub_v3_v3v3(v, end, start);
+	mul_v3_fl(v, 1.f/dtime);
+
+	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);
+		}
+	}
+	/* 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;
+	pressure =  fluid->stiffness_k * (p - fluid->rest_density);
+	pressure_near = fluid->stiffness_knear * pnear;
+
+	for(n=1; n<neighbours; n++) {
+		q = ptn[n].dist/radius;
+
+		D =  dtime * dtime * (pressure*(1-q) + pressure_near*(1-q)*(1-q))* 0.5f;
+		madd_v3_v3fl(end, ptn[n].co, -D * massfactor);
+	} 	
+
+	/* Artificial buoyancy force in negative gravity direction  */
+	if (fluid->buoyancy >= 0.f && psys_uses_gravity(sim)) {
+		float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
+		madd_v3_v3fl(end, sim->scene->physics_settings.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;}
+}
+
+static void apply_particle_fluidsim(ParticleSystem *psys, ParticleData *pa, ParticleSettings *part, ParticleSimulationData *sim, float dfra, float cfra){
+	ParticleTarget *pt;
+	float dtime = dfra*psys_get_timestep(sim);
+	float particle_mass = part->mass;
+
+	particle_fluidsim(psys, pa, part, sim, dfra, cfra, particle_mass);
+	
+	/*----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);
+
+		if(epsys)
+			particle_fluidsim(epsys, pa, epsys->part, sim, dfra, cfra, particle_mass);
+	}
+	/*----------------------------------------------------------------*/	 	 
+}
+
 /************************************************/
 /*			Newtonian physics					*/
 /************************************************/
@@ -2799,7 +2930,7 @@ static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, flo
 				deflections=max_deflections;
 			}
 			else {
-				float nor_vec[3], tan_vec[3], tan_vel[3], vel[3];
+				float nor_vec[3], tan_vec[3], tan_vel[3];
 				float damp, frict;
 				float inp, inp_v;
 				
@@ -3248,6 +3379,14 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 				psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
 		}
 	}
+	else if(part->phystype==PART_PHYS_FLUID){
+		ParticleTarget *pt = psys->targets.first;
+		psys_update_particle_tree(psys, cfra);
+		
+		for(; pt; pt=pt->next) {  /* Updating others systems particle tree for fluid-fluid interaction */
+			if(pt->ob) psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
+		}
+	}
 
 	/* main loop: calculate physics for all particles */
 	LOOP_SHOWN_PARTICLES {
@@ -3318,6 +3457,22 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 					}
 					break;
 				}
+				case PART_PHYS_FLUID:
+				{	
+					/* do global forces & effectors */
+					apply_particle_forces(sim, p, pa_dfra, cfra);
+
+					/* do fluid sim */
+					apply_particle_fluidsim(psys, pa, part, sim, pa_dfra, cfra);
+
+					/* deflection */
+ 					if(sim->colliders)
+						deflect_particle(sim, p, pa_dfra, cfra);
+ 					
+					/* rotations, SPH particles are not physical particles, just interpolation particles,  thus rotation has not a direct sense for them */	
+					rotate_particle(part, pa, pa_dfra, timestep);  
+ 					break;
+				} 
 			}
 
 			if(pa->alive == PARS_DYING){
@@ -3727,6 +3882,21 @@ void psys_check_boid_data(ParticleSystem *psys)
 				pa->boid = NULL;
 		}
 }
+
+static void fluid_default_settings(ParticleSettings *part){
+	SPHFluidSettings *fluid = part->fluid;
+
+	fluid->radius = 0.5f;
+	fluid->spring_k = 0.f;
+	fluid->rest_length = 0.5f;
+	fluid->viscosity_omega = 2.f;
+	fluid->viscosity_beta = 0.f;
+	fluid->stiffness_k = 0.1f;
+	fluid->stiffness_knear = 0.05f;
+	fluid->rest_density = 10.f;
+	fluid->buoyancy = 0.f;
+}
+
 static void psys_changed_physics(ParticleSimulationData *sim)
 {
 	ParticleSettings *part = sim->psys->part;
@@ -3756,6 +3926,10 @@ static void psys_changed_physics(ParticleSimulationData *sim)
 		state->flag |= BOIDSTATE_CURRENT;
 		BLI_addtail(&part->boids->states, state);
 	}
+	else if(part->phystype == PART_PHYS_FLUID && part->fluid == NULL) {
+		part->fluid = MEM_callocN(sizeof(SPHFluidSettings), "SPH Fluid Settings");
+		fluid_default_settings(part);
+	}
 
 	psys_check_boid_data(sim->psys);
 }

source/blender/blenloader/intern/readfile.c

@@ -3038,6 +3038,7 @@ static void direct_link_particlesettings(FileData *fd, ParticleSettings *part)
 	link_list(fd, &part->dupliweights);
 
 	part->boids= newdataadr(fd, part->boids);
+	part->fluid= newdataadr(fd, part->fluid);
 
 	if(part->boids) {
 		BoidState *state;

source/blender/blenloader/intern/writefile.c

@@ -652,6 +652,9 @@ static void write_particlesettings(WriteData *wd, ListBase *idbase)
 				for(; state; state=state->next)
 					write_boid_state(wd, state);
 			}
+			if(part->fluid && part->phystype == PART_PHYS_FLUID){
+				writestruct(wd, DATA, "SPHFluidSettings", 1, part->fluid); 
+			}
 		}
 		part= part->id.next;
 	}

source/blender/makesdna/DNA_particle_types.h

@@ -114,11 +114,20 @@ typedef struct ParticleData {
 	short alive;			/* the life state of a particle */
 } ParticleData;
 
+typedef struct SPHFluidSettings {
+	/*Particle Fluid*/
+	float spring_k, radius, rest_length;
+	float viscosity_omega, viscosity_beta;
+	float stiffness_k, stiffness_knear, rest_density;
+	float buoyancy;
+} SPHFluidSettings;
+
 typedef struct ParticleSettings {
 	ID id;
 	struct AnimData *adt;
 
 	struct BoidSettings *boids;
+	struct SPHFluidSettings *fluid;
 
 	struct EffectorWeights *effector_weights;
 
@@ -322,6 +331,7 @@ typedef struct ParticleSystem{				/* note, make sure all (runtime) are NULL's in
 #define PART_PHYS_NEWTON	1
 #define PART_PHYS_KEYED		2
 #define PART_PHYS_BOIDS		3
+#define PART_PHYS_FLUID		4
 
 /* part->kink */
 #define PART_KINK_NO		0

source/blender/makesrna/RNA_access.h

@@ -351,6 +351,7 @@ extern StructRNA RNA_ParticleHairKey;
 extern StructRNA RNA_ParticleInstanceModifier;
 extern StructRNA RNA_ParticleKey;
 extern StructRNA RNA_ParticleSettings;
+extern StructRNA RNA_SPHFluidSettings;
 extern StructRNA RNA_ParticleSystem;
 extern StructRNA RNA_ParticleSystemModifier;
 extern StructRNA RNA_ParticleTarget;

source/blender/makesrna/intern/rna_particle.c

@@ -837,6 +837,74 @@ static void rna_def_particle_dupliweight(BlenderRNA *brna)
 	RNA_def_property_update(prop, 0, "rna_Particle_redo");
 }
 
+static void rna_def_fluid_settings(BlenderRNA *brna)
+{
+	StructRNA *srna;
+	PropertyRNA *prop;
+	
+	srna = RNA_def_struct(brna, "SPHFluidSettings", NULL);
+	RNA_def_struct_ui_text(srna, "SPH Fluid Settings", "Settings for particle fluids physics");
+	
+	/* Fluid settings */
+ 	prop= RNA_def_property(srna, "spring_k", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "spring_k");
+ 	RNA_def_property_range(prop, 0.0f, 1.0f);
+ 	RNA_def_property_ui_text(prop, "Spring", "Spring force constant");
+        RNA_def_property_update(prop, 0, "rna_Particle_reset");
+  
+ 	prop= RNA_def_property(srna, "fluid_radius", PROP_FLOAT, PROP_NONE);
+ 	RNA_def_property_float_sdna(prop, NULL, "radius");
+ 	RNA_def_property_range(prop, 0.0f, 2.0f);
+	RNA_def_property_ui_text(prop, "Radius", "Fluid interaction Radius");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+ 	prop= RNA_def_property(srna, "rest_length", PROP_FLOAT, PROP_NONE);
+ 	RNA_def_property_float_sdna(prop, NULL, "rest_length");
+ 	RNA_def_property_range(prop, 0.0f, 1.0f);
+ 	RNA_def_property_ui_text(prop, "Rest Length", "The Spring Rest Length (factor of interaction radius)");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset"); 
+ 
+ 	/* Viscosity */
+ 	prop= RNA_def_property(srna, "viscosity_omega", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "viscosity_omega");
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Viscosity", "Linear viscosity");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+  
+ 	prop= RNA_def_property(srna, "viscosity_beta", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "viscosity_beta");
+	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Square viscosity", "Square viscosity factor");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	/* Double density relaxation */
+ 	prop= RNA_def_property(srna, "stiffness_k", PROP_FLOAT, PROP_NONE);
+ 	RNA_def_property_float_sdna(prop, NULL, "stiffness_k");
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Stiffness ", "Constant K - Stiffness");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+  
+ 	prop= RNA_def_property(srna, "stiffness_knear", PROP_FLOAT, PROP_NONE);
+ 	RNA_def_property_float_sdna(prop, NULL, "stiffness_knear");
+	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Repulsion", "Repulsion factor: stiffness_knear");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+  
+ 	prop= RNA_def_property(srna, "rest_density", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "rest_density");
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Rest Density", "Density");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+	
+	/* Buoyancy */
+	prop= RNA_def_property(srna, "buoyancy", PROP_FLOAT, PROP_NONE);
+ 	RNA_def_property_float_sdna(prop, NULL, "buoyancy");
+ 	RNA_def_property_range(prop, 0.0f, 1.0f);
+ 	RNA_def_property_ui_text(prop, "Buoyancy", "");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+	
+}
+
 static void rna_def_particle_settings(BlenderRNA *brna)
 {
 	StructRNA *srna;
@@ -861,6 +929,7 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 		{PART_PHYS_NEWTON, "NEWTON", 0, "Newtonian", ""},
 		{PART_PHYS_KEYED, "KEYED", 0, "Keyed", ""},
 		{PART_PHYS_BOIDS, "BOIDS", 0, "Boids", ""},
+		{PART_PHYS_FLUID, "FLUID", 0, "Fluid", ""},
 		{0, NULL, 0, NULL, NULL}
 	};
 
@@ -1773,9 +1842,14 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 	RNA_def_property_struct_type(prop, "BoidSettings");
 	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
 	RNA_def_property_ui_text(prop, "Boid Settings", "");
-
+	
+	/* Fluid particles */
+	prop= RNA_def_property(srna, "fluid", PROP_POINTER, PROP_NONE);
+	RNA_def_property_struct_type(prop, "SPHFluidSettings");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "SPH Fluid Settings", ""); 
+	
 	/* draw objects & groups */
-
 	prop= RNA_def_property(srna, "dupli_group", PROP_POINTER, PROP_NONE);
 	RNA_def_property_pointer_sdna(prop, NULL, "dup_group");
 	RNA_def_property_struct_type(prop, "Group");
@@ -2162,12 +2236,13 @@ void RNA_def_particle(BlenderRNA *brna)
 
 	rna_def_particle_hair_key(brna);
 	rna_def_particle_key(brna);
-
+	rna_def_fluid_settings(brna);
+	
 	rna_def_child_particle(brna);
 	rna_def_particle(brna);
 	rna_def_particle_dupliweight(brna);
 	rna_def_particle_system(brna);
-	rna_def_particle_settings(brna);
+	rna_def_particle_settings(brna);	
 }
 
 #endif