Basic solver result feedback from the mass-spring (cloth/hair) solver.

This returns a general status (success/no-convergence/other) along with
basic statistics (min/max/average) for the error value and the number
of iterations. It allows some general estimation of the simulation
quality and detection of critical settings that could become a problem.
Better visualization and extended feedback can follow later.

4 files changed, 75 insertions, 9 deletions

diff --git a/source/blender/physics/BPH_mass_spring.h b/source/blender/physics/BPH_mass_spring.h
index 14c03d8e2b0..a1d7af9c023 100644
--- a/source/blender/physics/BPH_mass_spring.h
+++ b/source/blender/physics/BPH_mass_spring.h
@@ -35,9 +35,18 @@ extern "C" {
 struct Implicit_Data;
 struct ClothModifierData;
 
+typedef enum eMassSpringSolverStatus {
+	BPH_SOLVER_SUCCESS              = (1 << 0),
+	BPH_SOLVER_NUMERICAL_ISSUE      = (1 << 1),
+	BPH_SOLVER_NO_CONVERGENCE       = (1 << 2),
+	BPH_SOLVER_INVALID_INPUT        = (1 << 3),
+} eMassSpringSolverStatus;
+
 struct Implicit_Data *BPH_mass_spring_solver_create(int numverts, int numsprings);
 void BPH_mass_spring_solver_free(struct Implicit_Data *id);
 
+const struct MassSpringSolverResult *BPH_mass_spring_solver_result(struct Implicit_Data *data);
+
 int BPH_cloth_solver_init(struct Object *ob, struct ClothModifierData *clmd);
 void BPH_cloth_solver_free(struct ClothModifierData *clmd);
 int BPH_cloth_solve(struct Object *ob, float frame, struct ClothModifierData *clmd, struct ListBase *effectors);
diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 6a64e17aee2..413089547c3 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -530,6 +530,46 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB
 	}
 }
 
