=trunk=

Recommitted eltopo collision code (but disabled by default)
with Genscher's permission.

To use, you need to install liblapack and libblas

1 files changed, 865 insertions, 9 deletions

diff --git a/source/blender/blenkernel/intern/collision.c b/source/blender/blenkernel/intern/collision.c
index ebdbbfcf7b4..064f805d046 100644
--- a/source/blender/blenkernel/intern/collision.c
+++ b/source/blender/blenkernel/intern/collision.c
@@ -48,6 +48,9 @@
 #include "BLI_math.h"
 #include "BLI_edgehash.h"
 #include "BLI_utildefines.h"
+#include "BLI_ghash.h"
+#include "BLI_memarena.h"
+#include "BLI_rand.h"
 
 #include "BKE_DerivedMesh.h"
 #include "BKE_global.h"
@@ -63,6 +66,10 @@
 #include "BLI_kdopbvh.h"
 #include "BKE_collision.h"
 
+#ifdef USE_ELTOPO
+#include "eltopo-capi.h"
+#endif
+
 
 /***********************************
 Collision modifier code start
@@ -486,7 +493,7 @@ DO_INLINE void collision_interpolateOnTriangle ( float to[3], float v1[3], float
 	VECADDMUL ( to, v3, w3 );
 }
 
-
+#ifndef USE_ELTOPO
 static int cloth_collision_response_static ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
 {
 	int result = 0;
@@ -601,12 +608,799 @@ static int cloth_collision_response_static ( ClothModifierData *clmd, CollisionM
 	}
 	return result;
 }
+#endif
+
+#ifdef USE_ELTOPO
+typedef struct edgepairkey {
+	int a1, a2, b1, b2;
+} edgepairkey;
+
+unsigned int edgepair_hash(void *vkey)
+{
+	edgepairkey *key = vkey;
+	int keys[4] = {key->a1, key->a2, key->b1, key->b2};
+	int i, j;
+	
+	for (i=0; i<4; i++) {
+		for (j=0; j<3; j++) {
+			if (keys[j] >= keys[j+1]) {
+				SWAP(int, keys[j], keys[j+1]);
+			}
+		}
+	}
+	
+	return keys[0]*101 + keys[1]*72 + keys[2]*53 + keys[3]*34;
+}
+
+int edgepair_cmp(const void *va, const void *vb)
+{
+	edgepairkey *a = va, *b = vb;
+	int keysa[4] = {a->a1, a->a2, a->b1, a->b2};
+	int keysb[4] = {b->a1, b->a2, b->b1, b->b2};
+	int i;
+	
+	for (i=0; i<4; i++) {
+		int j, ok=0;
+		for (j=0; j<4; j++) {
+			if (keysa[i] == keysa[j]) {
+				ok = 1;
+				break;
+			}
+		}
+		if (!ok)
+			return -1;
+	}
+	
+	return 0;
+}
+
+static void get_edgepairkey(edgepairkey *key, int a1, int a2, int b1, int b2)
+{
+	key->a1 = a1;
+	key->a2 = a2;
+	key->b1 = b1;
+	key->b2 = b2;
+}
+
+/*an immense amount of duplication goes on here. . .a major performance hit, I'm sure*/
+static CollPair* cloth_edge_collision ( ModifierData *md1, ModifierData *md2, 
+										BVHTreeOverlap *overlap, CollPair *collpair,
+										GHash *visithash, MemArena *arena)
+{
+	ClothModifierData *clmd = ( ClothModifierData * ) md1;
+	CollisionModifierData *collmd = ( CollisionModifierData * ) md2;
+	MFace *face1=NULL, *face2 = NULL;
+	ClothVertex *verts1 = clmd->clothObject->verts;
+	double distance = 0;
+	edgepairkey *key, tstkey;
+	float epsilon1 = clmd->coll_parms->epsilon;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
+	float no[3], uv[3], t, relnor;
+	int i, i1, i2, i3, i4, i5, i6;
+	Cloth *cloth = clmd->clothObject;
+	float n1[3], n2[3], off[3], v1[2][3], v2[2][3], v3[2][3], v4[2][3], v5[2][3], v6[2][3];
+	void **verts[] = {v1, v2, v3, v4, v5, v6};
+	int j, ret, bp1, bp2, bp3, ap1, ap2, ap3, table[6];
+	
+	face1 = & ( clmd->clothObject->mfaces[overlap->indexA] );
+	face2 = & ( collmd->mfaces[overlap->indexB] );
+
+	// check all 4 possible collisions
+	for ( i = 0; i < 4; i++ )
+	{
+		if ( i == 0 )
+		{
+			// fill faceA
+			ap1 = face1->v1;
+			ap2 = face1->v2;
+			ap3 = face1->v3;
+
+			// fill faceB
+			bp1 = face2->v1;
+			bp2 = face2->v2;
+			bp3 = face2->v3;
+		}
+		else if ( i == 1 )
+		{
+			if ( face1->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v3;
+				ap3 = face1->v4;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v2;
+				bp3 = face2->v3;
+			}
+			else {
+				continue;
+			}
+		}
+		if ( i == 2 )
+		{
+			if ( face2->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v2;
+				ap3 = face1->v3;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v3;
+				bp3 = face2->v4;
+			}
+			else {
+				continue;
+			}
+		}
+		else if ( i == 3 )
+		{
+			if ( face1->v4 && face2->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v3;
+				ap3 = face1->v4;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v3;
+				bp3 = face2->v4;
+			}
+			else {
+				continue;
+			}
+		}
+		
+		copy_v3_v3(v1[0], cloth->verts[ap1].txold); 
+		copy_v3_v3(v1[1], cloth->verts[ap1].tx);
+		copy_v3_v3(v2[0], cloth->verts[ap2].txold);
+		copy_v3_v3(v2[1], cloth->verts[ap2].tx);
+		copy_v3_v3(v3[0], cloth->verts[ap3].txold);
+		copy_v3_v3(v3[1], cloth->verts[ap3].tx);
+		
+		copy_v3_v3(v4[0], collmd->current_x[bp1].co);
+		copy_v3_v3(v4[1], collmd->current_xnew[bp1].co);
+		copy_v3_v3(v5[0], collmd->current_x[bp2].co);
+		copy_v3_v3(v5[1], collmd->current_xnew[bp2].co);
+		copy_v3_v3(v6[0], collmd->current_x[bp3].co);
+		copy_v3_v3(v6[1], collmd->current_xnew[bp3].co);
+		
+		normal_tri_v3(n2, v4[1], v5[1], v6[1]);
+
+		/*offset new positions a bit, to account for margins*/
+		i1 = ap1; i2 = ap2; i3 = ap3;
+		i4 = bp1; i5 = bp2; i6 = bp3;
+
+		for (j=0; j<3; j++) {
+			int collp1, collp2, k, j2 = (j+1)%3;
+			
+			table[0] = ap1; table[1] = ap2; table[2] = ap3;
+			table[3] = bp1; table[4] = bp2; table[5] = bp3;
+			for (k=0; k<3; k++) {
+				float p1[3], p2[3];
+				int k2 = (k+1)%3;
+				
+				get_edgepairkey(&tstkey, table[j], table[j2], table[k+3], table[k2+3]);
+				//if (BLI_ghash_haskey(visithash, &tstkey))
+				//	continue;
+				
+				key = BLI_memarena_alloc(arena, sizeof(edgepairkey));
+				*key = tstkey;
+				BLI_ghash_insert(visithash, key, NULL);
+
+				sub_v3_v3v3(p1, verts[j], verts[j2]);
+				sub_v3_v3v3(p2, verts[k+3], verts[k2+3]);
+				
+				cross_v3_v3v3(off, p1, p2);
+				normalize_v3(off);
+
+				if (dot_v3v3(n2, off) < 0.0)
+					negate_v3(off);
+				
+				mul_v3_fl(off,  epsilon1 + epsilon2 + ALMOST_ZERO);
+				copy_v3_v3(p1, verts[k+3]);
+				copy_v3_v3(p2, verts[k2+3]);
+				add_v3_v3(p1, off);
+				add_v3_v3(p2, off);
+				
+				ret = eltopo_line_line_moving_isect_v3v3_f(verts[j], table[j], verts[j2], table[j2], 
+														   p1, table[k+3], p2, table[k2+3], 
+														   no, uv, &t, &relnor);
+				/*cloth vert versus coll face*/
+				if (ret) {
+					collpair->ap1 = table[j]; collpair->ap2 = table[j2]; 
+					collpair->bp1 = table[k+3]; collpair->bp2 = table[k2+3];
+					
+					/*I'm not sure if this is correct, but hopefully it's 
+					  better then simply ignoring back edges*/
+					if (dot_v3v3(n2, no) < 0.0) {
+						negate_v3(no);
+					}
+					
+					copy_v3_v3(collpair->normal, no);
+					mul_v3_v3fl(collpair->vector, collpair->normal, relnor);
+					collpair->distance = relnor;
+					collpair->time = t;
+					
+					copy_v2_v2(collpair->bary, uv);
+					
+					collpair->flag = COLLISION_IS_EDGES;
+					collpair++;
+				}
+			}
+		}
+	}
+	
+	return collpair;
+}
+
+static int cloth_edge_collision_response_moving ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
+{
+	int result = 0;
+	Cloth *cloth1;
+	float w1, w2;
+	float v1[3], v2[3], relativeVelocity[3];
+	float magrelVel, pimpulse[3];
+
+	cloth1 = clmd->clothObject;
+
+	for ( ; collpair != collision_end; collpair++ )
+	{
+		if (!(collpair->flag & COLLISION_IS_EDGES))
+			continue;
+		
+		// was: txold
+		w1 = collpair->bary[0]; w2 = collpair->bary[1];			
+		
+		// Calculate relative "velocity".
+		VECADDFAC(v1, cloth1->verts[collpair->ap1].tv, cloth1->verts[collpair->ap2].tv, w1);
+		VECADDFAC(v2, collmd->current_v[collpair->bp1].co, collmd->current_v[collpair->bp2].co, w2);
+		
+		VECSUB ( relativeVelocity, v2, v1);
+		
+		// Calculate the normal component of the relative velocity (actually only the magnitude - the direction is stored in 'normal').
+		magrelVel = INPR ( relativeVelocity, collpair->normal );
+
+		// If v_n_mag < 0 the edges are approaching each other.
+		if ( magrelVel > ALMOST_ZERO )
+		{
+			// Calculate Impulse magnitude to stop all motion in normal direction.
+			float magtangent = 0, repulse = 0, d = 0;
+			double impulse = 0.0;
+			float vrel_t_pre[3];
+			float temp[3], spf;
+			
+			zero_v3(pimpulse);
+			
+			// calculate tangential velocity
+			VECCOPY ( temp, collpair->normal );
+			mul_v3_fl( temp, magrelVel );
+			VECSUB ( vrel_t_pre, relativeVelocity, temp );
+
+			// Decrease in magnitude of relative tangential velocity due to coulomb friction
+			// in original formula "magrelVel" should be the "change of relative velocity in normal direction"
+			magtangent = MIN2 ( clmd->coll_parms->friction * 0.01 * magrelVel,sqrt ( INPR ( vrel_t_pre,vrel_t_pre ) ) );
+
+			// Apply friction impulse.
+			if ( magtangent > ALMOST_ZERO )
+			{
+				normalize_v3( vrel_t_pre );
+
+				impulse = magtangent; 
+				VECADDMUL ( pimpulse, vrel_t_pre, impulse);
+			}
+
+			// Apply velocity stopping impulse
+			// I_c = m * v_N / 2.0
+			// no 2.0 * magrelVel normally, but looks nicer DG
+			impulse =  magrelVel;
+			
+			mul_v3_fl(collpair->normal, 0.5);
+			VECADDMUL ( pimpulse, collpair->normal, impulse);
+
+			// Apply repulse impulse if distance too short
+			// I_r = -min(dt*kd, m(0,1d/dt - v_n))
+			spf = (float)clmd->sim_parms->stepsPerFrame / clmd->sim_parms->timescale;
+
+			d = collpair->distance;
+			if ( ( magrelVel < 0.1*d*spf && ( d > ALMOST_ZERO ) ) )
+			{
+				repulse = MIN2 ( d*1.0/spf, 0.1*d*spf - magrelVel );
+
+				// stay on the safe side and clamp repulse
+				if ( impulse > ALMOST_ZERO )
+					repulse = MIN2 ( repulse, 5.0*impulse );
+				repulse = MAX2 ( impulse, repulse );
+
+				impulse = repulse / ( 5.0 ); // original 2.0 / 0.25
+				VECADDMUL ( pimpulse, collpair->normal, impulse);
+			}
+			
+			w2 = 1.0f-w1;
+			if (w1 < 0.5)
+				w1 *= 2.0;
+			else
+				w2 *= 2.0;
+			
+			VECADDFAC(cloth1->verts[collpair->ap1].impulse, cloth1->verts[collpair->ap1].impulse, pimpulse, w1*2.0);
+			VECADDFAC(cloth1->verts[collpair->ap2].impulse, cloth1->verts[collpair->ap2].impulse, pimpulse, w2*2.0);
+			
+			cloth1->verts[collpair->ap1].impulse_count++;
+			cloth1->verts[collpair->ap2].impulse_count++;
+			
+			result = 1;
+		}
+	} 
+	
+	return result;
+}
+
+static int cloth_collision_response_moving ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
+{
+	int result = 0;
+	Cloth *cloth1;
+	float w1, w2, w3, u1, u2, u3;
+	float v1[3], v2[3], relativeVelocity[3];
+	float magrelVel;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
+	
+	cloth1 = clmd->clothObject;
+
+	for ( ; collpair != collision_end; collpair++ )
+	{
+		if (collpair->flag & COLLISION_IS_EDGES)
+			continue;
+		
+		if ( collpair->flag & COLLISION_USE_COLLFACE ) {
+			// was: txold
+			w1 = collpair->bary[0]; w2 = collpair->bary[1]; w3 = collpair->bary[2];			
+
+			// Calculate relative "velocity".
+			collision_interpolateOnTriangle ( v1, collmd->current_v[collpair->bp1].co, collmd->current_v[collpair->bp2].co, collmd->current_v[collpair->bp3].co, w1, w2, w3);
+			
+			VECSUB ( relativeVelocity, v1, cloth1->verts[collpair->collp].tv);
+			
+			// Calculate the normal component of the relative velocity (actually only the magnitude - the direction is stored in 'normal').
+			magrelVel = INPR ( relativeVelocity, collpair->normal );
+	
+			// If v_n_mag < 0 the edges are approaching each other.
+			if ( magrelVel > ALMOST_ZERO )
+			{
+				// Calculate Impulse magnitude to stop all motion in normal direction.
+				float magtangent = 0, repulse = 0, d = 0;
+				double impulse = 0.0;
+				float vrel_t_pre[3];
+				float temp[3], spf;
+	
+				// calculate tangential velocity
+				VECCOPY ( temp, collpair->normal );
+				mul_v3_fl( temp, magrelVel );
+				VECSUB ( vrel_t_pre, relativeVelocity, temp );
+	
+				// Decrease in magnitude of relative tangential velocity due to coulomb friction
+				// in original formula "magrelVel" should be the "change of relative velocity in normal direction"
+				magtangent = MIN2 ( clmd->coll_parms->friction * 0.01 * magrelVel,sqrt ( INPR ( vrel_t_pre,vrel_t_pre ) ) );
+	
+				// Apply friction impulse.
+				if ( magtangent > ALMOST_ZERO )
+				{
+					normalize_v3( vrel_t_pre );
+	
+					impulse = magtangent; // 2.0 * 
+					VECADDMUL ( cloth1->verts[collpair->collp].impulse, vrel_t_pre, impulse);
+				}
+	
+				// Apply velocity stopping impulse
+				// I_c = m * v_N / 2.0
+				// no 2.0 * magrelVel normally, but looks nicer DG
+				impulse =  magrelVel/2.0;
+	
+				VECADDMUL ( cloth1->verts[collpair->collp].impulse, collpair->normal, impulse);
+				cloth1->verts[collpair->collp].impulse_count++;
+	
+				// Apply repulse impulse if distance too short
+				// I_r = -min(dt*kd, m(0,1d/dt - v_n))
+				spf = (float)clmd->sim_parms->stepsPerFrame / clmd->sim_parms->timescale;
+	
+				d = -collpair->distance;
+				if ( ( magrelVel < 0.1*d*spf ) && ( d > ALMOST_ZERO ) )
+				{
+					repulse = MIN2 ( d*1.0/spf, 0.1*d*spf - magrelVel );
+	
+					// stay on the safe side and clamp repulse
+					if ( impulse > ALMOST_ZERO )
+						repulse = MIN2 ( repulse, 5.0*impulse );
+					repulse = MAX2 ( impulse, repulse );
+	
+					impulse = repulse / ( 5.0 ); // original 2.0 / 0.25
+					VECADDMUL ( cloth1->verts[collpair->collp].impulse, collpair->normal, impulse);
+				}
+	
+				result = 1;
+			}
+		} else {	
+			w1 = collpair->bary[0]; w2 = collpair->bary[1]; w3 = collpair->bary[2];			
+
+			// Calculate relative "velocity".
+			collision_interpolateOnTriangle ( v1, cloth1->verts[collpair->ap1].tv, cloth1->verts[collpair->ap2].tv, cloth1->verts[collpair->ap3].tv, w1, w2, w3 );
+	
+			VECSUB ( relativeVelocity, collmd->current_v[collpair->collp].co, v1);
+			
+			// Calculate the normal component of the relative velocity (actually only the magnitude - the direction is stored in 'normal').
+			magrelVel = INPR ( relativeVelocity, collpair->normal );
+	
+			// If v_n_mag < 0 the edges are approaching each other.
+			if ( magrelVel > ALMOST_ZERO )
+			{
+				// Calculate Impulse magnitude to stop all motion in normal direction.
+				float magtangent = 0, repulse = 0, d = 0;
+				double impulse = 0.0;
+				float vrel_t_pre[3], pimpulse[3] = {0.0f, 0.0f, 0.0f};
+				float temp[3], spf;
+	
+				// calculate tangential velocity
+				VECCOPY ( temp, collpair->normal );
+				mul_v3_fl( temp, magrelVel );
+				VECSUB ( vrel_t_pre, relativeVelocity, temp );
+	
+				// Decrease in magnitude of relative tangential velocity due to coulomb friction
+				// in original formula "magrelVel" should be the "change of relative velocity in normal direction"
+				magtangent = MIN2 ( clmd->coll_parms->friction * 0.01 * magrelVel,sqrt ( INPR ( vrel_t_pre,vrel_t_pre ) ) );
+	
+				// Apply friction impulse.
+				if ( magtangent > ALMOST_ZERO )
+				{
+					normalize_v3( vrel_t_pre );
+	
+					impulse = magtangent; // 2.0 * 
+					VECADDMUL ( pimpulse, vrel_t_pre, impulse);
+				}
+	
+				// Apply velocity stopping impulse
+				// I_c = m * v_N / 2.0
+				// no 2.0 * magrelVel normally, but looks nicer DG
+				impulse =  magrelVel/2.0;
+	
+				VECADDMUL ( pimpulse, collpair->normal, impulse);
+	
+				// Apply repulse impulse if distance too short
+				// I_r = -min(dt*kd, m(0,1d/dt - v_n))
+				spf = (float)clmd->sim_parms->stepsPerFrame / clmd->sim_parms->timescale;
+	
+				d = -collpair->distance;
+				if ( ( magrelVel < 0.1*d*spf ) && ( d > ALMOST_ZERO ) )
+				{
+					repulse = MIN2 ( d*1.0/spf, 0.1*d*spf - magrelVel );
+	
+					// stay on the safe side and clamp repulse
+					if ( impulse > ALMOST_ZERO )
+						repulse = MIN2 ( repulse, 5.0*impulse );
+					repulse = MAX2 ( impulse, repulse );
+	
+					impulse = repulse / ( 2.0 ); // original 2.0 / 0.25
+					VECADDMUL ( pimpulse, collpair->normal, impulse);
+				}
+				
+				if (w1 < 0.5) w1 *= 2.0;
+				if (w2 < 0.5) w2 *= 2.0;
+				if (w3 < 0.5) w3 *= 2.0;
+				
+				VECADDMUL(cloth1->verts[collpair->ap1].impulse, pimpulse, w1*2.0);
+				VECADDMUL(cloth1->verts[collpair->ap2].impulse, pimpulse, w2*2.0);
+				VECADDMUL(cloth1->verts[collpair->ap3].impulse, pimpulse, w3*2.0);;
+				cloth1->verts[collpair->ap1].impulse_count++;
+				cloth1->verts[collpair->ap2].impulse_count++;
+				cloth1->verts[collpair->ap3].impulse_count++;
+				
+				result = 1;
+			}
+		}
+	} 
+	
+	return result;
+}
+
+
+typedef struct tripairkey {
+	int p, a1, a2, a3;
+} tripairkey;
+
+unsigned int tripair_hash(void *vkey)
+{
+	tripairkey *key = vkey;
+	int keys[4] = {key->p, key->a1, key->a2, key->a3};
+	int i, j;
+	
+	for (i=0; i<4; i++) {
+		for (j=0; j<3; j++) {
+			if (keys[j] >= keys[j+1]) {
+				SWAP(int, keys[j], keys[j+1]);
+			}
+		}
+	}
+	
+	return keys[0]*101 + keys[1]*72 + keys[2]*53 + keys[3]*34;
+}
+
+int tripair_cmp(const void *va, const void *vb)
+{
+	tripairkey *a = va, *b = vb;
+	int keysa[4] = {a->p, a->a1, a->a2, a->a3};
+	int keysb[4] = {b->p, b->a1, b->a2, b->a3};
+	int i;
+	
+	for (i=0; i<4; i++) {
+		int j, ok=0;
+		for (j=0; j<4; j++) {
+			if (keysa[i] == keysa[j]) {
+				ok = 1;
+				break;
+			}
+		}
+		if (!ok)
+			return -1;
+	}
+	
+	return 0;
+}
+
+static void get_tripairkey(tripairkey *key, int p, int a1, int a2, int a3)
+{
+	key->a1 = a1;
+	key->a2 = a2;
+	key->a3 = a3;
+	key->p = p;
+}
+
+static int checkvisit(MemArena *arena, GHash *gh, int p, int a1, int a2, int a3)
+{
+	tripairkey key, *key2;
+	
+	get_tripairkey(&key, p, a1, a2, a3);
+	if (BLI_ghash_haskey(gh, &key))
+		return 1;
+	
+	key2 = BLI_memarena_alloc(arena, sizeof(*key2));
+	*key2 = key;
+	BLI_ghash_insert(gh, key2, NULL);
+	
+	return 0;
+}
+
+int cloth_point_tri_moving_v3v3_f(float v1[2][3], int i1, float v2[2][3], int i2,
+                                   float v3[2][3],  int i3, float v4[2][3], int i4,
+                                   float normal[3], float bary[3], float *t, 
+								   float *relnor, GHash *gh, MemArena *arena)
+{
+	if (checkvisit(arena, gh, i1, i2, i3, i4))
+		return 0;
+	
+	return eltopo_point_tri_moving_v3v3_f(v1, i1, v2, i2, v3, i3, v4, i4, normal, bary, t, relnor);
+}
+
+static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTreeOverlap *overlap, 
+								   CollPair *collpair, double dt, GHash *gh, MemArena *arena)
+{
+	ClothModifierData *clmd = ( ClothModifierData * ) md1;
+	CollisionModifierData *collmd = ( CollisionModifierData * ) md2;
+	MFace *face1=NULL, *face2 = NULL;
+	ClothVertex *verts1 = clmd->clothObject->verts;
+	double distance = 0;
+	float epsilon1 = clmd->coll_parms->epsilon;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
+	float no[3], uv[3], t, relnor;
+	int i, i1, i2, i3, i4, i5, i6;
+	Cloth *cloth = clmd->clothObject;
+	float n1[3], sdis, p[3], l, n2[3], off[3], v1[2][3], v2[2][3], v3[2][3], v4[2][3], v5[2][3], v6[2][3];
+	int j, ret, bp1, bp2, bp3, ap1, ap2, ap3;
+	
+	face1 = & ( clmd->clothObject->mfaces[overlap->indexA] );
+	face2 = & ( collmd->mfaces[overlap->indexB] );
+
+	// check all 4 possible collisions
+	for ( i = 0; i < 4; i++ )
+	{
+		if ( i == 0 )
+		{
+			// fill faceA
+			ap1 = face1->v1;
+			ap2 = face1->v2;
+			ap3 = face1->v3;
+
+			// fill faceB
+			bp1 = face2->v1;
+			bp2 = face2->v2;
+			bp3 = face2->v3;
+		}
+		else if ( i == 1 )
+		{
+			if ( face1->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v3;
+				ap3 = face1->v4;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v2;
+				bp3 = face2->v3;
+			}
+			else {
+				continue;
+			}
+		}
+		if ( i == 2 )
+		{
+			if ( face2->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v2;
+				ap3 = face1->v3;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v3;
+				bp3 = face2->v4;
+			}
+			else {
+				continue;
+			}
+		}
+		else if ( i == 3 )
+		{
+			if ( face1->v4 && face2->v4 )
+			{
+				// fill faceA
+				ap1 = face1->v1;
+				ap2 = face1->v3;
+				ap3 = face1->v4;
+
+				// fill faceB
+				bp1 = face2->v1;
+				bp2 = face2->v3;
+				bp3 = face2->v4;
+			}
+			else {
+				continue;
+			}
+		}
+		
+		copy_v3_v3(v1[0], cloth->verts[ap1].txold); 
+		copy_v3_v3(v1[1], cloth->verts[ap1].tx);
+		copy_v3_v3(v2[0], cloth->verts[ap2].txold);
+		copy_v3_v3(v2[1], cloth->verts[ap2].tx);
+		copy_v3_v3(v3[0], cloth->verts[ap3].txold);
+		copy_v3_v3(v3[1], cloth->verts[ap3].tx);
+		
+		copy_v3_v3(v4[0], collmd->current_x[bp1].co);
+		copy_v3_v3(v4[1], collmd->current_xnew[bp1].co);
+		copy_v3_v3(v5[0], collmd->current_x[bp2].co);
+		copy_v3_v3(v5[1], collmd->current_xnew[bp2].co);
+		copy_v3_v3(v6[0], collmd->current_x[bp3].co);
+		copy_v3_v3(v6[1], collmd->current_xnew[bp3].co);
+		
+		normal_tri_v3(n2, v4[1], v5[1], v6[1]);
+		
+		sdis = clmd->coll_parms->distance_repel + epsilon2 + FLT_EPSILON;
+
+		/*apply a repulsion force, to help the solver along*/
+		copy_v3_v3(off, n2);
+		negate_v3(off);
+		if (isect_ray_plane_v3(v1[1], off, v4[1], v5[1], v6[1], &l, 0)) {
+			if (l >= 0.0 && l < sdis) {
+				mul_v3_fl(off, (l-sdis)*cloth->verts[ap1].mass*dt*clmd->coll_parms->repel_force*0.1);
+
+				add_v3_v3(cloth->verts[ap1].tv, off);
+				add_v3_v3(cloth->verts[ap2].tv, off);
+				add_v3_v3(cloth->verts[ap3].tv, off);
+			}
+		}
+
+		/*offset new positions a bit, to account for margins*/
+		copy_v3_v3(off, n2);
+		mul_v3_fl(off,  epsilon1 + epsilon2 + ALMOST_ZERO);
+		add_v3_v3(v4[1], off); add_v3_v3(v5[1], off); add_v3_v3(v6[1], off);
+		
+		i1 = ap1; i2 = ap2; i3 = ap3;
+		i4 = bp1+cloth->numverts; i5 = bp2+cloth->numverts; i6 = bp3+cloth->numverts;
+		
+		for (j=0; j<6; j++) {
+			int collp;
+
+			switch (j) {
+			case 0:
+				ret = cloth_point_tri_moving_v3v3_f(v1, i1, v4, i4, v5, i5, v6, i6, no, uv, &t, &relnor, gh, arena);
+				collp = ap1;
+				break;
+			case 1:
+				collp = ap2;
+				ret = cloth_point_tri_moving_v3v3_f(v2, i2, v4, i4, v5, i5, v6, i6, no, uv, &t, &relnor, gh, arena);
+				break;
+			case 2:
+				collp = ap3;
+				ret = cloth_point_tri_moving_v3v3_f(v3, i3, v4, i4, v5, i5, v6, i6, no, uv, &t, &relnor, gh, arena);
+				break;
+			case 3:
+				collp = bp1;
+				ret = cloth_point_tri_moving_v3v3_f(v4, i4, v1, i1, v2, i2, v3, i3, no, uv, &t, &relnor, gh, arena);
+				break;
+			case 4:
+				collp = bp2;				
+				ret = cloth_point_tri_moving_v3v3_f(v5, i5, v1, i1, v2, i2, v3, i3, no, uv, &t, &relnor, gh, arena);
+				break;
+			case 5:
+				collp = bp3;
+				ret = cloth_point_tri_moving_v3v3_f(v6, i6, v1, i1, v2, i2, v3, i3, no, uv, &t, &relnor, gh, arena);
+				break;
+			}
+			
+			/*cloth vert versus coll face*/
+			if (ret && j < 3) {
+				collpair->bp1 = bp1; collpair->bp2 = bp2; collpair->bp3 = bp3;
+				collpair->collp = collp;
+				
+				copy_v3_v3(collpair->normal, no);
+				mul_v3_v3fl(collpair->vector, collpair->normal, relnor);
+				collpair->distance = relnor;
+				collpair->time = t;
+				
+				copy_v3_v3(collpair->bary, uv);
+				
+				collpair->flag = COLLISION_USE_COLLFACE;
+				collpair++;
+			} else if (ret && j >= 3) { /*coll vert versus cloth face*/
+				collpair->ap1 = ap1; collpair->ap2 = ap2; collpair->ap3 = ap3;
+				collpair->collp = collp;
+				
+				copy_v3_v3(collpair->normal, no);
+				mul_v3_v3fl(collpair->vector, collpair->normal, relnor);
+				collpair->distance = relnor;
+				collpair->time = t;
+				
+				copy_v3_v3(collpair->bary, uv);
+
+				collpair->flag = 0;
+				collpair++;
+			}
+		}
+	}
+	
+	return collpair;
+}
+
+static void machine_epsilon_offset(Cloth *cloth) {
+	ClothVertex *cv;
+	int i, j;
+	
+	cv = cloth->verts;
+	for (i=0; i<cloth->numverts; i++, cv++) {
+		/*aggrevatingly enough, it's necassary to offset the coordinates
+		 by a multiple of the 32-bit floating point epsilon when switching
+		 into doubles*/
+		#define RNDSIGN (float)(-1*(BLI_rand()%2==0)|1)
+		for (j=0; j<3; j++) {
+			cv->tx[j] += FLT_EPSILON*30.0f*RNDSIGN;
+			cv->txold[j] += FLT_EPSILON*30.0f*RNDSIGN;
+			cv->tv[j] += FLT_EPSILON*30.0f*RNDSIGN;
+		}		
+	}
+}
+
+#else
 
 //Determines collisions on overlap, collisions are written to collpair[i] and collision+number_collision_found is returned
