Here's an idea of what the table of friction coefficients might look like (and sorry about the formatting)... each cell contains a pair of values... static and dynamic friction coefficients for given pair of materials... or just 0 for "unknown" - this is just me starting to collate some useful data, wanna help fill some of these fields?

Slippery conditions = something slippery is lubricating the surface(s), eg oil or water
Sticky conditions = something like race tyres on concrete

``.dataMyTableElement struct	_Static  real4 ?	_Dynamic real4 ?MyTableElement endsMyTableRow struct	Wood	MyTableElement <>	Stone	MyTableElement <>	Ice		MyTableElement <>	Glass	MyTableElement <>	Metal	MyTableElement <>	Rubber	MyTableElement <>	Velcro	MyTableElement <>	Meat	MyTableElement <>MyTableRow ends;--------------------------------------------------------------------------------------------------------------------------------------;				WOOD			STONE			ICE			GLASS			METAL			RUBBER			VELCRO			MEAT;--------------------------------------------------------------------------------------------------------------------------------------;-------------------;NORMAL CONDITIONS;-------------------;WOODMyTableRow	{<0.0f,0.0f>,	<0.5f,0.4f>,	<0.2f,0.1f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};STONEMyTableRow	{<0.5f,0.4f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<1.0f,0.8f>,	<0.0f,0.0f>,	<0.0f,0.0f>};ICEMyTableRow	{<0.2f,0.1f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.1f,0.03f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};GLASSMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.1f,0.03f>,	<0.95f,0.4f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};METALMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.6f,0.4f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};RUBBERMyTableRow	{<0.0f,0.0f>,	<1.0f, 0.8f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};VELCROMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<6.0f,4.0f>,	<0.0f,0.0f>};MEATMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};-------------------;STICKY CONDITIONS;-------------------;WOODMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};STONEMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<1.5f,1.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};ICEMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};GLASSMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};METALMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};RUBBERMyTableRow	{<0.0f,0.0f>,	<1.5f,1.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};VELCROMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};MEATMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};---------------------;SLIPPERY CONDITIONS;---------------------;WOODMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};STONEMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.7f,0.5f>,	<0.0f,0.0f>,	<0.0f,0.0f>};ICEMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};GLASSMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};METALMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.1f,0.05f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};RUBBERMyTableRow	{<0.0f,0.0f>,	<0.7f,0.5f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};VELCROMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>};MEATMyTableRow	{<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>,	<0.0f,0.0f>}.code``
Posted on 2008-06-02 22:44:26 by Homer

Here's a pretty picture of a collision thats about to occur.
Its showing some vectors with respect to the Center of Mass, but in fact we'll be interested in these vectors with respect to each Point of Collision.

The green line indicates the closing velocity of the Body.
The red line indicates the component of the closing velocity which is perpendicular with the collision normal.... our coefficient of restitution plays against the Red Line.
The yellow line indicates the component of the closing velocity which is tangent to the collision normal.... our coefficients of friction play against the Yellow Line.

I just wanted to visually show you what these vectors are because we'll be using them.

Posted on 2008-06-02 23:21:08 by Homer
RELATIVE VELOCITY
If our Body is colliding with some unmovable scenery in the World (such as our World Bounding Planes), then the Relative Velocity of the Collision is simply that of the Body.

RelativeVelocity = BodyA.CMVelocity

But if our Body is colliding with another Body, then the Relative Velocity is given as:

RelativeVelocity = BodyA.CMVelocity - BodyB.CMVelocity

COLLISION NORMAL
If our Body is colliding with some unmovable scenery in the World (such as our World Bounding Planes) then the Normal is that of the collision surface.. ie for collision with a flat floor plane, the Collision Normal would be <0.0, 1.0, 0.0>

CollisionNormal = SurfaceNormal

But if our Body is colliding with another Body, then the Collision Normal is taken from the relative velocities of the Bodies (at center of mass), and by convention, it is with respect to Body A:

CollisionNormal = Normalize (RelativeVelocity )

Note that the Collision Normal is a NORMAL (has a length of 1.0)

CLOSING VELOCITY
With the correct Collision Normal in hand, the closing velocity is given as:

ClosingVelocity = DotProduct (RelativeVelocity , CollisionNormal)

That works regardless of what kind of collision we have.
Note that this scalar value is with respect to the Center(s) of Mass, and does not tell us enough information alone to resolve the collision - it does not take rotation into account ("angular effects"). But its a good start.

In the next post, we'll begin to look at how we resolve a collision at the Particle level.
Then we'll expand on that to handle RigidBody collisions.

Posted on 2008-06-03 03:43:06 by Homer
Let us consider a single Point-collision between two Bodies (or a Body and a fixed surface).
We are handed the WorldSpace coordinate of the point of collision, and we are able to transform WorldSpace coordinates into any given BodySpace.
At the moment of impact, we are certain that if we transform the collision-point into the space of any Body involved in the collision, that Point will rest apon the Bounding Hull of that Body.
In other words, for a Body-Body collision, we can find the BodySpace coordinate of the collision point for EACH body's space.
I'm trying to reinforce something I mentioned earlier - that a Point collision actually involves TWO points - one on each of the colliding entities that happen to meet when the collision occurs.
These Points are the places on each entity where the collision-response Impulse should be applied :)

So we can begin to describe a general-purpose method for resolving point-based collisions...
Posted on 2008-06-03 05:55:53 by Homer
Here's an example function which resolves the collision of two PARTICLES.
Particles don't have any Orientation, so this code does NOT handle 'Angular Effects'.
Because of this, we will *NOT* be using this function.
It's simply an example.
Remember that ALL collisions can be expressed in terms of Particle collisions, so the solution for a Rigid Body is sure to be quite similar to this :)