+static void cloth_clear_result(ClothModifierData *clmd)
+{
+	ClothSolverResult *sres = clmd->solver_result;
+	
+	sres->status = 0;
+	sres->max_error = sres->min_error = sres->avg_error = 0.0f;
+	sres->max_iterations = sres->min_iterations = 0;
+	sres->avg_iterations = 0.0f;
+}
+
+static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *result, int steps)
+{
+	ClothSolverResult *sres = clmd->solver_result;
+	
+	if (sres->status) { /* already initialized ? */
+		/* error only makes sense for successful iterations */
+		if (result->status == BPH_SOLVER_SUCCESS) {
+			sres->min_error = min_ff(sres->min_error, result->error);
+			sres->max_error = max_ff(sres->max_error, result->error);
+			sres->avg_error += result->error / (float)steps;
+		}
+		
+		sres->min_iterations = min_ii(sres->min_iterations, result->iterations);
+		sres->max_iterations = max_ii(sres->max_iterations, result->iterations);
+		sres->avg_iterations += (float)result->iterations / (float)steps;
+	}
+	else {
+		/* error only makes sense for successful iterations */
+		if (result->status == BPH_SOLVER_SUCCESS) {
+			sres->min_error = sres->max_error = result->error;
+			sres->avg_error += result->error / (float)steps;
+		}
+		
+		sres->min_iterations = sres->max_iterations  = result->iterations;
+		sres->avg_iterations += (float)result->iterations / (float)steps;
+	}
+	
+	sres->status |= result->status;
+}
+
 int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors)
 {
 	unsigned int i=0;
@@ -546,6 +586,10 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *
 	
 	BKE_sim_debug_data_clear_category(clmd->debug_data, "collision");
 	
+	if (!clmd->solver_result)
+		clmd->solver_result = (ClothSolverResult *)MEM_callocN(sizeof(ClothSolverResult), "cloth solver result");
+	cloth_clear_result(clmd);
+	
 	if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) { /* do goal stuff */
 		for (i = 0; i < numverts; i++) {
 			// update velocities with constrained velocities from pinned verts
@@ -565,6 +609,7 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *
 	}
 	
 	while (step < tf) {
+		ImplicitSolverResult result;
 		
 		/* copy velocities for collision */
 		for (i = 0; i < numverts; i++) {
@@ -597,7 +642,8 @@ 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);
+		BPH_mass_spring_solve(id, dt, &result);
+		cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);
 		
 		BPH_mass_spring_apply_result(id);
 		
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index b11bc5d18a9..2bc491d4266 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -61,6 +61,13 @@ extern "C" {
 struct Implicit_Data;
 struct SimDebugData;
 
+typedef struct ImplicitSolverResult {
+	int status;
+	
+	int iterations;
+	float error;
+} ImplicitSolverResult;
+
 BLI_INLINE void implicit_print_matrix_elem(float v)
 {
     printf("%-8.3f", v);
@@ -118,7 +125,7 @@ void BPH_mass_spring_add_constraint_ndof0(struct Implicit_Data *data, int index,
 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);
+bool BPH_mass_spring_solve(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result);
 void BPH_mass_spring_apply_result(struct Implicit_Data *data);
 
 /* Clear the force vector at the beginning of the time step */
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index 425dcf9a384..b7c4345da37 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -834,7 +834,7 @@ static int  cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z
 }
 #endif
 
-static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S)
+static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, ImplicitSolverResult *result)
 {
 	// Solves for unknown X in equation AX=B
 	unsigned int conjgrad_loopcount=0, conjgrad_looplimit=100;
@@ -847,14 +847,15 @@ static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z,
 	lfVector *c = create_lfvector(numverts);
 	lfVector *q = create_lfvector(numverts);
 	lfVector *s = create_lfvector(numverts);
-	float delta_new, delta_old, delta_target, alpha;
+	float bnorm2, delta_new, delta_old, delta_target, alpha;
 	
 	cp_lfvector(ldV, z, numverts);
 	
 	/* d0 = filter(B)^T * P * filter(B) */
 	cp_lfvector(fB, lB, numverts);
 	filter(fB, S);
-	delta_target = conjgrad_epsilon*conjgrad_epsilon * dot_lfvector(fB, fB, numverts);
+	bnorm2 = dot_lfvector(fB, fB, numverts);
+	delta_target = conjgrad_epsilon*conjgrad_epsilon * bnorm2;
 	
 	/* r = filter(B - A * dV) */
 	mul_bfmatrix_lfvector(AdV, lA, ldV);
@@ -914,6 +915,10 @@ static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z,
 	del_lfvector(s);
 	// printf("W/O conjgrad_loopcount: %d\n", conjgrad_loopcount);
 
+	result->status = conjgrad_loopcount < conjgrad_looplimit ? BPH_SOLVER_SUCCESS : BPH_SOLVER_NO_CONVERGENCE;
+	result->iterations = conjgrad_loopcount;
+	result->error = bnorm2 > 0.0f ? sqrt(delta_new / bnorm2) : 0.0f;
+
 	return conjgrad_loopcount < conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
 }
 
@@ -1105,10 +1110,9 @@ static int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector
 }
 #endif
 
-bool BPH_mass_spring_solve(Implicit_Data *data, float dt)
+bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *result)
 {
 	unsigned int numverts = data->dFdV[0].vcount;
-	bool ok;
 
 	lfVector *dFdXmV = create_lfvector(numverts);
 	zero_lfvector(data->dV, numverts);
@@ -1123,7 +1127,7 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt)
 
 	// itstart();
 
-	ok = cg_filtered(data->dV, data->A, data->B, data->z, data->S); /* conjugate gradient algorithm to solve Ax=b */
+	cg_filtered(data->dV, data->A, data->B, data->z, data->S, result); /* conjugate gradient algorithm to solve Ax=b */
 	// cg_filtered_pre(id->dV, id->A, id->B, id->z, id->S, id->P, id->Pinv, id->bigI);
 
 	// itend();
@@ -1136,7 +1140,7 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt)
 
 	del_lfvector(dFdXmV);
 	
-	return ok;
+	return result->status == BPH_SOLVER_SUCCESS;
 }
 
 void BPH_mass_spring_apply_result(Implicit_Data *data)