Stable Joints

[raigan2]

hi,
In case anyone's interested in an alternate approach to joints, I've managed to get a particle-based method working which is rather rigid, and quite well-behaved. Check it out here: http://harveycartel.org/raigan/pbdTest16.swf

(click in the window to give the flashplayer focus)
a,s: change the joint angles
z: grab a particle w/ the mouse
1,2: change which particle you're grabbing
~: toggle the first "body" between pinned and free
spacebar: reset

This is based on the method presented in "Position Based Dynamics" by Muller, however it uses a different formulation for the angle joints (using atan instead of acos). After several failed attempts at doing it by hand, I cheated and used LiveMath to determine the formula for the Jacobian

I'm using only 4 iterations; things are much nicer/stiffer with more, but as it is I'm pushing the flashplayer.

raigan

[Jan Bender]

Do you use the same equations for the preview of the vertices as Mueller? Because this is physical incorrect since he first computes the new velocities using the external forces and then the positions with the new velocities (instead of using the old velocities and the forces). If a body falls down for one second (due to gravity g) this leads to a position error of 0.5*g which is quite big (in my opinion).

Or how do you implemented the preview?

Jan

[raigan2]

Integration is something I could stand to learn a bit more about; I'm just using good old semi-implicit Euler/"Newton-Stormer-Verlet"/whatever it's called, which I've always used, simply because it's always worked for me and I've never bothered to look for anything better:

vel += accel
pos += vel

I actually thought that using the new velocity to update position was a benefit, as it reduced the error due to "extrapolating" blindly compared to regular Euler. AFAIK many integrations schemes use future velocity to update current position -- if I understand correctly, that's the basis of predictor-corrector integrators.

I'm not sure if it's accurate, but I think it looks convincing: bodies in freefall accelerate downwards, bursts of particles have arc-like trajectories. Also, the only use I've ever had for "forces" is gravity, which should be the same regardless of whether you're evaluating it at the old or new position. In the demo I posted, g = 0.1, so 0.5*g is only 0.05 pixels; an error of 1 pixel per 20 frames isn't that bad, is it?

I'm certainly quite ignorant about the mathematical vagaries of numerical integration though.

[Dirk Gregorius]

If a body falls down for one second (due to gravity g) this leads to a position error of 0.5*g which is quite big (in my opinion).

Given a gravity of -10 and one timestep of 1s the error would be indeed 5m, but this is an unrealistic assumption. The position difference scales with the size of the timestep. So for a more realistic assumption of a timestep of 1/60 the absolute position difference is 83 mm after 1s. Also noting that all particles are subject to quite alot constraints this is really negligible. The only thing you want to make sure is that the solver sees the gravity directly and not with one frame delay like for the normal Euler. Yours and Muller's solver do this.

I have written a small sample app so you can verify this.

Code: Select all

#include <iostream>

const float g = -10.0f;
const float dt = 1.0f / 60.0f;


int main( void )
	{
	float x1 = 10.0f;
	float v1 =  0.0f;

	float x2 = 10.0f;
	float v2 =  0.0f;

	float x3 = 10.0f;
	float v3 =  0.0f;

	for ( int i = 0; i < 60; ++i )
		{
		//
		// Euler
		//
		// x(t+dt) = x(t) + v(t) * dt
		// v(t+dt) = v(t) + a(t) * dt
		//
		x1 += v1 * dt;
		v1 +=  g * dt;

		// 
		// Semi-implicit Euler
		//
		// v(t+dt) = v(t) + a(  t ) * dt
		// x(t+dt) = x(t) + v(t+dt) * dt
		// 
		v2 +=  g * dt;
		x2 += v2 * dt;

		//
		// 2. Order Euler
		//
		// x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt^2
		// v(t+dt) = v(t) + a(t) * dt
		x3 += v3 * dt + 0.5f * g * dt * dt;
		v3 +=  g * dt;

		printf( "E1: (%.3f, %.3f), E2: (%.3f, %.3f), E3: (%.3f, %.3f)\n", x1, v1, x2, v2, x3, v3 );
		}

	return 0;
	}

[Erin Catto]

You cannot simulate a simple mass-spring system with 1st order or 2nd order Euler. You have to add damping. A small time step won't help.

That is quite unrealistic.

[Dirk Gregorius]

This is no spring. It is a particle vertically falling down. I just wanted to show that the position difference between first and second order Euler is not 0.5 * g after one second for this particular case. Basically the point was if it makes a difference in the "Position Based Dynamics" paper to use this:

IntegrateAllParticlesUsingSemiImplictEuler();
SatisfyAllConstraints();

versus

IntegrateAllParticlesUsingSecondOrderEuler();
SatisfyAllConstraints();

Jan is correct for pointing this out and I guess you can use this, but I questioned that this makes such a big difference as he stated above. Another reason for the first solution might be that they use CCD what might make their life a little easier when dealing with the collision constraints. Just out of couriosity - has anybody success in simulating skirts with this method or any other more complicated contact setup?


Cheers,
-Dirk

[Jan Bender]

Given a gravity of -10 and one timestep of 1s the error would be indeed 5m, but this is an unrealistic assumption. The position difference scales with the size of the timestep. So for a more realistic assumption of a timestep of 1/60 the absolute position difference is 83 mm after 1s.

8.3 cm error during 1/60 s is quite much in my opinion but since you already know I am a fan of accuracy. In my PhD thesis I often work with a maximal allowed tolerance of 10^-6 m. Compared to this 8.3 cm is very much

Jan