``Method Particle.resolveVelocity,uses esi,pOtherParticle, pvContactNormal, duration:real8LOCAL separatingVelocity:real8,newSepVelocity:real8,restitution:real8, accCausedSepVelocity:real8, impulse:real8LOCAL deltaVelocity:real8,ftemp:real8LOCAL accCausedVelocity:Vec3, impulsePerIMass:Vec3	SetObject esi	; Find the velocity in the direction of the contact.	OCall esi.calculateSeparatingVelocity, pOtherParticle	fstp separatingVelocity		; Check whether it needs to be resolved.	.if separatingVelocity > 0		; The contact is either separating or stationary - there�s no impulse required.		ExitMethod	.endif		OCall esi.calculateRestitution, pOtherBody	fstp restitution	; Calculate the new separating velocity.	; newSepVelocity = -separatingVelocity * restitution;	fld separatingVelocity	fchs	fmul restitution	fstp newSepVelocity		; Check the velocity build-up due to acceleration only.	OCall esi.getLastFrameAcceleration	__StowFloat3 accCausedVelocity	 	.if pOtherParticle!=0		OCall pOtherBody::Particle.getLastFrameAcceleration		__SubFromFloat3 accCausedVelocity, accCausedVelocity	.endif		;accCausedSepVelocity = accCausedVelocity * contactNormal * duration	mov edx,pvContactNormal	DotProduct accCausedVelocity, .Vec3	fmul duration	fstp accCausedVelocity		; If we�ve got a closing velocity due to acceleration build-up,	; remove it from the new separating velocity.	.if accCausedSepVelocity < 0		;newSepVelocity += restitution * accCausedSepVelocity;		fld accCausedSepVelocity		fmul restitution		fadd newSepVelocity		fstp newSepVelocity		; Make sure we haven�t removed more than was there to remove.		.if (newSepVelocity < 0)			fmov newSepVelocity , 0.0		.endif	.endif		;deltaVelocity = newSepVelocity - separatingVelocity	fld newSepVelocity	fsub separatingVelocity	fstp deltaVelocity	; We apply the change in velocity to each object in proportion to	; its inverse mass (i.e., those with lower inverse mass  get less change in velocity).	fld .inverseMass		.if pOtherParticle!=0		mov edx,pOtherParticle		fadd .Particle.inverseMass		 totalInverseMass += particle->getInverseMass();	.endif	fstp totalInverseMass		; If all particles have infinite mass, then impulses have no effect.	.if totalInverseMass <= 0		ExitMethod	.endif		; Calculate the impulse to apply.	;impulse = deltaVelocity / totalInverseMass;	fld deltaVelocity	fdiv totalInverseMass	fstp impulse		; Find the amount of impulse per unit of inverse mass.	; impulsePerIMass = contactNormal * impulse	mov edx,pvContactNormal	CrossProduct .Vec3, impulse	__StowFloat3 impulsePerIMass		; Apply impulses in the direction of the contact, and proportional to the inverse mass.		;particle.velocity += impulsePerIMass * particle.inverseMass	__ScaleFloat3 impulsePerIMass, .inverseMass	__AddToFloat3 .Velocity, .Velocity	.if pOtherParticle!=0			; Particle 1 goes in the opposite direction - note how we flip the sign		;particle.velocity += impulsePerIMass * -particle.inverseMass		mov edx,pOtherParticle		fld .Particle.inverseMass		fchs		fstp ftemp				__ScaleFloat3 impulsePerIMass, ftemp		__AddToFloat3 .Particle.Velocity, .Particle.Velocity			.endif	MethodEnd``

Next we'll expand on this to handle angular effects for rigid bodies, and then we'll add isotropic friction.
Posted on 2008-06-03 22:49:48 by Homer
Before I go on, I'd like to talk about Contact Space.
Some of the math problems we're going to face (especially friction) can be rather difficult to solve if we attempt to solve them in the spatial systems we're familiar with (worldspace, bodyspace).... but can easily be solved if we move the problem into CONTACT space.
In this spatial coordinate system, one of the Axes (say, X) will be our Contact Normal.
The other two othogonal axes (Y, Z) are arbitrary - we'll look at one way to choose them.

In order to transform coordinates from worldspace to Contact space, we'll build a 3x3 transformation matrix. Since this matrix is actually a pure rotation matrix, we can use the Transpose Trick to transform coordinates back the other way (contact-to-world).

Now I'm going to show you a NEW TRICK.
This is a really efficient way to build a Rotation Matrix, given just one Vector, and the assumption that this Vector will always be the X Axis in the new coordinate system.

The 3x3 matrix we're building can be thought of as a set of 3 column-vectors.
Given a vector with components abc, the matrix will look like this:
a d g
b e h
c f i

So we'll need to find two vectors (def and ghi) which are orthogonal to abc and each other.
Orthogonal means that they're 90 degrees apart.
The most common ways to find these vectors involve finding the major and minor axes used by our input vector, then performing some cross-products.
I'm going to show you a cheating way which just makes the above a stupid idea.
Imagine if we have a 2D (XY) set of axes, if we swap the X and Y values, we just rotated the system by 90 degrees - yes?
The same is true in a 3D system - except that we ROTATE the axes.
abc = acb
def = bac
ghi = cba

So all we have to do is write the same vector component values out into each Column of our matrix, rotating the components as we go.
If the input vector was (0.0,0.86,0.5) , then the matrix would be:

0.0 0.5 0.86
0.86 0.0 0.5
0.5 0.86 0.0

We'll shove our Contact Normal in the first column, and fill out the rest as described.
Once we've shoved those orthogonal axis vectors in the two remaining columns, we're done, thats our transform matrix - how we FOUND these orthogonal axes is NOT RELEVANT - if we can find them "as a Given" and with no cost - then why shouldn't we use them?

Perhaps this solution isnt quite right - maybe we need to also toggle the Sign in places - who can contribute to solve this rather simple seeming problem?

Actually, I'm SURE that some of the signs are wrong, because that looks suspiciously like a General Rotation Matrix to me !!
However - we'll be transforming data into and out of the coordinate system described by our matrix - its only being used as a TRANSIENT space - we don't care about the DIRECTION of the two discovered axes - so the SIGN issue really is not an issue at all :)

Posted on 2008-06-04 02:27:19 by Homer

Anyway, just in case my little cheat is insufficient / Direction is important to us, the following macro was added to the Handy Math include (attached previously).
It creates a rotation matrix that transforms coordinates from 'Contact Space' into WorldSpace.
Since its a pure Rotation matrix, we can use the jolly Transpose trick to transform from WorldSpace to ContactSpace ("TransMult" function)...

``;Creates an Orthogonal set of Normal Vectors (axes) from a given Vector, as Columns of a Mat33;DO NOT USE EDXMat33_OrthoNormalBasis macro Mat,V	;The X column vector is taken directly from our collision Normal    m2m Mat.m00, V.x    m2m Mat.m10, V.y    m2m Mat.m20, V.z        ;Pick an axis which is orthogonal to the Major axis    fAbsMax V.x, V.y    fstpReg edx    .if edx==V.x            fld Mat.m00        fchs                fld Mat.m20                 ;s = 1 /sqrt(m20*m20 + m00*m00)        ; The new X-axis is at right angles to the world Y-axis     ;   m01 = m20*s     ;   m11 = 0     ;   m21 = m00*s        ; The new Y-axis is at right angles to the new X- and Z- axes      ;  m02 = m10*m01      ;  m12 = m20*m01 - m00*m21      ;  m22 = -m10*m01               fld1         fld Mat.m20         fmul st(0),st(0)         fld Mat.m00         fmul st(0),st(0)         fadd         fsqrt         fdiv                  fmul st(2),st(0)         fmul        mov  Mat.m11,0        fstp Mat.m01               fstp Mat.m21        fld  Mat.m10        fmul Mat.m01        fstp Mat.m02        fld  Mat.m20        fmul Mat.m01        fld  Mat.m00        fmul Mat.m21        fsub                fstp Mat.m12        fld  Mat.m10                fmul Mat.m01        fchs        fstp Mat.m22    .else               fld Mat.m10        fld Mat.m20        fchs                ;s = 1/sqrt(m20*m20 + m10*m10)        ; The new X-axis is at right angles to the world X-axis       ; m01 = 0       ; m11 = -m20*s       ; m21 = m10*s        ; The new Y-axis is at right angles to the new X- and Z- axes        ;m02 = m10*m21 - m20*m11        ;m12 = -m00*m21        ;m22 = m00*m11                 fld1         fld Mat.m20         fmul st(0),st(0)         fld Mat.m10         fmul st(0),st(0)         fadd         fsqrt         fdiv		 ;		 fmul st(2),st(0)		 fmul		 ;        mov  Mat.m01,0        fstp Mat.m11        fstp Mat.m21        fld  Mat.m10        fmul Mat.m21        fld  Mat.m20        fmul Mat.m11        fsub        fstp Mat.m02        fld  Mat.m00                fmul Mat.m21        fchs                fstp Mat.m12        fld  Mat.m00        fmul Mat.m11        fstp Mat.m22    .endif     endm``

