Contact Joint Summary

April 11, 2006

There are three DOF in the contact joint. A generic constraint equation is written as

J v = c (1)

The first DOF is in the normal direction. The LHS of Eq. (1) is

where r1 = cp – p1 and r2 = cp – p2. The RHS of Eq. (1) is

c = erp/h⋅depth

But c is limited by maxvel, so

if (c > maxvel) c = maxvel

Since cfm is associated with the constraint DOF, a specific cfm can be entered here. There are 3 possible cfm, which are set up by soft_cfm, slip1, and slip2, respectively.

cfm(1) = soft_cfm

The impact has a direct effect on c

pick the larger value between contact and impact

if(c_i > c) c = c_i

The lower and upper bounds are

lo(1) = 0
hi(1) = infinity

The second DOF is the first tangent direction. Let t1 and t2 be the tangent directions to normal n

c = motion1 (so motion1 and motion2 are velocities)
lo(2) = -mu (if mu=0, the friction force is practically set to zero)
hi(2) = mu
cfm(2) = slip1

The LHS of the constraint equation in the second tangent direction is

The RHS of the constraint equation is

c = motion2
lo(3) = -mu2
hi(3) = mu2
cfm(3) = slip2

Dual-space Expansion for Estimating Penetration depth

April 10, 2006

So I have been implementing Dual-space Expansion for Estimating Penetration depth (DEEP) for a while now. In the past few days I put it to work and it sucked. The amount of frustration from doing Contact Manifold & Signed Distance Map and getting Trouble With The Projection kept piling up to this point, and this DEEP blow just pushed the frustration to an astronomical level. This morning I finally had a revelation that I might just need to go back to make box-box collision work instead of finding new algorithms that add more unmanageable complexity without any realistic benefit.

But here is the twist of the story. After all, the time I have invested in this effort will not be wasted. An innocent bug was accidentally found and DEEP w/ contact manifold code finally works!! It turns out the problem is when only four Contact are used for every contact that is happenning. Somehow I still cannot duplicate the performance of test_crash, but with careful tunning of the contact parameters, the collision of the boxes will show some physically believable behaviors. To me it has been a major progress considering how many frustrating hours I have spent on coding this. For example, erp=0.8 never works and neither without using maxvel to restraint the correcting contact velocity.

Now a new problem emerges: DEEPGeom will break HashSpace …

Trouble With The Projection

March 13, 2006

The Gauss map on the right shows the overlay of a feature pair FV, where the intersection can be easily missed if this pair is projected onto a plane. There is definitely a chance I misunderstand how the projection works. Otherwise I have to workaround without the projection.

After some web searching, I learned that I was using orthographic projection, while I should use is gnomonic projection. The advantage for gnomonic projection is that the great circles will always be mapped to straight lines.

Contact Manifold & Signed Distance Map

January 19, 2006

A quick review of my work for the past two months:

Contact Manifold

As I have been struggling with the box-box collision detection, I started to get fascinated by the idea of constructing a contact manifold, which contains a number of selected contact points. My first attempt is to adapt Bullet code from Erwin Coumans, who claims to successfully build correct contact manifolds. Eventually I did not succeed to make the adapted Java code run correctly. Thus I stopped working on the contact manifold idea with frustration.

Signed Distance Map

Then I bought a book from Christer Ericson with the hope of getting some working box-box collision and/or contact manifold codes. It turned out there was nothing in the book I could use immediately. In the meantime, the idea of signed distance map came to me but I don’t understand how it really works. I tried to resolve this problem by adapting PhysicsEngine. So far I have not yet made the adapted Java code work and inevitably the frustration grew.

Intersection and Penetration Depth

November 28, 2005

So it has been two weeks since I started translating C++ SOLID code into Java code. Over the Thanksgiving weekend I intensely debugged my code such that it now runs reasonablely. One notable bug I have is that within the subroutine I cannot simply reassign the argument variable to a new object.

