svn merge -r 15988:16077 https://svn.blender.org/svnroot/bf-blender/trunk/blender

To have the 50% faster nearest_surface point.
Changed mesh_faces_nearest_point to return the face normal instead of collision normal

6 files changed, 975 insertions, 371 deletions

diff --git a/source/blender/blenkernel/intern/bvhutils.c b/source/blender/blenkernel/intern/bvhutils.c
index 10e92b8b705..ae449843d2a 100644
--- a/source/blender/blenkernel/intern/bvhutils.c
+++ b/source/blender/blenkernel/intern/bvhutils.c
@@ -48,9 +48,6 @@
 
 /* Math stuff for ray casting on mesh faces and for nearest surface */
 
-static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest);
-
-#define ISECT_EPSILON 1e-6
 static float ray_tri_intersection(const BVHTreeRay *ray, const float m_dist, const float *v0, const float *v1, const float *v2)
 {
 	float dist;
@@ -79,170 +76,324 @@ static float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, con
 	return FLT_MAX;
 }
 
+
 /*
- * This calculates the distance from point to the plane
- * Distance is negative if point is on the back side of plane
+ * Function adapted from David Eberly's distance tools (LGPL)
+ * http://www.geometrictools.com/LibFoundation/Distance/Distance.html
  */
-static float point_plane_distance(const float *point, const float *plane_point, const float *plane_normal)
+static float nearest_point_in_tri_surface(const float *v0,const float *v1,const float *v2,const float *p, int *v, int *e, float *nearest )
 {
-	float pp[3];
-	VECSUB(pp, point, plane_point);
-	return INPR(pp, plane_normal);
-}
-static float choose_nearest(const float v0[2], const float v1[2], const float point[2], float closest[2])
-{
-	float d[2][2], sdist[2];
-	VECSUB2D(d[0], v0, point);
-	VECSUB2D(d[1], v1, point);
-
-	sdist[0] = d[0][0]*d[0][0] + d[0][1]*d[0][1];
-	sdist[1] = d[1][0]*d[1][0] + d[1][1]*d[1][1];
-
-	if(sdist[0] < sdist[1])
+	float diff[3];
+	float e0[3];
+	float e1[3];
+	float A00;
+	float A01;
+	float A11;
+	float B0;
+	float B1;
+	float C;
+	float Det;
+	float S;
+	float T;
+	float sqrDist;
+	int lv = -1, le = -1;
+	
+	VECSUB(diff, v0, p);
+	VECSUB(e0, v1, v0);
+	VECSUB(e1, v2, v0);
+	
+	A00 = INPR ( e0, e0 );
+	A01 = INPR( e0, e1 );
+	A11 = INPR ( e1, e1 );
+	B0 = INPR( diff, e0 );
+	B1 = INPR( diff, e1 );
+	C = INPR( diff, diff );
+	Det = fabs( A00 * A11 - A01 * A01 );
+	S = A01 * B1 - A11 * B0;
+	T = A01 * B0 - A00 * B1;
+
+	if ( S + T <= Det )
 	{
-		if(closest)
-			VECCOPY2D(closest, v0);
-		return sdist[0];
+		if ( S < 0.0f )
+		{
+			if ( T < 0.0f )  // Region 4
+			{
+				if ( B0 < 0.0f )
+				{
+					T = 0.0f;
+					if ( -B0 >= A00 )
+					{
+						S = (float)1.0;
+						sqrDist = A00 + 2.0f * B0 + C;
+						lv = 1;
+					}
+					else
+					{
+						if(fabs(A00) > FLT_EPSILON)
+							S = -B0/A00;
+						else
+							S = 0.0f;
+						sqrDist = B0 * S + C;
+						le = 0;
+					}
+				}
+				else
+				{
+					S = 0.0f;
+					if ( B1 >= 0.0f )
+					{
+						T = 0.0f;
+						sqrDist = C;
+						lv = 0;
+					}
+					else if ( -B1 >= A11 )
+					{
+						T = 1.0f;
+						sqrDist = A11 + 2.0f * B1 + C;
+						lv = 2;
+					}
+					else
+					{
+						if(fabs(A11) > FLT_EPSILON)
+							T = -B1 / A11;
+						else
+							T = 0.0f;
+						sqrDist = B1 * T + C;
+						le = 1;
+					}
+				}
+			}
+			else  // Region 3
+			{
+				S = 0.0f;
+				if ( B1 >= 0.0f )
+				{
+					T = 0.0f;
+					sqrDist = C;
+					lv = 0;
+				}
+				else if ( -B1 >= A11 )
+				{
+					T = 1.0f;
+					sqrDist = A11 + 2.0f * B1 + C;
+					lv = 2;
+				}
+				else
+				{
+					if(fabs(A11) > FLT_EPSILON)
+						T = -B1 / A11;
+					else
+						T = 0.0;
+					sqrDist = B1 * T + C;
+					le = 1;
+				}
+			}
+		}
+		else if ( T < 0.0f )  // Region 5
+		{
+			T = 0.0f;
+			if ( B0 >= 0.0f )
+			{
+				S = 0.0f;
+				sqrDist = C;
+				lv = 0;
+			}
+			else if ( -B0 >= A00 )
+			{
+				S = 1.0f;
+				sqrDist = A00 + 2.0f * B0 + C;
+				lv = 1;
+			}
+			else
+			{
+				if(fabs(A00) > FLT_EPSILON)
+					S = -B0 / A00;
+				else
+					S = 0.0f;
+				sqrDist = B0 * S + C;
+				le = 0;
+			}
+		}
+		else  // Region 0
+		{
+            // Minimum at interior lv
+			float invDet;
+			if(fabs(Det) > FLT_EPSILON)
+				invDet = 1.0f / Det;
+			else
+				invDet = 0.0f;
+			S *= invDet;
+			T *= invDet;
+			sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0) +
+					T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
+		}
 	}
 	else
 	{
-		if(closest)
-			VECCOPY2D(closest, v1);
-		return sdist[1];
-	}
-}
-/*
- * calculates the closest point between point-tri (2D)
- * returns that tri must be right-handed
- * Returns square distance
- */
-static float closest_point_in_tri2D(const float point[2], /*const*/ float tri[3][2], float closest[2])
-{
-	float edge_di[2];
-	float v_point[2];
-	float proj[2];					//point projected over edge-dir, edge-normal (witouth normalized edge)
-	const float *v0 = tri[2], *v1;
-	float edge_slen, d;				//edge squared length
-	int i;
-	const float *nearest_vertex = NULL;
-
-
-	//for each edge
-	for(i=0, v0=tri[2], v1=tri[0]; i < 3; v0=tri[i++], v1=tri[i])
-	{
-		VECSUB2D(edge_di,    v1, v0);
-		VECSUB2D(v_point, point, v0);
-
-		proj[1] =  v_point[0]*edge_di[1] - v_point[1]*edge_di[0];	//dot product with edge normal
-
-		//point inside this edge
-		if(proj[1] < 0)
-			continue;
-
-		proj[0] = v_point[0]*edge_di[0] + v_point[1]*edge_di[1];
+		float tmp0, tmp1, numer, denom;
 
-		//closest to this edge is v0
-		if(proj[0] < 0)
+		if ( S < 0.0f )  // Region 2
 		{
-			if(nearest_vertex == NULL || nearest_vertex == v0)
-				nearest_vertex = v0;
+			tmp0 = A01 + B0;
+			tmp1 = A11 + B1;
+			if ( tmp1 > tmp0 )
+			{
+				numer = tmp1 - tmp0;
+				denom = A00 - 2.0f * A01 + A11;
+				if ( numer >= denom )
+				{
+					S = 1.0f;
+					T = 0.0f;
+					sqrDist = A00 + 2.0f * B0 + C;
+					lv = 1;
+				}
+				else
+				{
+					if(fabs(denom) > FLT_EPSILON)
+						S = numer / denom;
+					else
+						S = 0.0f;
+					T = 1.0f - S;
+					sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
+							T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
+					le = 2;
+				}
+			}
 			else
 			{
-				//choose nearest
-				return choose_nearest(nearest_vertex, v0, point, closest);
+				S = 0.0f;
+				if ( tmp1 <= 0.0f )
+				{
+					T = 1.0f;
+					sqrDist = A11 + 2.0f * B1 + C;
+					lv = 2;
+				}
+				else if ( B1 >= 0.0f )
+				{
+					T = 0.0f;
+					sqrDist = C;
+					lv = 0;
+				}
+				else
+				{
+					if(fabs(A11) > FLT_EPSILON)
+						T = -B1 / A11;
+					else
+						T = 0.0f;
+					sqrDist = B1 * T + C;
+					le = 1;
+				}
 			}
-			i++;	//We can skip next edge
-			continue;
 		}
-
-		edge_slen = edge_di[0]*edge_di[0] + edge_di[1]*edge_di[1];	//squared edge len
-		//closest to this edge is v1
-		if(proj[0] > edge_slen)
+		else if ( T < 0.0f )  // Region 6
 		{
-			if(nearest_vertex == NULL || nearest_vertex == v1)
-				nearest_vertex = v1;
+			tmp0 = A01 + B1;
+			tmp1 = A00 + B0;
+			if ( tmp1 > tmp0 )
+			{
+				numer = tmp1 - tmp0;
+				denom = A00 - 2.0f * A01 + A11;
+				if ( numer >= denom )
+				{
+					T = 1.0f;
+					S = 0.0f;
+					sqrDist = A11 + 2.0f * B1 + C;
+					lv = 2;
+				}
+				else
+				{
+					if(fabs(denom) > FLT_EPSILON)
+						T = numer / denom;
+					else
+						T = 0.0f;
+					S = 1.0f - T;
+					sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
+							T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
+					le = 2;
+				}
+			}
 			else
 			{
-				return choose_nearest(nearest_vertex, v1, point, closest);
+				T = 0.0f;
+				if ( tmp1 <= 0.0f )
+				{
+					S = 1.0f;
+					sqrDist = A00 + 2.0f * B0 + C;
+					lv = 1;
+				}
+				else if ( B0 >= 0.0f )
+				{
+					S = 0.0f;
+					sqrDist = C;
+					lv = 0;
+				}
+				else
+				{
+					if(fabs(A00) > FLT_EPSILON)
+						S = -B0 / A00;
+					else
+						S = 0.0f;
+					sqrDist = B0 * S + C;
+					le = 0;
+				}
+			}
+		}
+		else  // Region 1
+		{
+			numer = A11 + B1 - A01 - B0;
+			if ( numer <= 0.0f )
+			{
+				S = 0.0f;
+				T = 1.0f;
+				sqrDist = A11 + 2.0f * B1 + C;
+				lv = 2;
+			}
+			else
+			{
+				denom = A00 - 2.0f * A01 + A11;
+				if ( numer >= denom )
+				{
+					S = 1.0f;
+					T = 0.0f;
+					sqrDist = A00 + 2.0f * B0 + C;
+					lv = 1;
+				}
+				else
+				{
+					if(fabs(denom) > FLT_EPSILON)
+						S = numer / denom;
+					else
+						S = 0.0f;
+					T = 1.0f - S;
+					sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
+							T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
+					le = 2;
+				}
 			}
-			continue;
 		}
-
-		//nearest is on this edge
-		d= proj[1] / edge_slen;
-		closest[0] = point[0] - edge_di[1] * d;
-		closest[1] = point[1] + edge_di[0] * d;
-
-		return proj[1]*proj[1]/edge_slen;
-	}
-
-	if(nearest_vertex)
-	{
-		VECSUB2D(v_point, nearest_vertex, point);
-		VECCOPY2D(closest, nearest_vertex);
-		return v_point[0]*v_point[0] + v_point[1]*v_point[1];
-	}
-	else
-	{
-		VECCOPY(closest, point);	//point is already inside
-		return 0.0f;
 	}
-}
-
-/*
- * Returns the square of the minimum distance between the point and a triangle surface
- * If nearest is not NULL the nearest surface point is written on it
- */
-static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest)
-{
-	//Lets solve the 2D problem (closest point-tri)
-	float normal_dist, plane_sdist, plane_offset;
-	float du[3], dv[3], dw[3];	//orthogonal axis (du=(v0->v1), dw=plane normal)
-
-	float p_2d[2], tri_2d[3][2], nearest_2d[2];
-
-	CalcNormFloat((float*)v0, (float*)v1, (float*)v2, dw);
-
-	//point-plane distance and calculate axis
-	normal_dist = point_plane_distance(point, v0, dw);
-
-	// OPTIMIZATION
-	//	if we are only interested in nearest distance if its closer than some distance already found
-	//  we can:
-	//		if(normal_dist*normal_dist >= best_dist_so_far) return FLOAT_MAX;
-	//
 
-	VECSUB(du, v1, v0);
-	Normalize(du);
-	Crossf(dv, dw, du);
-	plane_offset = INPR(v0, dw);
-
-	//project stuff to 2d
-	tri_2d[0][0] = INPR(du, v0);
-	tri_2d[0][1] = INPR(dv, v0);
-
-	tri_2d[1][0] = INPR(du, v1);
-	tri_2d[1][1] = INPR(dv, v1);
-
-	tri_2d[2][0] = INPR(du, v2);
-	tri_2d[2][1] = INPR(dv, v2);
-
-	p_2d[0] = INPR(du, point);
-	p_2d[1] = INPR(dv, point);
-
-	//we always have a right-handed tri
-	//this should always happen because of the way normal is calculated
-	plane_sdist = closest_point_in_tri2D(p_2d, tri_2d, nearest_2d);
-
-	//project back to 3d
-	if(nearest)
+	// Account for numerical round-off error
+	if ( sqrDist < FLT_EPSILON )
+		sqrDist = 0.0f;
+	
 	{
-		nearest[0] = du[0]*nearest_2d[0] + dv[0] * nearest_2d[1] + dw[0] * plane_offset;
-		nearest[1] = du[1]*nearest_2d[0] + dv[1] * nearest_2d[1] + dw[1] * plane_offset;
-		nearest[2] = du[2]*nearest_2d[0] + dv[2] * nearest_2d[1] + dw[2] * plane_offset;
+		float w[3], x[3], y[3], z[3];
+		VECCOPY(w, v0);
+		VECCOPY(x, e0);
+		VecMulf(x, S);
+		VECCOPY(y, e1);
+		VecMulf(y, T);
+		VECADD(z, w, x);
+		VECADD(z, z, y);
+		//VECSUB(d, p, z);
+		VECCOPY(nearest, z);
+		// d = p - ( v0 + S * e0 + T * e1 );
 	}
+	*v = lv;
+	*e = le;
 
-	return plane_sdist + normal_dist*normal_dist;
+	return sqrDist;
 }
 
 