Posted on 2008-06-04 23:16:29 by Homer

Now that we're able to toss our collision Normal at a macro and spit out a matrix, we're ready to go on.

Our goal is to calculate the impulse at the collision, right?
And from there, we'd like to be able to calculate the change in velocity of each object, given that  impulse.
For the moment, we're NOT considering Friction - and for Frictionless collisions, the impulses generated at the contact are only applied along the collision normal.
Read that again - the IMPULSE is applied ONLY along the collision normal.
We�d like to end up with a simple number that tells us the change in velocity AT THE CONTACT POINT, in the direction of the contact normal, for each unit of impulse, something like a SIGNED SCALAR (signed float used as multiplier) would be nice, yes?
It's important to realize that the linear and angular stuff in our Bodies is separate stuff. We have to treat them separately, and then combine them at the end.

The linear component is very simple, noting that this is ONLY for the Linear component
The linear change in velocity for a unit impulse will be in the direction of the impulse, with a magnitude given by the inverse mass:
DeltaVelocity * Impulse = inverseMass
If theres two Bodies involved, inverseMass will be the Sum of their masses (add them).

The angular stuff is a bit more tricky, we have to do it in three parts.
``torquePerUnitImpulse =CrossProduct (relativeContactPosition , contactNormal)rotationPerUnitImpulse =inverseInertiaTensor * torquePerUnitImpulsevelocityPerUnitImpulse =CrossProduct (rotationPerUnitImpulse , relativeContactPosition)``

The first equation tells us the amount of impulsive torque generated from a unit of impulse.
The second one tells us the change in angular velocity for a unit of impulsive torque.
And the third tells us the velocity at the contact point which is due only to the body rotating...

It is a velocity in world space. When it comes to implementing Friction, we'll need to transform this Worldspace vector into contact coordinates. This would give us a vector of velocities that a unit impulse would cause. Note that since we DONT have Friction, we're only interested in the velocity in the direction of the contact normal.
In contact coordinates this is the X axis, so our value is the x component of the resulting vector:
velocityPerUnitImpulseContact = Mat33_transMult (ContactToWorld,velocityPerUnitImpulse)
angularComponent = velocityPerUnitImpulseContact.x

Of course we COULD do it that way.
But we already know a faster way to find the scalar component of one Vector along another:
angularComponent  = DotProduct (velocityPerUnitImpulse , contactNormal)

For frictionless collisions, the result will be the same, so for frictionless collisions,we don't actually NEED to use the ContactSpace matrix.
But I wanted to show it to you now so that it will 'click' in your mind when we do get around to adding Friction.

For each object in the collision, we're now able to calculate the change in both Linear and Angular velocity: at the contact point, and PER UNIT impulse. We now add those velocity values together to get one overall value (even if theres two bodies) which describes the total change in velocity at the collision point.

We now know the closing velocity at the contact point.
The desired separating velocity is the negative of that.
Since we now know the desired separating velocity, we can now calculate an impulse which would produce that velocity change.
Then we just need to Apply it..
Posted on 2008-06-05 00:10:50 by Homer
It's time to start putting together everything we've learned :)

``;This is the entrypoint for collision response.Method CollisionHull.CalculateResponse, uses esi, pTargetConfigA,pTargetConfigB,pvContactPoint,pvContactNormal, pBodyB, bWantFriction:BOOLLOCAL relativeContactPositionA:Vec3LOCAL relativeContactPositionB:Vec3LOCAL contactVelocity:Vec3LOCAL desiredDeltaVelocity:real8,  restitution:real8LOCAL velocityFromAcc:real4LOCAL impulsiveTorque:Vec3,impulse:Vec3 ;,rotationChange:Vec3LOCAL ContactToWorld:Mat33	SetObject esi    ; Calculate a set of axes at the contact point.    mov eax,pvContactNormal    Mat33_OrthoNormalBasis ContactToWorld, .Vec3    ; Find the relative position of the contact    mov eax,pvContactPoint    mov edx,pTargetConfigA    __SubFloat3 .Vec3, .configuration.CMPosition    __StowFloat3 relativeContactPositionA        fld .CoefficientOfRestitution    .if pBodyB!=0	    mov eax,pvContactPoint	    mov edx,pTargetConfigB	    __SubFloat3 .Vec3, .configuration.CMPosition	    __StowFloat3 relativeContactPositionB	    mov edx,pBodyB	    fadd .CollisionHull.CoefficientOfRestitution	    fmul r4_Half    .endif    fstp restitution    ; Find the relative velocity at the contact point.    OCall esi.calculateLocalVelocity, pTargetConfigA,addr relativeContactPositionA, addr ContactToWorld	__StowFloat3 contactVelocity    .if pBodyB!=0    	OCall pBodyB::CollisionHull.calculateLocalVelocity, pTargetConfigB, addr relativeContactPositionB, addr ContactToWorld    	__SubFromFloat3 contactVelocity, contactVelocity    .endif    ; Calculate the acceleration-induced velocity accumulated this frame    mov velocityFromAcc,0    .if .bIsAwake==TRUE    	mov eax,pTargetConfigA    	mov edx,pvContactNormal        DotProduct .Vec3, .configuration.CMVelocityChange        fstp velocityFromAcc    .endif	mov edx,pBodyB    .if edx!=0 && .CollisionHull.bIsAwake==TRUE     	    	mov eax,pTargetConfigB    	mov edx,pvContactNormal        DotProduct .Vec3, .configuration.CMVelocityChange        fchs        fadd velocityFromAcc        fstp velocityFromAcc    .endif     velocityLimit equ r4_Quarter     ; Calculate the desired change in velocity for resolution	; desiredDeltaVelocity = -contactVelocity.x -thisRestitution * (contactVelocity.x - velocityFromAcc)    ; If the velocity is very slow, limit the restitution    fld contactVelocity.x    fchs    fld contactVelocity.x    fabs    fsub velocityLimit    fstpReg eax    .ifBitSet eax,BIT31	    ;if abs(contactVelocity.x) < velocityLimit            ;thisRestitution = 0    .else    	fsub restitution    .endif    fld contactVelocity.x    fsub velocityFromAcc    fmul    fstp desiredDeltaVelocity``

Setting up stuff that we'll need to resolve the collision... pretty much everything we do will be with respect to the Contact point... we're solving each point-collision with respect to BOTH bodies (if theres two), which reduces the amount of calculations overall.

1. The first thing we do is set up our Contact-to-World transform.
2. Then we calculate the position of the contact point, relative to the Body or Bodies.
We also compute the average Coefficient of Restitution for the Body or Bodies.
3. Next, we calculate the absolute velocity of the contact point, relative to the Bodies, and in Contact coords.
4. And then we calculate the magnitude of the change in velocity due solely to acceleration during this Frame along the collision normal.
5. Finally, we calculate the change in velocity we'd LIKE to see along the collision normal.
This, combined with Step 4, will allow us to SCALE the output to exactly what we'd LIKE to see.

