Removed block matrix indices storage from ClothSpring struct.

This is not necessary: the implicit solver data can keep track instead
of how many off-diagonal matrix blocks are in use (provided the
allocation limit is calculated correctly). Every time a spring is
created it then simply increments this counter and uses the block index
locally - no need to store this persistently.

4 files changed, 33 insertions, 38 deletions

diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h
index f1b3f03037a..b3d0c46779d 100644
--- a/source/blender/blenkernel/BKE_cloth.h
+++ b/source/blender/blenkernel/BKE_cloth.h
@@ -134,10 +134,6 @@ typedef struct ClothSpring {
 	int	kl;		/* Pkl from the paper, one end of the spring.	*/
 	int mn;
 	float	restlen;	/* The original length of the spring.	*/
-	/* needed for implicit solver (fast lookup) */
-	int	matrix_ij_kl;
-	int	matrix_kl_mn;
-	int	matrix_ij_mn;
 	int	type;		/* types defined in BKE_cloth.h ("springType") */
 	int	flags; 		/* defined in BKE_cloth.h, e.g. deactivated due to tearing */
 	float dfdx[3][3];
diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 6a768f10d11..72639fbfc74 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -92,19 +92,6 @@ int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd)
 		BPH_mass_spring_set_vertex_mass(id, i, verts[i].mass);
 	}
 	
-	// init springs 
-	i = 0;
-	LinkNode *link = cloth->springs;
-	for (; link; link = link->next) {
-		ClothSpring *spring = (ClothSpring *)link->link;
-		
-		spring->matrix_ij_kl = BPH_mass_spring_init_spring(id, i++, spring->ij, spring->kl);
-		if (spring->type == CLOTH_SPRING_TYPE_BENDING_ANG) {
-			spring->matrix_kl_mn = BPH_mass_spring_init_spring(id, i++, spring->kl, spring->mn);
-			spring->matrix_ij_mn = BPH_mass_spring_init_spring(id, i++, spring->ij, spring->mn);
-		}
-	}
-	
 	for (i = 0; i < cloth->numverts; i++) {
 		BPH_mass_spring_set_motion_state(id, i, verts[i].x, ZERO);
 	}
@@ -379,10 +366,10 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		if (s->type & CLOTH_SPRING_TYPE_SEWING) {
 			// TODO: verify, half verified (couldn't see error)
 			// sewing springs usually have a large distance at first so clamp the force so we don't get tunnelling through colission objects
-			BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->matrix_ij_kl, s->restlen, k, parms->Cdis, no_compress, parms->max_sewing, s->f, s->dfdx, s->dfdv);
+			BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, parms->max_sewing, s->f, s->dfdx, s->dfdv);
 		}
 		else {
-			BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->matrix_ij_kl, s->restlen, k, parms->Cdis, no_compress, 0.0f, s->f, s->dfdx, s->dfdv);
+			BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, 0.0f, s->f, s->dfdx, s->dfdv);
 		}
 #endif
 	}
@@ -400,7 +387,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		scaling = parms->goalspring + s->stiffness * fabsf(parms->max_struct - parms->goalspring);
 		k = verts[s->ij].goal * scaling / (parms->avg_spring_len + FLT_EPSILON);
 		
-		BPH_mass_spring_force_spring_goal(data, s->ij, s->matrix_ij_kl, goal_x, goal_v, k, parms->goalfrict * 0.01f, s->f, s->dfdx, s->dfdv);
+		BPH_mass_spring_force_spring_goal(data, s->ij, goal_x, goal_v, k, parms->goalfrict * 0.01f, s->f, s->dfdx, s->dfdv);
 #endif
 	}
 	else if (s->type & CLOTH_SPRING_TYPE_BENDING) {  /* calculate force of bending springs */
@@ -415,7 +402,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		scaling = parms->bending_damping;
 		cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 		
-		BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->matrix_ij_kl, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);
+		BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);
 #endif
 	}
 	else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) {
@@ -431,7 +418,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 		
 		/* 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->matrix_ij_kl, s->matrix_kl_mn, s->matrix_ij_mn, s->target, kb, cb);
+		BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb);
 		
 		{
 			float x_kl[3], x_mn[3], v[3], d[3];
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 40b6174e250..ad714debe73 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -118,7 +118,7 @@ void BPH_mass_spring_set_velocity(struct Implicit_Data *data, int index, const f
 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);
 
-int BPH_mass_spring_init_spring(struct Implicit_Data *data, int index, int v1, int v2);
+int BPH_mass_spring_add_block(struct Implicit_Data *data, int v1, int v2);
 
 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]);
@@ -143,18 +143,18 @@ void BPH_mass_spring_force_face_wind(struct Implicit_Data *data, int v1, int v2,
 /* Wind force, acting on an edge */
 void BPH_mass_spring_force_edge_wind(struct Implicit_Data *data, int v1, int v2, const float (*winvec)[3]);
 /* Linear spring force between two points */
-bool BPH_mass_spring_force_spring_linear(struct Implicit_Data *data, int i, int j, int spring_index, float restlen,
+bool BPH_mass_spring_force_spring_linear(struct Implicit_Data *data, int i, int j, float restlen,
                                          float stiffness, float damping, bool no_compress, float clamp_force,
                                          float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3]);
 /* 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, int spring_index, float restlen,
+bool BPH_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int j, float restlen,
                                           float kb, float cb,
                                           float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3]);
 /* 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, int block_ij, int block_jk, int block_ik,
+bool BPH_mass_spring_force_spring_bending_angular(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, int spring_index, const float goal_x[3], const float goal_v[3],
+bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3],
                                        float stiffness, float damping,
                                        float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3]);
 
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index 51baf1a79ae..75dd80ee9c0 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -682,6 +682,7 @@ typedef struct Implicit_Data  {
 	fmatrix3x3 *M;				/* masses */
 	lfVector *F;				/* forces */
 	fmatrix3x3 *dFdV, *dFdX;	/* force jacobians */
+	int num_blocks;				/* number of off-diagonal blocks (springs) */
 	
 	/* motion state data */
 	lfVector *X, *Xnew;			/* positions */
@@ -1214,9 +1215,11 @@ void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
 	mul_m3_fl(data->M[index].m, mass);
 }
 