@@ -268,22 +419,15 @@ static void mesh_faces_nearest_point(void *userdata, int index, const float *co,
 	do
 	{	
 		float nearest_tmp[3], dist;
-		float vec[3][3];
+		int vertex, edge;
 		
-		// only insert valid triangles / quads with area > 0
-		VECSUB(vec[0], t2, t1);
-		VECSUB(vec[1], t0, t1);
-		Crossf(vec[2], vec[0], vec[1]);
-		if(INPR(vec[2], vec[2]) >= FLT_EPSILON)
+		dist = nearest_point_in_tri_surface(t0, t1, t2, co, &vertex, &edge, nearest_tmp);
+		if(dist < nearest->dist)
 		{
-			dist = nearest_point_in_tri_surface(co,t0, t1, t2, nearest_tmp);
-			if(dist < nearest->dist)
-			{
-				nearest->index = index;
-				nearest->dist = dist;
-				VECCOPY(nearest->co, nearest_tmp);
-				CalcNormFloat((float*)t0, (float*)t1, (float*)t2, nearest->no); //TODO.. (interpolate normals from the vertexs coordinates?
-			}
+			nearest->index = index;
+			nearest->dist = dist;
+			VECCOPY(nearest->co, nearest_tmp);
+			CalcNormFloat(t0, t1, t2, nearest->no);
 		}
 
 		t1 = t2;
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
index 1c4d5e323ee..d8535e10be5 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.c
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -405,7 +405,7 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 
 
 	//Now, everything is ready to project the vertexs!
-//#pragma omp parallel for private(i,hit) schedule(static)
+#pragma omp parallel for private(i,hit) schedule(static)
 	for(i = 0; i<calc->numVerts; ++i)
 	{
 		float *co = calc->vertexCos[i];
diff --git a/source/blender/blenkernel/intern/softbody.c b/source/blender/blenkernel/intern/softbody.c
index d5b5ab6d63e..d465c058d30 100644
--- a/source/blender/blenkernel/intern/softbody.c
+++ b/source/blender/blenkernel/intern/softbody.c
@@ -69,6 +69,8 @@ variables on the UI for now
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
 #include "BLI_ghash.h"
+#include "BLI_threads.h"
+
 #include "BKE_curve.h"
 #include "BKE_effect.h"
 #include "BKE_global.h"
@@ -118,6 +120,20 @@ typedef struct SBScratch {
 	float aabbmin[3],aabbmax[3];
 }SBScratch;
 
+typedef struct  SB_thread_context{
+        Object *ob;
+		float forcetime;
+		float timenow;
+		int ifirst;
+		int ilast;
+		ListBase *do_effector;
+		int do_deflector;
+		float fieldfactor;
+		float windfactor;
+		int nr;
+		int tot;
+}SB_thread_context;
+
 #define NLF_BUILD  1 
 #define NLF_SOLVE  2 
 
@@ -1514,17 +1530,15 @@ int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp
 
 
 
-void scan_for_ext_spring_forces(Object *ob,float timenow)
+void _scan_for_ext_spring_forces(Object *ob,float timenow,int ifirst,int ilast, struct ListBase *do_effector)
 {
 	SoftBody *sb = ob->soft;
-	ListBase *do_effector;
 	int a;
 	float damp; 
 	float feedback[3];
-	do_effector= pdInitEffectors(ob,NULL);
 
 	if (sb && sb->totspring){
-		for(a=0; a<sb->totspring; a++) {
+		for(a=ifirst; a<ilast; a++) {
 			BodySpring *bs = &sb->bspring[a];
 			bs->ext_force[0]=bs->ext_force[1]=bs->ext_force[2]=0.0f; 
 			feedback[0]=feedback[1]=feedback[2]=0.0f;
@@ -1584,9 +1598,88 @@ void scan_for_ext_spring_forces(Object *ob,float timenow)
 			}
 		}
 	}
-	if(do_effector)
-		pdEndEffectors(do_effector);
 }
+
+
+void scan_for_ext_spring_forces(Object *ob,float timenow)
+{
+  SoftBody *sb = ob->soft;
+  ListBase *do_effector= NULL; 
+  do_effector= pdInitEffectors(ob,NULL);
+  if (sb){
+  _scan_for_ext_spring_forces(ob,timenow,0,sb->totspring,do_effector);
+  }
+  if(do_effector)
+  pdEndEffectors(do_effector);
+}
+
+void *exec_scan_for_ext_spring_forces(void *data)
+{
+	SB_thread_context *pctx = (SB_thread_context*)data;
+	_scan_for_ext_spring_forces(pctx->ob,pctx->timenow,pctx->ifirst,pctx->ilast,pctx->do_effector);
+	return 0;
+} 
+
+void sb_sfesf_threads_run(struct Object *ob, float timenow,int totsprings,int *ptr_to_break_func())
+{
+    ListBase *do_effector = NULL; 
+	ListBase threads;
+	SB_thread_context *sb_threads;
+	int i, totthread,left,dec;
+	int lowsprings =10; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */
+
+	do_effector= pdInitEffectors(ob,NULL);
+
+	/* figure the number of threads while preventing pretty pointless threading overhead */
+	if(totsprings < lowsprings) {totthread=1;}
+	else{
+		if(G.scene->r.mode & R_FIXED_THREADS)
+			totthread= G.scene->r.threads;
+		else
+			totthread= BLI_system_thread_count();
+	}
+	/*left to do--> what if we got zillions of CPUs running but 'totsprings' tasks to spread*/
+
+	sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBSpringsThread");
+	memset(sb_threads, 0, sizeof(SB_thread_context)*totthread);
+	left = totsprings;
+	dec = totsprings/totthread +1;
+	for(i=0; i<totthread; i++) {
+		sb_threads[i].ob = ob; 
+		sb_threads[i].forcetime = 0.0; // not used here 
+		sb_threads[i].timenow = timenow; 
+		sb_threads[i].ilast   = left; 
+		left = left - dec;
+		if (left >0){
+			sb_threads[i].ifirst  = left;
+		}
+		else
+			sb_threads[i].ifirst  = 0; 
+        sb_threads[i].do_effector = do_effector;
+        sb_threads[i].do_deflector = 0;// not used here
+		sb_threads[i].fieldfactor = 0.0f;// not used here
+		sb_threads[i].windfactor  = 0.0f;// not used here
+		sb_threads[i].nr= i;
+		sb_threads[i].tot= totthread;
+	}
+	if(totthread > 1) {
+		BLI_init_threads(&threads, exec_scan_for_ext_spring_forces, totthread);
+
+		for(i=0; i<totthread; i++)
+			BLI_insert_thread(&threads, &sb_threads[i]);
+
+		BLI_end_threads(&threads);
+	}
+	else
+		exec_scan_for_ext_spring_forces(&sb_threads[0]);
+    /* clean up */
+	MEM_freeN(sb_threads);
+
+	  if(do_effector)
+  pdEndEffectors(do_effector);
+}
+
+
 /* --- the spring external section*/
 
 int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb,float*cc)
@@ -2023,109 +2116,72 @@ static void sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float fo
 }
 
 
-static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags)
+/* since this is definitely the most CPU consuming task here .. try to spread it */
+/* core function _softbody_calc_forces_slice_in_a_thread */
+/* result is int to be able to flag user break */
+int _softbody_calc_forces_slice_in_a_thread(Object *ob, float forcetime, float timenow,int ifirst,int ilast,int *ptr_to_break_func(),ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor)
 {
-/* rule we never alter free variables :bp->vec bp->pos in here ! 
- * this will ruin adaptive stepsize AKA heun! (BM) 
- */
+	float iks;
+	int bb,do_selfcollision,do_springcollision,do_aero;
+	int number_of_points_here = ilast - ifirst;
 	SoftBody *sb= ob->soft;	/* is supposed to be there */
 	BodyPoint  *bp;
-	BodyPoint *bproot;
-	BodySpring *bs;	
-	ListBase *do_effector;
-	float iks, ks, kd, gravity;
-	float fieldfactor = 1000.0f, windfactor  = 250.0f;   
-	float tune = sb->ballstiff;
-	int a, b,  do_deflector,do_selfcollision,do_springcollision,do_aero;
-
-
-/* jacobian
-	NLboolean success;
-
-	if(nl_flags){
-		nlBegin(NL_SYSTEM);
-		nlBegin(NL_MATRIX);
-	}
-*/
-	
-	
-	gravity = sb->grav * sb_grav_force_scale(ob);	
-	
+	/* intitialize */
+	if (sb) {
 	/* check conditions for various options */
-	do_deflector= query_external_colliders(ob);
+    /* +++ could be done on object level to squeeze out the last bits of it */
 	do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
 	do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
 	do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
-	
-	iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
-	bproot= sb->bpoint; /* need this for proper spring addressing */
-	
-	if (do_springcollision || do_aero)  scan_for_ext_spring_forces(ob,timenow);
-	/* after spring scan because it uses Effoctors too */
-	do_effector= pdInitEffectors(ob,NULL);
+    /* --- could be done on object level to squeeze out the last bits of it */
+	}
+	else {
+		printf("Error expected a SB here \n");
+		return (999);
+	}
 
-	if (do_deflector) {
-		float defforce[3];
-		do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
+/* debugerin */
+	if  (sb->totpoint < ifirst) {
+		printf("Aye 998");
+		return (998);
 	}
+/* debugerin */
 
-	for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+
+	bp = &sb->bpoint[ifirst]; 
+	for(bb=number_of_points_here; bb>0; bb--, bp++) {
 		/* clear forces  accumulator */
 		bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
-		if(nl_flags & NLF_BUILD){
-			//int ia =3*(sb->totpoint-a);
-			//int op =3*sb->totpoint;
-			/* dF/dV = v */ 
-			/* jacobioan
-			nlMatrixAdd(op+ia,ia,-forcetime);
-			nlMatrixAdd(op+ia+1,ia+1,-forcetime);
-			nlMatrixAdd(op+ia+2,ia+2,-forcetime);
-     		
-			nlMatrixAdd(ia,ia,1);
-			nlMatrixAdd(ia+1,ia+1,1);
-			nlMatrixAdd(ia+2,ia+2,1);
-
-			nlMatrixAdd(op+ia,op+ia,1);
-			nlMatrixAdd(op+ia+1,op+ia+1,1);
-			nlMatrixAdd(op+ia+2,op+ia+2,1);
-			*/
-
-
-		}
-
 		/* naive ball self collision */
 		/* needs to be done if goal snaps or not */
 		if(do_selfcollision){
 			 	int attached;
 				BodyPoint   *obp;
+				BodySpring *bs;	
 				int c,b;
 				float velcenter[3],dvel[3],def[3];
 				float distance;
 				float compare;
+     	        float bstune = sb->ballstiff;
 
-				for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) {
-					
-					//if ((bp->octantflag & obp->octantflag) == 0) continue;
-
+				for(c=sb->totpoint, obp= sb->bpoint; c>=ifirst+bb; c--, obp++) {
 					compare = (obp->colball + bp->colball);		
 					VecSubf(def, bp->pos, obp->pos);
-
 					/* rather check the AABBoxes before ever calulating the real distance */
 					/* mathematically it is completly nuts, but performace is pretty much (3) times faster */
 					if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
-
                     distance = Normalize(def);
 					if (distance < compare ){
 						/* exclude body points attached with a spring */
 						attached = 0;
 						for(b=obp->nofsprings;b>0;b--){
 							bs = sb->bspring + obp->springs[b-1];
-							if (( sb->totpoint-a == bs->v2)  || ( sb->totpoint-a == bs->v1)){
+							if (( ilast-bb == bs->v2)  || ( ilast-bb == bs->v1)){
 								attached=1;
 								continue;}
 						}
 						if (!attached){
-							float f = tune/(distance) + tune/(compare*compare)*distance - 2.0f*tune/compare ;
+							float f = bstune/(distance) + bstune/(compare*compare)*distance - 2.0f*bstune/compare ;
 
 							VecMidf(velcenter, bp->vec, obp->vec);
 							VecSubf(dvel,velcenter,bp->vec);
@@ -2134,38 +2190,12 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 							Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
 							Vec3PlusStVec(bp->force,sb->balldamp,dvel);
 
-							if(nl_flags & NLF_BUILD){
-								//int ia =3*(sb->totpoint-a);
-								//int ic =3*(sb->totpoint-c);
-								//int op =3*sb->totpoint;
-								//float mvel = forcetime*sb->nodemass*sb->balldamp;
-								//float mpos = forcetime*tune*(1.0f-sb->balldamp);
-								/*some quick and dirty entries to the jacobian*/
-								//dfdx_goal(ia,ia,op,mpos);
-								//dfdv_goal(ia,ia,mvel);
-								/* exploit force(a,b) == -force(b,a) part1/2 */
-								//dfdx_goal(ic,ic,op,mpos);
-								//dfdv_goal(ic,ic,mvel);
-								
-								
-								/*TODO sit down an X-out the true jacobian entries*/
-								/*well does not make to much sense because the eigenvalues
-								of the jacobian go negative; and negative eigenvalues
-								on a complex iterative system z(n+1)=A * z(n) 
-								give imaginary roots in the charcateristic polynom
-								--> solutions that to z(t)=u(t)* exp ( i omega t) --> oscilations we don't want here 
-								where u(t) is a unknown amplitude function (worst case rising fast)
-								*/ 
-							}
-
 							/* exploit force(a,b) == -force(b,a) part2/2 */
 							VecSubf(dvel,velcenter,obp->vec);
 							VecMulf(dvel,sb->nodemass);
 
 							Vec3PlusStVec(obp->force,sb->balldamp,dvel);
 							Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def);
-
-
 						}
 					}
 				}
@@ -2179,20 +2209,13 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 			/* do goal stuff */
 			if(ob->softflag & OB_SB_GOAL) {
 				/* true elastic goal */
+				float ks,kd;
 				VecSubf(auxvect,bp->pos,bp->origT);
 				ks  = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
 				bp->force[0]+= -ks*(auxvect[0]);
 				bp->force[1]+= -ks*(auxvect[1]);
 				bp->force[2]+= -ks*(auxvect[2]);
 
-				if(nl_flags & NLF_BUILD){
-					//int ia =3*(sb->totpoint-a);
-					//int op =3*(sb->totpoint);
-					/* depending on my pos */ 
-					//dfdx_goal(ia,ia,op,ks*forcetime);
-				}
-
-
 				/* calulate damping forces generated by goals*/
 				VecSubf(velgoal,bp->origS, bp->origE);
 				kd =  sb->goalfrict * sb_fric_force_scale(ob) ;
@@ -2202,13 +2225,6 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 					bp->force[0]-= kd * (auxvect[0]);
 					bp->force[1]-= kd * (auxvect[1]);
 					bp->force[2]-= kd * (auxvect[2]);
-					if(nl_flags & NLF_BUILD){
-						//int ia =3*(sb->totpoint-a);
-						Normalize(auxvect);
-						/* depending on my vel */ 
-						//dfdv_goal(ia,ia,kd*forcetime);
-					}
-
 				}
 				else {
 					bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
@@ -2218,14 +2234,15 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 			}
 			/* done goal stuff */
 			
-			
 			/* gravitation */
+			if (sb){ 
+			float gravity = sb->grav * sb_grav_force_scale(ob);	
 			bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
-			//bp->force[1]-= gravity*sb->nodemass; /* individual mass of node here */
-
+			}
 			
 			/* particle field & vortex */
 			if(do_effector) {
+				float kd;
 				float force[3]= {0.0f, 0.0f, 0.0f};
 				float speed[3]= {0.0f, 0.0f, 0.0f};
 				float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */
@@ -2246,21 +2263,12 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 			}
 			else {
 				/* BP friction in media (not) moving*/
-				kd= sb->mediafrict* sb_fric_force_scale(ob);  
+				float kd = sb->mediafrict* sb_fric_force_scale(ob);  
 				/* assume it to be proportional to actual velocity */
 				bp->force[0]-= bp->vec[0]*kd;
 				bp->force[1]-= bp->vec[1]*kd;
 				bp->force[2]-= bp->vec[2]*kd;
 				/* friction in media done */
-				if(nl_flags & NLF_BUILD){
-					//int ia =3*(sb->totpoint-a);
-					/* da/dv =  */ 
-
-//					nlMatrixAdd(ia,ia,forcetime*kd);
-//					nlMatrixAdd(ia+1,ia+1,forcetime*kd);
-//					nlMatrixAdd(ia+2,ia+2,forcetime*kd);
-				}
-
 			}
 			/* +++cached collision targets */
 			bp->choke = 0.0f;
@@ -2268,44 +2276,25 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 			bp->flag &= ~SBF_DOFUZZY;
 			if(do_deflector) {
 				float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion;
-				kd = 1.0f;
+				float kd = 1.0f;
 
 				if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
-					if ((!nl_flags)&&(intrusion < 0.0f)){
-						/*bjornmose:  uugh.. what an evil hack 
-						violation of the 'don't touch bp->pos in here' rule 
-						but works nice, like this-->
-						we predict the solution beeing out of the collider
-						in heun step No1 and leave the heun step No2 adapt to it
-						so we kind of introduced a implicit solver for this case 
-						*/
-						Vec3PlusStVec(bp->pos,-intrusion,facenormal);
-
-						sb->scratch->flag |= SBF_DOFUZZY;
-						bp->flag |= SBF_DOFUZZY;
-					    bp->choke = sb->choke*0.01f;
-					}
-					else{
+
 							VECSUB(cfforce,bp->vec,vel);
 							Vec3PlusStVec(bp->force,-cf*50.0f,cfforce);
-					}
-					Vec3PlusStVec(bp->force,kd,defforce);
-					if (nl_flags & NLF_BUILD){
-						// int ia =3*(sb->totpoint-a);
-						// int op =3*sb->totpoint;
-						//dfdx_goal(ia,ia,op,mpos); // don't do unless you know
-						//dfdv_goal(ia,ia,-cf);
-
-					}
-  
+					
+					Vec3PlusStVec(bp->force,kd,defforce);  
 				}
 
 			}
 			/* ---cached collision targets */
 
 			/* +++springs */
+			iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
 			if(ob->softflag & OB_SB_EDGES) {
 				if (sb->bspring){ /* spring list exists at all ? */
+					int b;
+					BodySpring *bs;	
 					for(b=bp->nofsprings;b>0;b--){
 						bs = sb->bspring + bp->springs[b-1];
 						if (do_springcollision || do_aero){
@@ -2315,90 +2304,513 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
 
 						}
                         // sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
-						sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,nl_flags);
+						sb_spring_force(ob,ilast-bb,bs,iks,forcetime,0);
 					}/* loop springs */
 				}/* existing spring list */ 
 			}/*any edges*/
 			/* ---springs */
 		}/*omit on snap	*/
 	}/*loop all bp's*/
