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

[Dirk Gregorius]

I am not too deep in FEM to answer this. FEM is only on my radar. The think the problem with constraint stabilizations is indeed that you want to solve such that C( x(t+dt) ) = 0. Baumgarte takes the error from the begining of the time step ( C(x(t) ) = 0 which is not what you want. Or as Uri Asher phrased it: Baumgarte doesn't necessarily stabilize on the best trajectory. There are other methods like Jan Bender, but personally I am not convinced of this approach.

I was lucky to have some private talks with Claude who is one of the authors in the mentioned papers above. Personally I would stick with his stabilization and regularization methods. This is very mathimatically sound as far as I can tell. When you have a good solver and not the crappy iterative Gauss-Seidel this might work very well in practice.


Cheers,
-Dirk

[FD]

Dirk,

I'm struggling to make it work on the crappy Gauss-Seidel. Besides that, I'm implementing my own technique, though it seems to be virtually the same as in the paper, but I still I've come to similar formulas in a substantially different way and I don't want to limit myself to the formulas in the paper.

I'm asking you not about FEM, I'm asking you about springs. Maybe you have tried to make spring joints between rigid bodies, not making them springy for stability but actually modelling spring-like behavior. The question is what methods can I use to stabilize such a joint, in the sense of providing a stable spring-like behavior, not making it stiff.

I've been experimaneting today with residual-based corrections. I was trying to solve a system in positions that would eliminate forces produced by residuals. It seems to be stable enough only if I don't account for masses, which makes the system rather unnatural..

[Dirk Gregorius]

You can derive springs using distance constraints (using primes for dots since I find it easier to read - it denotes the time derivative so C' = dC/dt).

C = | x2 - x1 | - L0
C' = ( x2 - x1 ) / | x2 - x1 | * ( v2 - v1 )

So the Jacobian becomes:

Define n = ( x2 - x1 ) / | x2 - x1 |
-> J = ( -n n )

The time derivetive has a special form:

C' = J * v

Instead of solving C' = 0 we solve (this is actually an index reduction but forget about the math for a second here):

C' + ERP * C / dt + CFM * lambda = 0

So you end up with this formula:

(J*W*JT + CFM) * lambda = -ERP * C / dt - J * v

Here you find a good explanation and how ERP and CFM match to spring constants:

http://www.ode.org/ode-0.5-userguide.html#sec_3_7_0
http://www.ode.org/ode-0.5-userguide.html#sec_3_8_0


Let me know if you have further questions!

-Dirk

[nisse]

The "How to make new joints in ODE" http://ode.org/joints.pdf has some even more details on this as well.
One thing that is not stressed enough in that text, I think, is that for the mapping of the spring stiffness and damping
coefficients to CMF and ERP to hold relies critically on that the Jacobian precisely matches the constraint function in the sense

J = dC / dx

This is important when you want to construct new constraints (e.g. stiff finite elements or a complicated bio-mechanical joint)
and give them a desired elasticity.

[Dirk Gregorius]

This is very true and could not be stressed enough!!! Don't start estimating the position error which is what a lot of people do for angular constraints (e.g. in generic 6DOF joints). This gives spongy results.

The velocity constraint has a special form: dC/dt = J*v + del_C / del_t

Here del_C / del_t is the *partial* time derivative. The good news is that it is always zero in our cases since the constraints don't depend explicetely on time. Here is a link why the derivative has a special form:

http://en.wikipedia.org/wiki/Total_derivative



HTH,
-Dirk

[FD]

OMG, Dirk!

The problem is not making the spring constarint. The formulas that you have given here are well-known (however, it's useful to have them here explicitly in our discussion). I wonder how to make this crap STABLE. You see, if the solver is iterative and produces a substantial error, the new positions on each iteration are likely to be different from the theoretical result for the spring, and hence it is likely that forces of increasingly higher magnitude will be produces on each iteration, threfore increasing the position error produced by the solver and further increasing the forces, and the system will diverge in some time. This is an importaint problem for very stiff springs. This is why Baumgarte stabilization is unstable for high stiffness. But it is possible to substitute Baumgarte stabilization with other methods like pseudoimpulses or some kind of position-based correction step. It is because the spongy behavior of Baumgarte-stabilized hard joints is somewhat undesireable, and we are happy to introduce a stabilization that would not produce such a behavior at all. A different story is for spring. For them, we WANT to observe the spongy behavior, but still we DON'T WANT the system to diverge and blow us into the face. So I am interested if there're stabilization techniques that increase stablity on tyhe one hand but allow us to keep the springs working correctly on the other hand.

Thanks.

[nisse]

FD,

Have you (or anyone else) tried preconditioning of the iterative solver? That is, you warm-start the solver with lagrange multipliers that is the best estimation of the solution you can do a priori - instead of using the initial iteration solve. The estimation could be a direct solve of a subsystem or starting from the lagrange multipliers from the last time step, or on a guess based on that solution. This is hopefully a better initial guess than what is obtained from the first (blind) iteration solve.

[Dirk Gregorius]

Honestly, you cannot expect to solve a system with more then 1000 unknowns solving each spring locally in an iterative way (aka Gauss-Seidel). In physics this is why we use projection methods as discussed here: http://www.gphysics.com/archives/35

Let me know if you find a solution for the problem!


-Dirk

[Dirk Gregorius]

Warmstarting is very common with this method. My observation is:

1) It works great with contacts
2) It doesn't work so great with joints when using Baumgarte (high oscillations)
3) Warmstarting joints with projection works very good, but doesn't recover well if the error gets big

