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

[FD]

Hi people!

I'm playing with FEM simulation of deformable bodies. Basically, no matter how you do rearrange the system, you arrive at the matrix form of Hook's law F = -K * delta_x. I'm interested in more or less stiff matrials, and it turned out to be challanging because of the well-known integration problem with stiff springs (it often diverges even when using very small timesteps).

How do you make stable FEM?

Those who haven't implamented FEM, maybe you know something related to the general problem of stiff springs instability?

[nisse]

The problem is not a FEM problem but rather a numerical-integrator problem and araises also for mass-spring systems and more. You must pick a numerical integrator that can handle your choise of the material stiffness and time-step size. A common choise is the (linearized semi-) implicit Euler. High stiffness implies that the object may have high frequency normal modes of oscillations. If any of these frequencies are higher than the integration frequency (1/time-step) most integrators become numerical unstable. This is the reason for modlling really strong couplings using constraints - these are blind to that high frecuncy dynamics.

You might wan't to take a look at the following paper which explains this problem and suggests a stable real-time simulation method for volumetric deformable materiall based on finite elements and constraint regularization:

Interactive Simulation of Elastic Deformable Materials
M. Servin, C. Lacoursiére and N. Melin
SIGRAD 2006. The Annual SIGRAD Conference, Special Theme: Computer Games
http://www.ep.liu.se/ecp/019/005/index.html

[FD]

Surprisingly, I've been implementing the same method as described in the article. I thought that I was the first to use such an approach, that's why I didn't go into details.

The way the system is decomposed doesn't fix the problem of oscilating springs by itself. In fact what we get is very similar to a common techniqe of simmulating springs through joints by means of making a specific CFM matrix. It's basically the same Hooke's law, slighltly rearranged, and it is vulnerable to all the troubles with the intergator.

The way they make their simmulation stable is called "beta" in the article. They way the get to it is described very briefly in the papaer and the way one shall calculate it is not given at all. It seems that it is chosen arbitrarily, and its function is dissipating the energy of the system.

For high values of beta, the system behaves better in the sense that it doesn't explode at highly stiff materials. However, the simulation becomes inaccurate resulting in bodies being not stiff enough. Still I can't model stiff elastic bodies, unless I make implractically small timesteps, impractically high precision of the solver, and maybe even switch to double-precision floats.

What is the best way of choosing beta?

Are there any ideas how one can make stiff materials behave stiff?

[nisse]

1) Exactly what equations do you integrate? You derive a stiffness-matrix it seems. That means that you linearize the system in some way. I take it that you integrate an equation of the form:

m (d^2 delta_x / dt^2) = -K * delta_x - D * (d delta_x/dt)

That should work fine with a good choise of integrator for small displacements. But for large displacements you might run into problems,
In general the right hand side is not invariant under rotational transformations - so you cannot even simulate a slowly rotating body (even though it is not deformed) with this method. These effetcs are sometimes called "ghost forces". One cure can be found in the paper

Interactive Virtual Materials
Matthias Müller, Markus H. Gross: Interactive Virtual Materials. Graphics Interface 2004: 239-246
www.matthiasmueller.info/publications/GI2004.pdf

Another cure is not to linearize the model but to integrate the original equations from solid mechanics. This can be integrated at large time-steps using large time-steps using a variational integrator, e.g.,

Geometric, variational integrators for computer animation
L. Kharevych et al
Proceedings of the 2006 ACM SIGGRAPH/Eurographics symposium on Computer animation
http://portal.acm.org/citation.cfm?id=1218064.1218071

or to use the method of "constraint regularization" (as in "Interactive Simulation of Elastic Deformable Materials"). This is pretty close to a variational integrator but stability and efficiency comes with the prize of some numerical dissipation (increasing with the stiffness).

2) The "beta-parameter" is not important (at all) for stability in the approach in the paper "Interactive Simulation of Elastic Deformable Materials". The stability comes from the stepper you end up with after constraint regularization - Eq. (20). At high stiffness the high frequency normal modes that may develop into numerical instabilities are "filtered out" by the stepper and appears as energy dissipation. You should be able to set beta = 0 and achieve stable simulation for any value of Young's modulus. (At really high stiffness, however, you may run into problems solving the linear system depending on wether you solve Eq. (20) or Eq. (21) and what solver method you use). The beta-parameter can be used for increasing the dissipation, if you want.