-static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTreeOverlap *overlap, CollPair *collpair )
+static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, 
+	BVHTreeOverlap *overlap, CollPair *collpair, float dt )
 {
 	ClothModifierData *clmd = ( ClothModifierData * ) md1;
 	CollisionModifierData *collmd = ( CollisionModifierData * ) md2;
+	Cloth *cloth = clmd->clothObject;
 	MFace *face1=NULL, *face2 = NULL;
 #ifdef USE_BULLET
 	ClothVertex *verts1 = clmd->clothObject->verts;
@@ -614,6 +1408,7 @@ static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTree
 	double distance = 0;
 	float epsilon1 = clmd->coll_parms->epsilon;
 	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
+	float n2[3], sdis, l;
 	int i;
 
 	face1 = & ( clmd->clothObject->mfaces[overlap->indexA] );
@@ -685,7 +1480,28 @@ static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTree
 			else
 				break;
 		}
+		
+		normal_tri_v3(n2, collmd->current_xnew[collpair->bp1].co, 
+			collmd->current_xnew[collpair->bp2].co, 
+			collmd->current_xnew[collpair->bp3].co);
+		
+		sdis = clmd->coll_parms->distance_repel + epsilon2 + FLT_EPSILON;
+		
+		/*apply a repulsion force, to help the solver along.
+          this is kindof crude, it only tests one vert of the triangle*/
+		if (isect_ray_plane_v3(cloth->verts[collpair->ap1].tx, n2, collmd->current_xnew[collpair->bp1].co, 
+			collmd->current_xnew[collpair->bp2].co,
+			collmd->current_xnew[collpair->bp3].co, &l, 0))
+		{
+			if (l >= 0.0 && l < sdis) {
+				mul_v3_fl(n2, (l-sdis)*cloth->verts[collpair->ap1].mass*dt*clmd->coll_parms->repel_force*0.1);
 
+				add_v3_v3(cloth->verts[collpair->ap1].tv, n2);
+				add_v3_v3(cloth->verts[collpair->ap2].tv, n2);
+				add_v3_v3(cloth->verts[collpair->ap3].tv, n2);
+			}
+		}
+		
 #ifdef USE_BULLET
 		// calc distance + normal
 		distance = plNearestPoints (
@@ -741,6 +1557,8 @@ static CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTree
 	}
 	return collpair;
 }
