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

1 files changed, 21 insertions, 4 deletions

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;
 			}
 		}
 	}