+return 0; /*done fine*/
+}
+
+void *exec_softbody_calc_forces(void *data)
+{
+	SB_thread_context *pctx = (SB_thread_context*)data;
+    _softbody_calc_forces_slice_in_a_thread(pctx->ob,pctx->forcetime,pctx->timenow,pctx->ifirst,pctx->ilast,NULL,pctx->do_effector,pctx->do_deflector,pctx->fieldfactor,pctx->windfactor);
+	return 0;
+} 
+
+void sb_cf_threads_run(struct Object *ob, float forcetime, float timenow,int totpoint,int *ptr_to_break_func(),struct ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor)
+{
+	ListBase threads;
+	SB_thread_context *sb_threads;
+	int i, totthread,left,dec;
+	int lowpoints =10; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */
+
+	/* figure the number of threads while preventing pretty pointless threading overhead */
+	if(totpoint < lowpoints) {totthread=1;}
+	else{
+		if(G.scene->r.mode & R_FIXED_THREADS)
+			totthread= G.scene->r.threads;
+		else
+			totthread= BLI_system_thread_count();
+	}
+	/*left to do--> what if we got zillions of CPUs running but 'totpoint' tasks to spread*/
+
+	sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBThread");
+	memset(sb_threads, 0, sizeof(SB_thread_context)*totthread);
+	left = totpoint;
+	dec = totpoint/totthread +1;
+	for(i=0; i<totthread; i++) {
+		sb_threads[i].ob = ob; 
+		sb_threads[i].forcetime = forcetime; 
+		sb_threads[i].timenow = timenow; 
+		sb_threads[i].ilast   = left; 
+		left = left - dec;
+		if (left >0){
+			sb_threads[i].ifirst  = left;
+		}
+		else
+			sb_threads[i].ifirst  = 0; 
+        sb_threads[i].do_effector = do_effector;
+        sb_threads[i].do_deflector = do_deflector;
+		sb_threads[i].fieldfactor = fieldfactor;
+		sb_threads[i].windfactor  = windfactor;
+		sb_threads[i].nr= i;
+		sb_threads[i].tot= totthread;
+	}
+
+
+	if(totthread > 1) {
+		BLI_init_threads(&threads, exec_softbody_calc_forces, totthread);
+
+		for(i=0; i<totthread; i++)
+			BLI_insert_thread(&threads, &sb_threads[i]);
+
+		BLI_end_threads(&threads);
+	}
+	else
+		exec_softbody_calc_forces(&sb_threads[0]);
+    /* clean up */
+	MEM_freeN(sb_threads);
+}
+
+static void softbody_calc_forcesEx(Object *ob, float forcetime, float timenow, int nl_flags)
+{
+/* rule we never alter free variables :bp->vec bp->pos in here ! 
+ * this will ruin adaptive stepsize AKA heun! (BM) 
+ */
+	SoftBody *sb= ob->soft;	/* is supposed to be there */
+	BodyPoint *bproot;
+	ListBase *do_effector;
+	float iks, gravity;
+	float fieldfactor = 1000.0f, windfactor  = 250.0f;   
+	int   do_deflector,do_selfcollision,do_springcollision,do_aero;
+	
+	gravity = sb->grav * sb_grav_force_scale(ob);	
+	
+	/* check conditions for various options */
+	do_deflector= query_external_colliders(ob);
+	do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
+	do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
+	do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
+	
+	iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
+	bproot= sb->bpoint; /* need this for proper spring addressing */
+	
+	if (do_springcollision || do_aero)  
+	sb_sfesf_threads_run(ob,timenow,sb->totspring,NULL);	
+	
+	/* after spring scan because it uses Effoctors too */
+	do_effector= pdInitEffectors(ob,NULL);
 
+	if (do_deflector) {
+		float defforce[3];
+		do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
+	}
+
+	sb_cf_threads_run(ob,forcetime,timenow,sb->totpoint,NULL,do_effector,do_deflector,fieldfactor,windfactor);
 
 	/* finally add forces caused by face collision */
 	if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
 	
 	/* finish matrix and solve */
-#if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be usefull .. so keep that BM
-	if(nl_flags & NLF_SOLVE){
-		//double sct,sst=PIL_check_seconds_timer();
-		for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
-			int iv =3*(sb->totpoint-a);
-			int ip =3*(2*sb->totpoint-a);
-			int n;
-			for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);}
-			for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);}
+	if(do_effector) pdEndEffectors(do_effector);
+}
+
+
+
+
+static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags)
+{
+	/* redirection to the new threaded Version */
+	if (G.rt !=16){ 
+		softbody_calc_forcesEx(ob, forcetime, timenow, nl_flags);
+		return;
+	}
+	else{
+		/* so the following will die  */
+		/* |||||||||||||||||||||||||| */
+		/* VVVVVVVVVVVVVVVVVVVVVVVVVV */
+
+		/* rule we never alter free variables :bp->vec bp->pos in here ! 
+		* this will ruin adaptive stepsize AKA heun! (BM) 
+		*/
+		SoftBody *sb= ob->soft;	/* is supposed to be there */
+		BodyPoint  *bp;
+		BodyPoint *bproot;
+		BodySpring *bs;	
+		ListBase *do_effector;
+		float iks, ks, kd, gravity;
+		float fieldfactor = 1000.0f, windfactor  = 250.0f;   
+		float tune = sb->ballstiff;
+		int a, b,  do_deflector,do_selfcollision,do_springcollision,do_aero;
+
+
+		/* jacobian
+		NLboolean success;
+
+		if(nl_flags){
+		nlBegin(NL_SYSTEM);
+		nlBegin(NL_MATRIX);
 		}
-		nlEnd(NL_MATRIX);
-		nlEnd(NL_SYSTEM);
+		*/
+
+
+		gravity = sb->grav * sb_grav_force_scale(ob);	
+
+		/* check conditions for various options */
+		do_deflector= query_external_colliders(ob);
+		do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
+		do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
+		do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
+
+		iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
+		bproot= sb->bpoint; /* need this for proper spring addressing */
 
-		if ((G.rt >0) && (nl_flags & NLF_BUILD))
-		{ 
-			printf("####MEE#####\n");
-			nlPrintMatrix();
+		if (do_springcollision || do_aero)  scan_for_ext_spring_forces(ob,timenow);
+		/* after spring scan because it uses Effoctors too */
+		do_effector= pdInitEffectors(ob,NULL);
+
+		if (do_deflector) {
+			float defforce[3];
+			do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
 		}
 
-		success= nlSolveAdvanced(NULL, 1);
+		for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+			/* clear forces  accumulator */
+			bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
+			if(nl_flags & NLF_BUILD){
+				//int ia =3*(sb->totpoint-a);
+				//int op =3*sb->totpoint;
+				/* dF/dV = v */ 
+				/* jacobioan
+				nlMatrixAdd(op+ia,ia,-forcetime);
+				nlMatrixAdd(op+ia+1,ia+1,-forcetime);
+				nlMatrixAdd(op+ia+2,ia+2,-forcetime);
+
+				nlMatrixAdd(ia,ia,1);
+				nlMatrixAdd(ia+1,ia+1,1);
+				nlMatrixAdd(ia+2,ia+2,1);
+
+				nlMatrixAdd(op+ia,op+ia,1);
+				nlMatrixAdd(op+ia+1,op+ia+1,1);
+				nlMatrixAdd(op+ia+2,op+ia+2,1);
+				*/
 
-		// nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */
-		if(success){
-			float f;
-			int index =0;
-			/* for debug purpose .. anyhow cropping B vector looks like working */
-			if (G.rt >0)
-				for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
-					f=nlGetVariable(0,index);
-					printf("(%f ",f);index++;
-					f=nlGetVariable(0,index);
-					printf("%f ",f);index++;
-					f=nlGetVariable(0,index);
-					printf("%f)",f);index++;
+
+			}
+
+			/* naive ball self collision */
+			/* needs to be done if goal snaps or not */
+			if(do_selfcollision){
+				int attached;
+				BodyPoint   *obp;
+				int c,b;
+				float velcenter[3],dvel[3],def[3];
+				float distance;
+				float compare;
+
+				for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) {
+
+					//if ((bp->octantflag & obp->octantflag) == 0) continue;
+
+					compare = (obp->colball + bp->colball);		
+					VecSubf(def, bp->pos, obp->pos);
+
+					/* rather check the AABBoxes before ever calulating the real distance */
+					/* mathematically it is completly nuts, but performace is pretty much (3) times faster */
+					if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
+
+					distance = Normalize(def);
+					if (distance < compare ){
+						/* exclude body points attached with a spring */
+						attached = 0;
+						for(b=obp->nofsprings;b>0;b--){
+							bs = sb->bspring + obp->springs[b-1];
+							if (( sb->totpoint-a == bs->v2)  || ( sb->totpoint-a == bs->v1)){
+								attached=1;
+								continue;}
+						}
+						if (!attached){
+							float f = tune/(distance) + tune/(compare*compare)*distance - 2.0f*tune/compare ;
+
+							VecMidf(velcenter, bp->vec, obp->vec);
+							VecSubf(dvel,velcenter,bp->vec);
+							VecMulf(dvel,sb->nodemass);
+
+							Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
+							Vec3PlusStVec(bp->force,sb->balldamp,dvel);
+
+							if(nl_flags & NLF_BUILD){
+								//int ia =3*(sb->totpoint-a);
+								//int ic =3*(sb->totpoint-c);
+								//int op =3*sb->totpoint;
+								//float mvel = forcetime*sb->nodemass*sb->balldamp;
+								//float mpos = forcetime*tune*(1.0f-sb->balldamp);
+								/*some quick and dirty entries to the jacobian*/
+								//dfdx_goal(ia,ia,op,mpos);
+								//dfdv_goal(ia,ia,mvel);
+								/* exploit force(a,b) == -force(b,a) part1/2 */
+								//dfdx_goal(ic,ic,op,mpos);
+								//dfdv_goal(ic,ic,mvel);
+
+
+								/*TODO sit down an X-out the true jacobian entries*/
+								/*well does not make to much sense because the eigenvalues
+								of the jacobian go negative; and negative eigenvalues
+								on a complex iterative system z(n+1)=A * z(n) 
+								give imaginary roots in the charcateristic polynom
+								--> solutions that to z(t)=u(t)* exp ( i omega t) --> oscilations we don't want here 
+								where u(t) is a unknown amplitude function (worst case rising fast)
+								*/ 
+							}
+
+							/* exploit force(a,b) == -force(b,a) part2/2 */
+							VecSubf(dvel,velcenter,obp->vec);
+							VecMulf(dvel,sb->nodemass);
+
+							Vec3PlusStVec(obp->force,sb->balldamp,dvel);
+							Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def);
+
+
+						}
+					}
 				}
+			}
+			/* naive ball self collision done */
 
-				index =0;
-				for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+			if(bp->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */
+				float auxvect[3];  
+				float velgoal[3];
+
+				/* do goal stuff */
+				if(ob->softflag & OB_SB_GOAL) {
+					/* true elastic goal */
+					VecSubf(auxvect,bp->pos,bp->origT);
+					ks  = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
+					bp->force[0]+= -ks*(auxvect[0]);
+					bp->force[1]+= -ks*(auxvect[1]);
+					bp->force[2]+= -ks*(auxvect[2]);
+
+					if(nl_flags & NLF_BUILD){
+						//int ia =3*(sb->totpoint-a);
+						//int op =3*(sb->totpoint);
+						/* depending on my pos */ 
+						//dfdx_goal(ia,ia,op,ks*forcetime);
+					}
+
+
+					/* calulate damping forces generated by goals*/
+					VecSubf(velgoal,bp->origS, bp->origE);
+					kd =  sb->goalfrict * sb_fric_force_scale(ob) ;
+					VecAddf(auxvect,velgoal,bp->vec);
+
+					if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */
+						bp->force[0]-= kd * (auxvect[0]);
+						bp->force[1]-= kd * (auxvect[1]);
+						bp->force[2]-= kd * (auxvect[2]);
+						if(nl_flags & NLF_BUILD){
+							//int ia =3*(sb->totpoint-a);
+							Normalize(auxvect);
+							/* depending on my vel */ 
+							//dfdv_goal(ia,ia,kd*forcetime);
+						}
+
+					}
+					else {
+						bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
+						bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
+						bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
+					}
+				}
+				/* done goal stuff */
+
+
+				/* gravitation */
+				bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
+				//bp->force[1]-= gravity*sb->nodemass; /* individual mass of node here */
+
+
+				/* particle field & vortex */
+				if(do_effector) {
+					float force[3]= {0.0f, 0.0f, 0.0f};
+					float speed[3]= {0.0f, 0.0f, 0.0f};
+					float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */
+
+					pdDoEffectors(do_effector, bp->pos, force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);
+
+					/* apply forcefield*/
+					VecMulf(force,fieldfactor* eval_sb_fric_force_scale); 
+					VECADD(bp->force, bp->force, force);
+
+					/* BP friction in moving media */	
+					kd= sb->mediafrict* eval_sb_fric_force_scale;  
+					bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
+					bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
+					bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
+					/* now we'll have nice centrifugal effect for vortex */
+
+				}
+				else {
+					/* BP friction in media (not) moving*/
+					kd= sb->mediafrict* sb_fric_force_scale(ob);  
+					/* assume it to be proportional to actual velocity */
+					bp->force[0]-= bp->vec[0]*kd;
+					bp->force[1]-= bp->vec[1]*kd;
+					bp->force[2]-= bp->vec[2]*kd;
+					/* friction in media done */
+					if(nl_flags & NLF_BUILD){
+						//int ia =3*(sb->totpoint-a);
+						/* da/dv =  */ 
+
+						//					nlMatrixAdd(ia,ia,forcetime*kd);
+						//					nlMatrixAdd(ia+1,ia+1,forcetime*kd);
+						//					nlMatrixAdd(ia+2,ia+2,forcetime*kd);
+					}
+
+				}
+				/* +++cached collision targets */
+				bp->choke = 0.0f;
+				bp->choke2 = 0.0f;
+				bp->flag &= ~SBF_DOFUZZY;
+				if(do_deflector) {
+					float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion;
+					kd = 1.0f;
+
+					if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
+						if ((!nl_flags)&&(intrusion < 0.0f)){
+							/*bjornmose:  uugh.. what an evil hack 
+							violation of the 'don't touch bp->pos in here' rule 
+							but works nice, like this-->
+							we predict the solution beeing out of the collider
+							in heun step No1 and leave the heun step No2 adapt to it
+							so we kind of introduced a implicit solver for this case 
+							*/
+							Vec3PlusStVec(bp->pos,-intrusion,facenormal);
+
+							sb->scratch->flag |= SBF_DOFUZZY;
+							bp->flag |= SBF_DOFUZZY;
+							bp->choke = sb->choke*0.01f;
+						}
+						else{
+							VECSUB(cfforce,bp->vec,vel);
+							Vec3PlusStVec(bp->force,-cf*50.0f,cfforce);
+						}
+						Vec3PlusStVec(bp->force,kd,defforce);
+						if (nl_flags & NLF_BUILD){
+							// int ia =3*(sb->totpoint-a);
+							// int op =3*sb->totpoint;
+							//dfdx_goal(ia,ia,op,mpos); // don't do unless you know
+							//dfdv_goal(ia,ia,-cf);
+
+						}
+
+					}
+
+				}
+				/* ---cached collision targets */
+
+				/* +++springs */
+				if(ob->softflag & OB_SB_EDGES) {
+					if (sb->bspring){ /* spring list exists at all ? */
+						for(b=bp->nofsprings;b>0;b--){
+							bs = sb->bspring + bp->springs[b-1];
+							if (do_springcollision || do_aero){
+								VecAddf(bp->force,bp->force,bs->ext_force);
+								if (bs->flag & BSF_INTERSECT)
+									bp->choke = bs->cf; 
+
+							}
+							// sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
+							// rather remove nl_falgs from code .. will make things a lot cleaner
+							sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,0);
+						}/* loop springs */
+					}/* existing spring list */ 
+				}/*any edges*/
+				/* ---springs */
+			}/*omit on snap	*/
+		}/*loop all bp's*/
+
+
+		/* finally add forces caused by face collision */
+		if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
+
+		/* finish matrix and solve */
+#if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be usefull .. so keep that BM
+		if(nl_flags & NLF_SOLVE){
+			//double sct,sst=PIL_check_seconds_timer();
+			for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+				int iv =3*(sb->totpoint-a);
+				int ip =3*(2*sb->totpoint-a);
+				int n;
+				for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);}
+				for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);}
+			}
+			nlEnd(NL_MATRIX);
+			nlEnd(NL_SYSTEM);
+
+			if ((G.rt == 32) && (nl_flags & NLF_BUILD))
+			{ 
+				printf("####MEE#####\n");
+				nlPrintMatrix();
+			}
+
+			success= nlSolveAdvanced(NULL, 1);
+
+			// nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */
+			if(success){
+				float f;
+				int index =0;
+				/* for debug purpose .. anyhow cropping B vector looks like working */
+				if (G.rt ==32)
+					for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+						f=nlGetVariable(0,index);
+						printf("(%f ",f);index++;
+						f=nlGetVariable(0,index);
+						printf("%f ",f);index++;
+						f=nlGetVariable(0,index);
+						printf("%f)",f);index++;
+					}
+
+					index =0;
+					for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
+						f=nlGetVariable(0,index);
+						bp->impdv[0] = f; index++;
+						f=nlGetVariable(0,index);
+						bp->impdv[1] = f; index++;
+						f=nlGetVariable(0,index);
+						bp->impdv[2] = f; index++;
+					}
+					/*
+					for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 					f=nlGetVariable(0,index);
-					bp->impdv[0] = f; index++;
+					bp->impdx[0] = f; index++;
 					f=nlGetVariable(0,index);
-					bp->impdv[1] = f; index++;
+					bp->impdx[1] = f; index++;
 					f=nlGetVariable(0,index);
-					bp->impdv[2] = f; index++;
-				}
-				/*
+					bp->impdx[2] = f; index++;
+					}
+					*/
+			}
+			else{
+				printf("Matrix inversion failed \n");
 				for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
-				f=nlGetVariable(0,index);
-				bp->impdx[0] = f; index++;
-				f=nlGetVariable(0,index);
-				bp->impdx[1] = f; index++;
-				f=nlGetVariable(0,index);
-				bp->impdx[2] = f; index++;
+					VECCOPY(bp->impdv,bp->force);
 				}
