Hi All,
on ODE manual, it says, "ODE has hard contacts. This means that a special non-penetration constraint is used whenever two bodies collide."
But in my simulation example, which basicly was a ball rolling on ground, somehow, I read out the contact point position and found the Z - direction (vertical)
coordinates were -0.005. Which means the contact point was underground. How did it happened? Since ODE has Hard contacts for each collision.
Any way to avoid such penetration? Thanks.
Ruby
Hi!
"hard" here means the constraint will exist in the whole time step. In ODE case, we:
1. Set up the Jacobi correspond to the constraint
2. Solve for the constraint force fc
Then because the constraint is "hard" so we will apply fc to the involved bodies through the whole timestep. That's actually not "exact" because the contact might be broken in-between the time step.
This problem is inherent to any force-acceleration level simulation and along with friction, there is singularity problem also.
If you change the constraint system into velocity-impulse level then you don't have "hard" or "soft" contact anymore because you treat all of them as impulses.
Binh Nguyen
Nguyen,
I would argue that the ODE solves on the velocity level. Basically you have an DAE (for simplicity for holonom constraints):
dv/dt = W * ( F_ext + JT * lambda )
dx/dt = v
C(p,q) = 0
We use the differential quotient for the accelerations:
dv/dt ~= ( v(t+dt) - v(t) ) / dt -> v(t+dt) = v(t) + W * ( F_ext + JT * lambda ) * dt
We can find the velocity constraints through differentiation of the position constraints
dC/dt = dC/dx * dx/dt = J * v(t) = 0
Since we use the velocities of the end of the timestep (Symplectic Euler) we can require J * v(t+dt) = 0
If you combine these equations you end up with the well know LS (or MLCP in the presence of non-holonom constraints)
J*W*JT*lambda = -( v(t)/dt + W * f_ext )
You could multiply the whole equation by dt > 0 and get
J*W*JT*(lambda*dt) = -( v(t) + W * f_ext * dt )
You clearly see the forces and impulses are linear related in this model. Also I argue that in every discrete model a constraint can't break. Of course a constraint force could become zero (through clamping), but the ODE can handle this as well.
So regarding the original question. Because of the discrete collision model you can end up in penetration. This has nothing todo with hard or soft constraints. The ERP parameter basically handles how much of the error is corrected each timestep. Also note that ODE has some value that allows for some penetration before the error correction jumps in. This is in order to deal with jitter and to improve coherence. So the 0.005 you see might be the default value for the allowed penetration.
HTH,
Dirk Gregorius
A hard constraint basically means that you have a "hard" algebraic
constraint equation instead of plugging in a "soft" spring like in a
penalty method. The ODE requires for each contact:
1) The penetration at the contact should vanish
C(p,q) = (p1 + r1 - p2 - r2) * n - d = 0
2) The relative velocity should vanish
dC/dt = (v1 + cross( omega1, r1) - v2 - cross(omega2, r2) * n
-> J = ( n, cross( r1, n ), -n, -cross( r2, n ) )
Since we apply only forces when we have an approaching relative velocity
at the contact we basically end up with
w = A*x + b >= 0 and x >= 0 and xT * w = 0 which is a classic LCP or
MLCP in presents of additonal holonome constraints.
Actually each position constraint imposes also a velocity- and
acceleration constraint onto the system. Think of a ball socket joint.
The position constraint requires that the anchors are at the same point
in space. In order to maintain this constraint the relative velocity and
acceleration must vanish at the anchor. Otherwise the anchors would get
apart. So *all* constraints (position, velocity and acceleration) must
always be satisfied. The ODE solves only at the velocity level. Note
that even when you fully satisfied the velocity constraints, you still
can end up with violated position constraints (e.g. because of high
angular velocities). The ODE deals with this using Baumgarte
stabilization. Erin Catto's has a nice presentaion about this at
http://www.gphysics.com.
Dirk Gregorius
But ODE still solve for constraint forces. That makes me refer to ODE as acceleration-force model.
The equations you showed are basically discretized version of Newton-Euler equation. And yes, ODE use velocity as variables there.
Because ODE solve for the contact forces so the constraints will be satisfied at velocity level.
But in a velocity-impulse contact model will solve for contact impulses so it will maintain the constraints *exactly* at position level and you don't need Baumgarte stabilization which adds stiffness to the system.
Nguyen Binh
How can this be?
Think of a pendulum. Your velocity constraint requires that the realtive
velocity is tangential to the pendulum. Now we satisfy this (velocity)
constraint. If you integrate now forward using the constraint satisfying
velocities you will by no means meet the position constraint. You will
end up outside the circle described by the pendulum
So I assume I misunderstand you. Can you please give an example how you
find impulses to project on the velocity constraint manifold and then
after some update of the positions I also atomatically satisfy the
position constraint?
Here is an example:
In an impulse solver (equivalent to the ODE) you would solve as follows:
C = |p2 - p1| - L0
dC/dt = (p2 - p1) / |p2 - p1| * (v2 - v1) -> J = ( -(p2 - p1) / |p2 -
p1| (p2 - p1) / |p2 - p1| )
1) Update velocities
2) Apply iteratively sequential impulses until all velocity constraints
are saisfied
lambda = - J*v / J*W*JT
P1 = -(p2 - p1) / |p2 - p1| * lambda
P2 = -P1 = (p2 - p1) / |p2 - p1| * lambda
v1 += P1 / m1
v2 += P2 / m2
3) Update positions using new velocities
How would this look in your model?
Dirk Gregorius
Good catch!
You are right about the fact that it will end up outside the circle. But that's because we linearized the position constraint. In the pendulum case, we have C(q) = q^2 - r^2 = 0 but we only enforce it in velocity level.
That's very bad if you use acceleration-force model because the end point will move away from original circle every time step. With velocity-impulse, the end point will move in a sets of line that circumscribed the circle, which is not so bad.
Velocity-Impulse:
+Bilateral constraint: C(q) = 0
+Taylor series expansion: C(q)[t+1] = C(q)[t] + dC/dq*delta(q) + dq/dt*h = 0; (1)
h is time step, C(q)[t+1] means the value of C(q) at time t+1.
Note that we only use the linear terms, you can use more but it will make your system non-linear.
+ Divide by h to have: C(q)[t]/h + J*v[t+1] + delta(q) = 0; (2)
Here the variable is v[t+1], means we need to find impulses such that v[t+1] satisfy (2).
If you apply this to the pendulum case, you can see that the end point will never sway too far from the circle. And if you decided expand Taylor series in (1) then it will satisfy C(q) exactly.
That's the constraint we need to satisfy. You can see that it's not exactly the same as your constraint.
All of this method has been discussed quite alot in papers by Jeff Trinkle. Anitescu, Potra.
Nguyen Binh
