Dirk, sorry for waking up this post from the tumb, but when two polygons are coplanar, did you use some test to avoid the issue of clipping only to a triangle (which will makes cubes bounce by the solver)?

I'm asking just for precaution, because in the coplanar case I'm getting 4 points in the case of the box-ground contact, and all the negative contact points penetrations (below the reference face) are the same, which means that is perfectly correct.

The specific problem is that I can't even stack more than 3 boxes with 10 iterations using SAT (and the SI solver), which looks incorrect. Using GJK and EPA with an incremental manifold seems to be more stable tough is more numerically instable.

I say this because since SAT is one shot, I think there is no contact persistency, only the feature caching correct? For me, in the E-E case, the contact ID is the index of the two edges, and on the FF they're the face indices. The == operators just compares each type and indices.

(Off-topic) This is how I'm warm-starting the contact constraint if there is a match in the SAT features.

http://pastebin.com/8zCMteLb

Do you mean coplanar faces (polygons) in a hull? E.g. each face of a box is a triangle and two (coplanar) triangles build a side of your box? In this case: Yes, I merge coplanar faces in my hull builder. It is very important and I spend huge effort on this. Actually I gave a talk about this as well :):

"Implementing Quickhull" (GDC 2014)

http://www.valvesoftware.com/company/publications.html

Otherwise the code looks code! Maybe give me another example if I misunderstood your problem!

HTH,

-Dirk

I've checked your implementation on QuickHull and definately will read this .

The problem is that the solver is quite instable.

I think when the penetration of each contact point gathered with the SAT gets solved, the other one bounces, and it will keep like this sequentially. I've debugged for an entire day, and could not find the problem.

I'm pretty sure that both Jacobians and Lambdas are correct (matched with Erin slides). I'll post here:

#include "CContact.h"
#include <assert.h>

void CContact::ComputeConstants(float _fInvDt) {
	float InvMassSum = Body0->InvMass() + Body1->InvMass();

	for (unsigned int I = 0; I < Manifold.TotalPoints; ++I) {
		CManifold::CONTACT_POINT& Point = Manifold.ContactPoints[I];
		
		//Calculate M^-1JT

		// Normal Mass
		{
			CVector3 R0xN = Point.R0.Cross(Manifold.Normal);
			CVector3 R1xN = Point.R1.Cross(Manifold.Normal);
			float NormalMass = InvMassSum;
			if (!Body0->IsStatic()) {
				NormalMass += R0xN.Dot(Body0->m_mWorldInvInertia * R0xN);
			}
			if (!Body1->IsStatic()) {
				NormalMass += R1xN.Dot(Body1->m_mWorldInvInertia * R1xN);
			}

			Point.NormalMass = 1.0f / NormalMass;
		}

		// Tangent Mass
		for (unsigned int J = 0; J < 2; ++J) {
			CVector3 R0xT = Point.R0.Cross(Manifold.Tangents[J]);
			CVector3 R1xT = Point.R1.Cross(Manifold.Tangents[J]);
			float TangentMass = InvMassSum;
			if (!Body0->IsStatic()) {
				TangentMass += R0xT.Dot(Body0->m_mWorldInvInertia * R0xT);
			}
			if (!Body1->IsStatic()) {
				TangentMass += R1xT.Dot(Body1->m_mWorldInvInertia * R1xT);
			}
			Point.TangentMass[J] = 1.0f / TangentMass;
		}

		Point.Bias = -BAUMGARTE * _fInvDt * CMathLib::Min(0.0f, Point.Penetration + PENETRATION_SLOP);
		
		//Warm Start (this is called per step, not per solve tentative!)
		CVector3 vP = Manifold.Normal * Point.LambdaNormal;
		for (unsigned int J = 0; J < 2; ++J) {
			vP += Manifold.Tangents[J] * Point.LambdaTangent[J];
		}

		if (!Body0->IsStatic()) {
			Body0->m_vLinVel -= vP * Body0->InvMass();
			Body0->m_vAngVel -= Body0->m_mWorldInvInertia * Point.R0.Cross(vP);
		}

		if (!Body1->IsStatic()) {
			Body1->m_vLinVel += vP * Body1->InvMass();
			Body1->m_vAngVel += Body1->m_mWorldInvInertia * Point.R1.Cross(vP);
		}
	}
}