-				*/
-		}
-		else{
-			printf("Matrix inversion failed \n");
-			for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
-				VECCOPY(bp->impdv,bp->force);
+
 			}
 
+			//sct=PIL_check_seconds_timer();
+			//if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name);
 		}
-
-		//sct=PIL_check_seconds_timer();
-		//if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name);
-	}
-	/* cleanup */
+		/* cleanup */
 #endif
-	if(do_effector) pdEndEffectors(do_effector);
+		if(do_effector) pdEndEffectors(do_effector);
+	}
 }
+
 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags)
 {
 	/* time evolution */
@@ -2458,7 +2870,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *
 			/* x(t + dt) = x(t) + v(t~) * dt */ 
 			VecMulf(dx,forcetime);
 
-			/* the freezer */
+			/* the freezer coming sooner or later */
 			/*
 			if  ((Inpf(dx,dx)<freezeloc )&&(Inpf(bp->force,bp->force)<freezeforce )){
 				bp->frozen /=2;
@@ -3529,6 +3941,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
 							 * we don't want to lock up the system if physics fail
 		*/
 		int loops =0 ; 
+		
 		SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */
 		if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
 		
@@ -3546,13 +3959,13 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
 			sb->scratch->flag &= ~SBF_DOFUZZY;
 			/* do predictive euler step */
 			softbody_calc_forces(ob, forcetime,timedone/dtime,0);
-			softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags);
 