So currently I use warmstarted contacts and non-warmstarted joints with a slightly higher Baumgarte for joints and regularization. This works well in practice.

You can look at some of the results with this method and test yourself here: http://www.dondickied.de/download/brian.zip

There is some problem with setting joint ERP/CFM after a reset. The GUI makes problems here. You have to Pause -> Reset -> Set Parameter -> Run

Compare e.g. the chain with its default setting *with projection* instead of Baumgarte against maybe:
Warmstart: OFF
Stabilization: Baumgarte

Joint ERP: 0.1 - 0.4
Joint CFM: 1.e-3 - 1.e-5

To see the issues with projection methods under large errors use the mouse.




And yes, I know this is of topic. We could potentially split this discussion into two threads if there is interest in the direction the discussion is currently going.

Cheers,
-Dirk

[FD]

nisse,

yes, Dirk is right, warmastring is very common, it's nice for contacts, and for joints as well but it doesn't help much if you use Baumgarte stabilization..

Dirk,

you CAN solve big systems in an iterative manner as well, Gauss-Seidel and PGS will converge; however you'd need too many iterations to get an appropriate result )))

Currently I'm balancing between unreasonably strongly damped systems and unresoanbly many iterations; both things are bad. My attemts to adapt projective methods to the correction of the estimated error of the solver seem to have failed.

[Dirk Gregorius]

Hi FD,

I am sure you know this, but if you use CFM and implement PGS in terms of sequential impulses (SI) the derivation goes a little bit different. Let P be the current accumulated impulses and dP the delta impulse we will compute. So in the derivation in my earlier post dP is lambda and P is the sum of all individual lambda.

M ( v2 - v1 ) = JT * dP
J * v2 + CFM * (P + dP) + ERP * C / dt = 0

( J*M^-1*JT + CFM )* dP = -ERP * C / dt - J * v1 - CFM * P


Note the additional term on the right hand side. This is sometimes overlooked when using pairwise relaxation (SI).


Cheers,
-Dirk


Edit: Fixed typo

[bone]

Currently I'm balancing between unreasonably strongly damped systems and unresoanbly many iterations; both things are bad. My attemts to adapt projective methods to the correction of the estimated error of the solver seem to have failed.

It sounds like you're attacking roughly the same problem as me. I haven't had perfect success yet, but I made a big jump in stability by increasing damping only when the frequency exceeds the sampling rate. To me this isn't really a problem - you can't get perfect results in that case no matter what you do.

However, I'm not solving the spring constraints in exactly the same form as the CFM/ERP method at the moment, so I don't know how helpful this advice is. (Note: the reason it's not the exact same form is that either due to an implementation mistake or an actual limitation of the CFM/ERP formulation, I don't seem to get the correct results for springs in parallel (EDIT: or was it in series? can't remember!)).

[FD]

bone,

how do you cakculate the frequency?

[h4tt3n]

Just out of curiosity, does anything prevent you from implementing a force-based Hooke's law spring and solve for velocity and position with an integration algorithm? I'm using an iterative velocity verlet method which allows me to implement undamped springs with ridiculously high K values without blowing up. This method also preserves energy very well, which I think iterative constraint solvers do not (correct me if I'm wrong).

[bone]

bone,

how do you cakculate the frequency?

The spring rate is an input to my system, and the K matrix (J*M^-1*J^T) is the reciprocal of the effective mass. I'm iteratively solving line-by-line, essentially, so the K matrix is just 1x1. The frequency then is just sqrt( spring_rate/mass).

Like I said, though, this still doesn't work perfectly and I'm still frustrated by the issue. No matter what you do, solver inaccuracies and integration are going to result in unexpected displacements in a complex system, which adds phantom energy to the spring. My dampers are perfectly stable, though, so I just add enough to get rid of a random amount of the phantom energy, which works in many cases.

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?