Main solver step for generating a divergence-free hair velocity field

on the grid.

This uses the Eigen conjugate-gradient solver to solve the implicit
Poisson equation for the pressure Laplacian:

    div(grad(p)) = div(v)

As described in "Detail Preserving Continuum Simulation of Straight Hair"
(McAdams, Selle, 2009).

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

4 files changed, 239 insertions, 36 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index dfc08b5574b..fbb1cb7cf2c 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -745,7 +745,7 @@ static void cloth_continuum_fill_grid(HairGrid *grid, Cloth *cloth)
 	BPH_hair_volume_normalize_vertex_grid(grid);

 }

 

-static void cloth_continuum_step(ClothModifierData *clmd)

+static void cloth_continuum_step(ClothModifierData *clmd, float dt)

 {

 	ClothSimSettings *parms = clmd->sim_parms;

 	Cloth *cloth = clmd->clothObject;

@@ -774,6 +774,9 @@ static void cloth_continuum_step(ClothModifierData *clmd)
 		

 		cloth_continuum_fill_grid(grid, cloth);

 		

+		/* main hair continuum solver */

+		BPH_hair_volume_solve_divergence(grid, dt);

+		

 		for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {

 			float x[3], v[3], nv[3];

 			

@@ -811,17 +814,19 @@ static void cloth_continuum_step(ClothModifierData *clmd)
 			BKE_sim_debug_data_clear_category(clmd->debug_data, "grid velocity");

 			for (j = 0; j < size; ++j) {

 				for (i = 0; i < size; ++i) {

-					float x[3], v[3], gvel[3], gdensity;

+					float x[3], v[3], gvel[3], gvel_smooth[3], gdensity;

 					

 					madd_v3_v3v3fl(x, offset, a, (float)i / (float)(size-1));

 					madd_v3_v3fl(x, b, (float)j / (float)(size-1));

 					zero_v3(v);

 					

-					BPH_hair_volume_grid_interpolate(grid, x, &gdensity, gvel, NULL, NULL);

+					BPH_hair_volume_grid_interpolate(grid, x, &gdensity, gvel, gvel_smooth, NULL, NULL);

 					

 //					BKE_sim_debug_data_add_circle(clmd->debug_data, x, gdensity, 0.7, 0.3, 1, "grid density", hash_int_2d(hash_int_2d(i, j), 3111));

-					if (!is_zero_v3(gvel))

+					if (!is_zero_v3(gvel) || !is_zero_v3(gvel_smooth)) {

 						BKE_sim_debug_data_add_vector(clmd->debug_data, x, gvel, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));

+						BKE_sim_debug_data_add_vector(clmd->debug_data, x, gvel_smooth, 0.6, 4, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3113));

+					}

 				}

 			}

 		}

@@ -950,7 +955,7 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *
 		BPH_mass_spring_solve_velocities(id, dt, &result);

 		cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);

 		

-		cloth_continuum_step(clmd);

+		cloth_continuum_step(clmd, dt);

 		

 		BPH_mass_spring_solve_positions(id, dt);

 		

diff --git a/source/blender/physics/intern/eigen_utils.h b/source/blender/physics/intern/eigen_utils.h
index 2112386781f..cd3bd5ee45d 100644
--- a/source/blender/physics/intern/eigen_utils.h
+++ b/source/blender/physics/intern/eigen_utils.h
@@ -53,21 +53,21 @@ typedef float Scalar;
 /* slightly extended Eigen vector class
  * with conversion to/from plain C float array
  */