void test() {
  Vector3d v = new Vector3d();
  System.out.println(v.x + " " + v.y + " " + v.z); // still (0.0, 0.0, 0.0)

void foo(Vector3d v) {
  Vector3d a = new Vector3d(1.0, 2.0, 3.0);
  v = a; // this has no effect on v

When checking the results from, it is found that the calculation of penetration depth is soemewhat unstable and slow. My to-do list will just have to expand:

  • To understand the GJK algorithm more deeply.
  • To thoroughly test my current implementation.

I think it is quite reasonable to assume only one contact point when two convex objects contact or intersect/impact. Thus the penetration depth and the impact normal can be conveniently calculated from the GJK subroutines. I implemented this idea by modifying today and the test results are satisfactory so far with minor observed intersections. The results certainly look better than the ODE implementation, which I don’t fully understand either. Let’s assume now the convex problems are resolved. Then I still need to investigate the complex/trimesh problems. More items on to-do list:

  • Check SOLID complex package and try to implement a corresponding one for my trimesh. Can RAPID code be integrated here?
  • Try two packaging and delivery methods – applet and web start.

Crank Model Now Works

November 1, 2005

Last Friday I expanded the piston assembly model to a full cranktrain model. The original Smith method worked immediately, but Baraff method did not. Today I finally have the bug fixed 🙂 My confidence in the Baraff method rises as the results match pretty well with those of the Smith method. Now the performace of the Baraff method is still a question mark. I don’t feel the program running any faster using the new method.

The bug fix is within buildTree() in

  TreeNode buildTree(Body b) {
    // get b = the next enabled, untagged body, and tag it
    if (b.tag!=0 || (b.flags & Body.Disabled)!=0) return null;
    b.tag = 1;

    LinkedList stack = new LinkedList();
    // tag all bodies and joints starting from b
    TreeNode root = new TreeNode(b, null);
    TreeNode node = root, childnode;
    while (stack.size() > 0) {
      node = (TreeNode)stack.removeFirst();  // pop body node off stack
      b = node.body;

      // traverse and tag all body's joints,
      // add untagged connected bodies to stack
      Iterator j = b.nodes.iterator();
      while(j.hasNext()) {
        JointNode jn = (JointNode);
        if(jn.joint.tag==0) {
          jn.joint.tag = 1;
          if(jn.body!=null && jn.body.tag==0) {
            childnode = new TreeNode(jn.joint, node);
            childnode = new TreeNode(jn.body, childnode);
            jn.body.tag = 1;
            aux.add(n.joint); // put boundary joints into auxiliary joint list
    return root;

Constraint Force Mixing

October 26, 2005

Today I learned the magic of CFM (constraint force mixing). If CFM is not introduced, the matrix A will somehow manages to get itself into singularity at some uncertain timesteps. I think physically the system just moves to a lock-up position when A becomes singular. The simplest solution used by ODE is to add a (diminutive) positive value (CFM) to the diagonal of A such that A will keep away from singularity. A quote from the ODE manual – A nonzero value of CFM allows the original constraint equation to be violated by an amount proportional to CFM times the restoring force lambda that is needed to enforce the constraint.

However, since the global CFM is an arbitrary value added to the diagonal of A, it is inherently unfair comparison when two different As are formulated and applied with the same value of CFM. This case happens when I try to compare the results of the Smith and Baraff formulations. Without the errors introduced by the CFM, the results match pretty well while subjected to some numerical roundoff errors. But the risk is either result may blow up under different circumstances without the help of the CFM. Adding the CFM allows the errors to accumulate and my guess is the errors will become substantial after many timesteps.

The system I used for this method comparison is of course my favorite piston assembly simulation. I may just try imposing the CFM to only one of the boundary joints to reduce the errors while retaining stability. After looking around for a while, it seems that only the CFM of the limit motor can be changed …

Linear-Time Dynamics using Lagrange Multipliers

October 26, 2005

I finally implemented Baraff’s algorithm from the paper, “Linear-Time Dynamics using Lagrange Multipliers”. I have compared the results of the Baraff method with those of the original method (Smith’s ODE) using the piston assembly simulation. The holy grail seems to be sparse LCP, but now I am content with the current implementation. According to the paper, the idea is to first solve for the primary constraints using the sparse LDLT decomposition, and then solve for the auxiliary constraints using the LCP method. The followings are the essential steps.

Given the external force F and the velocity vt at the current time step, the goal is to find the velocity at the next time step, vt + h, where h is the time step size. The system equation can be formulated as

where M is the system inertia matrix, J1 is the primary constraint Jacobian, J2 is the auxiliary constraint Jacobian, λ is the primary constraint impulse, and μ is the auxiliary constraint impulse.

First, solve

using sparse LDLT decomposition to yield

Next, formulate A as

A = J2 M – 1 K

where K is defined as

where B is solved from the following equation

Note that the LDLT factorization has been done for the matrix

in the previous step, so solving for B is rather efficient.

Finally, formulate and solve the following equation for μ using the LCP method.

The velocity vt + h is given by

Piston Assembly Simulation

September 27, 2005

I am not aware that it took me two weeks to finish the piston model. Now the piston model seems to work fine, and it is actually visually stunning. I can just stare satisfyingly at the running piston for a long period of time. My plan is to package and distribute it as a demo. My other thought is to painstakingly re-implement the whole thing in Matlab and then accustomize it with Baraff’s linear-time Lagrange constraint solver. Cross-checking the results in Java and Matlab should provide a solid foundation for this new methodology (well, 10 year old).
On second thought, I might just think it through and do it in Java and compare the results of the new method to those of the old method.

Implementing Piston Assembly Model

September 14, 2005

Very difficult running the C version of ODE and comparing it to my Java version. Instead of debugging the hinge joint, I continued to implement the slider joint. Then I moved all the test programs to the test package. Furthermore, I merged the rapid project with this TestHinge project. Now the TestHinge project includes Trimesh and I can actually continue to implement the piston project.