+			softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags);
 
 			/* crop new slope values to do averaged slope step */
 			softbody_calc_forces(ob, forcetime,timedone/dtime,0);
-			softbody_apply_forces(ob, forcetime, 2, &err,mid_flags);
 
+			softbody_apply_forces(ob, forcetime, 2, &err,mid_flags);
 			softbody_apply_goalsnap(ob);
 			
 			if (err > SoftHeunTol) { /* error needs to be scaled to some quantity */
@@ -3603,7 +4016,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
 		//				if(G.f & G_DEBUG){
 		if(sb->solverflags & SBSO_MONITOR ){
 			if (loops > HEUNWARNLIMIT) /* monitor high loop counts */
-				printf("\r       needed %d steps/frame ",loops);
+				printf("\r needed %d steps/frame",loops);
 		}
 		
 	}
@@ -3627,7 +4040,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
 
 	if(sb->solverflags & SBSO_MONITOR ){
 		sct=PIL_check_seconds_timer();
-		if (sct-sst > 0.5f) printf(" solver time %f %s \r",sct-sst,ob->id.name);
+		if (sct-sst > 0.5f) printf(" solver time %f sec %s \n",sct-sst,ob->id.name);
 	}
 }
 
diff --git a/source/blender/python/api2_2x/Texture.c b/source/blender/python/api2_2x/Texture.c
index 696c78f2bac..7ba8ad88ea6 100644
--- a/source/blender/python/api2_2x/Texture.c
+++ b/source/blender/python/api2_2x/Texture.c
@@ -105,6 +105,10 @@
 #define EXPP_TEX_LACUNARITY_MAX             6.0f
 #define EXPP_TEX_OCTS_MIN                   0.0f
 #define EXPP_TEX_OCTS_MAX                   8.0f
