Implementation of a target density feature for the hair simulation.

This allows setting a target density which the fluid simulation will
take into account as an additional term in the pressure Poisson
equation. Based on two papers
"Detail Preserving Continuum Simulation of Straight Hair" (McAdams et al. 2009)
and
"Two-way Coupled SPH and Particle Level Set Fluid Simulation" (Losasso et al. 2008)

Currently the target pressure is specified directly, but it will be
a lot more convenient to define this in terms of a geometric value such
as "number of hairs per area" (combined with hair "thickness").

Conflicts:
	source/blender/physics/intern/BPH_mass_spring.cpp

release/scripts/startup/bl_ui/properties_particle.py

@@ -325,10 +325,9 @@ class PARTICLE_PT_hair_dynamics(ParticleButtonsPanel, Panel):
         col.label(text="Volume")
         col.prop(cloth, "air_damping", text="Air Drag")
         col.prop(cloth, "internal_friction", slider=True)
-        sub = col.row(align=True)
-        # XXX disabled due to stability issues
-        #sub.prop(cloth, "pressure", slider=True, text="Pressure")
-        #sub.prop(cloth, "pressure_threshold", slider=True, text="Threshold")
+        sub = col.column(align=True)
+        sub.prop(cloth, "density_target", text="Density Target")
+        sub.prop(cloth, "density_strength", slider=True, text="Strength")
         col.prop(cloth, "voxel_cell_size")
 
         split.separator()

source/blender/makesdna/DNA_cloth_types.h

@@ -69,8 +69,8 @@ typedef struct ClothSimSettings {
 	float	goalspring;
 	float	goalfrict;
 	float	velocity_smooth; /* smoothing of velocities for hair */
-	float	pressure;			/* pressure factor from hair density */
-	float	pressure_threshold;	/* minimum density for hair pressure */
+	float	density_target;		/* minimum density for hair */
+	float	density_strength;	/* influence of hair density */
 	float	collider_friction; /* friction with colliders */
 	float	vel_damping; /* damp the velocity to speed up getting to the resting position */
 	float	shrink_min;  /* min amount to shrink cloth by 0.0f (no shrink) - 1.0f (shrink to nothing) */

source/blender/makesrna/intern/rna_cloth.c

@@ -385,16 +385,16 @@ static void rna_def_cloth_sim_settings(BlenderRNA *brna)
 	RNA_def_property_ui_text(prop, "Collider Friction", "");
 	RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-	prop = RNA_def_property(srna, "pressure", PROP_FLOAT, PROP_NONE);
-	RNA_def_property_float_sdna(prop, NULL, "pressure");
-	RNA_def_property_range(prop, 0.0f, 1000.0f);
-	RNA_def_property_ui_text(prop, "Internal Pressure", "Generate outward force based on hair density");
+	prop = RNA_def_property(srna, "density_target", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "density_target");
+	RNA_def_property_range(prop, 0.0f, 1000000.0f);
+	RNA_def_property_ui_text(prop, "Target Density", "Maximum density of hair");
 	RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-	prop = RNA_def_property(srna, "pressure_threshold", PROP_FLOAT, PROP_NONE);
-	RNA_def_property_float_sdna(prop, NULL, "pressure_threshold");
+	prop = RNA_def_property(srna, "density_strength", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "density_strength");
 	RNA_def_property_range(prop, 0.0f, 1000.0f);
-	RNA_def_property_ui_text(prop, "Internal Pressure Threshold", "No pressure force is generated below this pressure value");
+	RNA_def_property_ui_text(prop, "Target Density Strength", "Influence of target density on the simulation");
 	RNA_def_property_update(prop, 0, "rna_cloth_update");
 
 	/* mass */

source/blender/physics/intern/BPH_mass_spring.cpp

@@ -1017,8 +1017,8 @@ static void cloth_continuum_step(ClothModifierData *clmd, float dt)
 	
 	const float fluid_factor = 0.95f; /* blend between PIC and FLIP methods */
 	float smoothfac = parms->velocity_smooth;
-	float pressfac = parms->pressure;
-	float minpress = parms->pressure_threshold;
+	float denstarget = parms->density_target;
+	float densfac = parms->density_strength;
 	float gmin[3], gmax[3];
 	int i;
 	
@@ -1030,14 +1030,15 @@ static void cloth_continuum_step(ClothModifierData *clmd, float dt)
 	hair_get_boundbox(clmd, gmin, gmax);
 	
 	/* gather velocities & density */
-	if (smoothfac > 0.0f || pressfac > 0.0f) {
+	if (smoothfac > 0.0f || densfac > 0.0f) {
 		HairGrid *grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_cell_size, gmin, gmax);
 		BPH_hair_volume_set_debug_data(grid, clmd->debug_data);
 		
+//		BPH_hair_volume_set_debug_value(grid, (int)(spring2->ij == clmd->sim_parms->density_target));
 		cloth_continuum_fill_grid(grid, cloth);
 		
 		/* main hair continuum solver */
-		BPH_hair_volume_solve_divergence(grid, dt);
+		BPH_hair_volume_solve_divergence(grid, dt, denstarget, densfac);
 		
 		for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {
 			float x[3], v[3], nv[3];

source/blender/physics/intern/hair_volume.cpp

@@ -530,10 +530,22 @@ void BPH_hair_volume_normalize_vertex_grid(HairGrid *grid)
 	}
 }
 
-bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt)
+static const float density_threshold = 0.001f; /* cells with density below this are considered empty */
+
+/* Contribution of target density pressure to the laplacian in the pressure poisson equation.
+ * This is based on the model found in
+ * "Two-way Coupled SPH and Particle Level Set Fluid Simulation" (Losasso et al., 2008)
+ */
+BLI_INLINE float hair_volume_density_divergence(float density, float target_density, float strength)
+{
+	if (density > density_threshold && density > target_density)
+		return strength * density * logf(target_density / density);
+	else
+		return 0.0f;
+}
+
+bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_density, float target_strength)
 {
-	const float density_threshold = 0.001f; /* cells with density below this are considered empty */
-	
 	const float flowfac = grid->cellsize / dt;
 	const float inv_flowfac = dt / grid->cellsize;
 	
@@ -576,12 +588,17 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt)
 				float dy = vert_py->velocity[1] - v[1];
 				float dz = vert_pz->velocity[2] - v[2];
 				
+				float divergence = (dx + dy + dz) * flowfac;
+				
+				/* adjustment term for target density */
+				float target = hair_volume_density_divergence(vert->density, target_density, target_strength);
+				
 				/* B vector contains the finite difference approximation of the velocity divergence.
 				 * Note: according to the discretized Navier-Stokes equation the rhs vector
 				 * and resulting pressure gradient should be multiplied by the (inverse) density;
 				 * however, this is already included in the weighting of hair velocities on the grid!
 				 */
-				B[u] = (dx + dy + dz) * flowfac;
+				B[u] = divergence + target;
 			}
 		}
 	}

source/blender/physics/intern/implicit.h

@@ -187,7 +187,7 @@ void BPH_hair_volume_add_segment(struct HairGrid *grid,
 
 void BPH_hair_volume_normalize_vertex_grid(struct HairGrid *grid);
 
-bool BPH_hair_volume_solve_divergence(struct HairGrid *grid, float dt);
+bool BPH_hair_volume_solve_divergence(struct HairGrid *grid, float dt, float target_density, float target_strength);
 #if 0 /* XXX weighting is incorrect, disabled for now */
 void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size);
 #endif