Fluid physics for particles by Raul Fernandez Hernandez (Farsthary) and Stephen Swhit...
authorJanne Karhu <jhkarh@gmail.com>
Sun, 4 Apr 2010 12:29:06 +0000 (12:29 +0000)
committerJanne Karhu <jhkarh@gmail.com>
Sun, 4 Apr 2010 12:29:06 +0000 (12:29 +0000)
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
source/blender/blenkernel/intern/particle.c
source/blender/blenkernel/intern/particle_system.c
source/blender/blenloader/intern/readfile.c
source/blender/blenloader/intern/writefile.c
source/blender/makesdna/DNA_particle_types.h
source/blender/makesrna/RNA_access.h
source/blender/makesrna/intern/rna_particle.c

index dbd2b90fea2c64d3a6f413a1faf23c334434dad3..54ddb6da8a69098ffbe4c47fc0a561a434c010bc 100644 (file)
@@ -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):
index bd9e041dab4326b0b7f87e7b321dc213119a74fd..7ad65820bbed43f2dd788cd2f867c391f15d9b39 100644 (file)
@@ -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)
index a8446c0009ff73805ccffc04d0509b404e966cbb..aa48336247c56d2c8934562aab28e6e5d06b2845 100644 (file)
@@ -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);
 }
index a51ecbd2a1535817c940a2152e95e3d859fb0d8c..b23b0edb285828db2a3eb2f0c0b96995c3278ada 100644 (file)
@@ -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;
index f5470d1ecec97ba3504c527924fb36d3077e4972..2b32bbdf0c130ffc2a96089455b487638e4c67f6 100644 (file)
@@ -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;
        }
index 36c144a1cfad65379322a397196781d3dc2b72aa..8367875aa7aeebe81248f3839b4057c9befd7125 100644 (file)
@@ -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
index cb792a8d7d85a4946fff7248541cef97761a7af6..e72c5eb9a4da41f1cde8f8ce92109e0c0b426c07 100644 (file)
@@ -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;
index e2d684ecbfb49e2372f92e90599151362f18c8c3..4e8776eff6a6330091c5d9b42d8a6599cd658253 100644 (file)
@@ -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