3) You might have a bug in the implementation. The method should work for arbitrary stiffness. I advise that you work on an implement, of whatever method you choose, for a simple mass-spring system in 3D (so that you include the complications of rotational displacements) and get it working there first for the desired values of stiffness, mass and time-step. Then it is a simple matter (well, maybe not simple but managable) of extending this to finite elements.

[FD]

1) I solve the system very similar to eq. 21 in the paper "Interactive Simulation of Elastic Deformable Materials". There are three differences that I believe are not importaint
a)Instead of multiplying both sides of the equation on 1/dt to get away from (lambda * dt) to simply lambda, I do actually multiply S term on dt.
b)I do not diagonalize what they call "alfa" into Q^TDQ and do not rotate Jacobians by Q and I do not break up the volume term into two square roots and I do not distribute it between the jacobians - instead I use the non-diagonal alfa and the volume term is also accounted for in alfa matrix, not in the Jacobians. Alfa matrix still remains PD and must be ok.
c) I haven't read their formulas for Jacobians and I'm not sure if they are eqivalent to mine. I use the B-matrix as the Jacobian, which is a well-known interpolating functions derivatives matrix as it is defined in most engeneering books describing linear FEM method, e.g. it transforms the displacements of nodes into the linear part of the infinitesmall strain tensor in Voigt notation; quadratic part is omitted.

I do not use the "stiffness wrapping" yet, but I am only tryng to simmulate a stiff rod with one of its ends fixed to the world, so it must be ok, right? I though I would add the "stiffness wrapping" later if the method is proven to be stable in this case where there're no rotations.

I haven't studied the non-linearized methods that you give here. I beleve they are prohibitively slow, aren't they? In this case I don't want them.

2)The dissipative term DOES affect stability, because it does suppress oscilations that otherwise lead to divergence, but the system starts to behave unnaturaly.

From what I observe on my model, without the beta-term high young modulus leads to deivergence, which can be compensated by smaller timesteps (e.g. more precise integration) and more iterations of the solver (e.g. more precise solution). If the beta term is introduced, the system becomes stable even at big timetsteps and small number of solver iterations, but the bodies behave unnaturally.

I do not see how "the stability could come from the stepper I end up with after constraint regularization", and how could "at high stiffness the high frequency normal modes that may develop into numerical instabilities be "filtered out" by the stepper and appear as energy dissipation." Essentialy, what is the constrain regularizaton? It looks like an ordinary impulse-based solver arrangement like the one that is used for jointed rigid bodies systems.

3)The system I'm working on is itself very simple and it;s writte to test this algorithm. I think it shall be possible to find the bug, if any, in the system as it is.

How do you solve the system? Do you iterate untill you get a precise or rather precise solution? I'm using the typical approach where the solver terminates after a fixed nu,ber of iterations, and increasing this number of iterations leads to more stable systems, of course. Is the method from the paper supposed to be stable in case of substansial errors produced by the solver?

[raigan2]

This paper has a few details on how DMM works, it might be worth checking out: http://www.cs.berkeley.edu/b-cam/Papers ... -2009-RTD/

[bone]

Wow, very interesting paper on DMM! An absolute ton of work was put into that system.

[falcor]

Here's an early binary of the approach described in "Interactive Simulation of Elastic Deformable Materials" (Windows only): http://www.acc.umu.se/~melin/arbprov/Deformable.zip

[nisse]

1) I solve the system very similar to eq. 21 in the paper "Interactive Simulation of Elastic Deformable Materials". There are three differences that I believe are not importaint a)Instead of multiplying both sides of the equation on 1/dt to get away from (lambda * dt) to simply lambda, I do actually multiply S term on dt. b) I do not diagonalize what they call "alfa" into Q^TDQ and do not rotate Jacobians by Q and I do not break up the volume term into two square roots and I do not distribute it between the jacobians - instead I use the non-diagonal alfa and the volume term is also accounted for in alfa matrix, not in the Jacobians. Alfa matrix still remains PD and must be ok.
c) I haven't read their formulas for Jacobians and I'm not sure if they are eqivalent to mine. I use the B-matrix as the Jacobian, which is a well-known interpolating functions derivatives matrix as it is defined in most engeneering books describing linear FEM method, e.g. it transforms the displacements of nodes into the linear part of the infinitesmall strain tensor in Voigt notation; quadratic part is omitted.

Yes the dt-factors and diagonalization should be unimportant.

For large deformations it is, however, important to either include the nonlinear terms (that means using a sensible strain-tensor) or to do "stiffness warping" or equivalent. To achieve stability in the indefinitely stiff limit with the method in "Interactive Simulation of Elastic Deformable Materials" it is important write the deformation energy as a square of \phi