Stay tuned for the rest of this method!
Posted on 2008-06-05 21:07:48 by Homer
Here's the rest of that method.

``	; Calculate impulse in Contact coords	.if bWantFriction==FALSE		OCall esi.CalculateFrictionlessImpulse, addr impulse, pvContactNormal, addr ContactToWorld, addr relativeContactPositionA, addr relativeContactPositionB, pTargetConfigA, pBodyB, pTargetConfigB	.else		;Look up the Friction Coefficients based on the Materials of both Bodies		mov eax,sizeof MyTableElement		mul .dMaterialType		.if pBodyB!=0			mov edx,pBodyB			push eax			mov eax,sizeof MyTableRow			mul .CollisionHull.dMaterialType			pop edx			add eax,edx		.endif		lea edx,FrictionTable		add eax,edx					OCall esi.CalculateFrictionImpulse, addr impulse, addr ContactToWorld, addr contactVelocity,addr relativeContactPositionA,addr relativeContactPositionB,pTargetConfigA,pBodyB,pTargetConfigB, eax	.endif    ; Convert impulse from Contact to World coords    Vec3_mult_Mat33 ContactToWorld, impulse	__StowFloat3 impulse    ; Split the impulse into linear and angular components,    ; and apply them to the Body(s)     CrossProduct relativeContactPositionA, impulse    __StowFloat3 impulsiveTorque       mov eax,pTargetConfigA   	; This is how we would convert ImpulseTorque into Delta-Angular-Velocity	; BUT WE WILL NOT BOTHER because this Simulator loves Angular Momentum 	; and so it will extract the Angular Velocity directly from Angular Momentum   ; Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor,  impulsiveTorque   ; __StowFloat3 rotationChange        ; Here we're getting linear acceleration (delta-velocity) from the impulse by removing Mass	__ScaleFloat3 impulse, .OneOverMass    __AddToFloat3 .vAccumulated_Delta_Velocity,  .vAccumulated_Delta_Velocity    ; Here we're just collecting the change in angular momentum due to Torque    __LoadFloat3 impulsiveTorque	__AddToFloat3 .vAccumulated_Delta_AngMomentum, .vAccumulated_Delta_AngMomentum			inc .dNumSimultaneousCollisions    .if pBodyB!=0            ; Work out body one's linear and angular changes                ; NOTE :: Crossed members have been SWITCHED        ; We have FLIPPED THE SIGN - we want to apply the equal, opposite impulse		CrossProduct impulse, relativeContactPositionB		__StowFloat3 impulsiveTorque  		     	;	This is how we would convert ImpulseTorque into Delta-Angular-Velocity	;	BUT WE WILL NOT BOTHER because this Simulator loves Angular Momentum 	;	and so it will extract the Angular Velocity directly from Angular Momentum	;	mov eax,pTargetConfigB	;	Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor, impulsiveTorque	;	__StowFloat3 rotationChange				; Here we're getting linear acceleration (delta-velocity) from the impulse by removing Mass		; Note that this time we have FLIPPED THE SIGN - we want to apply the equal, opposite impulse		mov edx,pBodyB				__NegScaleFloat3 impulse, .CollisionHull.OneOverMass    	__AddToFloat3 .CollisionHull.vAccumulated_Delta_Velocity,  .CollisionHull.vAccumulated_Delta_Velocity    	; Collect ye olde momentum    	__LoadFloat3 impulsiveTorque    		__AddToFloat3 .CollisionHull.vAccumulated_Delta_AngMomentum, .CollisionHull.vAccumulated_Delta_AngMomentum		inc .CollisionHull.dNumSimultaneousCollisions    .endif	MethodEnd``

1. We calculate the Impulse using either with or without Friction.
2. The impulse we calculated is in Contact coords, so we convert it to World coords.
3. We split the impulse into Linear and Angular components, and apply them to the Body(s).

As you can see, I'm not applying them directly - I'm accumulating them so that each Body has the average effect of N impulses applied to it simultaneously - this isnt strictly correct, the real problem is a simultaneous equation, and the real solution involves 'Jacobians', something I'm not too familiar with... but the average of the change in momenta due to N impulses is pretty close and looks believable at a fraction of the cost of the Jacobian solution.

Although theres a lot going on here, the really interesting stuff is buried out of sight.. Now I suppose you want to see the new methods whose names were exposed in that source?
We'll do that in the next post :)

PS : Attached is an update of the Handy Math includefile, which contains some new stuff including the __NegScaleFloat3 macro, as well as 90% of the other odd-looking macros you've seen me using (the remainder can be found in OA32's fMath includefile)

Posted on 2008-06-06 00:41:45 by Homer

Here's the first new method we have to look at.
This one calculates the velocity of a Point on a moving Body, and expresses it in Contact coordinates , where the X axis is the Collision Normal, and the Y and Z axes are 'PLANAR'.

``;Calculate the Velocity at CollisionPoint (in Contact coordinates);Returns Vec3 on fpuMethod CollisionHull.calculateLocalVelocity,uses esi,pTargetConfig,pmatContactToWorld, pvrelativeContactPositionLOCAL velocity:Vec3LOCAL accVelocity:Vec3LOCAL vout:Vec3	SetObject esi	    ; Work out the velocity of the contact point.    mov eax,pvrelativeContactPosition    mov edx,pTargetConfig    CrossProduct .configuration.AngularVelocity ,.Vec3    __AddToFloat3 velocity, .configuration.CMVelocity    __StowFloat3 velocity    ; Turn the velocity into contact-coordinates.    mov eax,pmatContactToWorld    Vec3_transMult .Mat33, velocity    __StowFloat3 vout    ; Convert "the change in velocity that was due to acceleration In This Frame" into contact-coordinates.    mov eax,pmatContactToWorld    Vec3_transMult .Mat33, .configuration.CMVelocityChange    __StowFloat3 accVelocity    ; We ignore any component of acceleration in the contact normal direction,     ; we are only interested in 'planar' acceleration    mov accVelocity.x,0    ; Add the planar velocities - if there's enough friction they will be removed during velocity resolution	__AddFloat3  vout, accVelocityMethodEnd``

Posted on 2008-06-06 22:35:31 by Homer
Nobody noticed the little bug in there so I won't mention EBX.

Here's your daily dose of fiber :)
Please excuse the 'debug' macros, as this code was recently debugged and implemented.

