[Bf-blender-cvs] [c036c72] master: Proper implementation of angular bending springs including jacobian derivatives for stabilization.

Lukas Tönne noreply at git.blender.org
Tue Jan 20 09:51:18 CET 2015


Commit: c036c72284ea1e0447181274a65d1556d7c4a0d9
Author: Lukas Tönne
Date:   Mon Sep 22 19:27:41 2014 +0200
Branches: master
https://developer.blender.org/rBc036c72284ea1e0447181274a65d1556d7c4a0d9

Proper implementation of angular bending springs including jacobian
derivatives for stabilization.

The bending forces are based on a simplified torsion model where each
neighboring point of a vertex creates a force toward a local goal. This
can be extended later by defining the goals in a local curve frame, so
that natural hair shapes other than perfectly straight hair are
supported.

Calculating the jacobians for the bending forces analytically proved
quite difficult and doesn't work yet, so the fallback method for now
is a straightforward finite difference method. This works very well and
is not too costly. Even the original paper ("Artistic Simulation of
Curly Hair") suggests this approach.

===================================================================

M	source/blender/blenkernel/BKE_cloth.h
M	source/blender/physics/intern/BPH_mass_spring.cpp
M	source/blender/physics/intern/implicit.h
M	source/blender/physics/intern/implicit_blender.c

===================================================================

diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h
index b7be7c6..acda9e1 100644
--- a/source/blender/blenkernel/BKE_cloth.h
+++ b/source/blender/blenkernel/BKE_cloth.h
@@ -134,7 +134,10 @@ typedef struct ClothSpring {
 	int	kl;		/* Pkl from the paper, one end of the spring.	*/
 	int mn;
 	float	restlen;	/* The original length of the spring.	*/
-	int	matrix_index; 	/* needed for implicit solver (fast lookup) */
+	/* 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 4130895..fc7d5ef 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -51,26 +51,58 @@ extern "C" {
 #include "BPH_mass_spring.h"
 #include "implicit.h"
 
+/* Number of off-diagonal non-zero matrix blocks.
+ * Basically there is one of these for each vertex-vertex interaction.
+ */
+static int cloth_count_nondiag_blocks(Cloth *cloth)
+{
+	LinkNode *link;
+	int nondiag = 0;
+	
+	for (link = cloth->springs; link; link = link->next) {
+		ClothSpring *spring = (ClothSpring *)link->link;
+		switch (spring->type) {
+			case CLOTH_SPRING_TYPE_BENDING_ANG:
+				/* angular bending combines 3 vertices */
+				nondiag += 3;
+				break;
+				
+			default:
+				/* all other springs depend on 2 vertices only */
+				nondiag += 1;
+				break;
+		}
+	}
+	
+	return nondiag;
+}
+
 int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd)
 {
 	Cloth *cloth = clmd->clothObject;
 	ClothVertex *verts = cloth->verts;
 	const float ZERO[3] = {0.0f, 0.0f, 0.0f};
 	Implicit_Data *id;
-	unsigned int i;
+	unsigned int i, nondiag;
 	
-	cloth->implicit = id = BPH_mass_spring_solver_create(cloth->numverts, cloth->numsprings);
+	nondiag = cloth_count_nondiag_blocks(cloth);
+	cloth->implicit = id = BPH_mass_spring_solver_create(cloth->numverts, nondiag);
 	
 	for (i = 0; i < cloth->numverts; i++) {
 		BPH_mass_spring_set_vertex_mass(id, i, verts[i].mass);
 	}
 	
 	// init springs 
+	i = 0;
 	LinkNode *link = cloth->springs;
-	for (i = 0; link; link = link->next, ++i) {
+	for (; link; link = link->next) {
 		ClothSpring *spring = (ClothSpring *)link->link;
 		
-		spring->matrix_index = BPH_mass_spring_init_spring(id, i, spring->ij, spring->kl);
+		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++) {
@@ -349,10 +381,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_index, 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->matrix_ij_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_index, 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->matrix_ij_kl, s->restlen, k, parms->Cdis, no_compress, 0.0f, s->f, s->dfdx, s->dfdv);
 		}
 #endif
 	}
@@ -370,7 +402,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_index, goal_x, goal_v, k, parms->goalfrict * 0.01f, s->f, s->dfdx, s->dfdv);
+		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);
 #endif
 	}
 	else if (s->type & CLOTH_SPRING_TYPE_BENDING) {  /* calculate force of bending springs */
@@ -382,7 +414,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);
 		cb = kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 		
-		BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->matrix_index, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);
+		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);
 #endif
 	}
 	else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) {
@@ -394,7 +426,8 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);
 		cb = kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
 		
-		BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->matrix_index, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);
+		/* 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->restlen, s->restlen, kb, cb);
 #endif
 	}
 }
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 2bc491d..8c860d6 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -151,9 +151,8 @@ bool BPH_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int
                                           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 spring_index, float restlen,
-                                                  float stiffness, float damping,
-                                                  float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3]);
+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,
+                                                  float restlen_ij, float restlen_jk, 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],
                                        float stiffness, float damping,
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index b7c4345..7a64df4 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -517,6 +517,32 @@ BLI_INLINE void cross_v3_identity(float r[3][3], const float v[3])
 	r[0][1] = -v[2];	r[1][1] = 0.0f;		r[2][1] = v[0];
 	r[0][2] = v[1];		r[1][2] = -v[0];	r[2][2] = 0.0f;
 }
+
+BLI_INLINE void madd_m3_m3fl(float r[3][3], float m[3][3], float f)
+{
+	r[0][0] += m[0][0] * f;
+	r[0][1] += m[0][1] * f;
+	r[0][2] += m[0][2] * f;
+	r[1][0] += m[1][0] * f;
+	r[1][1] += m[1][1] * f;
+	r[1][2] += m[1][2] * f;
+	r[2][0] += m[2][0] * f;
+	r[2][1] += m[2][1] * f;
+	r[2][2] += m[2][2] * f;
+}
+
+BLI_INLINE void madd_m3_m3m3fl(float r[3][3], float a[3][3], float b[3][3], float f)
+{
+	r[0][0] = a[0][0] + b[0][0] * f;
+	r[0][1] = a[0][1] + b[0][1] * f;
+	r[0][2] = a[0][2] + b[0][2] * f;
+	r[1][0] = a[1][0] + b[1][0] * f;
+	r[1][1] = a[1][1] + b[1][1] * f;
+	r[1][2] = a[1][2] + b[1][2] * f;
+	r[2][0] = a[2][0] + b[2][0] * f;
+	r[2][1] = a[2][1] + b[2][1] * f;
+	r[2][2] = a[2][2] + b[2][2] * f;
+}
 /////////////////////////////////////////////////////////////////
 
 ///////////////////////////
@@ -1620,47 +1646,304 @@ bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, int
 	}
 }
 
-/* Angular spring that pulls the vertex toward the local target
- * See "Artistic Simulation of Curly Hair" (Pixar technical memo #12-03a)
+/* Jacobian of a direction vector.
+ * Basically the part of the differential orthogonal to the direction,
+ * inversely proportional to the length of the edge.
+ * 
+ * dD_ij/dx_i = -dD_ij/dx_j = (D_ij * D_ij^T - I) / len_ij
  */
-bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, int j, int k, int spring_index, float restlen,
-                                                  float stiffness, float damping,
-                                                  float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3])
+BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3], float dir[3], float grad_dir[3][3])
 {
-	float target[3], targetdir[3];
-	float extent[3], dir[3];
-	float dist[3], vel[3];
-	float f[3], dfdx[3][3], dfdv[3][3];
+	float length;
 	
-	sub_v3_v3v3(targetdir, data->X[j], data->X[i]);
-	normalize_v3(targetdir);
-	mul_v3_v3fl(target, targetdir, restlen);
-	{
-		float x[3], d[3];
-		root_to_world_v3(data, j, x, data->X[j]);
-		root_to_world_v3(data, j, d, target);
-		BKE_sim_debug_data_add_vector(data->debug_data, x, d, 1, 0, 0, "spoint", hash_vertex(935, j));
+	sub_v3_v3v3(edge, data->X[j], data->X[i]);
+	length = normalize_v3_v3(dir, edge);
+	
+	if (length > ALMOST_ZERO) {
+		outerproduct(grad_dir, dir, dir);
+		sub_m3_m3m3(grad_dir, I, grad_dir);
+		mul_m3_fl(grad_dir, 1.0f / length);
+	}
+	else {
+		zero_m3(grad_dir);
 	}
+}
+
+BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
+                    

@@ Diff output truncated at 10240 characters. @@




More information about the Bf-blender-cvs mailing list