U = 1/2 \phi^T D \phi

\phi is then treated as a constraint, but regularized using the material parameters in D. D is constant in time. For large deformations the volume changes and the volume-factor must therefore be included in \phi.

I guess you can combine the "constraint regularization" approach with a linear FEM model (including only the linear terms in the Green strain tensor) and "stiffness wrapping". It would be interesting to see. The non-linear terms does not affect the performance much. But it does affect your implementation time.

2)The dissipative term DOES affect stability, because it does suppress oscilations that otherwise lead to divergence, but the system starts to behave unnaturaly....

From what I observe on my model, without the beta-term high young modulus leads to deivergence, which can be compensated by smaller timesteps (e.g. more precise integration) and more iterations of the solver (e.g. more precise solution). If the beta term is introduced, the system becomes stable even at big timetsteps and small number of solver iterations, but the bodies behave unnaturally.

I do not see how "the stability could come from the stepper I end up with after constraint regularization", and how could "at high stiffness the high frequency normal modes that may develop into numerical instabilities be "filtered out" by the stepper and appear as energy dissipation." Essentialy, what is the constrain regularizaton? It looks like an ordinary impulse-based solver arrangement like the one that is used for jointed rigid bodies systems.

I understand that the beta-term improves the stability if you allready have a stability problem. But the method should be stable without the term. "Constraint regularization" is presented quite thoroughly in the thesis

C. Lacoursière
Ghosts and machines: regularized variational methods for interactive simulations of multibodies with dry frictional contacts
http://umu.diva-portal.org/smash/record ... 0361&rvn=1

It is heavy reading if you go through all pages but you can probably see from the titles what chapters are relevant. The key is to first understand why a variational/symplectic stepper can preserve invariants (such as energy) exactly at finite time steps. The stepper in "Interactive Simulation of Elastic Deformable Materials" is a simplified version that is guarantied to be dissipative (but not as severely dissipative as linear implicit Euler).

How do you solve the system? Do you iterate untill you get a precise or rather precise solution? I'm using the typical approach where the solver terminates after a fixed nu,ber of iterations, and increasing this number of iterations leads to more stable systems, of course. Is the method from the paper supposed to be stable in case of substansial errors produced by the solver?

Well here is a sigificant difference which. We use a direct solver. This can be made very time-efficient if you use sparse data-representation and find ways to optimize the vector and matrix operations (using high-level BLAS). What type of iterative solver do you use? Gauss-Seidel? Conjugate gradient? I expect that Jacobi-iterations will fail severely, especially in the stiff limit.

[nisse]

This paper has a few details on how DMM works, it might be worth checking out: http://www.cs.berkeley.edu/b-cam/Papers ... -2009-RTD/

Indeed, this is very good work. Very nice to see to see it working in a commercial game and that the physics is essential for the game play as well.

I am interested in finding out more about the PID-controller approach to modelling frictional contacts. Do you know of anymore references on that besides the reference to Sellars in the Parker & O'Brien paper?

I see there are issues with the integration with a rigid body engine. Is the integration necessary? Is it possible to simulate indefinetly stiff tetrahedras with this framework, e.g., make a vehicle with wheel axes and suspension using constraints? In the stiff limit it should be possibility to make some shortcuts in the DMM computations to achieve the efficiency of rigid body simulations and automatic integration between deformable and stiff elements. But I guess it is critical what contact method and solver (direct/iterative) is used.

[FD]

To achieve stability in the indefinitely stiff limit with the method in "Interactive Simulation of Elastic Deformable Materials" it is important write the deformation energy as a square of \phi

U = 1/2 \phi^T D \phi

\phi is then treated as a constraint, but regularized using the material parameters in D. D is constant in time. For large deformations the volume changes and the volume-factor must therefore be included in \phi.

Unfortunately, I don't understand you.

You say it's importaint to write the deformation as a square of phi, but I neither write the deformation energy nor use phi explicitly. At the same time, I still don't get the idea behind regularization, I don't understand if it is something beyond the usual Lagrangian form as opposed to penalty-pased methods.

