Finite Elements Method: are there proble,s with stability?

[h4tt3n]

I'd like to hear more about how h4tt3n avoids the unexpected displacements - do the ridiculously high K (spring rate) values still behave stably if you have a bunch of them in series, say a pendulum with a bunch of links and a heavy mass at the end?

Yes they would. The stability of a force-based simulation does not depend on the number of interconnected masses or on which geometric shapes you create with them. It is purely a question of balancing timestep, spring stiffnes, and particle mass.

[XperienS]

h4tt3n
Could you explain the method in more details or show us some links on what you're calling iterative velocity verlet scheme?

bone
Have you tried some post-projection after the main solver steps, which affects only position of spring nodes (to correct solver inaccuracy)?

all
The thing about correction is that currently these methods used for absolutely stiff joints, which can be corrected even with Baumgarte, but for spring-like joints ERP/CFM is a part of its normal functionality. So it seems like you should apply ERP/CFM scheme for spring behaviour and post-stabilisation like pseudovelocities or more accurate for correcting that solver inncauracies and lower the ghost forces. But how shall we calculate error in spring joint position, since every position in world is valid for spring? (we can't say that for stiff joints like, say, ball-socket joint; we can always say if there is some position drift) We should somehow analyze new state of the spring according to the old one? For Gauss-Seidel solver residual vector b - A*x should help us a bit to determine the magnitude of ghost forces introduced by the solver.

[h4tt3n]

h4tt3n
Could you explain the method in more details or show us some links on what you're calling iterative velocity verlet scheme?

Certainly, it's almost ridiculously simple. The idea is as follows. Let's say we want to model a chain of particles held together by stiff constraints. Normally you would do this by implementing some iterative impulse based solver which spends, let's say ten iterations per game loop. Now the idea is that it will cost more or less the same in terms of cpu-cycles if we dropped the impulse based constraints and instead implemented force based springs between the particles and ran the integrator 10 times per game loop with 1/10 of a timestep. I realise that this method is only conditionally stable, but it's still doing a really good job compared to poor constraint solvers.

I start out with the velocity verlet scheme:

Code: Select all

dt = timestep
x += v * dt + 0.5 * a * dt * dt  // update position
v += 0.5 * a * dt                // update half velocity
a = f / m                        // call function that updates all forces based on new position and velocity
v += 0.5 * a * dt                // update half velocity

where x = position, v = velocity, a = acceleration, f = force, m = mass, dt = timestep. Vectors are shown in bold. Then I iterate through the scheme as many times as desired, updating the force vector for each particle per iteration:

Code: Select all

N = 10
dt = timestep / N
for i = 1 to N
  x += v * dt + 0.5 * a * dt * dt  // update position
  v += 0.5 * a * dt                // update half velocity
  a = f / m                        // call function that updates all forces based on new position and velocity
  v += 0.5 * a * dt                // update half velocity
Next i

where N is the number of iterations. Even with moderate N values this method can handle spring stiffnes values (k) ranging form 100.000 to 1.000.000 for particles with m = 1 and timestep = 0.01.

Here's a link for a little softbody toy application I made using this method. The .zip file includes a code::blocks project, full source, and a windows .exe (the actual integration algorithm is in softbody.cpp somewhere around lines 440 to 480). Enjoy

Cheers,
Mike

[bone]

bone
Have you tried some post-projection after the main solver steps, which affects only position of spring nodes (to correct solver inaccuracy)?

Yes, but the question is: what is the correct position? If I get rid of any position change that I didn't expect at the beginning of the timestep, I lose some valid effects.

Consider a mass swinging around a center point, attached by a spring. Say that the starting velocity is perpendicular to the spring, and that the spring is at its resting length. Initially, then, there should be no correction impulse. But at the end of the timestep, I see that there is an error, but not all of the error is necessarily invalid. If I correct all of it, then the mass will continue swinging around the center point at the spring's rest length, which is wrong.

Now I *did* try some fancier ideas here, for example looking at the error and saying "well, if this spring had known this would have happened, it would have reacted like this". But all of them just tended to increase stability problems in complex systems. My brain might not be big enough for this problem.

[raigan2]

instead implemented force based springs between the particles and ran the integrator 10 times per game loop with 1/10 of a timestep.

This reminded me of this paper: http://leri.univ-reims.fr/~nocent/papers/kacic03.pdf

Also I've spoken with the programmer who made Gish/Bridge Builder/etc. and they use a similar approach: force-based model with many substeps.

raigan

[h4tt3n]

This reminded me of this paper: http://leri.univ-reims.fr/~nocent/papers/kacic03.pdf

Sorry, the link is broken, and google doesn't know of any "kacic03.pdf". Have you got another link?

Also I've spoken with the programmer who made Gish/Bridge Builder/etc. and they use a similar approach: force-based model with many substeps.

Cool, I've always liked Edmund McMillen's work. Didn't know he originally made bridge builder?

Cheers,
Mike

[raigan2]

Cool, I've always liked Edmund McMillen's work. Didn't know he originally made bridge builder?

AFAIK Edmund only did the graphics for Gish; the programmer for all these games ( http://crypticsea.com/ ) is Alex Austin. I've been trying to get him to write up a presentation on his simulator for years

[XperienS]

Yes, but the question is: what is the correct position?

Well, as I wrote above - error in spring goes from solver (and a part of it from the integrator due to discretization), so basically if we know how wrong is our solver, we can determine what error contains the spring. The problem is what should we do with that error.

h4tt3n
I believe that mentioned paper is "A practical dynamics system", Zoran Kačić-Alesić

[h4tt3n]

@ExperienS

Thanks a ton for the link. The solver described is very similar to mine, except they seem to use the leap-frog method where I use velocity verlet. I found chapter 4.2 about stability and the size of the time step very interesting. For quite a while I've been looking for a method to theoretically determine the highest allowable timestep, and there it was, plane and simple. I'll be implementing it in my engine asap, which should make it de-facto unconditionally stable. Wohoo!

Cheers,
Mike

[raigan2]

Thanks for correcting the link, sorry about that! It took a while just to find a link to the actual paper and not some stupid ACM page, I failed to test it

[FD]

I feel I must say that after we solved the system with Conjugate Gradients method, the stability turned out to be quite OK, so the method works perfectly if an appropriate solver is used. However, I know no way to combine CG with inequality constarints, so the main goal - integration of FEM into the system of contrainted particles and bodies - is not achieved since contact joints lead to LCP and a system with cointacts can't be therefore solved with CG.

[Dirk Gregorius]

Did you look at the Pixelux presentation from the GDC? The use some kind of special element and still can solve using CG. Also if you use CG you might use some kind of constraint satisfaction technique as Baraff in his cloth paper (Large Steps in Cloth Simulation).


HTH,
-Dirk

[FD]

Hi Dirk!

Thx for the links.. The Baraff's approach looks like a dismal piece of crap, this is, however, an option. I'd prefer to know what Pixelux guys are doing, don't you know?

[Dirk Gregorius]

I visited their session, but cannot remember exactly. There is a Siggraph paper by O'Brian which talks about Pixelux. Maybe you can also find the GDC presentation on the GDC website.


Cheers,
-Dirk

[job]

I visited their session, but cannot remember exactly. There is a Siggraph paper by O'Brian which talks about Pixelux. Maybe you can also find the GDC presentation on the GDC website.


Cheers,
-Dirk

SCA 09 paper on DMM:
http://www.cs.berkeley.edu/b-cam/Papers ... -2009-RTD/

SIGGRAPH 09 paper that discusses a similar problem and uses a Newmark integrator:
http://www.cs.berkeley.edu/b-cam/Papers ... -2009-ISN/

GDC does not really have official papers to go with the talks, but the details you're interested in should be in the SCA paper.

If you're having issues with not getting the expected stability improvements with an implicit integrator, it's probably a good idea to check the Jacobian matrix carefully and then double check it. An error computing the matrix will often still give you a working simulation, but you won't get the stability improvements you were expecting.