Perform grid-based velocity smoothing for hair outside of the implicit

solver step.

Calculating forces and jacobians from linearly interpolated grid values
is problematic due to discontinuities at the grid boundaries. The new
approach of modifying velocities after the backward euler solver step
was suggested in a newer paper

"Detail Preserving Continuum Simulation of Straight Hair"
(McAdams, Selle 2009)

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

5 files changed, 218 insertions, 66 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 153993dc9b8..821b08ce334 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -453,54 +453,6 @@ static void hair_get_boundbox(ClothModifierData *clmd, float gmin[3], float gmax
 	}
 }
 
-static void cloth_calc_volume_force(ClothModifierData *clmd)
-{
-	ClothSimSettings *parms = clmd->sim_parms;
-	Cloth *cloth = clmd->clothObject;
-	Implicit_Data *data = cloth->implicit;
-	int numverts = cloth->numverts;
-	
-	/* 2.0f is an experimental value that seems to give good results */
-	float smoothfac = 2.0f * parms->velocity_smooth;
-	// float collfac = 2.0f * parms->collider_friction;
-	float pressfac = parms->pressure;
-	float minpress = parms->pressure_threshold;
-	float gmin[3], gmax[3];
-	int i;
-	
-	hair_get_boundbox(clmd, gmin, gmax);
-	
-	/* gather velocities & density */
-	if (smoothfac > 0.0f || pressfac > 0.0f) {
-		HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);
-		
-		for (i = 0; i < numverts; i++) {
-			float x[3], v[3];
-			
-			BPH_mass_spring_get_motion_state(data, i, x, v);
-			BPH_hair_volume_add_vertex(vertex_grid, x, v);
-		}
-		BPH_hair_volume_normalize_vertex_grid(vertex_grid);
-		
-#if 0
-		/* apply velocity filter */
-		BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);
-#endif
-		
-		for (i = 0; i < numverts; i++) {
-			float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3];
-			
-			/* calculate volumetric forces */
-			BPH_mass_spring_get_motion_state(data, i, x, v);
-			BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv);
-			/* apply on hair data */
-			BPH_mass_spring_force_extern(data, i, f, dfdx, dfdv);
-		}
-		
-		BPH_hair_volume_free_vertex_grid(vertex_grid);
-	}
-}
-
 static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time)
 {
 	/* Collect forces and derivatives:  F, dFdX, dFdV */
@@ -525,7 +477,7 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB
 	}
 #endif
 
-	cloth_calc_volume_force(clmd);
+	/* cloth_calc_volume_force(clmd); */
 
 #ifdef CLOTH_FORCE_DRAG
 	BPH_mass_spring_force_drag(data, drag);
@@ -570,6 +522,125 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB
 	}
 }
 
