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

3 files changed, 27 insertions, 9 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index a20f1617505..d911a431b23 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/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];

diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp
index 5dadd45c154..4d10c042f46 100644
--- a/source/blender/physics/intern/hair_volume.cpp
+++ b/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;
 			}
 		}
 	}
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 203e6f3c235..246b916e353 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/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