Cloth: Implement angular bending springs

This implements angular bending springs for cloth simulation. This also
adds shearing springs for n-gons.

This angular spring implementation does not include Jacobian matrices,
as the springs can exist between polygons of different vertex counts,
rendering their relationships asymmetrical, and thus impossible to solve
with the current implementation. This means that the bending component
is solved explicitly. However, this is usually not a big problem, as
bending springs contribute less to instability than structural springs.

The the old linear bending model can still be used, and is the default for
existing files, to keep compatibility. However, the new angular bending
model is the default for any new simulation.

This commit makes small breaking changes, in that shearing springs are
now created on n-gons (also in linear bending mode), while n-gons were
previously ignored.

Reviewed By: brecht

Differential Revision: http://developer.blender.org/D3662

1 files changed, 123 insertions, 15 deletions

diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index 677e566ff39..0d897893ff5 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -1664,6 +1664,114 @@ bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, flo
 	}
 }
 
+BLI_INLINE void poly_avg(lfVector *data, int *inds, int len, float r_avg[3])
+{
+	float fact = 1.0f / (float)len;
+
+	zero_v3(r_avg);
+
+	for (int i = 0; i < len; i++) {
+		madd_v3_v3fl(r_avg, data[inds[i]], fact);
+	}
+}
+
+BLI_INLINE void poly_norm(lfVector *data, int i, int j, int *inds, int len, float r_dir[3])
+{
+	float mid[3];
+
+	poly_avg(data, inds, len, mid);
+
+	normal_tri_v3(r_dir, data[i], data[j], mid);
+}
+
+BLI_INLINE void edge_avg(lfVector *data, int i, int j, float r_avg[3])
+{
+	r_avg[0] = (data[i][0] + data[j][0]) * 0.5f;
+	r_avg[1] = (data[i][1] + data[j][1]) * 0.5f;
+	r_avg[2] = (data[i][2] + data[j][2]) * 0.5f;
+}
+
+BLI_INLINE void edge_norm(lfVector *data, int i, int j, float r_dir[3])
+{
+	sub_v3_v3v3(r_dir, data[i], data[j]);
+	normalize_v3(r_dir);
+}
+
+BLI_INLINE float bend_angle(float dir_a[3], float dir_b[3], float dir_e[3])
+{
+	float cos, sin;
+	float tmp[3];
+
+	cos = dot_v3v3(dir_a, dir_b);
+
+	cross_v3_v3v3(tmp, dir_a, dir_b);
+	sin = dot_v3v3(tmp, dir_e);
+
+	return atan2f(sin, cos);
+}
+
+BLI_INLINE void spring_angle(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b,
+                             float r_dir_a[3], float r_dir_b[3],
+                             float *r_angle, float r_vel_a[3], float r_vel_b[3])
+{
+	float dir_e[3], vel_e[3];
+
+	poly_norm(data->X, j, i, i_a, len_a, r_dir_a);
+	poly_norm(data->X, i, j, i_b, len_b, r_dir_b);
+
+	edge_norm(data->X, i, j, dir_e);
+
+	*r_angle = bend_angle(r_dir_a, r_dir_b, dir_e);
+
+	poly_avg(data->V, i_a, len_a, r_vel_a);
+	poly_avg(data->V, i_b, len_b, r_vel_b);
+
+	edge_avg(data->V, i, j, vel_e);
+
+	sub_v3_v3(r_vel_a, vel_e);
+	sub_v3_v3(r_vel_b, vel_e);
+}
+
+/* Angular springs roughly based on the bending model proposed by Baraff and Witkin in "Large Steps in Cloth Simulation". */
+bool BPH_mass_spring_force_spring_angular(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b,
+                                          float restang, float stiffness, float damping)
+{
+	float angle, dir_a[3], dir_b[3], vel_a[3], vel_b[3];
+	float f_a[3], f_b[3], f_e[3];
+	float force;
+	int x;
+
+	spring_angle(data, i, j, i_a, i_b, len_a, len_b,
+	             dir_a, dir_b, &angle, vel_a, vel_b);
+
+	/* spring force */
+	force = stiffness * (angle - restang);
+
+	/* damping force */
+	force += -damping * (dot_v3v3(vel_a, dir_a) + dot_v3v3(vel_b, dir_b));
+
+	mul_v3_v3fl(f_a, dir_a, force / len_a);
+	mul_v3_v3fl(f_b, dir_b, force / len_b);
+
+	for (x = 0; x < len_a; x++) {
+		add_v3_v3(data->F[i_a[x]], f_a);
+	}
+
+	for (x = 0; x < len_b; x++) {
+		add_v3_v3(data->F[i_b[x]], f_b);
+	}
+
+	mul_v3_v3fl(f_a, dir_a, force * 0.5f);
+	mul_v3_v3fl(f_b, dir_b, force * 0.5f);
+
+	add_v3_v3v3(f_e, f_a, f_b);
+
+	sub_v3_v3(data->F[i], f_e);
+	sub_v3_v3(data->F[j], f_e);
+
+	return true;
+}
+
 /* Jacobian of a direction vector.
  * Basically the part of the differential orthogonal to the direction,
  * inversely proportional to the length of the edge.
@@ -1687,7 +1795,7 @@ BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3]
 	}
 }
 
-BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
+BLI_INLINE void spring_hairbend_forces(Implicit_Data *data, int i, int j, int k,
                                       const float goal[3],
                                       float stiffness, float damping,
                                       int q, const float dx[3], const float dv[3],
@@ -1736,7 +1844,7 @@ BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
 }
 
 /* Finite Differences method for estimating the jacobian of the force */
-BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k,
+BLI_INLINE void spring_hairbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k,
                                              const float goal[3],
                                              float stiffness, float damping,
                                              int q, float dfdx[3][3])
@@ -1755,11 +1863,11 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j,
 	/* XXX TODO offset targets to account for position dependency */
 
 	for (a = 0; a < 3; ++a) {
-		spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
+		spring_hairbend_forces(data, i, j, k, goal, stiffness, damping,
 		                      q, dvec_pos[a], dvec_null[a], f);
 		copy_v3_v3(dfdx[a], f);
 
-		spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
+		spring_hairbend_forces(data, i, j, k, goal, stiffness, damping,
 		                      q, dvec_neg[a], dvec_null[a], f);
 		sub_v3_v3(dfdx[a], f);
 
@@ -1770,7 +1878,7 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j,
 }
 
 /* Finite Differences method for estimating the jacobian of the force */
-BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k,
+BLI_INLINE void spring_hairbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k,
                                              const float goal[3],
                                              float stiffness, float damping,
                                              int q, float dfdv[3][3])
@@ -1789,11 +1897,11 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j,
 	/* XXX TODO offset targets to account for position dependency */
 
 	for (a = 0; a < 3; ++a) {
-		spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
+		spring_hairbend_forces(data, i, j, k, goal, stiffness, damping,
 		                      q, dvec_null[a], dvec_pos[a], f);
 		copy_v3_v3(dfdv[a], f);
 
-		spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
+		spring_hairbend_forces(data, i, j, k, goal, stiffness, damping,
 		                      q, dvec_null[a], dvec_neg[a], f);
 		sub_v3_v3(dfdv[a], f);
 
@@ -1806,7 +1914,7 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j,
 /* Angular spring that pulls the vertex toward the local target
  * See "Artistic Simulation of Curly Hair" (Pixar technical memo #12-03a)
  */
-bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, int j, int k,
+bool BPH_mass_spring_force_spring_bending_hair(Implicit_Data *data, int i, int j, int k,
                                                   const float target[3], float stiffness, float damping)
 {
 	float goal[3];
@@ -1822,18 +1930,18 @@ bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, in
 
 	world_to_root_v3(data, j, goal, target);
 
-	spring_angbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk);
+	spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk);
 	negate_v3_v3(fj, fk); /* counterforce */
 
-	spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi);
-	spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj);
-	spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk);
+	spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi);
+	spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj);
+	spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk);
 	copy_m3_m3(dfj_dxi, dfk_dxi); negate_m3(dfj_dxi);
 	copy_m3_m3(dfj_dxj, dfk_dxj); negate_m3(dfj_dxj);
 
-	spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi);
-	spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj);
-	spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk);
+	spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi);
+	spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj);
+	spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk);
 	copy_m3_m3(dfj_dvi, dfk_dvi); negate_m3(dfj_dvi);
 	copy_m3_m3(dfj_dvj, dfk_dvj); negate_m3(dfj_dvj);