**Looks like the eqn tag died, tried to make readable "ascii math"**

I tried asking this question on MathOverflow but it doesn't seem to be gaining any speed there... This forum seems much better suited. So I saw some papers on laproscopic surgery and simulating thread and thought "That would make some wicked cool rope to play with" Something to shove together with my Razer Hydra and Oculus. The most current of the papers is this which in turn references this other paper.

In Müller's paper he talks about constraint functions being C : R^3 -> R Which makes sense because using the constraint solver you solve attempt to get each constraint function either equal to 0 or greater than or equal to 0.

However if you look at fratarcangeli's paper he gives the contact constraint function as

C(p) = [p - (p_n0 + p_v)]

Where p must be some vertex of the rope, p_v is the penetration vector and p_n0 is "the current position of point p. This is where things stop making sense for me. Because it appears that fratarcangeli's constraint equation is in R^3 and not R}. Perhaps I'm miss-understanding the equation?

My second issue with his constraint function is

p_v = (|p_n0 - p_n1| - r) dot (p_n0 - p_n1)/|p_n0 - p_n1|

and he gives a very loose definition of what p, p_n0, p_n1 are.

Perhaps someone can explain what his constraint function is supposed to be?

---------------------------

I attempted to figure out what the constraint function should be.

I assumed if two segments were colliding I would have to apply a contact constraint to all 4 points. Since I'm applying the constraint to both sides I only need to move each mass point halfway out of the collision.

When a collision occurs between two line segments p_1 -> q_1 and p_2 -> q_2 And p = p1. I get the two closest points c_1 and c_2 on those segments respectively.Most of the time c_1 != p. So I have to define my constraint function with that in mind.

I call the collision normal n = (c_1 - c_2) / |c_1 - c_2| and an offset o = (p - c_1) dot n which is the offset along the collision normal of p down to the contact point. This handles when c_1 != p Note: o is calculated once at the beginning of a collision, it's expected to stay constant.

The goal is to have the constraint function equal to 0 when p has moved halfway to resolve the collision (the other segment will move the other half) and > 0 when the point has moved further than halfway.

So I define

C(p) = -(2r - ((p - c_2) dot n - o))/2

Here's a poorly drawn diagram to illustrate my thoughts.

Which when I punch that through the method described in Müller's paper. I get delta_p = -2C(p)n However when I plug this into my simulation it's stable up until I tie a knot which given the nature of the papers means I've done something wrong. Can anyone elaborate on where I'm going wrong?