I'll try to say in a few words how I see the problem. There is a very common approach to (rigid body) dynamics that implies that constarint forces are reprsented in the form Fc = J^t * lambda, where J is (dg/dxi) matrix of the derivatives of all constraints given in the form g(x) = 0 with respect to each of the coordiantes xi of each body of the system. At the smae time, J * v + C * lambda = rhs must be satisfied. This yields the system
M*a - J^t * lambda = Fext;
J * v + C * lambda = rhs;
After approximating Ma = M(v - v0) /dt, this yields
M * v - J^t * lambda * dt = M*v0 + dt * Fext;
J * v + C * lambda = rhs;

This can be solved by evaluating v from the upper equation and plugging it into the lower equation and solving for lambda. No headache about energy potentials, no need to get deeper into the justification of this method.

Since it can be proved that F = K*(x - x0), where K is the stiffness matrix and F is all the force including inertia that affect the nodes of the FEM-represented body, we denote Felast = - K(x-x0). Since for linear finite elements the derivatives of shape functions are constant, K can be written as K = B^t * D * B * V, where D is the matrix that expresses the relation between a strain tensor and a stress tensor for this material and represents the Hooke's law and V is the volume of the finite element. B matrix tranansforms the offsets of the nodes to the infintesmall strain tensor, omitting the rotational part of the tensor and representing the tensor in Voigt notation. If we approximate xnew = x + dt * v and put the Felat into the Lagrangian dynamics problem as given above, we arrive with -dt * B^t * D * B * V * ( v * dt + x - x0) term in the upper part of the system. We can denote lambda1 = - D * B * V * ( v * dt + x - x0), so that the term transforms into B^t * lambda1 * dt, which is very similar to the usual form used to express constraints. Similarly, lambda1 = - D * B * V * ( v * dt + x - x0) can be rearranged B * v + D^-1 * 1/V * 1/dt * lambda1 = -1/dt * B(x - x0). It is essentialy the same velocity-based constraint equation with CFM = D^-1 * 1/V * 1/dt and RHS = -1/dt * B(x - x0). I've got here without thinking about the deformation potential energy at all.

So here is the question: where is the place for "regularization" in my logic? Does the form of the Lagrangian dynamics problem given above solve the problem by itself or not? I'm used to it and it seems to me that all the talking about regularization is actually concentrated on advocating this particular form of the problem, am I right? Or is there anything beyond it?

Second, you say that for large deformations volume changes must be included in phi. But phi is nothing more than a deformation tensor! Of course since there're no second-orper terms it will be incorrect for large deformations, but stiffness wrapping basically must solve this problem more or less. In either case, I see no place for volume in phi since volume in K appears from the integration of B^t*D*B (constant) over the volume of the finite element, and K remains constant both for small deformations approach and stiffness wrapping (it's wrapped into rotation matrices, but not recalculated). Why shall the volume term be in phi and in what way?

Third, I think that a direct solver must be a remedy, since using very high number of iterations makes the system almost unconditionally stable. The challange is however to combine FEM with the regular rigid body approach, so it would be nice if the system didn't explode at an approximate solution of the linear system. As far as I understand, an approximate solution leads to high forces on the next timeframe and the system explodes very quickly.

[nisse]

It seems the use of iterative v.s. direct solver might explain the difference in stability. But I am not sure. I use iterative solver sometimes as well and don't recall any particular difficulties with stability in the stiff limit. The problem usually shows up as the material being overly elastic when the error is big rather than blowing up in your face. I still think the regularization procedure might be part of the difference in stability as well.

I understand your questioning of the regularization involning energy potentials. It seems to be just an overhead of theory and computations ending up in almost the same system. But that "almost" is plays a huge difference sometimes, e.g., in the stiff limit and for over-constrained systems. I will try to explain the point of it.

You say you solve the system

M*a - J^t * lambda = Fext;
J * v + C * lambda = rhs;

This is regularized form of the constraints and the C-parameter is regularization (relaxation) parameter. Small C -> stiff constraint and large C -> soft constraint. But there is a family of equations of this form that have different stability properties. The difference show up in the terms in the right hand side (rhs) and comes from how you or if you do "index reduction" of the original set of differential algebraic equations (DAE). What is your right hand side? Does it involve the constraint phi (or called g or c or something)? If it doesn't you are probably using the classical "acceleration method" that can be found in a number of places. This is not flawed by constraint drifting. This is cured by Baumgarte stabilization - you add to lambda two terms "alpha*phi + beta*d_phi/d_t" proportional to the constraint and its time derivative. This is not stable for any alpha and beta and teh tuning of these parameters is difiicult and virtually impossible to map to real material parameters (lie Young's modulus of elasticity). The regularization term C improves the stability.

