Hi, I've been reading through Muller's paper, and I'm interested in applying PBD to 2D rigid bodies. Does anyone have a suggestion in regard to updating angular velocity and applying friction / restitution?
Thanks!
Position Based Dynamics (PBD)
-
- Posts: 10
- Joined: Mon Apr 16, 2007 7:11 pm
- Location: Washington, DC
-
- Posts: 861
- Joined: Sun Jul 03, 2005 4:06 pm
- Location: Kirkland, WA
Re: Position Based Dynamics (PBD)
For angular constraints you check the current relative angle and set it back to the initial relative angle. Then as for the linear velocities you set the angular velocity as the difference between the old and the new angle divided by the time step.
I am not aware of any good solutions for friction or restitution in position based dynamics since you don't have any geometric (position) constraints for this. Several cheats exist. You can search the forum for some suggestions. I think Mueller also suggests something for friction in his paper.
-Dirk
I am not aware of any good solutions for friction or restitution in position based dynamics since you don't have any geometric (position) constraints for this. Several cheats exist. You can search the forum for some suggestions. I think Mueller also suggests something for friction in his paper.
-Dirk
-
- Posts: 197
- Joined: Sat Aug 19, 2006 11:52 pm
Re: Position Based Dynamics (PBD)
We have a 2D rigid-body simulator based on Jakobsen/Muller, however each body modeled as is a pair of particles connected by a stick rather than a true 2D rigid bodies (i.e 2D position + 1D orientation). You should definitely be able to apply the same ideas (least-squares fitting of the system state to the constraint/error functions via Gauss-Newton) to proper 2D rigid bodies.
Our "stick-based" approach was chosen in order to simplify the calculation of the Jacobians; this is similar to Jakobsen's idea of using a triangle of particles, but simpler and faster. You can of course embed constraint geometry (collision shapes, joint positions or directions/axes) in each body's space.
To support velocity constraints, there are a couple options:
1) Explicitly solve at both position and velocity level. You can just plug in velocities instead of positions into the solver, for the most part the Jacobian for positions is the same as for velocities. For example to constrain the position of two points p0,p1 along an axis parallel to n you would have C = Dot(n,p1 - p0), to constrain the velocity of the points you just use C = Dot(n,v1 - v0). In this case the velocity-level Jacobian is the same as for the position constraints, except that you would replace p with v in the formula.
I'm not sure if it's best to solve positions and then velocities, or vice-versa. We used velocity-before-position because this makes it like semi-implicit Euler integration, and also Box2D uses this order so it must be good
2) Implicitly solve velocity constraints at the position level: since you know that the velocity v will be inferred as: v = p_post - p_pre (where p_pre is the old/before-solving position and p_post is after solving), you can infer how to alter p_post in order to satisfy a velocity constraint. For the above example, you would have C = Dot(n, (p1_post - p1_pre) - (p0_post - p0_pre)), where the _pre values are taken as fixed/constant and you solve for the _post values.
Of course, this will violate some position constraints (since you're moving the position of things), but that's no different than any other situation -- through iterative relaxation the solver finds an amicable solution which violates all constraints the least amount possible. The benefit of this approach is that you're only solving at the position level so it's less work. When we first thought of it I expected problems since, for example with friction, you'll be moving the contact point along the surface tangent, which will move it out of collision if the surface isn't flat (i.e if it's two circles colliding for instance). In practice we haven't noticed this to be a problem, however we're not doing really hard-core accurate simulation.
Both approaches should let you model friction.
For constraints related to orientation we define directions as vectors and use atan2, i.e C = atan2(d1) - atan2(d0) would constrain d0 and d1 to be parallel; angular velocity is inferred from pre-solver state as with linear velocities, i.e C = (atan2(d1_post) - atan2(d1_pre)) - (atan2(d0_post) - atan2(d0_pre)) would constrain d1 and d0 to rotate at the same rate.
A more interesting/complicated problem is how to properly model motors in a position-based simulation. We've tried two different methods, one is to just copy what Box2D does (i.e motors are constraints whose summed response is clamped to within some upper/lower limit) and another is much weirder; neither is perfect so far unfortunately.
Collision is also a pain in comparison to velocity-level solvers such as Box2D; since the solver can move positions, you really should run collision detection at each solver iteration, but this is quite expensive.
Our "stick-based" approach was chosen in order to simplify the calculation of the Jacobians; this is similar to Jakobsen's idea of using a triangle of particles, but simpler and faster. You can of course embed constraint geometry (collision shapes, joint positions or directions/axes) in each body's space.
To support velocity constraints, there are a couple options:
1) Explicitly solve at both position and velocity level. You can just plug in velocities instead of positions into the solver, for the most part the Jacobian for positions is the same as for velocities. For example to constrain the position of two points p0,p1 along an axis parallel to n you would have C = Dot(n,p1 - p0), to constrain the velocity of the points you just use C = Dot(n,v1 - v0). In this case the velocity-level Jacobian is the same as for the position constraints, except that you would replace p with v in the formula.
I'm not sure if it's best to solve positions and then velocities, or vice-versa. We used velocity-before-position because this makes it like semi-implicit Euler integration, and also Box2D uses this order so it must be good
2) Implicitly solve velocity constraints at the position level: since you know that the velocity v will be inferred as: v = p_post - p_pre (where p_pre is the old/before-solving position and p_post is after solving), you can infer how to alter p_post in order to satisfy a velocity constraint. For the above example, you would have C = Dot(n, (p1_post - p1_pre) - (p0_post - p0_pre)), where the _pre values are taken as fixed/constant and you solve for the _post values.
Of course, this will violate some position constraints (since you're moving the position of things), but that's no different than any other situation -- through iterative relaxation the solver finds an amicable solution which violates all constraints the least amount possible. The benefit of this approach is that you're only solving at the position level so it's less work. When we first thought of it I expected problems since, for example with friction, you'll be moving the contact point along the surface tangent, which will move it out of collision if the surface isn't flat (i.e if it's two circles colliding for instance). In practice we haven't noticed this to be a problem, however we're not doing really hard-core accurate simulation.
Both approaches should let you model friction.
For constraints related to orientation we define directions as vectors and use atan2, i.e C = atan2(d1) - atan2(d0) would constrain d0 and d1 to be parallel; angular velocity is inferred from pre-solver state as with linear velocities, i.e C = (atan2(d1_post) - atan2(d1_pre)) - (atan2(d0_post) - atan2(d0_pre)) would constrain d1 and d0 to rotate at the same rate.
A more interesting/complicated problem is how to properly model motors in a position-based simulation. We've tried two different methods, one is to just copy what Box2D does (i.e motors are constraints whose summed response is clamped to within some upper/lower limit) and another is much weirder; neither is perfect so far unfortunately.
Collision is also a pain in comparison to velocity-level solvers such as Box2D; since the solver can move positions, you really should run collision detection at each solver iteration, but this is quite expensive.
-
- Posts: 10
- Joined: Mon Apr 16, 2007 7:11 pm
- Location: Washington, DC
Re: Position Based Dynamics (PBD)
Many thanks for the great insight! Using a stick model is a very interesting concept, although I'm attempting to take the traditional RB approach before I try anything fancy
p and the surface normal (n) at this position. Then I use the constraint, C = (p - q) dot n, as shown in section 3.4, to solve for dP and dQ. This actually gives me very good results when bodies have no angular momentum; I can stack many objects (boxes) with very good stability - visually I detect absolutely no bounce....
However, once I induce angular momentum and solve for change in orientation, and update angular velocity, resting contacts bounce and vibrate all over the place - there doesn't seem to be any loss in kinetic energy. Is this where friction should come into play?
Since the solver is updating position, why also constrain linear velocity? Why not update linear velocity as show on line 13 in the algorithm presented in section 3.1, and also do the same for angular velocity once you solve for the change in orientation?raigan2 wrote: I'm not sure if it's best to solve positions and then velocities, or vice-versa.
This is an interesting idea. I also simply measured the angle between the point of contact before and after position correction, which is less work computationally than solving equation (9) with the gradient of C = atan2(d1) - atan2(d0).... Is this not a good method?raigan2 wrote: For constraints related to orientation we define directions as vectors and use atan2, i.e C = atan2(d1) - atan2(d0) would constrain d0 and d1 to be parallel
I've been using MPR (http://code.google.com/p/mpr2d/), to try something new, which gives me the deepest point of penetration (p), the surface point (q) which is closest toraigan2 wrote: Collision is also a pain
p and the surface normal (n) at this position. Then I use the constraint, C = (p - q) dot n, as shown in section 3.4, to solve for dP and dQ. This actually gives me very good results when bodies have no angular momentum; I can stack many objects (boxes) with very good stability - visually I detect absolutely no bounce....
However, once I induce angular momentum and solve for change in orientation, and update angular velocity, resting contacts bounce and vibrate all over the place - there doesn't seem to be any loss in kinetic energy. Is this where friction should come into play?
-
- Posts: 197
- Joined: Sat Aug 19, 2006 11:52 pm
Re: Position Based Dynamics (PBD)
I simply meant that if you wish to, you can solve at both the position and velocity level in two separate phases, as in Box2D -- this is one possible method to add velocity-level constraints to PBD. This will probably end up being identical to Box2D with NGS position solving.zzzzrrr wrote:Since the solver is updating position, why also constrain linear velocity? Why not update linear velocity as show on line 13 in the algorithm presented in section 3.1, and also do the same for angular velocity once you solve for the change in orientation?
I don't understand what you mean by "measured the angle between the point of contact"; is this for friction? Do you mean than, rather than decomposing collision constraints into two axial/linear constraints -- one along the surface normal and the other along the tangent -- you're instead using two polar constraints, i.e distance and angle?zzzzrrr wrote: This is an interesting idea. I also simply measured the angle between the point of contact before and after position correction, which is less work computationally than solving equation (9) with the gradient of C = atan2(d1) - atan2(d0).... Is this not a good method?
The main problem with collisions from my experience is that, since position and velocity are coupled, this can add energy -- when you separate penetrating objects you introduce some outward velocity. Since the solver is always moving objects, this means that you can't simply generate collision constraints at the start of the frame (as in Box2D for instance) because later on during solving some new collisions may have been created.zzzzrrr wrote: I've been using MPR (http://code.google.com/p/mpr2d/), to try something new, which gives me the deepest point of penetration (p), the surface point (q) which is closest to
p and the surface normal (n) at this position. Then I use the constraint, C = (p - q) dot n, as shown in section 3.4, to solve for dP and dQ. This actually gives me very good results when bodies have no angular momentum; I can stack many objects (boxes) with very good stability - visually I detect absolutely no bounce....
For instance, you'd like your simulation to be something like this:
Code: Select all
Step_Simulation()
{
step world forward in time
find all colliding pairs A,B
for(n iterations)
{
Solve_Joints()
for(each pair): Solve_Collision(A,B)
}
}
Something like this will work better:
Code: Select all
Step_Simulation()
{
step world forward in time
for(n iterations)
{
Solve_Joints()
for(each body A)
for(each body B)
if(A collides with B): Solve_Collision(A,B)
}
}
a) you're now doing a full pass of collision detection at each solver iteration
b) it's still not a perfect solution, because collisions solved later may violate collisions solved earlier (i.e there may still be some residual penetration at the start of the next frame).
c) you can't use any "all-pairs" collision detection method (which are typically faster than pair-by-pair) since you're moving each pair
c) it's not actually possible to do this (as far as I can tell) if you want to avoid O(n^2) collision tests; since A can be moved after the first B, you basically need to re-query your structure every time you move A.. at least, I haven't found a nice way to use broadphase collision systems. Such a way might exist though, and I'm just missing something.
Currently we're doing this:
Code: Select all
Step_Simulation()
{
step world forward in time
find all *potentially* colliding pairs A,B
for(n iterations)
{
Solve_Joints()
for(each pair):
if(A collides with B): Solve_Collision(A,B)
}
}
If you're not seeing problems of this sort, I'd be _very_ interesting in learning specifically how you're processing collisions, since this has so far been the #1 pain of position-based methods.
I'm not sure -- in our case stacks were stable without friction (they would eventually slide/slip apart though), but we're not using proper rigid bodies. I suspect that friction would only help hide whatever the fundamental problem is in this case. Are you modelling bounce/restitution?zzzzrrr wrote: However, once I induce angular momentum and solve for change in orientation, and update angular velocity, resting contacts bounce and vibrate all over the place - there doesn't seem to be any loss in kinetic energy. Is this where friction should come into play?
We simply decided to forget about bounce and only model totally inelastic collisions to start with, since bounce can always be modelled via gameplay logic (i.e when collision is detected, create a velocity constraint for next frame); in our game the vast majority of objects shouldn't bounce anyway.
-
- Posts: 10
- Joined: Mon Apr 16, 2007 7:11 pm
- Location: Washington, DC
Re: Position Based Dynamics (PBD)
Here's the velocity constraint from Erin Catto's 2009 GDC slides: Here's equation (9) from the PBD paper: Would it be feasible to solve for position, or velocity correction, by using the above constraint equation? From what you are suggesting it would be ok to substitute position for velocity and orientation for angular velocity...?raigan2 wrote:... to constrain the position of two points p0,p1 along an axis parallel to n you would have C = Dot(n,p1 - p0), to constrain the velocity of the points you just use C = Dot(n,v1 - v0). In this case the velocity-level Jacobian is the same as for the position constraints, except that you would replace p with v in the formula.
EDIT:
Using the excellent (and free!) SAGE math tool, I get the following for angular velocity:
C.diff(w2) = (p.y - x2.y)*n.y + (p.x - x2.x)*n.x
C.diff(w1) = -(p.x - x1.y)*n.y - (p.x - x1.x)*n.x
I'm wondering if this would be the same gradient for orientation... e.g. is it safe to assume this:
C.diff(w2) == C.diff(q2)
C.diff(w1) == C.diff(q1)
With q1 and q2 representing orientation.
You do not have the required permissions to view the files attached to this post.
-
- Posts: 861
- Joined: Sun Jul 03, 2005 4:06 pm
- Location: Kirkland, WA
Re: Position Based Dynamics (PBD)
You actually end up with a very similar formula for the velocity and the position constraint (using primes instead of dots for simpler reading and typing ):
J*W*JT * lambda = -C
J*W*JT * lambda = -C'
In your example this would be:
C = x2 + r2 - x1 - r1 where r = p - x
C' = J*v = v2 + cross( omega2, r2 ) - v1 - cross( omega1, r1 )
The gotcha is to realize that when solving for the position error and applying the correction (pseudo) impulse doesn't necessarily imply that the constraint is fully satisfied by this since this is a *non-linear* problem. In other words there still can be an error! Remember that you used Taylor expansion:
C(x + dx) ~= C(x) + J * dx
You could take this actually a step further and solve a quadratic formulation (which is unpractical in my opinion, but maybe helps understanding the problem):
C(x + dx) ~= C(x) + J * dx + transpose( dx ) * H * dx
Here H is the Hessian.
Cheers,
-Dirk
J*W*JT * lambda = -C
J*W*JT * lambda = -C'
In your example this would be:
C = x2 + r2 - x1 - r1 where r = p - x
C' = J*v = v2 + cross( omega2, r2 ) - v1 - cross( omega1, r1 )
The gotcha is to realize that when solving for the position error and applying the correction (pseudo) impulse doesn't necessarily imply that the constraint is fully satisfied by this since this is a *non-linear* problem. In other words there still can be an error! Remember that you used Taylor expansion:
C(x + dx) ~= C(x) + J * dx
You could take this actually a step further and solve a quadratic formulation (which is unpractical in my opinion, but maybe helps understanding the problem):
C(x + dx) ~= C(x) + J * dx + transpose( dx ) * H * dx
Here H is the Hessian.
Cheers,
-Dirk