+#endif
+
 
 #if 0
 static int cloth_collision_response_moving( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
@@ -1446,17 +2264,45 @@ void free_collider_cache(ListBase **colliders)
 	}
 }
 
-static void cloth_bvh_objcollisions_nearcheck ( ClothModifierData * clmd, CollisionModifierData *collmd, CollPair **collisions, CollPair **collisions_index, int numresult, BVHTreeOverlap *overlap)
+
+static void cloth_bvh_objcollisions_nearcheck ( ClothModifierData * clmd, CollisionModifierData *collmd,
+	CollPair **collisions, CollPair **collisions_index, int numresult, BVHTreeOverlap *overlap, double dt)
 {
 	int i;
+#ifdef USE_ELTOPO
+	GHash *visithash = BLI_ghash_new(edgepair_hash, edgepair_cmp, "visthash, collision.c");
+	GHash *tri_visithash = BLI_ghash_new(tripair_hash, tripair_cmp, "tri_visthash, collision.c");
+	MemArena *arena = BLI_memarena_new(1<<16, "edge hash arena, collision.c");
+#endif
 	
-	*collisions = ( CollPair* ) MEM_mallocN ( sizeof ( CollPair ) * numresult * 4, "collision array" ); //*4 since cloth_collision_static can return more than 1 collision
+	*collisions = ( CollPair* ) MEM_mallocN ( sizeof ( CollPair ) * numresult * 64, "collision array" ); //*4 since cloth_collision_static can return more than 1 collision
 	*collisions_index = *collisions;
+	
+#ifdef USE_ELTOPO
+	machine_epsilon_offset(clmd->clothObject);
 
 	for ( i = 0; i < numresult; i++ )
 	{
-		*collisions_index = cloth_collision ( ( ModifierData * ) clmd, ( ModifierData * ) collmd, overlap+i, *collisions_index );
+		*collisions_index = cloth_collision ( ( ModifierData * ) clmd, ( ModifierData * ) collmd,
+											  overlap+i, *collisions_index, dt, tri_visithash, arena );
 	}
+
+	for ( i = 0; i < numresult; i++ )
+	{
+		*collisions_index = cloth_edge_collision ( ( ModifierData * ) clmd, ( ModifierData * ) collmd,
+												   overlap+i, *collisions_index, visithash, arena );
+	}
+	BLI_ghash_free(visithash, NULL, NULL);
+	BLI_ghash_free(tri_visithash, NULL, NULL);
+	BLI_memarena_free(arena);
+#else
+	for ( i = 0; i < numresult; i++ )
+	{
+		*collisions_index = cloth_collision ( ( ModifierData * ) clmd, ( ModifierData * ) collmd,
+											  overlap+i, *collisions_index, dt );
+	}
+#endif	
+
 }
 
 static int cloth_bvh_objcollisions_resolve ( ClothModifierData * clmd, CollisionModifierData *collmd, CollPair *collisions, CollPair *collisions_index)