void CContact::Solve() {
	for (unsigned int I = 0; I < Manifold.TotalPoints; ++I) {
		CManifold::CONTACT_POINT& Point = Manifold.ContactPoints[I];

		// Tangent
		{
			for (unsigned int J = 0; J < 2; ++J) {
				CVector3 vTangent = Manifold.Tangents[J];

				CVector3 vRVCP = Body1->m_vLinVel + Body1->m_vAngVel.Cross(Point.R1) - Body0->m_vLinVel - Body0->m_vAngVel.Cross(Point.R0);
				float fRVN = vRVCP.Dot(vTangent);
				float fLambda = Point.TangentMass[J] * -fRVN;

				float fOldLambda = Point.LambdaTangent[J];
				float fMaxLambda = 0.5f * Point.LambdaNormal;
				Point.LambdaTangent[J] = CMathLib::Clamp(fOldLambda + fLambda, -fMaxLambda, fMaxLambda);
				fLambda = Point.LambdaTangent[J] - fOldLambda;

				CVector3 vPT = vTangent * fLambda;

				if (!Body0->IsStatic()) {
					Body0->m_vLinVel -= vPT * Body0->InvMass();
					Body0->m_vAngVel -= Body0->m_mWorldInvInertia * Point.R0.Cross(vPT);
				}

				if (!Body1->IsStatic()) {
					Body1->m_vLinVel += vPT * Body1->InvMass();
					Body1->m_vAngVel += Body1->m_mWorldInvInertia * Point.R1.Cross(vPT);
				}
			}
		}

		{
			// Normal

			CVector3 vRVCP = Body1->m_vLinVel + Body1->m_vAngVel.Cross(Point.R1) - Body0->m_vLinVel - Body0->m_vAngVel.Cross(Point.R0);
			float fRVN = vRVCP.Dot(Manifold.Normal);
			float fLambda = Point.NormalMass * (-fRVN + Point.Bias);

			float fOldLambda = Point.LambdaNormal;
			Point.LambdaNormal = CMathLib::Max(fOldLambda + fLambda, 0.0f);
			fLambda = Point.LambdaNormal - fOldLambda;

			CVector3 vPN = Manifold.Normal * fLambda;

			if (!Body0->IsStatic()) {
				Body0->m_vLinVel -= vPN * Body0->InvMass();
				Body0->m_vAngVel -= Body0->m_mWorldInvInertia * Point.R0.Cross(vPN);
			}

			if (!Body1->IsStatic()) {
				Body1->m_vLinVel += vPN * Body1->InvMass();
				Body1->m_vAngVel += Body1->m_mWorldInvInertia * Point.R1.Cross(vPN);
			}
		}
	}
}

Baumgarte is 0.1 (10%), and the slop is 0.01 cm, but no good results. I use two friction directions, calculating a basis when the rel. velocity at the contact point vanishes.

Do you have any idea of what could be a common mistake on this?

This looks correct! BAUMGARTE should be 0.1 - 0.2.

Here are some more debug tips:

- Disable penetration recovery (maybe you compute the penetration wrong). Just set BAUMGARTE to zero. Simple test.

- Color your contact points. E.g. green for points you were able to wamstart and red for points that you cannot match to old points. If you have a blinking 'Christmas' tree you know that this is the problem

How do you handle more than four points? Check that your reduction code is working. You need to project all contact points onto the reference plane for culling.

Again, I think it is very likely some simple mistake and not a fundamental flaw in your math or implementation. I say this, because I always think that and then it is often just a stupid typo.