+#define EXPP_TEX_OFST_MIN                   0.0f
+#define EXPP_TEX_OFST_MAX                   6.0f
+#define EXPP_TEX_GAIN_MIN                   0.0f
+#define EXPP_TEX_GAIN_MAX                   6.0f
 #define EXPP_TEX_ISCALE_MIN                 0.0f
 #define EXPP_TEX_ISCALE_MAX                 10.0f
 #define EXPP_TEX_EXP_MIN                    0.010f
@@ -430,6 +434,8 @@ GETFUNC( getNoiseDepth );
 GETFUNC( getNoiseSize );
 GETFUNC( getNoiseType );
 GETFUNC( getOcts );
+GETFUNC( getOffset );
+GETFUNC( getGain );
 GETFUNC( getRepeat );
 GETFUNC( getRGBCol );
 GETFUNC( getSType );
@@ -478,6 +484,8 @@ SETFUNC( setNoiseDepth );
 SETFUNC( setNoiseSize );
 SETFUNC( setNoiseType );
 SETFUNC( setOcts );
+SETFUNC( setOffset );
+SETFUNC( setGain );
 SETFUNC( setRepeat );
 SETFUNC( setRGBCol );
 SETFUNC( setSType );
@@ -646,6 +654,14 @@ static PyGetSetDef BPy_Texture_getseters[] = {
 	 (getter)Texture_getLacunarity, (setter)Texture_setLacunarity,
 	 "Gap between succesive frequencies (for Musgrave textures)",
 	 NULL},
+	{"offset",
+	 (getter)Texture_getOffset, (setter)Texture_setOffset,
+	 "Fractal offset (for Musgrave textures)",
+	 NULL},
+	{"gain",
+	 (getter)Texture_getGain, (setter)Texture_setGain,
+	 "Gain multiplier (for Musgrave textures)",
+	 NULL},
 	{"noiseBasis",
 	 (getter)Texture_getNoiseBasis, (setter)Texture_setNoiseBasis,
 	 "Noise basis type (wood, stucci, marble, clouds, Musgrave, distorted noise)",
@@ -1837,6 +1853,20 @@ static int Texture_setLacunarity( BPy_Texture * self, PyObject * value )
 								EXPP_TEX_LACUNARITY_MAX );
 }
 
+static int Texture_setOffset( BPy_Texture * self, PyObject * value )
+{
+	return EXPP_setFloatClamped ( value, &self->texture->mg_offset,
+								EXPP_TEX_OFST_MIN,
+								EXPP_TEX_OFST_MAX );
+}
+
+static int Texture_setGain( BPy_Texture * self, PyObject * value )
+{
+	return EXPP_setFloatClamped ( value, &self->texture->mg_gain,
+								EXPP_TEX_GAIN_MIN,
+								EXPP_TEX_GAIN_MAX );
+}
+
 static int Texture_setOcts( BPy_Texture * self, PyObject * value )
 {
 	return EXPP_setFloatClamped ( value, &self->texture->mg_octaves,
@@ -2168,6 +2198,16 @@ static PyObject *Texture_getOcts( BPy_Texture *self )
 	return PyFloat_FromDouble( self->texture->mg_octaves );
 }
 
+static PyObject *Texture_getOffset( BPy_Texture *self )
+{
+	return PyFloat_FromDouble( self->texture->mg_offset );
+}
+
+static PyObject *Texture_getGain( BPy_Texture *self )
+{
+	return PyFloat_FromDouble( self->texture->mg_gain );
+}
+
 static PyObject *Texture_getRepeat( BPy_Texture *self )
 {
 	return Py_BuildValue( "(i,i)", self->texture->xrepeat,
diff --git a/source/blender/python/api2_2x/doc/Texture.py b/source/blender/python/api2_2x/doc/Texture.py
index 7f42d06240b..cebb7de7011 100644
--- a/source/blender/python/api2_2x/doc/Texture.py
+++ b/source/blender/python/api2_2x/doc/Texture.py
@@ -344,6 +344,12 @@ class Texture:
 	@ivar octs:  Number of frequencies (for Musgrave textures).
 	Value is clamped to the range [0.0,8.0].
 	@type octs:  float
+	@ivar offset:  Fractal offset (for hetero terrain and multifractal Musgrave textures).
+	Value is clamped to the range [0.0,6.0].
+	@type offset:  float
+	@ivar gain:  Gain multiplier (for multifractal Musgrave textures).
+	Value is clamped to the range [0.0,6.0].
+	@type gain:  float
 	@ivar repeat:  Repetition multiplier (for image textures).
 	@type repeat:  tuple of 2 ints
 	@ivar rgbCol:  RGB color tuple.
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index 139943ae5ba..f7606884f29 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -6720,3 +6720,4 @@ void editing_panels()
 	}
 	uiClearButLock();
 }
+