``;desiredDeltaVelocity is passed silently as input on FPUMethod CollisionHull.CalculateFrictionlessImpulse,uses esi,pvimpulseContactOut,pvContactNormal, pmatContactBasis,pvRelativeContactPositionA, pvRelativeContactPositionB, pconfigA,pBodyB,pconfigBLOCAL deltaVelWorld:Vec3LOCAL deltaVelocity:real8LOCAL desiredDeltaVelocity:real8	fstp desiredDeltaVelocity		;<--	DbgFloat desiredDeltaVelocity		SetObject esi	    ; Build a vector that shows the change in velocity (in world space)    ; for a UNIT impulse in the direction of the contact normal.           ;Calculate the Torque per Unit of Impulse    mov eax,pvContactNormal				    mov ecx,pvRelativeContactPositionA	    CrossProduct .Vec3, .Vec3 	    __StowFloat3 deltaVelWorld    			;torquePerUnitImpulse = CrossProduct(relativeContactPosition,contactNormal)        lea edx,deltaVelWorld    DbgVec3 edx,"torquePerUnitImpulse"        ;Use InverseInertiaTensor to convert Torque into rotationPerUnitImpulse    mov ecx,pconfigA    Vec3_mult_Mat33  .configuration.InverseWorldInertiaTensor, deltaVelWorld    __StowFloat3 deltaVelWorld     lea edx,deltaVelWorld    DbgVec3 edx,"rotationPerUnitImpulse"	;Find the Linear Velocity that is due to rotation and that is along the collision normal    mov ecx,pvRelativeContactPositionA		;velocityPerUnitImpulse = CrossProduct(rotationPerUnitImpulse, relativeContactPosition)	CrossProduct deltaVelWorld,  .Vec3	__StowFloat3 deltaVelWorld	    ; Work out the change in velocity in contact coordinates.    ;ImpulseDenominator = Body.OneOverMass + DotProduct(velocityPerUnitImpulse,CollisionNormal)		    mov eax,pvContactNormal    DotProduct deltaVelWorld, .Vec3    fadd .OneOverMass    fstp deltaVelocity    ; Check if we need the second body's data    .if pBodyB!=0        	int 3        ; Go through the same transformation sequence again	    mov eax,pvContactNormal	    mov ecx,pvRelativeContactPositionB	    CrossProduct .Vec3, .Vec3 	    __StowFloat3 deltaVelWorld   	    mov ecx,pconfigB	    Vec3_mult_Mat33  .configuration.InverseWorldInertiaTensor, deltaVelWorld		__StowFloat3 deltaVelWorld 	    mov ecx,pvRelativeContactPositionB		CrossProduct deltaVelWorld,  .Vec3		__StowFloat3 deltaVelWorld		        ; Add the change in velocity due to rotation        mov eax,pvContactNormal        DotProduct deltaVelWorld, .Vec3        fadd deltaVelocity        mov ecx,pBodyB        fadd .CollisionHull.OneOverMass        fstp deltaVelocity    .endif    ; Calculate the required size of the impulse    mov edx,pvimpulseContactOut    fld desiredDeltaVelocity    fdiv deltaVelocity    fstp .Vec3.x    mov  .Vec3.y,0     mov  .Vec3.z,0    DbgVec3 edx,"IMPULSE contact space"    MethodEnd``

The goal of this method is to produce a vector that contains an Impulse, given in ContactSpace coordinates. Remember that in Contact space, the X axis points along our Collision Normal, so X will contain the impulse along the collision normal, the one we really want. The Y and Z axes, which lay along the collision surface, are totally not used in a Frictionless collision, so we just set them to Zero.
You can take a look back to the Response method if you're interested in what happens to the returned impulse vector.

One method to go, and I believe we're done here.
Next we'll be looking at the FRICTIONAL impulse calculation.
Take the time to study the above method, as the Frictional version is quite similar, although a little more complex (we'll be using those Y and Z fields of the contactspace impulse).

After we've covered the last method, you'll realize that THIS method could have been implemented without any mention of 'contact space' at all.
The reason we've wasted our time converting to and from contact space is to give you the foundation knowledge to understand the Frictional version, where this WILL be required.
Posted on 2008-06-08 21:43:57 by Homer
Here's the final piece of the puzzle.
With this method in place, we have a fully-functional collision detection and response system that accurately simulates collisions between multiple rigid bodies comprised of materials with known physical properties.

``;desiredDeltaVelocity is passed silently as input on FPUMethod CollisionHull.CalculateFrictionImpulse,uses esi,pvimpulseContactOut, pmatContactBasis,pvContactVelocity,pvrelativeContactPositionA,pvrelativeContactPositionB,pconfigA,pBodyB,pconfigB, LP_frictionsLOCAL impulseToTorque:Mat33, tempmat:Mat33LOCAL deltaVelWorld:Mat33, deltaVelWorld2:Mat33LOCAL deltaVelocity:Mat33, impulseMatrix:Mat33LOCAL velKill:Vec3LOCAL inverseMass:real8, planarImpulse:real8LOCAL desiredDeltaVelocity:real8	fstp desiredDeltaVelocity		SetObject esi		fld .OneOverMass	fstp inverseMass    ; The equivalent of a cross product in matrices is multiplication    ; by a skew symmetric matrix - we build the matrix for converting    ; between linear and angular quantities.    mov eax,pvrelativeContactPositionA    Mat33_star impulseToTorque, .Vec3    DbgMat33 impulseToTorque, "impulseToTorque matrix"        ; Build the matrix to convert contact impulse --> change in velocity in world coordinates.    mov eax,pconfigA    invoke Mat33_multiply,addr tempmat, 		  addr impulseToTorque, addr .configuration.InverseWorldInertiaTensor    invoke Mat33_multiply,addr deltaVelWorld, 	  addr tempmat ,		addr impulseToTorque     Mat33_neg deltaVelWorld	DbgMat33 deltaVelWorld,"impulseToDeltaVelocity WORLDSPACE"    .if pBodyB!=0    	mov eax,pvrelativeContactPositionA    	Mat33_star impulseToTorque, .Vec3	    invoke Mat33_multiply,addr tempmat, 		  addr impulseToTorque, addr .configuration.InverseWorldInertiaTensor	    invoke Mat33_multiply,addr deltaVelWorld2, 	  addr tempmat ,		addr impulseToTorque 		Mat33_neg deltaVelWorld2        ; Add to the total delta velocity.        ;deltaVelWorld += deltaVelWorld2;        lea eax,deltaVelWorld        lea edx,deltaVelWorld2        Mat33_addto eax,edx        fld inverseMass        mov eax,pBodyB        fadd .CollisionHull.OneOverMass        fstp inverseMass    .endif        ; Do a change-of-basis to convert into contact coordinates.    mov eax,pmatContactBasis    Mat33_Transpose deltaVelocity, .Mat33	;Contact-To-World transform    	DbgMat33 deltaVelocity,"Transpose of ContactToWorld = WorldToContact"    invoke Mat33_multiply,addr tempmat,		  addr deltaVelocity,addr deltaVelWorld    invoke Mat33_multiply,addr deltaVelocity, addr tempmat, 	 pmatContactBasis    DbgMat33 deltaVelocity,"impulseToDeltaVelocity CONTACTSPACE"    ; Add in the linear velocity change    fld inverseMass    fadd deltaVelocity.m00    fstp deltaVelocity.m00    fld inverseMass    fadd deltaVelocity.m11    fstp deltaVelocity.m11    fld inverseMass    fadd deltaVelocity.m22	fstp deltaVelocity.m22	DbgMat33 deltaVelocity,"impulseMatrix"    ; Invert to get the impulse needed per unit velocity    invoke Mat33_inv,addr impulseMatrix ,addr deltaVelocity	DbgMat33 impulseMatrix,"impulseMatrix inverted"    ; Find the target velocities to kill    mov edx,pvContactVelocity    fld desiredDeltaVelocity    fstp velKill.x    fld .Vec3.y    fchs    fstp velKill.y    fld .Vec3.z    fchs    fstp velKill.z    ; Find the impulse to kill target velocities    Vec3_mult_Mat33 impulseMatrix, velKill    mov eax,pvimpulseContactOut    __StowFloat3 .Vec3    ; Check if theres enough PLANAR IMPULSE to overcome STATIC FRICTION    fld .Vec3.y    fmul st(0),st(0)    fld .Vec3.z    fmul st(0),st(0)    fadd    fsqrt    fstp planarImpulse    ;if (planarImpulse > impulseContact.x * STATIC friction)    fld .Vec3.x    mov edx,LP_frictions    fmul .MyTableElement._Static    fsub planarImpulse    fstpReg eax    .ifBitSet eax,BIT31        ; We need to use DYNAMIC friction    	mov eax,pvimpulseContactOut    	                fld .Vec3.y        fdiv planarImpulse        fstp .Vec3.y        fld .Vec3.z        fdiv planarImpulse        fstp .Vec3.z        ;impulseContact.x =  deltaVelocity.data +        ;    			    deltaVelocity.data*friction*impulseContact.y +        ;    				deltaVelocity.data*friction*impulseContact.z;        ;impulseContact.x = desiredDeltaVelocity / impulseContact.x;        ;impulseContact.y *= friction * impulseContact.x;        ;impulseContact.z *= friction * impulseContact.x;        fld desiredDeltaVelocity        fld deltaVelocity.m11        fmul .MyTableElement._Dynamic        fmul .Vec3.y         fadd deltaVelocity.m00        fld deltaVelocity.m22        fmul .MyTableElement._Dynamic        fmul .Vec3.z        fadd        fdiv        fstp .Vec3.x                fld  .Vec3.y        fmul .MyTableElement._Dynamic        fstp .Vec3.y        fld  .Vec3.z        fmul .MyTableElement._Dynamic        fstp .Vec3.z         .endifMethodEnd``

Posted on 2008-06-09 00:47:40 by Homer
Well, perhaps not the final piece.
Having implemented all this stuff in a demo, I found that my demo was 'mostly stable' - but that under certain conditions, bodies were able to penetrate the world bounding planes and cause the simulator to enter into an infinite loop.
The above system was supposed to make penetrations an impossible thing, but nonetheless under rare circumstances the system is failing, and the immediate problem is unexpected (inter)penetration of one or more bodies.
So, we're going to need to resolve these penetrations when they occur.
One way we could do that is simply to 'move' the bodies out of penetration along the collision normal, but that is not accurate.. it assumes that only LINEAR motion was responsible for the penetration, and indeed it causes error in our body's state since we didn't correct the velocity(s) and other stuff.. so what are we gonna do?

To resolve a penetration, we wish to 'go backwards in Time' until the Penetration becomes a Collision, resolve the collision, and then 'go forwards in Time' so that our Body is back in sync with the rest of the simulation. But our goal won't be to calculate how much negative Time we need...
The reason that we didn't detect this Penetration BEFORE it occurred (as we planned to) is that there's a relatively high velocity at the point of collision, so the Time involved is so small that it becomes problematic in terms of numerical precision... it escaped our temporal search , it must be really small, so we can't use Time.
We'll need to directly manipulate the physical states of our body(s), and work with the Velocities instead.

Let us assume that our system can at least tell us how deep the penetration is..
The penetration was caused by the motion of the body (or bodies).
We must recognize that this motion has two components - linear and angular.
Each component is partly responsible for the collision.
Can we determine HOW MUCH each is responsible? Yep, we sure can.
If we review the calculation of a collision response impulse, the first thing we did was calculate the total velocity at the collisionpoint due to both angular and linear motion.
We were able to determine how much velocity was due to Angular motion alone, yes?.
That is what we need now.
We need to find two Scaling values which represent the contributions of the linear and angular motions to the total penetration depth, and from them, calculate the translation and rotation required to bring that penetration depth to Zero.
I think we'll tackle this one in stages, its not huge, but its mathematically hairy - even though it only uses stuff we've already seen.

Posted on 2008-06-12 00:45:48 by Homer

We need to separate the angular and linear parts of this problem... we're trying to find the amount of the penetration due to EACH component, based on their 'relative intensities'.
Here's the most important method we'll require to solve this problem - it tells us the momentum of a given Body, along a given direction :)
If we wanted the Velocity, we could remove the Mass component (momemtum=velocity * mass).