+#if 0
+static void cloth_calc_volume_force(ClothModifierData *clmd)
+{
+	ClothSimSettings *parms = clmd->sim_parms;
+	Cloth *cloth = clmd->clothObject;
+	Implicit_Data *data = cloth->implicit;
+	int numverts = cloth->numverts;
+	ClothVertex *vert;
+	
+	/* 2.0f is an experimental value that seems to give good results */
+	float smoothfac = 2.0f * parms->velocity_smooth;
+	float collfac = 2.0f * parms->collider_friction;
+	float pressfac = parms->pressure;
+	float minpress = parms->pressure_threshold;
+	float gmin[3], gmax[3];
+	int i;
+	
+	hair_get_boundbox(clmd, gmin, gmax);
+	
+	/* gather velocities & density */
+	if (smoothfac > 0.0f || pressfac > 0.0f) {
+		HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);
+		
+		vert = cloth->verts;
+		for (i = 0; i < numverts; i++, vert++) {
+			float x[3], v[3];
+			
+			if (vert->solver_index < 0) {
+				copy_v3_v3(x, vert->x);
+				copy_v3_v3(v, vert->v);
+			}
+			else {
+				BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);
+			}
+			BPH_hair_volume_add_vertex(vertex_grid, x, v);
+		}
+		BPH_hair_volume_normalize_vertex_grid(vertex_grid);
+		
+#if 0
+		/* apply velocity filter */
+		BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);
+#endif
+		
+		vert = cloth->verts;
+		for (i = 0; i < numverts; i++, vert++) {
+			float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3];
+			
+			if (vert->solver_index < 0)
+				continue;
+			
+			/* calculate volumetric forces */
+			BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);
+			BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv);
+			/* apply on hair data */
+			BPH_mass_spring_force_extern(data, vert->solver_index, f, dfdx, dfdv);
+		}
+		
+		BPH_hair_volume_free_vertex_grid(vertex_grid);
+	}
+}
+#endif
+
+static void cloth_continuum_step(ClothModifierData *clmd)
+{
+	ClothSimSettings *parms = clmd->sim_parms;
+	Cloth *cloth = clmd->clothObject;
+	Implicit_Data *data = cloth->implicit;
+	int numverts = cloth->numverts;
+	ClothVertex *vert;
+	
+	const float fluid_factor = 0.95f; /* blend between PIC and FLIP methods */
+	/* 2.0f is an experimental value that seems to give good results */
+	float smoothfac = 2.0f * parms->velocity_smooth;
+	float collfac = 2.0f * parms->collider_friction;
+	float pressfac = parms->pressure;
+	float minpress = parms->pressure_threshold;
+	float gmin[3], gmax[3];
+	int i;
+	
+	hair_get_boundbox(clmd, gmin, gmax);
+	
+	/* gather velocities & density */
+	if (smoothfac > 0.0f || pressfac > 0.0f) {
+		HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);
+		
+		vert = cloth->verts;
+		for (i = 0; i < numverts; i++, vert++) {
+			float x[3], v[3];
+			
+			BPH_mass_spring_get_motion_state(data, i, x, v);
+			BPH_hair_volume_add_vertex(vertex_grid, x, v);
+		}
+		BPH_hair_volume_normalize_vertex_grid(vertex_grid);
+		
+#if 0
+		/* apply velocity filter */
+		BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);
+#endif
+		
+		vert = cloth->verts;
+		for (i = 0; i < numverts; i++, vert++) {
+			float x[3], v[3], nv[3];
+			
+			/* calculate volumetric velocity influence */
+			BPH_mass_spring_get_position(data, i, x);
+			BPH_mass_spring_get_new_velocity(data, i, v);
+			
+			BPH_hair_volume_grid_velocity(vertex_grid, x, v, fluid_factor, nv);
+			
+			interp_v3_v3v3(nv, v, nv, smoothfac);
+			
+			/* apply on hair data */
+			BPH_mass_spring_set_new_velocity(data, i, nv);
+		}
+		
+		BPH_hair_volume_free_vertex_grid(vertex_grid);
+	}
+}
+
 static void cloth_clear_result(ClothModifierData *clmd)
 {
 	ClothSolverResult *sres = clmd->solver_result;
@@ -686,9 +757,13 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *
 		cloth_calc_force(clmd, frame, effectors, step);
 		
 		// calculate new velocity and position
-		BPH_mass_spring_solve(id, dt, &result);
+		BPH_mass_spring_solve_velocities(id, dt, &result);
 		cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);
 		
+		cloth_continuum_step(clmd);
+		
+		BPH_mass_spring_solve_positions(id, dt);
+		
 		BPH_mass_spring_apply_result(id);
 		
 		/* move pinned verts to correct position */
diff --git a/source/blender/physics/intern/hair_volume.c b/source/blender/physics/intern/hair_volume.c
index e72df493b79..87f202bbb2c 100644
--- a/source/blender/physics/intern/hair_volume.c
+++ b/source/blender/physics/intern/hair_volume.c
@@ -220,6 +220,18 @@ void BPH_hair_volume_vertex_grid_forces(HairVertexGrid *grid, const float x[3],
 	mul_m3_fl(dfdv, smoothfac);
 }
 
+void BPH_hair_volume_grid_velocity(HairVertexGrid *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];
+	
+	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->scale, x, &gdensity, gvelocity, ggrad, gvelgrad);
+	
+	/* XXX TODO implement FLIP method and use fluid_factor to blend between FLIP and PIC */
+	copy_v3_v3(r_v, gvelocity);
+}
+
 BLI_INLINE bool hair_grid_point_valid(const float vec[3], float gmin[3], float gmax[3])
 {
 	return !(vec[0] < gmin[0] || vec[1] < gmin[1] || vec[2] < gmin[2] ||
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 4678b273dcf..012125b4c87 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -111,20 +111,26 @@ BLI_INLINE int hash_collpair(int type, CollPair *collpair)
 
 void BPH_mass_spring_solver_debug_data(struct Implicit_Data *id, struct SimDebugData *debug_data);
 
+void BPH_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, float mass);
 void BPH_mass_spring_set_rest_transform(struct Implicit_Data *data, int index, float rot[3][3]);
 
 void BPH_mass_spring_set_motion_state(struct Implicit_Data *data, int index, const float x[3], const float v[3]);
 void BPH_mass_spring_set_position(struct Implicit_Data *data, int index, const float x[3]);
 void BPH_mass_spring_set_velocity(struct Implicit_Data *data, int index, const float v[3]);
 void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, float x[3], float v[3]);
-void BPH_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, float mass);
+void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3]);
+
+/* modified velocity during solver step */
+void BPH_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3]);
+void BPH_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3]);
 
 void BPH_mass_spring_clear_constraints(struct Implicit_Data *data);
 void BPH_mass_spring_add_constraint_ndof0(struct Implicit_Data *data, int index, const float dV[3]);
 void BPH_mass_spring_add_constraint_ndof1(struct Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3]);
 void BPH_mass_spring_add_constraint_ndof2(struct Implicit_Data *data, int index, const float c1[3], const float dV[3]);
 
-bool BPH_mass_spring_solve(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result);
+bool BPH_mass_spring_solve_velocities(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result);
+bool BPH_mass_spring_solve_positions(struct Implicit_Data *data, float dt);
 void BPH_mass_spring_apply_result(struct Implicit_Data *data);
 
 /* Clear the force vector at the beginning of the time step */
@@ -170,6 +176,20 @@ void BPH_hair_volume_normalize_vertex_grid(struct HairVertexGrid *grid);
 #if 0 /* XXX weighting is incorrect, disabled for now */
 void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size);
 #endif
+
+/* Effect of fluid simulation grid on velocities.
+ * fluid_factor controls blending between PIC (Particle-in-Cell)
+ *     and FLIP (Fluid-Implicit-Particle) methods (0 = only PIC, 1 = only FLIP)
+ */
+void BPH_hair_volume_grid_velocity(struct HairVertexGrid *grid, const float x[3], const float v[3],
+                                   float fluid_factor,
+                                   float r_v[3]);
+/* XXX Warning: expressing grid effects on velocity as a force is not very stable,
+ * due to discontinuities in interpolated values!
+ * Better use hybrid approaches such as described in
+ * "Detail Preserving Continuum Simulation of Straight Hair"
+ * (McAdams, Selle 2009)
+ */
 void BPH_hair_volume_vertex_grid_forces(struct HairVertexGrid *grid, const float x[3], const float v[3],
                                         float smoothfac, float pressurefac, float minpressure,
                                         float f[3], float dfdx[3][3], float dfdv[3][3]);
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index 6b82dd3f736..865f100e7e6 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -1138,7 +1138,7 @@ static int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector
 }
 #endif
 
-bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *result)
+bool BPH_mass_spring_solve_velocities(Implicit_Data *data, float dt, ImplicitSolverResult *result)
 {
 	unsigned int numverts = data->dFdV[0].vcount;
 
@@ -1163,14 +1163,22 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *
 
 	// advance velocities
 	add_lfvector_lfvector(data->Vnew, data->V, data->dV, numverts);
-	// advance positions
-	add_lfvector_lfvectorS(data->Xnew, data->X, data->Vnew, dt, numverts);
 
 	del_lfvector(dFdXmV);
 	
 	return result->status == BPH_SOLVER_SUCCESS;
 }
 