-class fVector : public Eigen::Vector3f {
+class Vector3 : public Eigen::Vector3f {
 public:
 	typedef float *ctype;
 	
-	fVector()
+	Vector3()
 	{
 	}
 	
-	fVector(const ctype &v)
+	Vector3(const ctype &v)
 	{
 		for (int k = 0; k < 3; ++k)
 			coeffRef(k) = v[k];
 	}
 	
-	fVector& operator = (const ctype &v)
+	Vector3& operator = (const ctype &v)
 	{
 		for (int k = 0; k < 3; ++k)
 			coeffRef(k) = v[k];
@@ -83,22 +83,22 @@ public:
 /* slightly extended Eigen matrix class
  * with conversion to/from plain C float array
  */
-class fMatrix : public Eigen::Matrix3f {
+class Matrix3 : public Eigen::Matrix3f {
 public:
 	typedef float (*ctype)[3];
 	
-	fMatrix()
+	Matrix3()
 	{
 	}
 	
-	fMatrix(const ctype &v)
+	Matrix3(const ctype &v)
 	{
 		for (int k = 0; k < 3; ++k)
 			for (int l = 0; l < 3; ++l)
 				coeffRef(l, k) = v[k][l];
 	}
 	
-	fMatrix& operator = (const ctype &v)
+	Matrix3& operator = (const ctype &v)
 	{
 		for (int k = 0; k < 3; ++k)
 			for (int l = 0; l < 3; ++l)
@@ -112,19 +112,21 @@ public:
 	}
 };
 
+typedef Eigen::VectorXf lVector;
+
 /* Extension of dense Eigen vectors,
  * providing 3-float block access for blenlib math functions
  */
-class lVector : public Eigen::VectorXf {
+class lVector3f : public Eigen::VectorXf {
 public:
 	typedef Eigen::VectorXf base_t;
 	
-	lVector()
+	lVector3f()
 	{
 	}
 	
 	template <typename T>
-	lVector& operator = (T rhs)
+	lVector3f& operator = (T rhs)
 	{
 		base_t::operator=(rhs);
 		return *this;
@@ -151,8 +153,8 @@ typedef Eigen::SparseMatrix<Scalar> lMatrix;
  * This should be used for building lMatrix instead of writing to such lMatrix directly (which is very inefficient).
  * After all elements have been defined using the set() method, the actual matrix can be filled using construct().
  */
-struct lMatrixCtor {
-	lMatrixCtor()
+struct lMatrix3fCtor {
+	lMatrix3fCtor()
 	{
 	}
 	
@@ -167,7 +169,7 @@ struct lMatrixCtor {
 		m_trips.reserve(numverts * 9);
 	}
 	
-	void add(int i, int j, const fMatrix &m)
+	void add(int i, int j, const Matrix3 &m)
 	{
 		i *= 3;
 		j *= 3;
@@ -176,7 +178,7 @@ struct lMatrixCtor {
 				m_trips.push_back(Triplet(i + k, j + l, m.coeff(l, k)));
 	}
 	
-	void sub(int i, int j, const fMatrix &m)
+	void sub(int i, int j, const Matrix3 &m)
 	{
 		i *= 3;
 		j *= 3;
@@ -199,7 +201,7 @@ typedef Eigen::ConjugateGradient<lMatrix, Eigen::Lower, Eigen::DiagonalPrecondit
 
 using Eigen::ComputationInfo;
 
-BLI_INLINE void print_lvector(const lVector &v)
+BLI_INLINE void print_lvector(const lVector3f &v)
 {
 	for (int i = 0; i < v.rows(); ++i) {
 		if (i > 0 && i % 3 == 0)
diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp
index cda56946345..cd4e6503965 100644
--- a/source/blender/physics/intern/hair_volume.cpp
+++ b/source/blender/physics/intern/hair_volume.cpp
@@ -39,6 +39,7 @@
 #include "BKE_effect.h"
 
 #include "implicit.h"
+#include "eigen_utils.h"
 
 /* ================ Volumetric Hair Interaction ================
  * adapted from
@@ -113,7 +114,7 @@ BLI_INLINE int hair_grid_interp_weights(const int res[3], const float gmin[3], f
 }
 
 BLI_INLINE void hair_grid_interpolate(const HairGridVert *grid, const int res[3], const float gmin[3], float scale, const float vec[3],
-                                      float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3])
+                                      float *density, float velocity[3], float vel_smooth[3], float density_gradient[3], float velocity_gradient[3][3])
 {
 	HairGridVert data[8];
 	float uvw[3], muvw[3];
@@ -151,6 +152,16 @@ BLI_INLINE void hair_grid_interpolate(const HairGridVert *grid, const int res[3]
 		}
 	}
 	
+	if (vel_smooth) {
+		int k;
+		for (k = 0; k < 3; ++k) {
+			vel_smooth[k] = muvw[2]*( muvw[1]*( muvw[0]*data[0].velocity_smooth[k] + uvw[0]*data[1].velocity_smooth[k] )   +
+			                           uvw[1]*( muvw[0]*data[2].velocity_smooth[k] + uvw[0]*data[3].velocity_smooth[k] ) ) +
+			                 uvw[2]*( muvw[1]*( muvw[0]*data[4].velocity_smooth[k] + uvw[0]*data[5].velocity_smooth[k] )   +
+			                           uvw[1]*( muvw[0]*data[6].velocity_smooth[k] + uvw[0]*data[7].velocity_smooth[k] ) );
+		}
+	}
+	
 	if (density_gradient) {
 		density_gradient[0] = muvw[1] * muvw[2] * ( data[0].density - data[1].density ) +
 		                       uvw[1] * muvw[2] * ( data[2].density - data[3].density ) +
@@ -180,7 +191,7 @@ void BPH_hair_volume_vertex_grid_forces(HairGrid *grid, const float x[3], const
 {
 	float gdensity, gvelocity[3], ggrad[3], gvelgrad[3][3], gradlen;
 	
-	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, ggrad, gvelgrad);
+	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, NULL, ggrad, gvelgrad);
 	
 	zero_v3(f);
 	sub_v3_v3(gvelocity, v);
@@ -199,21 +210,32 @@ void BPH_hair_volume_vertex_grid_forces(HairGrid *grid, const float x[3], const
 }
 
 void BPH_hair_volume_grid_interpolate(HairGrid *grid, const float x[3],
-                                      float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3])
+                                      float *density, float velocity[3], float velocity_smooth[3], float density_gradient[3], float velocity_gradient[3][3])
 {
-	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, density, velocity, density_gradient, velocity_gradient);
+	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, density, velocity, velocity_smooth, density_gradient, velocity_gradient);
 }
 
 void BPH_hair_volume_grid_velocity(HairGrid *grid, const float x[3], const float v[3],
                                    float fluid_factor,
                                    float r_v[3])
 {
-	float gdensity, gvelocity[3], ggrad[3], gvelgrad[3][3];
+	float gdensity, gvelocity[3], gvel_smooth[3], ggrad[3], gvelgrad[3][3];
 	
-	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, ggrad, gvelgrad);
+	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, gvel_smooth, ggrad, gvelgrad);
 	
 	/* XXX TODO implement FLIP method and use fluid_factor to blend between FLIP and PIC */
-	copy_v3_v3(r_v, gvelocity);
+	copy_v3_v3(r_v, gvel_smooth);
+}
+
+void BPH_hair_volume_grid_clear(HairGrid *grid)
+{
+	const int size = hair_grid_size(grid->res);
+	int i;
+	for (i = 0; i < size; ++i) {
+		zero_v3(grid->verts[i].velocity);
+		zero_v3(grid->verts[i].velocity_smooth);
+		grid->verts[i].density = 0.0f;
+	}
 }
 
 BLI_INLINE bool hair_grid_point_valid(const float vec[3], float gmin[3], float gmax[3])
@@ -464,6 +486,182 @@ void BPH_hair_volume_normalize_vertex_grid(HairGrid *grid)
 	}
 }
 
+bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt)
+{
+	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;
+	
+	/*const int num_cells = hair_grid_size(grid->res);*/
+	const int inner_res[3] = { grid->res[0] - 2, grid->res[1] - 2, grid->res[2] - 2 };
+	
+	const int stride0 = 1;
+	const int stride1 = grid->res[0];
+	const int stride2 = grid->res[1] * grid->res[0];
+	
+	const int strideA0 = 1;
+	const int strideA1 = grid->res[0]-2;
+	const int strideA2 = (grid->res[1]-2) * (grid->res[0]-2);
+	
+	/* NB: to avoid many boundary checks, we only solve the system
+	 * for the inner vertices, excluding a 1-cell margin.
+	 */
+	const int inner_cells_start = stride0 + stride1 + stride2;
+	const int num_inner_cells = inner_res[0] * inner_res[1] * inner_res[2];
+	
+	HairGridVert *vert;
+	int i, j, k, u;
+	
+	BLI_assert(num_inner_cells >= 1);
+	
+	/* Calculate divergence */
+	lVector B(num_inner_cells);
+	for (k = 0; k < inner_res[2]; ++k) {
+		for (j = 0; j < inner_res[1]; ++j) {
+			for (i = 0; i < inner_res[0]; ++i) {
+				u = i * strideA0 + j * strideA1 + k * strideA2;
+				vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2;
+				
+				HairGridVert *vert_px = vert + stride0;
+				HairGridVert *vert_py = vert + stride1;
+				HairGridVert *vert_pz = vert + stride2;
+				
+				const float *v = vert->velocity;
+				float dx = vert_px->velocity[0] - v[0];
+				float dy = vert_py->velocity[1] - v[1];
+				float dz = vert_pz->velocity[2] - v[2];
+				
+				/* 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;
+			}
+		}
+	}
+	
+	/* Main Poisson equation system:
+	 * This is derived from the discretezation of the Poisson equation
+	 *   div(grad(p)) = div(v)
+	 * 
+	 * The finite difference approximation yields the linear equation system described here:
+	 * http://en.wikipedia.org/wiki/Discrete_Poisson_equation
+	 */
+	lMatrix A(num_inner_cells, num_inner_cells);
+	/* Reserve space for the base equation system (without boundary conditions).
+	 * Each column contains a factor 6 on the diagonal
+	 * and up to 6 factors -1 on other places.
+	 */
+	A.reserve(Eigen::VectorXi::Constant(num_inner_cells, 7));
+	
+	for (k = 0; k < inner_res[2]; ++k) {
+		for (j = 0; j < inner_res[1]; ++j) {
+			for (i = 0; i < inner_res[0]; ++i) {
+				u = i * strideA0 + j * strideA1 + k * strideA2;
+				vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2;
+				
+				if (vert->density > density_threshold) {
+					int neighbors_lo = 0;
+					int neighbors_hi = 0;
+					int non_solid_neighbors = 0;
+					int neighbor_lo_index[3];
+					int neighbor_hi_index[3];
+					int n;
+					
+					/* check for upper bounds in advance
+					 * to get the correct number of neighbors,
+					 * needed for the diagonal element
+					 */
+					if (k >= 1) {
+						if ((vert - stride2)->density > density_threshold)
+							neighbor_lo_index[neighbors_lo++] = u - strideA2;
+					}
+					if (j >= 1) {
+						if ((vert - stride1)->density > density_threshold)
+							neighbor_lo_index[neighbors_lo++] = u - strideA1;
+					}
+					if (i >= 1) {
+						if ((vert - stride0)->density > density_threshold)
+							neighbor_lo_index[neighbors_lo++] = u - strideA0;
+					}
+					if (i < inner_res[0] - 1) {
+						if ((vert + stride0)->density > density_threshold)
+							neighbor_hi_index[neighbors_hi++] = u + strideA0;
+					}
+					if (j < inner_res[1] - 1) {
+						if ((vert + stride1)->density > density_threshold)
+							neighbor_hi_index[neighbors_hi++] = u + strideA1;
+					}
+					if (k < inner_res[2] - 1) {
+						if ((vert + stride2)->density > density_threshold)
+							neighbor_hi_index[neighbors_hi++] = u + strideA2;
+					}
+					
+					/*int liquid_neighbors = neighbors_lo + neighbors_hi;*/
+					non_solid_neighbors = 6;
+					
+					for (n = 0; n < neighbors_lo; ++n)
+						A.insert(neighbor_lo_index[n], u) = -1.0f;
+					A.insert(u, u) = (float)non_solid_neighbors;
+					for (n = 0; n < neighbors_hi; ++n)
+						A.insert(neighbor_hi_index[n], u) = -1.0f;
+				}
+				else {
+					A.insert(u, u) = 1.0f;
+				}
+			}
+		}
+	}
+	
+	ConjugateGradient cg;
+	cg.setMaxIterations(100);
+	cg.setTolerance(0.01f);
+	
+	cg.compute(A);
+	
+	lVector p = cg.solve(B);
+	
+	if (cg.info() == Eigen::Success) {
+		/* Calculate velocity = grad(p) */
+		for (k = 0; k < inner_res[2]; ++k) {
+			for (j = 0; j < inner_res[1]; ++j) {
+				for (i = 0; i < inner_res[0]; ++i) {
+					u = i * strideA0 + j * strideA1 + k * strideA2;
+					vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2;
+					
+					if (vert->density > density_threshold) {
+						float p0 = p[u];
+						
+						/* finite difference estimate of pressure gradient */
+						float grad_p[3];
+						grad_p[0] = i >= 1 ? p0 - p[u - strideA0] : 0.0f;
+						grad_p[1] = j >= 1 ? p0 - p[u - strideA1] : 0.0f;
+						grad_p[2] = k >= 1 ? p0 - p[u - strideA2] : 0.0f;
+						
+						/* pressure gradient describes velocity delta */
+						madd_v3_v3v3fl(vert->velocity_smooth, vert->velocity, grad_p, inv_flowfac);
+					}
+					else {
+						zero_v3(vert->velocity_smooth);
+					}
+				}
+			}
+		}
+		
+		return true;
+	}
+	else {
+		/* Clear result in case of error */
+		for (i = inner_cells_start, vert = grid->verts + inner_cells_start; i < num_inner_cells; ++i, ++vert) {
+			zero_v3(vert->velocity_smooth);
+		}
+		
+		return false;
+	}
+}
+
 #if 0 /* XXX weighting is incorrect, disabled for now */
 /* Velocity filter kernel
  * See http://en.wikipedia.org/wiki/Filter_%28large_eddy_simulation%29
@@ -584,14 +782,8 @@ HairGrid *BPH_hair_volume_create_vertex_grid(float cellsize, const float gmin[3]
 	copy_v3_v3(grid->gmax, gmax_margin);
 	grid->cellsize = cellsize;
 	grid->inv_cellsize = scale;
-	grid->verts = (HairGridVert *)MEM_mallocN(sizeof(HairGridVert) * size, "hair voxel data");
+	grid->verts = (HairGridVert *)MEM_callocN(sizeof(HairGridVert) * size, "hair voxel data");
 
-	/* initialize grid */
-	for (i = 0; i < size; ++i) {
-		zero_v3(grid->verts[i].velocity);
-		grid->verts[i].density = 0.0f;
-	}
-	
 	return grid;
 }
 
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 373a73cf53f..58883a26c85 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -177,6 +177,7 @@ void BPH_hair_volume_free_vertex_grid(struct HairGrid *grid);
 void BPH_hair_volume_set_debug_data(struct HairGrid *grid, struct SimDebugData *debug_data);
 void BPH_hair_volume_grid_geometry(struct HairGrid *grid, float *cellsize, int res[3], float gmin[3], float gmax[3]);
 
+void BPH_hair_volume_grid_clear(struct HairGrid *grid);
 void BPH_hair_volume_add_vertex(struct HairGrid *grid, const float x[3], const float v[3]);
 void BPH_hair_volume_add_segment(struct HairGrid *grid,
                                  const float x1[3], const float v1[3], const float x2[3], const float v2[3],
@@ -184,12 +185,15 @@ void BPH_hair_volume_add_segment(struct HairGrid *grid,
                                  const float dir1[3], const float dir2[3], const float dir3[3]);
 
 void BPH_hair_volume_normalize_vertex_grid(struct HairGrid *grid);
+
+bool BPH_hair_volume_solve_divergence(struct HairGrid *grid, float dt);
 #if 0 /* XXX weighting is incorrect, disabled for now */
 void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size);
 #endif
 
 void BPH_hair_volume_grid_interpolate(struct HairGrid *grid, const float x[3],
-                                      float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3]);
+                                      float *density, float velocity[3], float velocity_smooth[3],
+                                      float density_gradient[3], float velocity_gradient[3][3]);
 
 /* Effect of fluid simulation grid on velocities.
  * fluid_factor controls blending between PIC (Particle-in-Cell)