• Advertisement
Sign in to follow this  

Sphere Collision Response with Verlet

This topic is 4409 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I'm updating the physics engine I'm working on from Euler to Verlet integration. In the old system, if two spheres overlapped I'd put a soft spring between them. The result was a soft repulsive effect that reminded me of colliding marshmallows. While there's nothing preventing me from doing the same in my new Verlet method, verlet is the sort of system that is supposed to handle this sort of thing quite easily. This is the code I'm using:
Vector4 ab = pos2 - pos1;
float mindist = radius1 + radius2;
mindist *= mindist;
            
if(float currdist = LengthSquared3(ab) < mindist)
{
  //relaxation collision technique
  //from http://www.gamasutra.com/resource_guide/20030121/jacobson_03.shtml
  
  ab *= (mindist / (mindist + currdist));
  
  //ab is the distance vector that is being overlapped
                
  pos1 -= ab * .5;
  pos2 += ab * .5;
}
The result is incredibly bouncy collisions that I seriously doubt are conserving momentum or kinetic energy. Also, the above doesn't take into account relative masses. All the collision response routines that I find are either of the version I used with the springs or kinematic versions using Euler integration (specifically, they change a velocity vector). Ideally I'd have a gradation the user could set between perfectly elastic and perfectly inelastic collisions, but even a single working version somewhere between the two extremes would make me content. How exactly do I respond to collisions kinematically without a velocity vector?

Share this post


Link to post
Share on other sites
Advertisement
What happens if you replace:
ab *= (mindist / (mindist + currdist));
with
ab *= sqrt(mindist - currdist);
which should give you extremely non-elastic collisions (albeit inefficiently)?

Share this post


Link to post
Share on other sites
It's even more 'bouncy'.

I'm sure it has something to do with the fact that opos isn't being modified (since velocity is implicitly given by the opos vector) but more than that I'm not sure.

Share this post


Link to post
Share on other sites
Well, there was no 'opos' in your original post ;)
I got my equation wrong... should have been
ab *= sqrt(mindist) - sqrt(currdist)
which is subtle as a brick but should be very non bouncy.

Share this post


Link to post
Share on other sites
Oops, well, I guess I sort of figured opos would be implied by verlet.

I'll look see if your idea works.

edit:

It's still orders of magnitude too 'bouncy'. To give you an idea of how far off it is,

pos1 -= ab * .5 * .0001;
pos2 += ab * .5 * .0001;

looks like collisions with those rubber bouncy balls.

Share this post


Link to post
Share on other sites
OK. Maybe if I post equations after my morning cup of tea we will get further ;)

I'm going to rename all your variables, and introduce some new ones, because using the same variable for two different meanings makes it hard to read. We can put them back how they were once it works :p

Vec gap = pos2 - pos1;
float distsqr = LengthSquared(gap);
float mindist = radius2 + radius1;

if (distsqr < mindist * mindist)
{
Vec gapdir = normal(gap);
float depth = mindist - sqrt(distsqr)
pos2 += depth * gapdir * 0.5;
pos1 -= depth * gapdir * 0.5;
}

Then they have an approximation to an optimization of the above, based on:
gapdir * depth == (gap / dist) * depth
depth == mindist - dist
so:
gapdir * depth = gap * (mindist - dist) / dist


also, just to make sure we aren't being really dumb, could you try changing
if(float currdist = LengthSquared3(ab) < mindist)
to
if((float currdist = LengthSquared3(ab)) < mindist)
?



Share this post


Link to post
Share on other sites
Ripped straight from the paper that started the verlet craze: Advanced Character Physics by Thomas Jakobsen

// Pseudo-code to satisfy (C2)
delta = x2-x1;
deltalength = sqrt(delta*delta);
diff = (deltalength-restlength) /(deltalength*(invmass1+invmass2));
x1 -= invmass1*delta*diff;
x2 += invmass2*delta*diff;

If your spheres have a non-zero radius (of course they do), then you probably need to subtract the sum of the radii from deltalength after the sqrt...

Cheers,
Mark

Share this post


Link to post
Share on other sites
After significant playing around, this is what I finally ended up with:

ab *= sqrt(mindist) / sqrt(currdist) - 1.0f;

pos1 += ab * 0.5f;
pos2 -= ab * 0.5f;

Which seems to work just fine.

I also switched from regular verlet (I don't know what it's called) to velocity verlet, which has an explicit velocity term. Using the velocity term it provides it isn't too hard to adapt most common collision response algorithms.

My next question is about the elastic/inelastic gradient I want to set. The only collision response algorithm I can find that provides such a capability is in my reference book "Tips of the Windows Game Programming Gurus".

It calls such a variable the "coefficient of restitution" and it's apparently given by:

sphere2's final velocity along normal - sphere1's final velocity along normal
sphere2's initial velocity along normal - sphere1's initial velocity along normal

Which makes sense. 0 would be totally inelastic and 1 would be perfectly elastic.

Unfortunately, the algorithm it gives doesn't use a vector class, and the variable names aren't all that descriptive, and it's for the 2D case. I could just copy and paste the code and use it that way, but I'd really prefer to construct the general case so if/when I move into 3D, I don't have to go rewriting the algorithm, or rewrite it only slightly.

Basically, the algorithm looks like:

1. Compute the normal and tangental unit vectors for the collision. (normal = pos2 - pos1 normalized).

2. Find the velocity of the two spheres in terms of the normal and tangental vectors.

3. Compute the new collision velocities using an equation that takes into account the coefficient of restitution.

4. The results are in terms of the normal and tangental vectors, so transform them back into the original x,y coordinate system.

Does this look familiar to anyone? How do I extend that to the general 3D case? I can elaborate more if anyone needs it, but I'm mostly hoping there's an article online somewhere that does all the work for me.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement