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

3 files changed, 152 insertions, 26 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 668e40e71cc..6ea2eeca6f8 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -67,7 +67,7 @@ static int cloth_count_nondiag_blocks(Cloth *cloth)
 	for (link = cloth->springs; link; link = link->next) {
 		ClothSpring *spring = (ClothSpring *)link->link;
 		switch (spring->type) {
-			case CLOTH_SPRING_TYPE_BENDING_ANG:
+			case CLOTH_SPRING_TYPE_BENDING_HAIR:
 				/* angular bending combines 3 vertices */
 				nondiag += 3;
 				break;
@@ -346,14 +346,29 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 
 	s->flags &= ~CLOTH_SPRING_FLAG_NEEDED;
 
-	// calculate force of structural + shear springs
-	if ((s->type & CLOTH_SPRING_TYPE_STRUCTURAL) || (s->type & CLOTH_SPRING_TYPE_SEWING)) {
+	/* Calculate force of bending springs. */
+	if (s->type & CLOTH_SPRING_TYPE_BENDING) {
+#ifdef CLOTH_FORCE_SPRING_BEND
+		float k, scaling;
+
+		s->flags |= CLOTH_SPRING_FLAG_NEEDED;
+
+		scaling = parms->bending + s->ang_stiffness * fabsf(parms->max_bend - parms->bending);
+		k = scaling * s->restlen * 0.1f; /* Multiplying by 0.1, just to scale the forces to more reasonable values. */
+
+		BPH_mass_spring_force_spring_angular(data, s->ij, s->kl, s->pa, s->pb, s->la, s->lb,
+		                                     s->restang, k, parms->bending_damping);
+#endif
+	}
+
+	/* Calculate force of structural + shear springs. */
+	if (s->type & (CLOTH_SPRING_TYPE_STRUCTURAL | CLOTH_SPRING_TYPE_SEWING)) {
 #ifdef CLOTH_FORCE_SPRING_STRUCTURAL
 		float k_tension, scaling_tension;
 
 		s->flags |= CLOTH_SPRING_FLAG_NEEDED;
 
-		scaling_tension = parms->tension + s->stiffness * fabsf(parms->max_tension - parms->tension);
+		scaling_tension = parms->tension + s->lin_stiffness * fabsf(parms->max_tension - parms->tension);
 		k_tension = scaling_tension / (parms->avg_spring_len + FLT_EPSILON);
 
 		if (s->type & CLOTH_SPRING_TYPE_SEWING) {
@@ -365,7 +380,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 		}
 		else {
 			float k_compression, scaling_compression;
-			scaling_compression = parms->compression + s->stiffness * fabsf(parms->max_compression - parms->compression);
+			scaling_compression = parms->compression + s->lin_stiffness * fabsf(parms->max_compression - parms->compression);
 			k_compression = scaling_compression / (parms->avg_spring_len + FLT_EPSILON);
 
 			BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen,
@@ -381,7 +396,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 
 		s->flags |= CLOTH_SPRING_FLAG_NEEDED;
 
-		scaling = parms->shear + s->stiffness * fabsf(parms->max_shear - parms->shear);
+		scaling = parms->shear + s->lin_stiffness * fabsf(parms->max_shear - parms->shear);
 		k = scaling / (parms->avg_spring_len + FLT_EPSILON);
 
 		BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->shear_damp,
@@ -394,7 +409,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 
 		s->flags |= CLOTH_SPRING_FLAG_NEEDED;
 
-		scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);
+		scaling = parms->bending + s->lin_stiffness * fabsf(parms->max_bend - parms->bending);
 		kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 
 		// Fix for [#45084] for cloth stiffness must have cb proportional to kb
@@ -403,7 +418,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 		BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->restlen, kb, cb);
 #endif
 	}
-	else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) {
+	else if (s->type & CLOTH_SPRING_TYPE_BENDING_HAIR) {
 #ifdef CLOTH_FORCE_SPRING_BEND
 		float kb, cb, scaling;
 
@@ -413,14 +428,14 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s)
 		 * this is crap, but needed due to cloth/hair mixing ...
 		 * max_bend factor is not even used for hair, so ...
 		 */
-		scaling = s->stiffness * parms->bending;
+		scaling = s->lin_stiffness * parms->bending;
 		kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 
 		// Fix for [#45084] for cloth stiffness must have cb proportional to kb
 		cb = kb * parms->bending_damping;
 
 		/* XXX assuming same restlen for ij and jk segments here, this can be done correctly for hair later */
-		BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb);
+		BPH_mass_spring_force_spring_bending_hair(data, s->ij, s->kl, s->mn, s->target, kb, cb);
 
 #if 0
 		{
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index ffe9dbbec04..f99812a8aa9 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -118,10 +118,13 @@ bool BPH_mass_spring_force_spring_linear(struct Implicit_Data *data, int i, int
                                          float stiffness_tension, float damping_tension,
                                          float stiffness_compression, float damping_compression,
                                          bool resist_compress, bool new_compress, float clamp_force);
+/* Angular spring force between two polygons */
+bool BPH_mass_spring_force_spring_angular(struct 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);
 /* Bending force, forming a triangle at the base of two structural springs */
 bool BPH_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int j, float restlen, float kb, float cb);
 /* Angular bending force based on local target vectors */
-bool BPH_mass_spring_force_spring_bending_angular(struct Implicit_Data *data, int i, int j, int k,
+bool BPH_mass_spring_force_spring_bending_hair(struct Implicit_Data *data, int i, int j, int k,
                                                   const float target[3], float stiffness, float damping);
 /* Global goal spring */
 bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3],
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);