``;**Returns float on FPU;Given a Collision Normal and a Contact Point, ;this method calculates the Linear Momentum;that is due ONLY to rotation, and which is in;the direction of the Collision Normal.Method CollisionHull.Calculate_AngularMomentum_along_Normal,uses esi,pvContactNormal, pvrelativeContactPosition, pconfigLOCAL angularInertiaWorld:Vec3	SetObject esi	;Calculate Torque per unit impulse at contact point	mov edx,pvContactNormal	mov eax,pvrelativeContactPosition	CrossProduct .Vec3, .Vec3			;torque = R cross N = change in angular momentum over time	__StowFloat3 angularInertiaWorld	mov edx,pconfig	;Convert Torque into Angular Momentum	Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor, angularInertiaWorld		;rotation = !I * torque = change in angular velocity		__StowFloat3 angularInertiaWorld	;Find the Linear part of that Angular motion	mov edx,pvrelativeContactPosition	CrossProduct angularInertiaWorld , .Vec3	__StowFloat3 angularInertiaWorld	;Find the Linear part that is along the Normal	mov edx,pvContactNormal	DotProduct angularInertiaWorld ,  .Vec3		MethodEnd``

This method returns a floating point value representing the given Body's momentum in the given direction... since its a direction, you MUST provide a Normal vector. We'll be passing the Collision Normal, but we COULD use this method to find the Body's momentum in ANY given direction.
It's telling us how much the Rotation Motion is responsible for total linear motion at a given Point on the Body.
Our solver will use this value to determine how much the Linear Motion is responsible for total linear motion at that Point, so we then know how much EACH is responsible.
From there we can solve our problem.
Posted on 2008-06-12 02:39:56 by Homer
The entrypoint for our penetration-solver is called Resolve_Penetration (oddly enough).
Here's the first few lines from that method.

``;Corrects a Penetration of two Bodies, but does not resolve the CollisionMethod CollisionHull.Resolve_Penetration, uses esi, pvContactNormal, pvrelativeContactPosition, pOtherBody, pSourceA,pTargetA,pSourceB,pTargetB,fPenetrationDepth:real8LOCAL angularInertiaWorld:Vec3LOCAL angularMoveA:real4, angularMoveB:real4LOCAL linearMoveA:real4, linearMoveB:real4LOCAL linearChangeA:Vec3, linearChangeB:Vec3LOCAL angularChangeA:Vec3, angularChangeB:Vec3LOCAL totalMoveA:real4,totalMoveB:real4local maxMagnitude:real8local totalInertia:real8local angularMomentumA:real8, angularMomentumB:real8LOCAL targetAngularDirectionA:Vec3,targetAngularDirectionB:Vec3LOCAL projection:Vec3LOCAL ftemp	SetObject esi    ; We need to find the momentum of each object in the direction    ; of the contact normal which is due to rotation only... 	OCall 						esi.Calculate_AngularMomentum_along_Normal, pvContactNormal, pvrelativeContactPosition, pSourceA 	fstp  angularMomentumA 	OCall pOtherBody::CollisionHull.Calculate_AngularMomentum_along_Normal, pvContactNormal, pvrelativeContactPosition, pSourceB 	fst  angularMomentumB 	 	;Now calculate the TOTAL momentum with respect to the collision normal 	fadd  angularMomentumA			 	 	fadd .OneOverMass					;Add the Linear component of momentum 	mov edx,pOtherBody 	fadd .CollisionHull.OneOverMass 	fstp  totalInertia``

