float mB = vc.invMassB;
float iA = vc.invIA;
float iB = vc.invIB;
int pointCount = vc.pointCount;
Vec2 vA = m_velocities[indexA].v;
float wA = m_velocities[indexA].w;
Vec2 vB = m_velocities[indexB].v;
float wB = m_velocities[indexB].w;
Vec2 normal = vc.normal;
tangent.x = 1.0f * vc.normal.y;
tangent.y = -1.0f * vc.normal.x;
final float friction = vc.friction;
assert (pointCount == 1 || pointCount == 2);
// Solve tangent constraints
for (int j = 0; j < pointCount; ++j) {
final VelocityConstraintPoint vcp = vc.points[j];
final Vec2 a = vcp.rA;
float dvx = -wB * vcp.rB.y + vB.x - vA.x + wA * a.y;
float dvy = wB * vcp.rB.x + vB.y - vA.y - wA * a.x;
// Compute tangent force
final float vt = dvx * tangent.x + dvy * tangent.y - vc.tangentSpeed;
float lambda = vcp.tangentMass * (-vt);
// Clamp the accumulated force
final float maxFriction = friction * vcp.normalImpulse;
final float newImpulse =
MathUtils.clamp(vcp.tangentImpulse + lambda, -maxFriction, maxFriction);
lambda = newImpulse - vcp.tangentImpulse;
vcp.tangentImpulse = newImpulse;
// Apply contact impulse
// Vec2 P = lambda * tangent;
final float Px = tangent.x * lambda;
final float Py = tangent.y * lambda;
// vA -= invMassA * P;
vA.x -= Px * mA;
vA.y -= Py * mA;
wA -= iA * (vcp.rA.x * Py - vcp.rA.y * Px);
// vB += invMassB * P;
vB.x += Px * mB;
vB.y += Py * mB;
wB += iB * (vcp.rB.x * Py - vcp.rB.y * Px);
}
// Solve normal constraints
if (vc.pointCount == 1) {
final VelocityConstraintPoint vcp = vc.points[0];
// Relative velocity at contact
// Vec2 dv = vB + Cross(wB, vcp.rB) - vA - Cross(wA, vcp.rA);
float dvx = -wB * vcp.rB.y + vB.x - vA.x + wA * vcp.rA.y;
float dvy = wB * vcp.rB.x + vB.y - vA.y - wA * vcp.rA.x;
// Compute normal impulse
final float vn = dvx * normal.x + dvy * normal.y;
float lambda = -vcp.normalMass * (vn - vcp.velocityBias);
// Clamp the accumulated impulse
float a = vcp.normalImpulse + lambda;
final float newImpulse = (a > 0.0f ? a : 0.0f);
lambda = newImpulse - vcp.normalImpulse;
vcp.normalImpulse = newImpulse;
// Apply contact impulse
float Px = normal.x * lambda;
float Py = normal.y * lambda;
// vA -= invMassA * P;
vA.x -= Px * mA;
vA.y -= Py * mA;
wA -= iA * (vcp.rA.x * Py - vcp.rA.y * Px);
// vB += invMassB * P;
vB.x += Px * mB;
vB.y += Py * mB;
wB += iB * (vcp.rB.x * Py - vcp.rB.y * Px);
} else {
// Block solver developed in collaboration with Dirk Gregorius (back in 01/07 on
// Box2D_Lite).
// Build the mini LCP for this contact patch
//
// vn = A * x + b, vn >= 0, , vn >= 0, x >= 0 and vn_i * x_i = 0 with i = 1..2
//
// A = J * W * JT and J = ( -n, -r1 x n, n, r2 x n )
// b = vn_0 - velocityBias
//
// The system is solved using the "Total enumeration method" (s. Murty). The complementary
// constraint vn_i * x_i
// implies that we must have in any solution either vn_i = 0 or x_i = 0. So for the 2D
// contact problem the cases
// vn1 = 0 and vn2 = 0, x1 = 0 and x2 = 0, x1 = 0 and vn2 = 0, x2 = 0 and vn1 = 0 need to be
// tested. The first valid
// solution that satisfies the problem is chosen.
//
// In order to account of the accumulated impulse 'a' (because of the iterative nature of
// the solver which only requires
// that the accumulated impulse is clamped and not the incremental impulse) we change the
// impulse variable (x_i).
//
// Substitute:
//
// x = a + d
//
// a := old total impulse
// x := new total impulse
// d := incremental impulse
//
// For the current iteration we extend the formula for the incremental impulse
// to compute the new total impulse:
//
// vn = A * d + b
// = A * (x - a) + b
// = A * x + b - A * a
// = A * x + b'
// b' = b - A * a;
final VelocityConstraintPoint cp1 = vc.points[0];
final VelocityConstraintPoint cp2 = vc.points[1];
a.x = cp1.normalImpulse;
a.y = cp2.normalImpulse;
assert (a.x >= 0.0f && a.y >= 0.0f);
// Relative velocity at contact
// Vec2 dv1 = vB + Cross(wB, cp1.rB) - vA - Cross(wA, cp1.rA);
dv1.x = -wB * cp1.rB.y + vB.x - vA.x + wA * cp1.rA.y;
dv1.y = wB * cp1.rB.x + vB.y - vA.y - wA * cp1.rA.x;
// Vec2 dv2 = vB + Cross(wB, cp2.rB) - vA - Cross(wA, cp2.rA);
dv2.x = -wB * cp2.rB.y + vB.x - vA.x + wA * cp2.rA.y;
dv2.y = wB * cp2.rB.x + vB.y - vA.y - wA * cp2.rA.x;
// Compute normal velocity
float vn1 = dv1.x * normal.x + dv1.y * normal.y;
float vn2 = dv2.x * normal.x + dv2.y * normal.y;
b.x = vn1 - cp1.velocityBias;
b.y = vn2 - cp2.velocityBias;
// System.out.println("b is " + b.x + "," + b.y);
// Compute b'
Mat22 R = vc.K;
b.x -= R.ex.x * a.x + R.ey.x * a.y;
b.y -= R.ex.y * a.x + R.ey.y * a.y;
// System.out.println("b' is " + b.x + "," + b.y);
// final float k_errorTol = 1e-3f;
// B2_NOT_USED(k_errorTol);
for (;;) {
//
// Case 1: vn = 0
//
// 0 = A * x' + b'
//
// Solve for x':
//
// x' = - inv(A) * b'
//
// Vec2 x = - Mul(c.normalMass, b);
Mat22.mulToOutUnsafe(vc.normalMass, b, x);
x.x *= -1;
x.y *= -1;
if (x.x >= 0.0f && x.y >= 0.0f) {
// System.out.println("case 1");
// Get the incremental impulse
// Vec2 d = x - a;
d.set(x).subLocal(a);
// Apply incremental impulse
// Vec2 P1 = d.x * normal;
// Vec2 P2 = d.y * normal;
P1.set(normal).mulLocal(d.x);
P2.set(normal).mulLocal(d.y);
/*
* vA -= invMassA * (P1 + P2); wA -= invIA * (Cross(cp1.rA, P1) + Cross(cp2.rA, P2));
*
* vB += invMassB * (P1 + P2); wB += invIB * (Cross(cp1.rB, P1) + Cross(cp2.rB, P2));
*/
temp1.set(P1).addLocal(P2);
temp2.set(temp1).mulLocal(mA);
vA.subLocal(temp2);
temp2.set(temp1).mulLocal(mB);
vB.addLocal(temp2);
wA -= iA * (Vec2.cross(cp1.rA, P1) + Vec2.cross(cp2.rA, P2));
wB += iB * (Vec2.cross(cp1.rB, P1) + Vec2.cross(cp2.rB, P2));
// Accumulate
cp1.normalImpulse = x.x;
cp2.normalImpulse = x.y;
/*
* #if B2_DEBUG_SOLVER == 1 // Postconditions dv1 = vB + Cross(wB, cp1.rB) - vA -
* Cross(wA, cp1.rA); dv2 = vB + Cross(wB, cp2.rB) - vA - Cross(wA, cp2.rA);
*
* // Compute normal velocity vn1 = Dot(dv1, normal); vn2 = Dot(dv2, normal);
*
* assert(Abs(vn1 - cp1.velocityBias) < k_errorTol); assert(Abs(vn2 - cp2.velocityBias)
* < k_errorTol); #endif
*/
if (DEBUG_SOLVER) {
// Postconditions
Vec2 dv1 =
vB.add(Vec2.cross(wB, cp1.rB).subLocal(vA).subLocal(Vec2.cross(wA, cp1.rA)));
Vec2 dv2 =
vB.add(Vec2.cross(wB, cp2.rB).subLocal(vA).subLocal(Vec2.cross(wA, cp2.rA)));
// Compute normal velocity
vn1 = Vec2.dot(dv1, normal);
vn2 = Vec2.dot(dv2, normal);
assert (MathUtils.abs(vn1 - cp1.velocityBias) < k_errorTol);
assert (MathUtils.abs(vn2 - cp2.velocityBias) < k_errorTol);
}
break;
}
//
// Case 2: vn1 = 0 and x2 = 0
//
// 0 = a11 * x1' + a12 * 0 + b1'
// vn2 = a21 * x1' + a22 * 0 + '
//
x.x = -cp1.normalMass * b.x;
x.y = 0.0f;
vn1 = 0.0f;
vn2 = vc.K.ex.y * x.x + b.y;
if (x.x >= 0.0f && vn2 >= 0.0f) {
// System.out.println("case 2");
// Get the incremental impulse
d.set(x).subLocal(a);
// Apply incremental impulse
// Vec2 P1 = d.x * normal;
// Vec2 P2 = d.y * normal;
P1.set(normal).mulLocal(d.x);
P2.set(normal).mulLocal(d.y);
/*
* Vec2 P1 = d.x * normal; Vec2 P2 = d.y * normal; vA -= invMassA * (P1 + P2); wA -=
* invIA * (Cross(cp1.rA, P1) + Cross(cp2.rA, P2));
*
* vB += invMassB * (P1 + P2); wB += invIB * (Cross(cp1.rB, P1) + Cross(cp2.rB, P2));
*/
temp1.set(P1).addLocal(P2);
temp2.set(temp1).mulLocal(mA);
vA.subLocal(temp2);
temp2.set(temp1).mulLocal(mB);
vB.addLocal(temp2);
wA -= iA * (Vec2.cross(cp1.rA, P1) + Vec2.cross(cp2.rA, P2));
wB += iB * (Vec2.cross(cp1.rB, P1) + Vec2.cross(cp2.rB, P2));
// Accumulate
cp1.normalImpulse = x.x;
cp2.normalImpulse = x.y;
/*
* #if B2_DEBUG_SOLVER == 1 // Postconditions dv1 = vB + Cross(wB, cp1.rB) - vA -
* Cross(wA, cp1.rA);
*
* // Compute normal velocity vn1 = Dot(dv1, normal);
*
* assert(Abs(vn1 - cp1.velocityBias) < k_errorTol); #endif
*/
if (DEBUG_SOLVER) {
// Postconditions
Vec2 dv1 =
vB.add(Vec2.cross(wB, cp1.rB).subLocal(vA).subLocal(Vec2.cross(wA, cp1.rA)));
// Compute normal velocity
vn1 = Vec2.dot(dv1, normal);
assert (MathUtils.abs(vn1 - cp1.velocityBias) < k_errorTol);
}
break;
}
//
// Case 3: wB = 0 and x1 = 0
//
// vn1 = a11 * 0 + a12 * x2' + b1'
// 0 = a21 * 0 + a22 * x2' + '
//
x.x = 0.0f;
x.y = -cp2.normalMass * b.y;
vn1 = vc.K.ey.x * x.y + b.x;
vn2 = 0.0f;
if (x.y >= 0.0f && vn1 >= 0.0f) {
// System.out.println("case 3");
// Resubstitute for the incremental impulse
d.set(x).subLocal(a);
// Apply incremental impulse
/*
* Vec2 P1 = d.x * normal; Vec2 P2 = d.y * normal; vA -= invMassA * (P1 + P2); wA -=
* invIA * (Cross(cp1.rA, P1) + Cross(cp2.rA, P2));
*
* vB += invMassB * (P1 + P2); wB += invIB * (Cross(cp1.rB, P1) + Cross(cp2.rB, P2));
*/
P1.set(normal).mulLocal(d.x);
P2.set(normal).mulLocal(d.y);
temp1.set(P1).addLocal(P2);
temp2.set(temp1).mulLocal(mA);
vA.subLocal(temp2);
temp2.set(temp1).mulLocal(mB);
vB.addLocal(temp2);
wA -= iA * (Vec2.cross(cp1.rA, P1) + Vec2.cross(cp2.rA, P2));
wB += iB * (Vec2.cross(cp1.rB, P1) + Vec2.cross(cp2.rB, P2));
// Accumulate
cp1.normalImpulse = x.x;
cp2.normalImpulse = x.y;
/*
* #if B2_DEBUG_SOLVER == 1 // Postconditions dv2 = vB + Cross(wB, cp2.rB) - vA -
* Cross(wA, cp2.rA);
*
* // Compute normal velocity vn2 = Dot(dv2, normal);
*
* assert(Abs(vn2 - cp2.velocityBias) < k_errorTol); #endif
*/
if (DEBUG_SOLVER) {
// Postconditions
Vec2 dv2 =
vB.add(Vec2.cross(wB, cp2.rB).subLocal(vA).subLocal(Vec2.cross(wA, cp2.rA)));
// Compute normal velocity
vn2 = Vec2.dot(dv2, normal);
assert (MathUtils.abs(vn2 - cp2.velocityBias) < k_errorTol);