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.

11 files changed, 157 insertions, 9 deletions

diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h
index 838388f428c..b7be7c60af5 100644
--- a/source/blender/blenkernel/BKE_cloth.h
+++ b/source/blender/blenkernel/BKE_cloth.h
@@ -64,6 +64,14 @@ typedef struct ClothHairRoot {
 	float rot[3][3];
 } ClothHairRoot;
 
+typedef struct ClothSolverResult {
+	int status;
+	
+	int max_iterations, min_iterations;
+	float avg_iterations;
+	float max_error, min_error, avg_error;
+} ClothSolverResult;
+
 /**
  * This structure describes a cloth object against which the
  * simulation can run.
diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c
index 53ceaa00703..709032b03eb 100644
--- a/source/blender/blenloader/intern/readfile.c
+++ b/source/blender/blenloader/intern/readfile.c
@@ -4683,6 +4683,7 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 				}
 			}
 			
+			clmd->solver_result = NULL;
 			clmd->debug_data = NULL;
 		}
 		else if (md->type == eModifierType_Fluidsim) {
diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h
index 3dd9ebb5f37..6c6d24cb7a7 100644
--- a/source/blender/makesdna/DNA_modifier_types.h
+++ b/source/blender/makesdna/DNA_modifier_types.h
@@ -567,6 +567,8 @@ typedef struct ClothModifierData {
 	/* XXX nasty hack, remove once hair can be separated from cloth modifier data */
 	struct ClothHairRoot *roots;
 	
+	struct ClothSolverResult *solver_result;
+	
 	struct SimDebugData *debug_data;      /* debug info */
 } ClothModifierData;
 