You can see we're calculating the linear momentum along the collision normal which is solely due to the Rotation of each Body, as well as the Total linear momentum (again, along the collision normal). Don't get confused about the terms Inertia and Momentum ... the only real difference is that one is the signed opposite of the other... although not strictly correct, I use these terms interchangeably in my variable names.
Now let's see how we use those results to calculate that part of PenetrationDepth that is attributable to the Linear and to the Angular motions:
``	;Determine how much the Linear and Angular motions each contributed to the Penetration,	;and thereby, how much each must contribute to the Correction.    fld angularMomentumA    fdiv totalInertia    fmul fPenetrationDepth    fstp angularMoveA    ;    fld .OneOverMass    fdiv totalInertia    fmul fPenetrationDepth    fstp linearMoveA    ;    fld angularMomentumB    fdiv totalInertia    fmul fPenetrationDepth    fchs    fstp angularMoveB    ;    mov edx,pOtherBody    fld .CollisionHull.OneOverMass    fdiv totalInertia    fmul fPenetrationDepth    fchs    fstp linearMoveB``

Now we are in business  8)
From these two partial-PenDepth values we can calculate the changes in linear and angular position that are required to balance the equations such that the PenDepth 'disappears' (is equal to zero), and then adjust the linear and angular Velocities to suit.
We'll see how to do that next time :)

Posted on 2008-06-12 08:17:28 by Homer
I've added a little 'floor grid' to my demo so that we can better visualize the collisions of the demo body and the floor plane.

Here's a wishlist of other stuff I basically forgot about.
Chances are this list will grow!

-Improved support for Edge and Face collisions (reduced number of Point collision responses)
-Some basic gui menu controls for manipulating and resetting the demo simulation
-()

And now we return you to your feature program.
Posted on 2008-06-13 02:53:48 by Homer
OK lets get some closure on the current issue.
Here's the entire procedure, with small recent changes.
Take note of the 'angularLimit' threshold, which is used to limit the amount of Rotational correction, which places more load on the Linear correction.
This threshold is mostly for body/body collisions, where the bodies may vary greatly in mass, causing acute projections for the body with smaller mass.

``;Corrects a Penetration but does not resolve the CollisionMethod CollisionHull.Resolve_Penetration, uses esi, pvContactNormal, pvrelativeContactPositionA, pvrelativeContactPositionB,pOtherBody, pSourceA,pTargetA,pSourceB,pTargetB,fPenetrationDepth:real8LOCAL angularInertiaWorld:Vec3LOCAL angularMoveA:real4, angularMoveB:real4LOCAL linearMoveA:real4, linearMoveB:real4LOCAL linearChangeA:Vec3, linearChangeB:Vec3LOCAL angularChangeA:Vec3, angularChangeB:Vec3LOCAL totalMoveA:real4,totalMoveB:real4local maxMagnitude:real8local totalInertia:real8local angularMomentumA:real8, angularMomentumB:real8LOCAL targetAngularDirectionA:Vec3,targetAngularDirectionB:Vec3LOCAL projection:Vec3LOCAL ftemp	SetObject esi		DbgFloat fPenetrationDepth	    ; We need to find the momentum of each object in the direction    ; of the contact normal which is due to rotation only... 	OCall esi.Calculate_AngularMomentum_along_Normal, pvContactNormal, pvrelativeContactPositionA, pSourceA 	fabs 	DbgFST angularMomentumA 	.if pOtherBody!=0	 	OCall pOtherBody::CollisionHull.Calculate_AngularMomentum_along_Normal, pvContactNormal, pvrelativeContactPositionB, pSourceB	 	fabs	 	DbgFST angularMomentumB 	 	fadd		 	mov edx,pOtherBody 		fadd .CollisionHull.OneOverMass	.endif		 	;Now calculate the TOTAL momentum with respect to the collision normal 	fadd .OneOverMass					;Add the Linear component of momentum 	DbgFSTP totalInertia	;Determine how much the Linear and Angular motions each contributed to the Penetration,	;and thereby, how much each must contribute to the Correction.    fld angularMomentumA    fdiv totalInertia    fmul fPenetrationDepth    DbgFSTP angularMoveA    ;whatevers left over is the linear correction    fld fPenetrationDepth    fsub angularMoveA    DbgFSTP linearMoveA    ;	    ; To avoid angular projections that are too great     ; (when mass is large but inertia tensor is small)     ; we will limit the angular move.    mov edx,pvContactNormal    mov eax,pvrelativeContactPositionA    DotProduct .Vec3,.Vec3    DbgFSTP ftemp    __ScaleFloat3 .Vec3, ftemp    __AddToFloat3 projection, .Vec3        ; Use the small angle approximation for the sine of the angle     ; i.e. the magnitude would be sine(angularLimit) * projection.magnitude    ; but we approximate sine(angularLimit) to angularLimit    ; maxMagnitude = angularLimit * projection.magnitude    DotProduct projection,projection    fsqrt    fmul angularLimit    fst maxMagnitude    fchs    DbgFSTP ftemp            fMin angularMoveA, ftemp				;if (angularMoveA < -maxMagnitude)    fstpReg eax    .if eax==angularMoveA            	DbgWarning "FLOOR CLAMP on angular correction"        fld linearMoveA        fadd angularMoveA        DbgFSTP totalMoveA        fld maxMagnitude        fchs        DbgFSTP angularMoveA        fld totalMoveA        fsub angularMoveA        DbgFSTP linearMoveA        .else									    fMax angularMoveA, maxMagnitude  	;elseif (angularMoveA > maxMagnitude)	    fstpReg eax	    .if eax==angularMoveA 		    DbgWarning "CEILING CLAMP on angular correction"       			fld linearMoveA			fadd angularMoveA			DbgFSTP totalMoveA		     		    fld maxMagnitude		    DbgFSTP angularMoveA	        fld totalMoveA	        fsub angularMoveA	        DbgFSTP linearMoveA	    .endif 	.endif    ; We have the linear amount of movement required by ROTATING the rigid body (in angularMoveA).     ; We now need to calculate the desired rotation to achieve that.     .if (angularMoveA == 0 || angularMoveA==80000000h)        		; Easy case - no angular movement means no rotation.		fldz		fst angularChangeA.x		fst angularChangeA.y		fstp angularChangeA.z           	.else        		; Work out the direction we'd like to rotate in.		; targetAngularDirection = CrossProduct (relativeContactPositionA,contactNormal)		mov edx,pvrelativeContactPositionA		mov eax,pvContactNormal		CrossProduct .Vec3, .Vec3		;torque		__StowFloat3 targetAngularDirectionA	;Its a Direction because we used a Normal vector		; Convert the 'Torque' into 'change in Angular Velocity'		mov eax,pTargetA		Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor, targetAngularDirectionA	;= desired angular velocity		;Scale the Angular Velocity by our DESIRED angular velocity		fld angularMoveA		fdiv angularMomentumA		fmul st(3),st(0)		fmul st(2),st(0)		fmul				__StowFloat3 angularChangeA		lea edx,angularChangeA		DbgVec3 edx,"angularChangeA"		;The result is the change in angular momentum required to produce the desired change in angular velocity	.endif		; Linear change is easier - it is just the linear movement along the contact normal.	mov edx,pvContactNormal	__ScaleFloat3 .Vec3, linearMoveA	__StowFloat3 linearChangeA			lea edx,linearChangeA	DbgVec3 edx,"linearChangeA"	; Now we can apply the values we've calculated.	OCall esi.Apply_Correction, addr linearChangeA, addr angularChangeA, pSourceA, pTargetA	;------------------------------------------------------------------------------------------------------		 	;Now we go through almost the exact same sequence, for the Other Body 	.if pOtherBody!=0		fld angularMomentumB		fdiv totalInertia		fmul fPenetrationDepth		fchs		DbgFSTP angularMoveB		mov edx,pOtherBody		fld .CollisionHull.OneOverMass		fdiv totalInertia		fmul fPenetrationDepth		fchs		DbgFSTP linearMoveB	    mov edx,pvContactNormal	    mov eax,pvrelativeContactPositionB	    DotProduct .Vec3,.Vec3	    DbgFSTP ftemp	    __ScaleFloat3 .Vec3, ftemp	    __AddToFloat3 projection, .Vec3	    	    ; Use the small angle approximation for the sine of the angle 	    ; i.e. the magnitude would be sine(angularLimit) * projection.magnitude	    ; but we approximate sine(angularLimit) to angularLimit	    ; maxMagnitude = angularLimit * projection.magnitude	    DotProduct projection,projection	    fsqrt	    fmul angularLimit	    fst maxMagnitude	    fchs	    DbgFSTP ftemp	    fMin angularMoveB, ftemp	    fstpReg eax	    .if eax==angularMoveB	        fld linearMoveB	        fadd angularMoveB	        DbgFSTP totalMoveB	        fld maxMagnitude	        fchs	        DbgFSTP angularMoveB	        fld totalMoveB	        fsub angularMoveB	        DbgFSTP linearMoveB	    .else		    fMax angularMoveB, maxMagnitude		    fstpReg eax		    .if eax==angularMoveB				fld linearMoveB				fadd angularMoveB				DbgFSTP totalMoveB	     			    fld maxMagnitude			    DbgFSTP angularMoveB		        fld totalMoveB		        fsub angularMoveB		        DbgFSTP linearMoveB		    .endif	 	.endif 		    .if (angularMoveB == 0 || angularMoveB==80000000h)        			fldz			fst angularChangeB.x			fst angularChangeB.y			fstp angularChangeB.z           		.else        			mov edx,pvContactNormal			mov eax,pvrelativeContactPositionB			CrossProduct .Vec3, .Vec3			__StowFloat3 targetAngularDirectionB			mov eax,pTargetB			Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor, targetAngularDirectionB			fld angularMoveB			fdiv angularMomentumB			fmul st(3),st(0)			fmul st(2),st(0)			fmul			__StowFloat3 angularChangeB			lea edx,angularChangeB			DbgVec3 edx,"angularChangeB"		.endif		mov edx,pvContactNormal		__ScaleFloat3 .Vec3, linearMoveB		__StowFloat3 linearChangeB		OCall pOtherBody::CollisionHull.Apply_Correction, addr linearChangeB, addr angularChangeB, pSourceB, pTargetB   	.endif		        		;Now we COULD move Forwards in time if we were really worried about time-syncing our bodies...	;...OR we can simply live with the tiny time-glitch	MethodEnd``