@@ -1481,11 +2327,19 @@ static int cloth_bvh_objcollisions_resolve ( ClothModifierData * clmd, Collision
 
 		if ( collmd->bvhtree )
 		{
+#ifdef USE_ELTOPO
+			result += cloth_collision_response_moving(clmd, collmd, collisions, collisions_index);
+			result += cloth_edge_collision_response_moving(clmd, collmd, collisions, collisions_index);
+#else
 			result += cloth_collision_response_static ( clmd, collmd, collisions, collisions_index );
-
+#endif
+#ifdef USE_ELTOPO
+			{
+#else
 			// apply impulses in parallel
 			if ( result )
 			{
+#endif
 				for ( i = 0; i < numverts; i++ )
 				{
 					// calculate "velocities" (just xnew = xold + v; no dt in v)
@@ -1518,7 +2372,7 @@ int cloth_bvh_objcollision (Object *ob, ClothModifierData * clmd, float step, fl
 
 	if ((clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ) || cloth_bvh==NULL)
 		return 0;
-
+	
 	verts = cloth->verts;
 	numfaces = cloth->numfaces;
 	numverts = cloth->numverts;
@@ -1557,6 +2411,7 @@ int cloth_bvh_objcollision (Object *ob, ClothModifierData * clmd, float step, fl
 				continue;
 			
 			/* move object to position (step) in time */
+			
 			collision_move_object ( collmd, step + dt, step );
 			
 			/* search for overlapping collision pairs */
@@ -1565,7 +2420,8 @@ int cloth_bvh_objcollision (Object *ob, ClothModifierData * clmd, float step, fl
 			// go to next object if no overlap is there
 			if( result && overlap ) {
 				/* check if collisions really happen (costly near check) */
-				cloth_bvh_objcollisions_nearcheck ( clmd, collmd, &collisions[i], &collisions_index[i], result, overlap);
+				cloth_bvh_objcollisions_nearcheck ( clmd, collmd, &collisions[i], 
+					&collisions_index[i], result, overlap, dt/(float)clmd->coll_parms->loop_count);
 			
 				// resolve nearby collisions
 				ret += cloth_bvh_objcollisions_resolve ( clmd, collmd, collisions[i],  collisions_index[i]);
@@ -1721,5 +2577,5 @@ int cloth_bvh_objcollision (Object *ob, ClothModifierData * clmd, float step, fl
 	if(collobjs)
 		MEM_freeN(collobjs);
 
-	return MIN2 ( ret, 1 );
+	return 1|MIN2 ( ret, 1 );
 }