diff --git a/source/blender/makesrna/intern/CMakeLists.txt b/source/blender/makesrna/intern/CMakeLists.txt
index e64743c5536..0778b14ed25 100644
--- a/source/blender/makesrna/intern/CMakeLists.txt
+++ b/source/blender/makesrna/intern/CMakeLists.txt
@@ -293,6 +293,7 @@ blender_include_dirs(
 	../../ikplugin
 	../../makesdna
 	../../nodes/
+	../../physics
 	../../windowmanager
 	../../editors/include
 	../../render/extern/include
diff --git a/source/blender/makesrna/intern/rna_cloth.c b/source/blender/makesrna/intern/rna_cloth.c
index 154ab7db52d..57dc3ea293d 100644
--- a/source/blender/makesrna/intern/rna_cloth.c
+++ b/source/blender/makesrna/intern/rna_cloth.c
@@ -38,6 +38,8 @@
 #include "BKE_cloth.h"
 #include "BKE_modifier.h"
 
+#include "BPH_mass_spring.h"
+
 #include "WM_api.h"
 #include "WM_types.h"
 
@@ -268,6 +270,64 @@ static char *rna_ClothCollisionSettings_path(PointerRNA *ptr)
 
 #else
 
+static void rna_def_cloth_solver_result(BlenderRNA *brna)
+{
+	StructRNA *srna;
+	PropertyRNA *prop;
+	
+	static EnumPropertyItem status_items[] = {
+	    {BPH_SOLVER_SUCCESS, "SUCCESS", 0, "Success", "Computation was successful"},
+	    {BPH_SOLVER_NUMERICAL_ISSUE, "NUMERICAL_ISSUE", 0, "Numerical Issue", "The provided data did not satisfy the prerequisites"},
+	    {BPH_SOLVER_NO_CONVERGENCE, "NO_CONVERGENCE", 0, "No Convergence", "Iterative procedure did not converge"},
+	    {BPH_SOLVER_INVALID_INPUT, "INVALID_INPUT", 0, "Invalid Input", "The inputs are invalid, or the algorithm has been improperly called"},
+		{0, NULL, 0, NULL, NULL}
+	};
+	
+	srna = RNA_def_struct(brna, "ClothSolverResult", NULL);
+	RNA_def_struct_ui_text(srna, "Solver Result", "Result of cloth solver iteration");
+	
+	RNA_define_verify_sdna(0);
+	
+	prop = RNA_def_property(srna, "status", PROP_ENUM, PROP_NONE);
+	RNA_def_property_enum_items(prop, status_items);
+	RNA_def_property_enum_sdna(prop, NULL, "status");
+	RNA_def_property_flag(prop, PROP_ENUM_FLAG);
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Status", "Status of the solver iteration");
+	
+	prop = RNA_def_property(srna, "max_error", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "max_error");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Maximum Error", "Maximum error during substeps");
+	
+	prop = RNA_def_property(srna, "min_error", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "min_error");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Minimum Error", "Minimum error during substeps");
+	
+	prop = RNA_def_property(srna, "avg_error", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "avg_error");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Average Error", "Average error during substeps");
+	
+	prop = RNA_def_property(srna, "max_iterations", PROP_INT, PROP_NONE);
+	RNA_def_property_int_sdna(prop, NULL, "max_iterations");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Maximum Iterations", "Maximum iterations during substeps");
+	
+	prop = RNA_def_property(srna, "min_iterations", PROP_INT, PROP_NONE);
+	RNA_def_property_int_sdna(prop, NULL, "min_iterations");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Minimum Iterations", "Minimum iterations during substeps");
+	
+	prop = RNA_def_property(srna, "avg_iterations", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "avg_iterations");
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Average Iterations", "Average iterations during substeps");
+	
+	RNA_define_verify_sdna(1);
+}
+
 static void rna_def_cloth_sim_settings(BlenderRNA *brna)
 {
 	StructRNA *srna;
@@ -660,6 +720,7 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
 
 void RNA_def_cloth(BlenderRNA *brna)
 {
+	rna_def_cloth_solver_result(brna);
 	rna_def_cloth_sim_settings(brna);
 	rna_def_cloth_collision_settings(brna);
 }
diff --git a/source/blender/makesrna/intern/rna_modifier.c b/source/blender/makesrna/intern/rna_modifier.c
index 1ba00056ea5..8e7a1d84d9a 100644
--- a/source/blender/makesrna/intern/rna_modifier.c
+++ b/source/blender/makesrna/intern/rna_modifier.c
@@ -2488,6 +2488,11 @@ static void rna_def_modifier_cloth(BlenderRNA *brna)
 	RNA_def_property_pointer_sdna(prop, NULL, "coll_parms");
 	RNA_def_property_ui_text(prop, "Cloth Collision Settings", "");
 	
+	prop = RNA_def_property(srna, "solver_result", PROP_POINTER, PROP_NONE);
+	RNA_def_property_struct_type(prop, "ClothSolverResult");
+	RNA_def_property_pointer_sdna(prop, NULL, "solver_result");
+	RNA_def_property_ui_text(prop, "Solver Result", "");
+	
 	prop = RNA_def_property(srna, "point_cache", PROP_POINTER, PROP_NONE);
 	RNA_def_property_flag(prop, PROP_NEVER_NULL);
 	RNA_def_property_ui_text(prop, "Point Cache", "");
diff --git a/source/blender/modifiers/intern/MOD_cloth.c b/source/blender/modifiers/intern/MOD_cloth.c
index 5c5e38f3444..1b513dd5727 100644
--- a/source/blender/modifiers/intern/MOD_cloth.c
+++ b/source/blender/modifiers/intern/MOD_cloth.c
@@ -179,6 +179,7 @@ static void copyData(ModifierData *md, ModifierData *target)
 	tclmd->point_cache->step = 1;
 	tclmd->clothObject = NULL;
 	tclmd->roots = NULL;
+	tclmd->solver_result = NULL;
 	tclmd->debug_data = NULL;
 }
 
@@ -211,6 +212,9 @@ static void freeData(ModifierData *md)
 		if (clmd->roots)
 			MEM_freeN(clmd->roots);
 		
+		if (clmd->solver_result)
+			MEM_freeN(clmd->solver_result);
+		
 		BKE_sim_debug_data_free(clmd->debug_data);
 	}
 }
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)