I like that you solve the friction constraints before non-penetrations since the friction constraints can violate the non-penetration. Nevertheless, you get better results from my experience the other way around. Until you don't have CCD and continuous physics you don't need to solve friction first.

How do you handle more than four points? Check that your reduction code is working. You need to project all contact points onto the reference plane for culling.

I'm using a maximum of 8 points, no hull approximation yet. The penetration depth of each contact point is the projected clipped verts on the reference face (those that are negative, below the ref. face).

I like that you solve the friction constraints before non-penetrations since the friction constraints can violate the non-penetration. Nevertheless, you get better results from my experience the other way around. Until you don't have CCD and continuous physics you don't need to solve friction first.

This makes completely sense.

Just disabled the penetration depth, and looks like that's the problem. Only with 2 iterations I have now a almost stable stack of 5 boxes.

Now, each contact point gets an penetration depth, and looking into SAT, all points will have the same separation depth for a box on plane scenario. I looked here and the separation for each contact point is the same separation value of the face-face query! Which is correct, since all points lies on the same plane, and is not different from the query itself.

Shoudn't the penetration depth be averaged for a single manifold? For me, this makes no sense, since is the same of solving each contact point, but maybe this could be considerable.

Make sure that you use the normal correctly. Usually you flip after all computations are done! You use the reference face normal to build the manifold. The flipping is only for the solver!

Here what you need to do:

- Get the clip polygon in world space

- Get the clip planes in world space

- Clip the polygon against each plane

- Get the reference plane in world space

- For each vertex of the clipped polygon compute the distance to the reference plane: float Distance = dot( Plane.Normal, Polygon.Vertex ) - Plane.Offset

- Keep points below reference plane ( e.g. Distance <= 0 - there is no need to test against FLT_EPSILON here since the allowed penetration handles this)

- Make sure your plane normals are normalized. Otherwise the distance will be wrong!

Here is how I transform a plane (not sure if yours is equivalent):

RN_FORCEINLINE rnPlane operator*( const rnTransform& Transform, const rnPlane& Plane )
{
    rnVector3 Normal = Transform.Rotation * Plane.Normal;  // Transform.Rotation is a 3x3 matrix
    return rnPlane( Normal, Plane.Offset + rnDot( Normal, Transform.Translation ) );  // The Plane constructor takes a normal and offset 
}

All the computations are in world space.

I made a video testing the solver, with 10 iterations at 60Hz:

As you see, this is not very stable. One thing I've changed was calculate tangents for each contact point (not per manifold, as the normal). Even this way it seems to not converge. The points are kind of blinking because the one shot manifold generation way, so, since that there is some hidden problem, it will compute different ones each step.

Disabling the angular impulses looks more correct:

According with Erin Iterative Dynamics this doesn't looks correct. At 10 iterations this should be stable, shouldn't it?

Looks like I need to use the old sphere-plane test, and step through this .

Not, sure! Can you stack 5 boxes? Another test is if the solver converges. E.g. if you run 20 or 100 iterations is the box stack stable. The iteration count is a function of the stacking height. You get some amortization due to warmstarting, but 10 iterations and 10 boxes is at the border. I would also add another test. E.g. a pyramid. They are usually stable so you should be able to do a pyramid with height 10.

I think Rubikon can do 10 and 10, but there are differences to your solver. E.g. I use position projection and not Baumgarte. I also use a different friction model which improves convergence as it decouples the friction and non-penetration constraints a bit. I also don't loose warmstart information if you slide to another edge. The stack you show is quite difficult if you are not careful. I would try something with 45 degrees rotated boxes since you will not loose warmstart information this way.

With 100 iterations a stack of 10 boxes bounces a little bit, but they do not fall on the ground.

In my inspection Lambda is approaching to zero or < 0 (since negative ones are allowed).

The solver isn't losing warm-start information since all contact IDs are matching (since there are only face-face contacts in this example).

The non-penetration lambdas must be greater zero, This can be checked easily for a plane. If this fails you have a bug. Friction lambdas can be negative though.