+bool BPH_mass_spring_solve_positions(Implicit_Data *data, float dt)
+{
+	int numverts = data->M[0].vcount;
+	
+	// advance positions
+	add_lfvector_lfvectorS(data->Xnew, data->X, data->Vnew, dt, numverts);
+	
+	return true;
+}
+
 void BPH_mass_spring_apply_result(Implicit_Data *data)
 {
 	int numverts = data->M[0].vcount;
@@ -1178,6 +1186,12 @@ void BPH_mass_spring_apply_result(Implicit_Data *data)
 	cp_lfvector(data->V, data->Vnew, numverts);
 }
 
+void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
+{
+	unit_m3(data->M[index].m);
+	mul_m3_fl(data->M[index].m, mass);
+}
+
 void BPH_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3])
 {
 #ifdef CLOTH_ROOT_FRAME
@@ -1210,12 +1224,23 @@ void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, flo
 	if (v) root_to_world_v3(data, index, v, data->V[index]);
 }
 
-void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
+void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
 {
-	unit_m3(data->M[index].m);
-	mul_m3_fl(data->M[index].m, mass);
+	root_to_world_v3(data, index, x, data->X[index]);
 }
 
+void BPH_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3])
+{
+	root_to_world_v3(data, index, v, data->Vnew[index]);
+}
+
+void BPH_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3])
+{
+	world_to_root_v3(data, index, data->Vnew[index], v);
+}
+
+/* -------------------------------- */
+
 static int BPH_mass_spring_add_block(Implicit_Data *data, int v1, int v2)
 {
 	int s = data->M[0].vcount + data->num_blocks; /* index from array start */
diff --git a/source/blender/physics/intern/implicit_eigen.cpp b/source/blender/physics/intern/implicit_eigen.cpp
index 4543cc0a3d3..ee95e94a843 100644
--- a/source/blender/physics/intern/implicit_eigen.cpp
+++ b/source/blender/physics/intern/implicit_eigen.cpp
@@ -503,7 +503,7 @@ BLI_INLINE void root_to_world_m3(Implicit_Data *data, int index, float r[3][3],
 
 /* ================================ */
 
-bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *result)
+bool BPH_mass_spring_solve_velocities(Implicit_Data *data, float dt, ImplicitSolverResult *result)
 {
 #ifdef USE_EIGEN_CORE
 	typedef ConjugateGradient solver_t;
@@ -555,7 +555,6 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *
 #endif
 	
 	data->Vnew = data->V + data->dV;
-	data->Xnew = data->X + data->Vnew * dt;
 	
 	switch (cg.info()) {
 		case Eigen::Success:        result->status = BPH_SOLVER_SUCCESS;         break;
@@ -567,7 +566,13 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *
 	result->iterations = cg.iterations();
 	result->error = cg.error();
 	
-	return cg.info() != Eigen::Success;
+	return cg.info() == Eigen::Success;
+}
+
+bool BPH_mass_spring_solve_positions(Implicit_Data *data, float dt)
+{
+	data->Xnew = data->X + data->Vnew * dt;
+	return true;
 }
 
 /* ================================ */
@@ -578,6 +583,14 @@ void BPH_mass_spring_apply_result(Implicit_Data *data)
 	data->V = data->Vnew;
 }
 
+void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
+{
+	float m[3][3];
+	copy_m3_m3(m, I);
+	mul_m3_fl(m, mass);
+	data->iM.add(index, index, m);
+}
+
 void BPH_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3])
 {
 #ifdef CLOTH_ROOT_FRAME
@@ -610,12 +623,19 @@ void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, flo
 	if (v) root_to_world_v3(data, index, v, data->V.v3(index));
 }
 
-void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
+void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
 {
-	float m[3][3];
-	copy_m3_m3(m, I);
-	mul_m3_fl(m, mass);
-	data->iM.add(index, index, m);
+	root_to_world_v3(data, index, x, data->X.v3(index));
+}
+
+void BPH_mass_spring_get_new_velocity(Implicit_Data *data, int index, float v[3])
+{
+	root_to_world_v3(data, index, v, data->V.v3(index));
+}
+
+void BPH_mass_spring_set_new_velocity(Implicit_Data *data, int index, const float v[3])
+{
+	world_to_root_v3(data, index, data->V.v3(index), v);
 }
 
 void BPH_mass_spring_clear_constraints(Implicit_Data *data)