And you'll also be interested in this:
``;Called by Resolve_PenetrationMethod CollisionHull.Apply_Correction,uses esi, pvLinearChange, pvAngularChange, pSource, pTargetLOCAL matTemp:Mat33LOCAL matTemp2:Mat33	SetObject esi		; Apply the change in position     	mov edx,pTarget	mov eax,pvLinearChange	__AddFloat3				.configuration.CMPosition,.Vec3	__StowFloat3 			.configuration.CMPosition         		; And the change in orientation 	mov edx,pTarget	mov eax,pvAngularChange	__AddFloat3 			.configuration.AngularMomentum,.Vec3	__StowFloat3 			.configuration.AngularMomentum		;We've changed the POSITION and ANGULAR MOMENTUM	;We need to fix up the Linear and Angular Velocities	;AngVelocity from AngMomentum as seen in INTEGRATE...	;Make a new Orientation matrix based on the Angular Change	mov eax,pvAngularChange	Mat33_star matTemp, .Vec3	mov edx,pSource	invoke Mat33_multiply,addr matTemp2, addr .configuration.Orientation, addr matTemp	mov eax,pSource	lea edx,.configuration.Orientation ;Source.Orientation	mov ebx,pTarget	lea ebx,.configuration.Orientation	;Target.Orientation	lea ecx,matTemp2						;Offset.Orientation	Mat33_Add  ebx, ecx, edx  	mov edx,pTarget	invoke OrthonormalizeOrientation,addr .configuration.Orientation		;Fix up the Inertia Tensor for the new Orientation	mov edx,pTarget	Mat33_SimilarityTransform .configuration.InverseWorldInertiaTensor, matTemp, .InverseBodyInertiaTensor, .configuration.Orientation	; Extract new Angular Velocity from the new Tensor and Momentum	; Target.AngularVelocity = Target.InverseWorldInertiaTensor * Target.AngularMomentum				mov eax,edx	Vec3_mult_Mat33 .configuration.InverseWorldInertiaTensor, .configuration.AngularMomentum	__StowFloat3 	.configuration.AngularVelocity		mov edx,pSource	;Ensure that our worldspace vertex array is up to date	OCall esi.TransformVertices,addr .configuration.CMPosition, addr .configuration.Orientation	;Scale the LINEAR Velocity as follows:	;CMVelocity and CMVelocityChange from CHANGE IN POSITION	;S						 = (new target.CMPosition - old target.CMPosition) / (old target.CMPosition-src.CMPosition)	;target.CMVelocity		 = target.CMVelocity * S	;target.CMVelocityChange = target.CMVelocityChange * S		MethodEnd``

Please note that I am NOT correcting the Linear Velocity - the last few lines of code are incomplete.
Also, note that the "angularLimit" threshold is completely arbitrary - you pretty much have to experiment to find a nice value but I find that 0.2 works nicely for me.

When you've had a chance to look that code over, I'll post an updated version of the simulator's CheckCollisions method, which calls this new penetration-correction stuff.
Posted on 2008-06-16 02:21:24 by Homer
Theres a reason why I did not bother to adjust the linear velocity.
And this will lead to a separate discussion of the solver, where we will see that positional correction is NOT NEEDED.
But I digress.

For now, please note that when Penetration occurs, it is ALWAYS SHALLOW - its limited to the TIME EPSILON of the "FindExactCollisionTime" method... so when we advanced the system into the penetration state, we did so by the smallest possible margin.
The positional error, for which we corrected, was small.
The velocity correction is small too.
So small.
Important?
Conservation of energy.

We corrected the position, and not the velocity.

What if we DIDNT correct the position - and used the penetration point for collision resolution?
The velocity issue would disappear - and we'd have a little angular error instead..
Would you not say that is more acceptable than wearing, or having to correct, the error in TOTAL momentum?
Posted on 2008-06-16 04:59:07 by Homer