It was proved long ago (RUBIN and UNGAR, 1957) that a constraint can be seen as the limit of a penalty force whith the stiffness going to infinity. By penalty force I mean a force of the type

F = - grad U

where U is a quadratic potential U = 1/2 * k * phi(x)^2 of some stiffness k and grad is the spatial derivative (of x). We call this energy function or potential energy. Observe that that the simple spring force have this form. If d is the spring displacement, the potential is U = 1/2 * k * d^2 with spring constant k and the spring force becomes F = - k * d. At high stiffness the potiential forces the system towards (or rather oscillating about) the constraint surface phi(x) = 0.

The numerical solvers in the references I cited previously in this post utilize this connection between constraints and energy functions, which usually contain the material parameters - the spring coefficient k can be measured experimentalyy or looked up in tables for different materials. As a result you end up with particular terms and coefficients C and on the right hand side. It is important that the Jacobian is precisely the gradient of the constraint occuring in rhs and that the coefficients are also important. You can make a similar approach to dissipative terms (using Rayleigh dissipation functions rather than energy potentials). There are two advantages of the approach:

1) It is very stable. Lacoursière has a proof of linear stability. The stability is due to the lagrangian/variational approach tayloring the integrator to preserve the symmetries of the original physical system.

2) You automatically regularize using parameters that are real material parameters that can be found from experiments/tables, e.g., elasticity and energy dissipation rate. You don't have to tune blindly.

I hope this helps. If you have a simple spring-version of the method implemented you should be able to see these differences. And for a simple spring you can study how stability depends on the number of iterations and benchmark with the direct solver
(the solve becomes a simple division for the spring system).

[Dirk Gregorius]

What you are describing here is exactly the ODE approach where the regularization is basically called CFM. Neither is the ODE infinitely stiff or stable! The same holds for Bullet and all other physics engines I am aware of. I think this is mainly due to the iterative solver which is usually run with 4-16 iterations. Also note that index reduction always results in drift. The regularization is getting some attention lately though it is around for quite a long time. As you say it makes the constraints soft and in my opinion this is not what people want (e.g. stretchy ragdolls).

Actually Erin showed lately that a seperate projection sweep over the postions constraints gives much better quality and stability than the index 1 form you are suggesting here. For a discussion see here:http://www.gphysics.com/archives/35

You might get very good results with a direct solver with the above system, but I think that this is unfeasible for games.

-Dirk

[nisse]

...it makes the constraints soft and in my opinion this is not what people want (e.g. stretchy ragdolls).-Dirk

Well, ideal is to be able to parametrize the joint stiffnes ranging from soft to indefinetly stiff AND use conventional material parameters for this (e.g. Young's modulus etc). In this thread we where speaking of deformable materials - but I think there are situations where you also wan't well defined elasticity in the joints between rigid bodies as well (robots, ragdolls, ropes - they are typically stiff but under large enough tension they are elastic). There are several ways to do regularization. Starting from energy potentials give you the link to these parameters. I have not seen this been shown for the CMF-approach (I have not ruled out that it is possible to make the connection on the other hand). Please correct me on this, I am very interested in finding such material!

You might get very good results with a direct solver with the above system, but I think that this is unfeasible for games.

We can probably have a long debate wether direct solvers will ever be feasible for games Regardless of this, there is good experience with stiffness and stability with this approach also using iterative methods, e.g., the physics engine in Phun/Algodoo (http://www.phunland.com). It is for the 3D-deformable material I lack stability data for iterative solver.

Instabilities can still occur if you strain the objects enough and there are degrees of freedoms left that can start oscillate, e.g., hanging a massive (>1000 kg) rigid body in chain of light bodies (1 kg). This instability occurs both for iterative and direct solvers.

There are several cures for this, e.g., projection methods which you mention (thanks for the link, haven't seen it before). This is what makes the Jacobsen method successful. A drawback of projection methods is that it is difficult to combine with predictable joint elasticity, especially if you want to be able to read off and use the deformation forces (maybe not a common thing to do in a game, though). I think for games, the Position Based Dynamics approach (Müller & Co) is a very good one as well.

[FD]

Dirk,

I wonder if there're stabilization techniques for springs. If so, they might be applicable for FEM as well. You see, all the techniques listed in the abovementioned post by Erin imply that we do know the correct posiitons and are able to estimate the position error in each joint. However, in case of a spring any position might be valid and will produce the adequate response force. For me it seems not clear at all how one can estimate and correct the position error in a spring while keeping the spring behavior.