-int BPH_mass_spring_init_spring(Implicit_Data *data, int index, int v1, int v2)
+int BPH_mass_spring_add_block(Implicit_Data *data, int v1, int v2)
 {
-	int s = data->M[0].vcount + index; /* index from array start */
+	int s = data->M[0].vcount + data->num_blocks; /* index from array start */
+	BLI_assert(s < data->M[0].vcount + data->M[0].scount);
+	++data->num_blocks;
 	
 	/* tfm and S don't have spring entries (diagonal blocks only) */
 	init_fmatrix(data->bigI + s, v1, v2);
@@ -1287,6 +1290,8 @@ void BPH_mass_spring_clear_forces(Implicit_Data *data)
 	zero_lfvector(data->F, numverts);
 	init_bfmatrix(data->dFdX, ZERO);
 	init_bfmatrix(data->dFdV, ZERO);
+	
+	data->num_blocks = 0;
 }
 
 void BPH_mass_spring_force_reference_frame(Implicit_Data *data, int index, const float acceleration[3], const float omega[3], const float domega_dt[3])
@@ -1552,21 +1557,23 @@ BLI_INLINE bool spring_length(Implicit_Data *data, int i, int j, float r_extent[
 	return true;
 }
 
-BLI_INLINE void apply_spring(Implicit_Data *data, int i, int j, int spring_index, const float f[3], float dfdx[3][3], float dfdv[3][3])
+BLI_INLINE void apply_spring(Implicit_Data *data, int i, int j, const float f[3], float dfdx[3][3], float dfdv[3][3])
 {
+	int block_ij = BPH_mass_spring_add_block(data, i, j);
+	
 	add_v3_v3(data->F[i], f);
 	sub_v3_v3(data->F[j], f);
 	
 	add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, dfdx);
 	add_m3_m3m3(data->dFdX[j].m, data->dFdX[j].m, dfdx);
-	sub_m3_m3m3(data->dFdX[spring_index].m, data->dFdX[spring_index].m, dfdx);
+	sub_m3_m3m3(data->dFdX[block_ij].m, data->dFdX[block_ij].m, dfdx);
 	
 	add_m3_m3m3(data->dFdV[i].m, data->dFdV[i].m, dfdv);
 	add_m3_m3m3(data->dFdV[j].m, data->dFdV[j].m, dfdv);
-	sub_m3_m3m3(data->dFdV[spring_index].m, data->dFdV[spring_index].m, dfdv);
+	sub_m3_m3m3(data->dFdV[block_ij].m, data->dFdV[block_ij].m, dfdv);
 }
 
-bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, int spring_index, float restlen,
+bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, float restlen,
                                          float stiffness, float damping, bool no_compress, float clamp_force,
                                          float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3])
 {
@@ -1591,7 +1598,7 @@ bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, int
 		dfdx_spring(dfdx, dir, length, restlen, stiffness);
 		dfdv_damp(dfdv, dir, damping);
 		
-		apply_spring(data, i, j, spring_index, f, dfdx, dfdv);
+		apply_spring(data, i, j, f, dfdx, dfdv);
 		
 		if (r_f) copy_v3_v3(r_f, f);
 		if (r_dfdx) copy_m3_m3(r_dfdx, dfdx);
@@ -1609,7 +1616,7 @@ bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, int
 }
 
 /* See "Stable but Responsive Cloth" (Choi, Ko 2005) */
-bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, int spring_index, float restlen,
+bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, float restlen,
                                           float kb, float cb,
                                           float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3])
 {
@@ -1629,7 +1636,7 @@ bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, int
 		/* XXX damping not supported */
 		zero_m3(dfdv);
 		
-		apply_spring(data, i, j, spring_index, f, dfdx, dfdv);
+		apply_spring(data, i, j, f, dfdx, dfdv);
 		
 		if (r_f) copy_v3_v3(r_f, f);
 		if (r_dfdx) copy_m3_m3(r_dfdx, dfdx);
@@ -1788,7 +1795,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, int block_ij, int block_jk, int block_ik,
+bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, int j, int k,
                                                   const float target[3], float stiffness, float damping)
 {
 	float goal[3];
@@ -1798,6 +1805,10 @@ bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, in
 	
 	const float vecnull[3] = {0.0f, 0.0f, 0.0f};
 	
+	int block_ij = BPH_mass_spring_add_block(data, i, j);
+	int block_jk = BPH_mass_spring_add_block(data, j, k);
+	int block_ik = BPH_mass_spring_add_block(data, i, k);
+	
 	world_to_root_v3(data, j, goal, target);
 	
 	spring_angbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk);
@@ -1816,6 +1827,7 @@ bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, in
 	copy_m3_m3(dfj_dvj, dfk_dvj); negate_m3(dfj_dvj);
 	
 	/* add forces and jacobians to the solver data */
+	
 	add_v3_v3(data->F[j], fj);
 	add_v3_v3(data->F[k], fk);
 	
@@ -1924,7 +1936,7 @@ bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, in
 	return true;
 }
 
-bool BPH_mass_spring_force_spring_goal(Implicit_Data *data, int i, int UNUSED(spring_index), const float goal_x[3], const float goal_v[3],
+bool BPH_mass_spring_force_spring_goal(Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3],
                                        float stiffness, float damping,
                                        float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3])
 {